# Anisotropy of Yu-Shiba-Rusinov states in NbSe2

Mateo Uldemolins, Freek Massee, Tristan Cren, Andrej Mesaros, Pascal Simon

► To cite this version:

Mateo Uldemolins, Freek Massee, Tristan Cren, Andrej Mesaros, Pascal Simon. Anisotropy of Yu-Shiba-Rusinov states in NbSe2. Physical Review B, 2024, 110 (22), pp.224519. 10.1103/Phys-RevB.110.224519. hal-04870834

HAL Id: hal-04870834

https://hal.science/hal-04870834v1

Submitted on 27 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Anisotropy of Yu-Shiba-Rusinov states in NbSe$_2$

Mateo Uldemolins,$^{1,2,*}$ Freek Massee,$^2$ Tristan Cren,$^3$ Andrej Mesaros,$^2$ and Pascal Simon$^2$

$^1$Donostia International Physics Center (DIPC), 20018, Donostia-San Sebastián, Spain

$^2$Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405, Orsay, France

$^3$Sorbonne Université, CNRS, Institut des Nanosciences de Paris, UMR7588, F-75252 Paris, France

(Dated: December 24, 2024)

The spatial structure of in-gap Yu-Shiba-Rusinov (YSR) bound states induced by a magnetic impurity in a superconductor is the essential ingredient for the possibility of engineering collective impurity states. Recently, a saddle-point approximation [Phys. Rev. B 105, 144503] revealed how the spatial form of a YSR state is controlled by an anisotropic exponential decay length, and an anisotropic prefactor, which depends on the Fermi velocity and Fermi-surface curvature. Here we analyze scanning tunnel microscope (STM) data on YSR states in NbSe$_2$, focusing on the key issue that the exponential decay length predicted theoretically from the small superconducting gap is much larger than the observed extent of YSR states. We confirm that the exponential decay can be neglected in the analysis of the anisotropy. Instead, we extract the anisotropic prefactor directly from the data, matching it to the theoretical prediction, and we establish that the theoretical expression for the prefactor alone captures the characteristic flower-like shape of the YSR state. Surprisingly, we find that up to linear order in the superconducting gap the anisotropic prefactor that determines the shape of YSR states is the same as the anisotropic response to the impurity in the underlying normal metal. Our work points out the correct way to analyze STM data on impurities in small-gap superconductors, and reveals the importance of the normal band structure's curvature and Fermi velocity in designing multi-impurity in-gap states in superconductors.

## I. INTRODUCTION

The local response to single magnetic impurities in conventional superconductors appears most directly in the form of pairs of in-gap excitations localized on the impurity known as Yu-Shiba-Rusinov (YSR) states [1-3]. YSR states are nowadays routinely probed by scanning tunneling microscopy (STM) experiments [4] and, while an old problem, they have attracted a renewed wave of interest owing to their potential to extract information on the properties of the host [5-7] and study magnetism [8, 9] and spin-transport phenomena [10-13] at the atomic scale. At the same time, much experimental and theoretical effort has been devoted to designing complex structures in real space (e.g., impurity chains or lattices) [14-27], which are suggested to realize topological states of matter harboring Majorana zero modes [28]. Understanding the spatial structure of YSR states is key for developing this program.

In an isotropic $d$-dimensional system, the local density of states (LDOS) of YSR states far away from the impurity behaves according to $\sim e^{-r/\xi}/r^{d-1}$, with the superconducting coherence length controlling the exponential decay, $\xi \sim \frac{v_F h}{\Delta} \sim \xi_{SC}$. The dependence of the power-law on the dimensionality of the substrate has been evidenced by STM measurements on two-dimensional superconductors ($d=2$), e.g., NbSe$_2$, where the in-gap LDOS extends over several nanometers [29-31], thereby rendering this class of materials a promising platform for constructing collective impurity states [32-36]. In general, the measured in-gap LDOS oftentimes displays a

marked anisotropy on length scales largely surpassing the typical impurity size and hence attributed to the lower symmetry of the band structure [30, 37]. The theory of spatial structure of YSR states is hence intricate, requiring a blend of information about the substrate's band structure, pairing function, and impurity coupling.

Recently, the theory of the quasiparticle focusing effect (QFE) [38, 39] was extended to superconducting substrates by some of the authors [40], leading to the following insights in the limit of weak pairing, $\Delta/\epsilon_F \ll 1$, and at long distances ($rk_F \gg 1$):

1. The QFE persists, i.e., the response to the impurity in a given real-space direction is controlled by momenta (denoted *critical*) on the Fermi surface at which the Fermi velocity lies parallel or antiparallel to that real-space direction. This selectivity in reciprocal space controls the spatial anisotropy of the YSR states; in particular, YSR states are focused in directions perpendicular to flatter sections of the Fermi surface.

2. Specifically to gapped superconductors, the in-gap YSR LDOS has an additional exponentially decaying factor with angularly dependent exponential decay length, apart from a prefactor and oscillations, all related to the geometry of the Fermi surface through a simple analytical expression. Importantly, the prefactor and exponential decay length are essentially always in phase, i.e. they agree on the direction along which the LDOS is focused.

This theory describes very accurately the anisotropy of the characteristic exponential decaying length of YSR states simulated on a lattice [40] and successfully predicts the main anisotropy directions of YSR states ob-

<footnote>* mateo.uldemolins@dipc.org</footnote>
served in NbSe$_2$. However, the tight-binding description of the Nb-derived bands at the Fermi level [41] in combination with the small superconducting gap ($\Delta \sim 1$ meV [42, 43]), yields an in-gap LDOS with large exponential decay length ($\xi_{TB} \sim 100$ nm) [29, 33, 35, 36] that widely surpasses the extent of the observed YSR states, as well as the superconducting coherence length reported in this material ($\xi_{SC}^{\text{exp}} \sim 10$ nm [44–46]). Therefore, there is a problem in the natural interpretation that the anisotropy of the coherence length in NbSe$_2$ fundamentally controls the shape of the YSR LDOS, and theoretical radial fitting of an exponential decay is questionable.

In this work, we show that the angular dependence of the LDOS prefactor accounts for the anisotropy of YSR states in this material, while the radial fitting of the LDOS can be circumvented, i.e., the (anisotropic) exponential decay can be neglected. This evokes the QFE in normal metals [38, 39] where standard Friedel oscillations decay algebraically and lack the exponential factor. In fact, our analytical expression of the prefactor has the same functional form, in terms of the Fermi velocity and curvature of the Fermi surface, as the LDOS in the case of a normal metal. A difference in principle arises when this expression is evaluated at the critical momenta, which vary as the band structure of the metal changes into the one of the superconductor. Surprisingly, we find that up to linear order in $\Delta/\epsilon_F \ll 1$ the prefactor of the impurity-induced LDOS in the superconductor exactly equals the one in the parent normal metal. We show that this prefactor can be reliably extracted from STM data on NbSe$_2$ and that the characteristic flower-shaped YSR LDOS is faithfully modeled by the prefactor alone. We also find that introducing the effect of the charge density wave (CDW) on the bandstructure does not play a significant role in the YSR LDOS.

The rest of the paper is organized as follows: we start by presenting the STM data in Sec. II. Then, in Sec. II A, we briefly recall the main analytical expressions derived in [40] and we analyze the spatial anisotropy of the LDOS prefactor. This analysis is complemented in Sec. II B, where we study the Fourier transform of the YSR states. We conclude with a short discussion in Sec. III. Additional details concerning the experiment and the theory expressions are presented in the Appendixes.

## II. RESULTS

We consider STM data on 2H-NbSe$_2$ crystals containing a few tens of ppm of magnetic atoms measured at 300 mK with a metallic tip (for additional details concerning the fabrication and measurement conditions, see Appendix A). YSR bound states manifest as an in-gap peak in the differential-conductance spectrum at the impurity location (see Appendix A). In Fig. 1, we present atomically resolved data around two independent impurities on the same sample. The constant-current image [Figs. 1(a) and 1(c)] shows the triangular geometry of

<figure><img src="image_2.png" /><figcaption>FIG. 1. Spatial data sets 1 and 2. (a) Atomically resolved constant-current image of the Se surface of 2H-NbSe$_2$, $V_{set} = 4$ mV and $I_{set} = 160$ pA. (b) In-gap current recorded at $V = -0.6$ mV showing a sixfold YSR state. (c),(d) The same as (a),(b) for another YSR state, with $V_{set} = 4.2$ mV, $I_{set} = 100$ pA, and $V = -0.6$ mV, respectively. The lines at the bottom-left corner indicate the crystallographic axes of 2H-NbSe$_2$. Their length corresponds to 2 nm in both data sets.</figcaption></figure>

the Se lattice without any visible imperfection, implying that the defect giving rise to the YSR state is most likely an atom substitution in the Nb layer. The corresponding current images recorded before the onset of the superconducting coherence peak [Figs. 1(b) and 1(d)] show a marked six-fold anisotropy that constitutes an excellent example of the QFE on a superconductor.

### A. Characterization of the LDOS prefactor in real space

Let us begin by briefly recalling the theoretical framework to describe the QFE on superconductors. We model the impurity-substrate system as a classical spin on an s-wave superconductor.

The Hamiltonian reads,

$$H = H_0 + H_{\text{imp}}, \quad (1a)$$

$$H_0 = \sum_{\mathbf{k}\sigma} \varepsilon_{\mathbf{k}\sigma} c_{\mathbf{k}\sigma}^{\dagger} c_{\mathbf{k}\sigma} + \Delta \sum_{\mathbf{k}} c_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger} + \text{h.c.}, \quad (1b)$$

$$H_{\text{imp}} = -J \left( c_{\mathbf{r}_0 \uparrow}^{\dagger} c_{\mathbf{r}_0 \uparrow} - c_{\mathbf{r}_0 \downarrow}^{\dagger} c_{\mathbf{r}_0 \downarrow} \right), \quad (1c)$$

where $\varepsilon_k$ is the spin-degenerate normal energy dispersion of the substrate, $\Delta$ is the BCS superconducting gap and $J$ is the amplitude of the magnetic exchange coupling
between the superconducting electrons and the spin of the impurity at $\mathbf{r}_0$. Specifically, $\varepsilon_k$ is an effective fifth-nearest neighbors tight-binding energy dispersion on a triangular lattice that describes one of the two lowest-lying Nb 4d bands of 2H-NbSe$_2$ [41] (see Appendix C). This description neglects interlayer coupling of the bulk crystal and assumes that the magnetic impurity couples to one band only [29, 33]. It also neglects spin-orbit coupling [47]; however, it provides a faithful description of the geometry of the Fermi surface and band structure of the substrate, which is the crucial element to characterize the QFE.

An important fact about the band structure of niobium diselenide is that it develops CDW at $T_{CDW} = 33$ K, which coexists with superconductivity ($T_c = 7.2$ K) at low temperatures [48–50]. Although density functional theory indicated a quite large deformation of the Fermi surface due to the CDW [51], for the particular case of 2H-NbSe$_2$ the ARPES data show that the shape and size of the Fermi contours are largely preserved well below the CDW transition temperature [41]. Therefore, neglecting the effect of CDW on the long-range spatial properties of the YSR LDOS should be a reasonable simplification as well. For completeness, in Appendix D we discuss in detail the effect of the CDW within our approach. We hence confirm that although the CDW introduces interesting effects in principle, for the effective parameter values derived from ARPES there is no change to any of our main conclusions. Hence for simplicity we neglect the presence of CDW in the rest of the main text.

In addition, Eq. (1c) assumes a fully local and isotropic scatterer, thereby neglecting the magnetic ion's d-orbital structure and the adsorption site's symmetry. These add to the YSR-state wave function some degree of anisotropy [8, 52]; however, the observed six-fold anisotropy of the LDOS [Figs. 1(b) and 1(d)] occurs on length scales of several nanometers that largely surpass the typical orbital size; therefore, it can be safely ascribed to the QFE and justifies the model's simplification (see also Appendix D).

The LDOS at a distance $\mathbf{r}$ from the impurity and at the bound-state energy $E_S$ is computed by means of a saddle-point approximation on the bare propagator valid in the far-field limit ($rk_{F,min} \gg 1$, with $k_{F,min}$ the minimum Fermi wave vector of the Fermi contour). By further taking the small-gap limit ($\Delta/\varepsilon_F \ll 1$), the LDOS takes the approximate form

$$ \rho(\mathbf{r}) \sim \frac{1}{r} \sum_{j,j'} \Gamma_{j,j'}(\theta_r) e^{-r/\varepsilon_{j,j'}(\theta_r)} f[\mathbf{k}_{j,j'}^{\pm}(\theta_r) \cdot \mathbf{r}], \quad (2) $$

where the angular-dependent prefactor reads

$$ \Gamma_{j,j'}(\theta_r) = \frac{1}{|\nabla \varepsilon_{\tilde{k}_j(\theta_r)}| \sqrt{\kappa_{\tilde{k}_j(\theta_r)}} |\nabla \varepsilon_{\tilde{k}_{j'}(\theta_r)}| \sqrt{\kappa_{\tilde{k}_{j'}(\theta_r)}}}, \quad (3) $$

with $|\nabla \varepsilon_{\tilde{k}_j(\theta_r)}|$ the Fermi velocity and $\kappa_{\tilde{k}_j(\theta_r)}$ the curvature of the Fermi contour evaluated at $\tilde{k}_j(\theta_r)$. The critical momenta $\tilde{k}_j(\theta_r)$ are solutions to the saddle-point

equations in the small-gap limit, namely, points on the Fermi contour for which the Fermi velocity is parallel to $\mathbf{r}$,

$$ \frac{\nabla \varepsilon_{\tilde{k}_j(\theta_r)}}{|\nabla \varepsilon_{\tilde{k}_j(\theta_r)}|} = \hat{\mathbf{r}}, \quad (4) $$

hence their dependence on the polar angle $\theta_r$ which controls the spatial anisotropy of the LDOS. We recall that the gapped nature of the superconductor yields complex solutions to the saddle-point equations [40]. However, up to and including the linear order in $\Delta/\varepsilon_F$, these solutions [i.e. the vectors $\tilde{k}_j(\theta_r)$] are real and positioned on the Fermi contour of the underlying metal. Therefore, the prefactor $\Gamma_{j,j'}(\theta_r)$ in Eq. (3) in the small-gap limit is simply determined by the band structure of the normal metal (see Appendix B). The exponential decay length $\xi_{j,j'}(\theta_r)$ depends on the Fermi velocity at $\tilde{k}_j(\theta_r)$ and $\tilde{k}_{j'}(\theta_r)$, and $f[\mathbf{k}_{j,j'}^{\pm}(\theta_r) \cdot \mathbf{r}]$ is a linear combination of oscillating functions with typical characteristic inverse lengths $\mathbf{k}_{j,j'}^{\pm}(\theta_r) = \tilde{k}_j(\theta_r) \pm \tilde{k}_{j'}(\theta_r)$. The indices $j, j' = 1, \dots, N$ label the various solutions to Eq. (4). The form in Eq. (2), derived in Appendix B, is more compact than the form first derived in Ref. [40], as it pairs up the two contributions from the wavevectors whose respective Fermi velocities are parallel and antiparallel to the observation direction [the latter satisfying Eq. (4)] with $\hat{\mathbf{r}} \rightarrow -\hat{\mathbf{r}}$. The Fermi contour of NbSe$_2$ consists of three disconnected pockets centered at $\Gamma, K, \text{ and } K'$ and, therefore, $N=3$.

To analyze the spatial anisotropy of the STM data we consider radial averaging in annuli around the impurity and hence introduce the following quantity:

$$ I_{\text{int}}(\theta_r) = \frac{1}{M} \sum_{\mathbf{R}_m \in M(\theta_r)} |\mathbf{R}_m| \cdot I(\mathbf{R}_m), \quad (5) $$

where $I(\mathbf{R}_m)$ is the measured current at point $\mathbf{R}_m$ sitting at a distance $|\mathbf{R}_m|$ from the impurity, while $M(\theta_r)$ is the set of $M$ points belonging to an angular sector centered on the impurity, having infinitesimal angular width $d\theta$, an inner radius $R_0$, and outer radius $R_0 + \delta$. As long as $R_0$ is sufficiently large to make the asymptotic approximation reasonable ($R_0 \gg 1/k_F^{\text{min}} \sim 0.2$ nm for the Fermi contour under consideration), we can relate $I_{\text{int}}(\theta_r)$ to the radially integrated LDOS,

$$ I_{\text{int}}(\theta_r) \sim \rho_{\text{int}}^{R_0, \delta}(\theta_r) = \frac{1}{\delta} \int_{R_0}^{R_0+\delta} dr r \rho_{\text{max}}(\mathbf{r}), \quad (6) $$

where $\rho_{\text{max}}(\mathbf{r})$ is the term in Eq. (2) with dominant prefactor, $\Gamma_{\text{max}}(\theta_r) = \max\{\Gamma_{j,j'}(\theta_r)\}$. This simplification is legitimate as the various prefactors $\Gamma_{j,j'}(\theta_r)$ are essentially in phase (see Appendix C for further details). We can further simplify Eq. (6) by making the following assumptions:

* The integration range $\delta$ is larger than the typical oscillation length; therefore, $f[\mathbf{k}_{j,j'}^{\pm}(\theta_r) \cdot \mathbf{r}]$ is averaged out and can be safely ignored.
* The exponential decay length is sufficiently large so that $e^{-r/\xi(\theta)} \sim 1$ over the integration range. This demands $R_0 + \delta \le \xi_{SC}^{\text{exp}} \sim 10$ nm $\ll \min\{\xi_{j,j'}(\theta_r)\}$. This approximation also relies on the fact that the angular dependence of the prefactor and the exponential decay length are in phase and therefore there are not any competition effects.

Under these assumptions $\rho_{\text{int}}^{R_0,\delta}(\theta_r)$ becomes independent of $\delta$ and, further, we have

$$ \rho_{\text{int}}^{R_0,\delta}(\theta_r) \sim \Gamma_{\text{max}}(\theta_r). \qquad (7) $$

<figure><img src="image_3.png" /><figcaption>FIG. 2. Angular dependence of the LDOS prefactor. Panels (a) and (b) correspond to data sets 1 and 2, respectively. Purple-to-yellow curves represent experimental data points $|I_{\text{int}}(\theta)|$ for $\delta$ values specified in the legend and $d\theta = 2\pi/100$. Legend inset shows the real-space current data with a lattice vector aligned with the horizontal axis ($\theta_r = 0$). We excluded a range of data points in data set 1 due to interference with another YSR state (partially visible on the bottom right of the image).</figcaption></figure>

In Fig. 2 we present the result of our analysis in comparison with the model curve

$$ \Gamma_{\text{model}}(\theta_r) = a\Gamma_{\text{max}}(\theta_r) + a_0. \qquad (8) $$

Here $a$ is a free parameter in the theory, and $a_0$ is constant shift on the order of the background signal. We estimate the latter by taking the signal average far away from the impurity times the typical integration length $\bar{R}_n \sim R_0 + \frac{\delta_{\text{min}}+\delta_{\text{max}}}{2}$.

Experimental data points $\rho_{\text{int}}^{R_0,\delta}(\theta)$ above the background signal collapse onto a universal curve, thereby supporting the assumptions stated in the previous paragraph.

The relative angular position of the model curve with respect to the data points is not a degree of freedom, but it is determined by the orientation of the lattice vectors (see legend insets in Fig. 2), the latter being extracted from the constant-current image of the atomic surface. We emphasize the agreement between the positions of the maxima and minima of the model and the experimental curves. While the six-fold nature of the YSR state may follow trivially from the symmetry of the band structure, its orientation with respect to the lattice does not. One could generally expect YSR states to be oriented according to high-symmetry lines; however, it is a priori equally reasonable for the “petals” of the LDOS to be along the lattice vectors or in between. The theory shows how this is, in fact, determined by the geometry of the Fermi surface and Fermi velocity.

We note that the periodicity of the data points with $\theta_r$ exhibits certain higher harmonics, namely, the signal becomes thinner at the crests and wider at the valleys, in other words, the “petals” are narrow. This nontrivial feature is also present in the model curve.

## B. Fourier-space analysis

To further confirm the validity of our interpretation, we study the Fourier transform (FT) of the in-gap recorded current [Fig. 3 (a)], and we compare it against the FT of the asymptotic expression of the LDOS, which contains the prefactor, the exponential decay, and the oscillations,

$$ \tilde{\rho}(\mathbf{q}) = \sum_{j,j'} \int_{0}^{2\pi} d\theta_{\mathbf{r}} \Gamma_{j,j'}(\theta_{\mathbf{r}}) \tilde{F}_{j,j'}(\mathbf{q}, \theta_{\mathbf{r}}), \qquad (9) $$

with $\tilde{F}_{j,j'}(\mathbf{q}, \theta_r) \equiv \int d\mathbf{r}e^{-r/\xi_{j,j'}(\theta_r)} f[(\mathbf{k}_{j,j'}^{\pm}(\theta_r) - \mathbf{q}) \cdot \mathbf{r}]$ the integral over the radial coordinate $r$ (see Appendix E for details).

The experimental FT of the YSR states exhibits two salient features that are qualitatively captured by the model [cf. Figs. 3(a), 3(d) and 3(b), 3(e)]. The first comprises the linear segments parallel to the edge of the hexagonal first Brillouin zone (FBZ), corresponding to the short wavelength oscillations of the YSR state. The second is the six-fold star at small $\mathbf{q}$ whose orientation is rotated by $\pi/6$ with respect to the YSR state in real-space (cf. Fig. 1) as expected from the fact that a focused signal in real-space gives a focused feature in Fourier space that is in orthogonal direction. The small star in momentum space is the manifestation of the QFE at large length scales. Finally, we corroborate that the long-range anisotropy of the YSR state, i.e., the small-$\mathbf{q}$ star feature in Fourier space, is controlled by the prefactor in the asymptotic expression of the LDOS: we set
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Fourier transform's modulus of the in-gap current shown in Fig. 1 (b). The blue hexagon indicates the FBZ and the blue cross marks one of the Bragg peaks. (b) Fourier transform's modulus of the analytical expression of the LDOS [Eq. (9)] evaluated on the FBZ. (c) Same as (b) with constant prefactor. (d), (e), (f) Zoom around the origin delimited by the dashed-line square in panels (a), (b), and (c), respectively.</figcaption></figure>

$\Gamma_{j,j'}(\theta_r) = 1$ in Eq. (9) and obtain a FT map which differs from the full $\tilde{\rho}(q)$ exactly in that it lacks the small-$q$ sixfold star [Figs. 3(c), 3(f)].

ACKNOWLEDGMENTS

M.U., A.M., and P.S. acknowledge the support of the French Agence Nationale de la Recherche (ANR), under Grant No. ANR-22-CE30-0037. M.U. acknowledges funding by the European Union NextGenerationEU/PRTR-C17.I1 and the Department of Education of the Basque Government (IKUR Strategy). F.M. acknowledges L. Cario for providing the samples and funding from H2020 Marie Skłodowska-Curie Actions (Grant No. 659247) and the ANR (Grant No. ANR-16-ACHN-0018-01).

III. CONCLUSIONS

We show that the theory of the quasiparticle focusing effect in two-dimensional superconductors [40] can be further simplified in application to STM data of superconductors such as 2H-NbSe₂, in which the coherence length derived from the band structure exceeds the spatial extent of YSR states. Our conclusions are established both by a practical analysis of annuli in real space, and by analysis of features in Fourier space.

The key physical quantity for all practical purposes becomes the prefactor $\Gamma(\theta_r)$ of the LDOS expression, which faithfully captures the symmetry, orientation, and angular profile of YSR states observed in STM experiments. The prefactor has a simple analytical expression depending on the Fermi velocity and curvature of the Fermi surface only. In particular, the LDOS around a magnetic impurity is preferentially focused in directions perpendicular to flatter sections of the Fermi contour, and in the small-gap limit it matches the QFE in the underlying normal metal [38, 39]. It would be interesting to test this matching in future measurements of local response to impurities in the normal metallic state of NbSe₂.

In conclusion, this work demonstrates the applicability of asymptotic approximations in the spirit of [38–40] for describing the local response to point defects, and underlines the importance of the band structure geometry for designing collective impurity states.
## Appendix A: Experimental details and additional data

2H-NbSe$_2$ single crystals were grown using an iodine transport method and were unintentionally doped by a few tens of ppm of magnetic species (Fe, Cr, Mn) contained in the niobium precursor; see Ref. [29] for more details. The crystals were mechanically cleaved in cryogenic vacuum at T ~ 20 K and directly inserted into the STM head [53] at 4.2 K. An etched atomically sharp and stable tungsten tip was used for all measurements, which were performed at T = 0.3 K. In Fig. 4 we present the $dI/dV$ spectra for the two impurities considered in this work and in Figs. 5(b), 5(d) we analyze the anisotropy of the signal recorded at positive bias to complement the data presented in the main text. We note that the spatial maps presented in this work show I (as opposed to $dI/dV$), but recorded in the fully gapped region between 0 mV and the coherence peak, such that only the signal from the sub-gap state is present.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Tunneling spectrum taken at the core of the YSR state in Figs. 1(a), 1(b) (a) and Figs. 1(c), 1(d) (b) showing the strongest in-gap peak at $V_S = -0.1$ meV and $V_S = 0.25$ meV, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Same analysis of the LDOS prefactor as in Fig. 2 of the main text for the in-gap current recorded at positive bias, V = 0.6 mV. Panels (a) and (b) correspond to data sets 1 and 2, respectively. The phenomenological parameters in the model curve [Eq. (8)] for panels [Figs. 2(a), 2(b), and Figs. 5(a), 5(b)] are: $a_0 = (-3, -3, -2.5, -2.1) \cdot 10^{-12}$ pA · nm and $a/a_c = (4.45, 4.45, 3.2, 4.45) \cdot 10^{-7}$ pA · nm · meV$^{1/2}$, with the model curve computed using Eqs. (C1) to (C3b) and “Band 2” parameters in Table I.</figcaption></figure>

## Appendix B: Derivation of the LDOS expression

We derive the expression for the in-gap LDOS in Eq. (2), which is equal in value to the main result from [40], but is further transformed to make explicit that LDOS is real, and that the doubling of critical momenta in the superconductor with respect to the normal metal does not essentially alter the prefactor. We further show that the prefactor up to linear order in $\Delta/\epsilon_F$ matches the one of the normal metal.

The bare propagator at energy $E < \Delta$ from $r_a = r$ to $r_b = 0$ on an s-wave superconductor (described by the Hamiltonian in Eq. (1b) of the main text) in Nambu space can be approximated in the far-field limit ($rk_{F,min} \gg 1$, with $k_{F,min}$ the minimum Fermi wave vector of the Fermi contour) as
$$ \hat{G}_0(\mathbf{r}, 0; E) \sim \frac{1}{\sqrt{r}} \sum_{j, \epsilon = \pm} \Gamma_{j, \epsilon}(\theta_r) \exp \left\{ -\frac{r}{\xi_{j, \epsilon}(\theta_r)} + i \left[ \mathbf{r} \cdot \text{Re}[\mathbf{k}_{j, \epsilon}(\theta_r)] + \epsilon \frac{\pi}{4} \right] \right\} \begin{pmatrix} E + i\epsilon\omega & \Delta \\ \Delta & E - i\epsilon\omega \end{pmatrix}, \quad (B1) $$

where $\theta_r$ is the polar angle on the plane defined by the vector $\mathbf{r}$ from the impurity to the point where we seek to evaluate the LDOS, $\omega = \sqrt{\Delta^2 - E^2}$ and

$$ \Gamma_{j, \epsilon}(\theta_r) = \frac{1}{|\nabla \epsilon_{\mathbf{k}_{j, \epsilon}(\theta_r)}| \sqrt{\kappa_{\mathbf{k}_{j, \epsilon}(\theta_r)}}}, \quad (B2a) $$

$$ \xi_{j, \epsilon}(\theta_r) = \frac{1}{\text{Im}[\mathbf{k}_{j, \epsilon}(\theta_r)] \cdot \hat{\mathbf{r}}}, \quad (B2b) $$

with $|\nabla\epsilon_{\mathbf{k}_{j, \epsilon}(\theta_r)}|$ norm of the gradient of the energy dispersion $\epsilon_{\mathbf{k}}$, and $\kappa_{\mathbf{k}_{j, \epsilon}(\theta_r)}$ the curvature of the Fermi contour $\epsilon_{\mathbf{k}} = 0$ evaluated at $\mathbf{k}_{j, \epsilon}(\theta_r)$. The points $\mathbf{k}_{j, \epsilon}(\theta_r)$ are the solutions to the saddle-point equations, namely complex momenta satisfying

$$ \epsilon_{\mathbf{k}_{j, \epsilon}(\theta_r)} = i\epsilon\omega, \quad (B3a) $$

$$ \frac{\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}}{|\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}|} = \epsilon \hat{\mathbf{r}}, \quad (B3b) $$

where $j = 1, \dots, N$ denote the solutions from multiple disconnected pockets of the Fermi contour, $\epsilon = \pm$, and $\hat{\mathbf{r}}$ is the unit vector defined by $\theta_r$. For our present purpose of analyzing a small-gap superconductor, we consider a perturbative expansion of the saddle-point equations in $\Delta/\epsilon_F$, noting that the normal metal is recovered as the pairing vanishes. This forces two simultaneous limits: (1) $\Delta/\epsilon_F \to 0$, and (2) $E_S \to \Delta$, which due to the only in-gap state being at $E = E_S$ effectively means that in Eq. (B3a) one has $\omega \to 0$. In Ref. [40] it was shown that to linear order in $\Delta/\epsilon_F$, the angle-dependent quantities in the exponent of the propagator (B1) take the following form:

$$ \text{Re}[\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)] = \tilde{\mathbf{k}}_{j, \epsilon}(\theta_r), \quad (B4a) $$

$$ \xi_{j, \epsilon}(\theta_r) = \frac{|\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}|}{\omega}, \quad (B4b) $$

with $\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)$ the solutions to the saddle-point equations in the $\omega \to 0$ limit,

$$ \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)} = 0, \quad (B5a) $$

$$ \frac{\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}}{|\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}|} = \epsilon \hat{\mathbf{r}}, \quad (B5b) $$

that is, real momenta positioned on the Fermi contour. To zeroth order in $\Delta/\epsilon_F$, it is obvious that the prefactor is real and takes the same form as in the normal-metal propagator, i.e., Eq. (B2a) evaluated at $\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)$,

$$ \Gamma_{j, \epsilon}(\theta_r) = \frac{1}{|\nabla \epsilon_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}| \sqrt{\kappa_{\tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)}}}. \quad (B6) $$

The question arises whether there is any non-trivial correction to next-leading order, i.e., linear in $\Delta/\epsilon_F$. To show that there is no correction, we consider an expansion of the prefactor in terms of a small complex quantity $z$ with $\mathbf{k}_{j, \epsilon}(\theta_r) = \tilde{\mathbf{k}}_{j, \epsilon}(\theta_r) + z$, knowing that $z \to 0$ as $\Delta/\epsilon_F \to 0$. The key point now is that the angular profile of the measured YSR state depends on the amplitude of the prefactor, so that we write $\Gamma_{j, \epsilon}(\theta_r) \equiv |\Gamma| \exp(i\Phi)$ and consider the expansion of the real function evaluated for a complex value of its variable, $|\Gamma(z)| \approx \sqrt{[\Gamma(k) + \frac{1}{2}\Gamma''(k)\text{Re}(z^2)]^2 + [\Gamma'(k)\text{Im}(z) + \frac{1}{2}\Gamma''(k)\text{Im}(z^2)]^2} \approx \Gamma(k) + O(\Delta^2)$, where $k \equiv \tilde{\mathbf{k}}_{j, \epsilon}(\theta_r)$ is the real solution in the metal [Eq. (B5)]. We used that $\Gamma(k)$ is real, and further assumed that the correction $\text{Im}(z)$ is at least linear in $\Delta$, while the $\text{Re}(z)$ is at least quadratic, which is reasonable due to the structure of Eq. (B3a). We note that these properties of $z$ can be confirmed in the particular example of an elliptic Fermi contour worked out in Ref. [40]. For the second step in the expansion we also assumed $\Gamma(k) \neq 0$, which is a condition for the entire saddle-point approach to be valid anyway. Hence the correction of the amplitude of the prefactor is at most second order in $\Delta/\epsilon_F$, so that up to and including first order the prefactor is calculated as for the normal metal [Eq. (B6)]. We note that the prefactor's phase $\Phi$ can in general be linear in $\Delta/\epsilon_F$, thereby contributing an additional $\theta_r$-dependent phase to the bare propagator (and hence the LDOS) whose analysis is beyond the scope of this work.
Continuing now the derivation of the LDOS, to all orders in perturbation theory in the strength of the impurity potential $\hat{V}$, the modification of the LDOS due to a point-like scatterer at $\mathbf{r} = \mathbf{0}$ is given by

$$ \rho(\mathbf{r}, E) = -\frac{1}{\pi} \operatorname{Im} \operatorname{Tr} \left[ \hat{M} \hat{G}_0(\mathbf{r}, \mathbf{0}; E) \hat{T}(E) \hat{G}_0(\mathbf{0}, \mathbf{r}; E) \right], \quad (B7) $$

where $\hat{T}(E) = \hat{V} [1 - \hat{G}_0(\mathbf{0}, \mathbf{0}; E) \hat{V}]^{-1}$ is the T-matrix and $\hat{M} = (\tau_0 + \tau_z)/2$ projects out the electron-electron component. Here and in the following $\tau_{0,x,y,z}$ denote the Pauli matrices in Nambu space. Note that the counter-propagator satisfies $\hat{G}_0(\mathbf{0}, \mathbf{r}; E) = \hat{G}_0^*(\mathbf{r}, \mathbf{0}; E)$. Since the anisotropy of the LDOS is encoded in the bare propagator (B1), we simplify by considering a magnetic scatterer only, without non-magnetic scattering amplitude, i.e. $\hat{V} = -J\tau_0$. Further, we neglect the anisotropy of the band structure when evaluating $\hat{G}_0(\mathbf{0}, \mathbf{0}; E)$ in the T-matrix,

$$ \hat{G}_0(\mathbf{0}, \mathbf{0}; E) \equiv \frac{1}{V} \sum_{\mathbf{k}} \hat{G}_0(\mathbf{k}, E) \sim -\frac{\pi \nu_0}{\sqrt{\Delta^2 - E^2}} \begin{pmatrix} E & \Delta \\ \Delta & E \end{pmatrix}, \quad (B8) $$

where $\nu_0$ is the density of states at the Fermi level. Under these assumptions and in the small-gap limit we have

$$ \rho(\mathbf{r}) \sim -\frac{1}{\pi} \operatorname{Im} \frac{1}{r} \sum_{j,j',\epsilon,\epsilon'} \Gamma^{\epsilon,\epsilon'}_{j,j'}(\theta_{\mathbf{r}}) \exp \left\{ -\frac{r}{\xi^{\epsilon,\epsilon'}_{j,j'}(\theta_{\mathbf{r}})} + i \left[ \mathbf{r} \cdot (\tilde{\mathbf{k}}_{j,\epsilon}(\theta_{\mathbf{r}}) - \tilde{\mathbf{k}}_{j',\epsilon'}(\theta_{\mathbf{r}})) + (\epsilon - \epsilon') \frac{\pi}{4} \right] \right\} \frac{N_{\epsilon,\epsilon'}}{D}, \quad (B9) $$

where

$$ \Gamma^{\epsilon, \epsilon'}_{j,j'}(\theta_{\mathbf{r}}) = \Gamma_{j,\epsilon}(\theta_{\mathbf{r}}) \Gamma_{j',\epsilon'}(\theta_{\mathbf{r}}), \quad (B10a) $$

$$ \xi^{\epsilon, \epsilon'}_{j,j'}(\theta_{\mathbf{r}}) = \left( \frac{1}{\xi_{j,\epsilon}(\theta_{\mathbf{r}})} + \frac{1}{\xi_{j',\epsilon'}(\theta_{\mathbf{r}})} \right)^{-1}, \quad (B10b) $$

$$ N_{\epsilon, \epsilon'} = t_0(\Delta^2 + E^2) + 2E\Delta t_x + t_0\omega^2\epsilon\epsilon' + i\omega(Et_0 + \Delta t_x)(\epsilon - \epsilon'), \quad (B10c) $$

$$ D = t_0^2 - t_x^2, \quad (B10d) $$

with $t_0 = 1 - \alpha \frac{E}{\omega}$, $t_x = \alpha \frac{\Delta}{\omega}$ and $\alpha = \nu_0 \pi J$, and $\Gamma_{j,\epsilon}(\theta_{\mathbf{r}})$ and $\xi_{j,\epsilon}(\theta_{\mathbf{r}})$ as defined in Eqs. (B6) and (B4b), respectively. The energy of the in-gap bound states are given by the poles of the T-matrix, i.e. $D = 0$, and takes the well known form,

$$ E_S^{\pm} = \pm \Delta \frac{1 - \alpha^2}{1 + \alpha^2}. \quad (B11) $$

To obtain the LDOS of the YSR state [54], we evaluate $N_{\epsilon, \epsilon'}$ at $E = E_S^+$ and expand $D$ to first order around $E = E_S^+$ to obtain up to a constant prefactor,

$$ \rho(\mathbf{r}) \sim -\frac{1}{\pi} \operatorname{Im} \frac{1}{r} \sum_{j,j',\epsilon,\epsilon'} \Gamma^{\epsilon,\epsilon'}_{j,j'}(\theta_{\mathbf{r}}) \exp \left\{ -\frac{r}{\xi^{\epsilon,\epsilon'}_{j,j'}(\theta_{\mathbf{r}})} + i\mathbf{r} \cdot (\tilde{\mathbf{k}}_{j,\epsilon}(\theta_{\mathbf{r}}) - \tilde{\mathbf{k}}_{j',\epsilon'}(\theta_{\mathbf{r}})) \right\} \frac{X_{\epsilon,\epsilon'}}{E - E_S}, \quad (B12) $$

where $X_{\epsilon, \epsilon'} = (1 + \alpha^2)e^{i\varphi_{\epsilon, \epsilon'}}$ with

$$ \varphi_{\epsilon, \epsilon'} = \begin{cases} 0 & \text{if } \epsilon = \epsilon', \\ \epsilon \left[ \frac{\pi}{2} + \arg(1 - \alpha^2 + i2\alpha) \right] & \text{if } \epsilon \neq \epsilon'. \end{cases} \quad (B13) $$

We note that including a nonmagnetic component to the scattering potential yields a different expression for $X_{\epsilon, \epsilon'}$, but crucially with the same dependence on $\epsilon, \epsilon'$.

The expression in Eq. (B12) can be written in a more compact form by noticing that, owing to the evenness of the energy dispersion $\epsilon_k$, for each point $\tilde{k}_{j,+}(\theta_{\mathbf{r}})$, there exists a point $\tilde{k}_{l,-}(\theta_{\mathbf{r}})$ such that $\tilde{k}_{l,-}(\theta_{\mathbf{r}}) = -\tilde{k}_{j,+}(\theta_{\mathbf{r}})$. If the j pocket is centered around the $\Gamma$ point, then $j = l$, but this is not the case in general. The prefactor and characteristic decay length satisfy that $\Gamma_{j,+}(\theta_{\mathbf{r}}) = \Gamma_{l,-}(\theta_{\mathbf{r}}) \equiv \Gamma_j(\theta_{\mathbf{r}})$ and $\xi_{j,+}(\theta_{\mathbf{r}}) = \xi_{l,-}(\theta_{\mathbf{r}}) \equiv \xi_j(\theta_{\mathbf{r}})$, therefore, by relabeling the
critical points as {$\epsilon \tilde{k}_j(\theta_r)}_{j=1...N}$, it becomes obvious that the only contribution to the LDOS comes from the poles of the T-matrix,

$$ \rho(\mathbf{r}) \sim \frac{1}{r} \sum_{j,j'} \Gamma_{j,j'}(\theta_r) e^{-r/\xi_{j,j'}(\theta_r)} \sum_{\epsilon, \epsilon'} \exp \left\{ i \left[ \mathbf{r} \cdot (\epsilon \tilde{k}_j(\theta_r) - \epsilon' \tilde{k}_{j'}(\theta_r)) + \varphi_{\epsilon, \epsilon'} \right] \right\} \delta(E - E_S), \quad (B14) $$

where Eq. (2) in the main text is recovered by introducing

$$ f[\mathbf{k}_{j,j'}^{\pm}(\theta_r) \cdot \mathbf{r}] \equiv \sum_{\epsilon, \epsilon'} \exp \left\{ i \left[ \mathbf{r} \cdot (\epsilon \tilde{\mathbf{k}}_j(\theta_r) - \epsilon' \tilde{\mathbf{k}}_{j'}(\theta_r)) + \varphi_{\epsilon, \epsilon'} \right] \right\}, \quad (B15) $$

that is a real-valued function. With the new notation the $\epsilon = -$ solutions to Eq. (C4a) are constructed from the $\epsilon = +$ solutions; therefore, Eq. (4) in the main text considers the parallel solutions only. Note, however, that the antiparallel solutions are encoded in $f[\mathbf{k}_{j,j'}^{\pm}(\theta_r) \cdot \mathbf{r}]$.

## Appendix C: Analysis of the tight-binding energy dispersion

The starting point to compute the model curve for the prefactor (Fig. 2 in the main text) is a tight-binding energy dispersion that describes the two lowest-lying Nb 4d bands of 2H-NbSe$_2$,

$$ \varepsilon_{\mathbf{k}} = \mu + \sum_{i=1}^{5} t_i g_i(\mathbf{k}), \quad (C1) $$

where $t_i$ is the hopping amplitude to the i-th nearest neighbor and

$$ g_1(\mathbf{k}) = 2 \cos \zeta \cos \eta + \cos 2\zeta, \quad (C2a) $$

$$ g_2(\mathbf{k}) = 2 \cos 3\zeta \cos \eta + \cos 2\eta, \quad (C2b) $$

$$ g_3(\mathbf{k}) = 2 \cos 2\zeta \cos 2\eta + \cos 4\zeta, \quad (C2c) $$

$$ g_4(\mathbf{k}) = \cos \zeta \cos 3\eta + \cos 5\zeta \cos \eta + \cos 4\zeta \cos 2\eta, \quad (C2d) $$

$$ g_5(\mathbf{k}) = 2 \cos 3\eta \cos 3\zeta + \cos 6\zeta, \quad (C2e) $$

with $\zeta = \frac{1}{2}k_x a_c$ and $\zeta = \frac{\sqrt{3}}{2}k_y a_c$, and $a_c$ is the lattice constant. The hopping amplitudes are presented in Table I. We follow Refs. [29, 33] and assume that the adatom primarily couples to one of the bands only. Specifically, we take Band 2, which we treat as a spin-degenerate energy dispersion.

<table><caption>TABLE I. Fitted tight-binding parameters for the two lowest-lying Nb 4d bands of 2H-NbSe$_2$. All parameters in meV, for $a_c = 3.444$ Å. Extracted from Ref. [41]</caption><thead><tr><th></th><th>μ</th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>Band 1</td><td>10.9</td><td>86.8</td><td>139.9</td><td>29.6</td><td>3.5</td><td>3.3</td></tr><tr><td>Band 2</td><td>203.0</td><td>46.0</td><td>257.5</td><td>4.4</td><td>-15.0</td><td>6.0</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Examples of critical points for various values of $\theta_r$. The blue lines indicate the Fermi contour for the chosen 2H-NbSe$_2$ model, i.e. Eq. (C1) with parameters of Band 2. The black arrow represents $\hat{r}$, i.e., a vector in real space from the impurity location to an arbitrary position where we evaluate the LDOS. The orange dots on the Fermi contour indicate the corresponding critical points $\epsilon \tilde{k}_j(\theta_r)$, i.e., the solutions to Eqs. (C4a and C4b) ($\epsilon = +$) and their spatially inverted counterparts ($\epsilon = -$). The orange arrows indicate the normalized gradient of the energy dispersion at that point, parallel ($\epsilon = +$) and antiparallel ($\epsilon = -$) to $\hat{r}$. Panels (a), (b), and (c) correspond to $\theta_r = \pi/6, \pi/3, 3\pi/2$, respectively.</figcaption></figure>

The Fermi velocity and the curvature of the Fermi contour that control the anisotropy of the LDOS are related to the energy dispersion as follows,

$$|\nabla \epsilon_k| = \sqrt{(\partial_{k_x} \epsilon_k)^2 + (\partial_{k_y} \epsilon_k)^2}, \qquad (C3a)$$

$$\kappa_k = \frac{(\partial_{k_x} \epsilon_k)^2 \partial_{k_y}^2 \epsilon_k + (\partial_{k_y} \epsilon_k)^2 \partial_{k_x}^2 \epsilon_k - 2 \partial_{k_x} \epsilon_k \partial_{k_y} \epsilon_k \partial_{k_x, k_y}^2 \epsilon_k}{|\nabla \epsilon_k|^3}. \qquad (C3b)$$

To compute the model curve for the prefactor we first obtain the critical points in the FBZ by numerically solving

$$\epsilon_{\tilde{k}_j}(\theta_r) = 0, \qquad (C4a)$$

$$\frac{\nabla \epsilon_{\tilde{k}_j}(\theta_r)}{|\nabla \epsilon_{\tilde{k}_j}(\theta_r)|} = \hat{r}, \qquad (C4b)$$

for a discrete set of $\theta_r$ values. As the Fermi contour of the energy dispersion in Eq. (C1) has three non-equivalent pockets, we obtain six critical points $\{\epsilon \tilde{k}_j(\theta_r)\}_{j=1,2,3}$ for a given direction in real space determined by $\theta_r$ (see Fig. 6).

This yields a set of six distinct $\Gamma_{j,j'}(\theta_r)$ curves that are essentially in phase [see Fig. 7 (a)]. Crucially, the prefactor curves are in phase with respect to the characteristic decay length curves as well [Fig. 7 (b)], which justifies the simplification of taking the dominant prefactor $\Gamma_{\max}(\theta_r) = \max\{\Gamma_{j,j'}(\theta_r)\}$ in the model curve.

## Appendix D: Effect of Charge Density Wave

The presence of charge density wave (CDW) in principle affects two aspects of the YSR LDOS, namely, there is a deformation of the underlying normal state Fermi surface, and there is a modification of the coupling of the impurity to the substrate.

Concerning the effect on the impurity first, Liebhaber and coworkers established in [33] that the relative position of the magnetic adatom with respect to the CDW ordering can significantly modify the properties (e.g., reduce the spatial symmetry) of the Yu-Shiba-Rusinov wave function. Impurities studied in this work are subsurface; therefore, it is challenging to determine their exact location with respect to the CDW modulation. However, since the observed LDOS is sixfold symmetric at long-range distances, we can safely assume that the corresponding magnetic impurities are situated on a maximum or a minimum of the CDW, or, alternatively, that the local CDW effect is negligible for buried defects. Hence, we will neglect any influence on the impurity.

Concerning the effect on the band structure, the Fermi surface deformed by the CDW may change the prefactor that is essential to the QFE theory. The conclusion revealed in the rest of this Appendix is that the CDW introduces
<figure><img src="image_8.png" /><figcaption>FIG. 7. (a) Angular dependence of the LDOS prefactor. (b) Angular dependence of LDOS exponential decay length with $\omega = 1$ meV. Each curve denotes an interference term ($j, j'$) with $j, j' = 1, 2, 3$ for pockets centered at $\Gamma = (0, 0)$, $K = (0, 4\pi/3a)$ and $K' = -K$. Certain curves are three-fold periodic only reflecting the reduced symmetry of the $K, K'$ pockets [e.g. (2, 2)], yet, the six-fold symmetry of the band structure is recovered when taking into account all the interference terms. The curve presented in Fig. 2 in the main text corresponds to curve (1, 1) in panel (a). Compare with Fig. 6 to observe that both the prefactor and the exponential decay length are maximal when $\dot{r}$ is perpendicular to the flattest sections of the Fermi contour.</figcaption></figure>

qualitative changes which may manifest strongly in other compounds, but the magnitude of these effects is negligible in 2H-NbSe$_2$.

We borrow the model for CDW from [33] where they consider an effective mean-field description of the electron-phonon coupling and select the phonon modes corresponding to the lowest Fourier harmonics only. Under these assumptions, the quasi-(3 × 3)-commensurate CDW ordering is described by the Hamiltonian

$$H_V = \sum_{\mathbf{r}} V(\mathbf{r}) c_{\mathbf{r},\sigma}^{\dagger} c_{\mathbf{r},\sigma}, \qquad (\text{D1})$$

where $V(\mathbf{r}) = V_0 \sum_{i=1}^3 \cos(\mathbf{Q}_i \cdot \mathbf{r} + \phi_i)$ is a periodic potential with amplitude $V_0 = -30.0$ meV,

$$\mathbf{Q}_1 = q(1 - \delta) \left(\frac{\sqrt{3}}{2}, -\frac{1}{2}\right), \qquad (\text{D2a})$$

$$\mathbf{Q}_2 = q(1 - \delta)(0, 1), \qquad (\text{D2b})$$

$$\mathbf{Q}_3 = -\mathbf{Q}_1 - \mathbf{Q}_2, \qquad (\text{D2c})$$

the CDW wavevectors, and the phases chosen to be $\phi_1 = \phi_2 = 0$ and $\phi_3 = -\pi/3$. In Eqs. (D2), $q = \frac{4\pi}{3\sqrt{3}a_c}$ and $\delta \ll 1$ accounts for the small incommensurability of the CDW with the Nb lattice that we shall neglect. In practice, we apply the method described in App. C to each of the six eigenbands at the Fermi energy of the following 9 × 9 Hamiltonian in the basis $(c_k, c_{k+\mathbf{Q}_1}, c_{k-\mathbf{Q}_1}, c_{k+\mathbf{Q}_2}, c_{k-\mathbf{Q}_2}, c_{k+\mathbf{Q}_3}, c_{k-\mathbf{Q}_3}, c_{k+\mathbf{Q}_4}, c_{k-\mathbf{Q}_4})^T$,

$$H_{\text{CDW}} = \begin{pmatrix}
\varepsilon(k) & V_1 & V_1^* & V_2 & V_2^* & V_3 & V_3^* & \cdot & \cdot \\
\varepsilon(k+\mathbf{Q}_1) & V_1 & \cdot & V_3 & \cdot & V_2 & V_2^* & V_3^* & V_3 \\
\varepsilon(k-\mathbf{Q}_1) & V_3^* & \cdot & V_2^* & \cdot & V_3 & V_2 & V_2^* & V_2 \\
\varepsilon(k+\mathbf{Q}_2) & V_2 & \cdot & V_1 & V_1^* & V_3^* & V_1^* & V_3 & V_3 \\
\varepsilon(k-\mathbf{Q}_2) & V_1^* & \cdot & V_1 & V_1^* & V_3 & V_2^* & V_1^* & V_3 \\
\varepsilon(k+\mathbf{Q}_3) & V_3 & V_3^* & V_2 & V_2^* & V_3 & V_1^* & V_2^* & V_1 \\
\varepsilon(k-\mathbf{Q}_3) & V_2 & V_2^* & V_1 & V_1^* & V_2 & V_1^* & V_1 & V_1 \\
\varepsilon(k+\mathbf{Q}_4) & V_3 & V_3^* & V_2 & V_2^* & V_3 & V_1^* & V_2^* & V_1^* \\
\varepsilon(k-\mathbf{Q}_4) & V_2 & V_2^* & V_1 & V_1^* & V_2 & V_1^* & V_1 & V_1
\end{pmatrix} \quad (\text{D3})$$

where $\mathbf{Q}_4 = \mathbf{Q}_1 - \mathbf{Q}_2$, $\delta = 0$ in Eqs. (D2), $\varepsilon(\mathbf{k})$ is the energy dispersion defined in Eq. (C1), $V_i = \frac{V_0}{18} e^{i\phi_i}$, and $\cdot = 0$, while “*” and “c.c.” denote complex conjugation.
<figure><img src="image_9.png" /><figcaption>FIG. 8. Effect of CDW on the focusing effect in 2H-NbSe$_2$. (a) Fermi contours with a finite CDW potential. Each color represents a band of the Hamiltonian in Eq. (D3). The circular markers indicate the inflection points (IPs). The black hexagon delimits the reduced FBZ, which is nine times smaller than the FBZ in Fig. 6 owing to band folding. (b) Angle of the Fermi velocity $\theta_{v_F}$ as a function of the polar angle in momentum space $\theta_k$ for each of the Fermi contours hosting IPs, which are indicated by the circular markers. The small arrows mark the corresponding value $\theta_{v_F}$ of the IPs (including those outside the shown $\theta_k$ range). The horizontal dashed black lines indicate the high-symmetry points on which the maxima of the prefactor curve in the absence of CDW potential lie. The horizontal axes span the range delimited by the dashed black radii in panel (a). (c) Angular dependence of the LDOS prefactor. The blue (red) curve corresponds to the calculation with finite (vanishing) CDW potential $V_0$ (the red curve is equivalent to the model curve in Fig. 2, with $a_0 = 0$ and $a = 1$). The black shaded rectangles indicate the $\theta_r$ regions for which the corresponding critical momenta fall on the vicinity of an IP, and, therefore, the QFE theory is not applicable. In addition to the fundamental break down of the theory, the calculation of Eq. (3) in these regions is numerically unstable owing to the vanishing second derivatives, hence the prefactor curve appears to be very noisy.</figcaption></figure>

If the amplitude of the CDW potential $V_0$ vanishes, the Hamiltonian in Eq. (D3) is diagonal and follows from a redefinition of the unit cell. The resulting band-folded Fermi contours naturally give rise to the same prefactor as the unfolded band structure. A finite (but small) CDW potential opens small gaps at the intersection of the band-folded Fermi contours [51]. Due to this hybridization, a pair of inflection points (IPs) with vanishing curvature $\kappa_k$ appears in the vicinity of the contour-touching points. The pair of IPs emerges in one of the contours involved in the crossing only, namely, the contour which goes from concave to convex and to concave again [see, for instance, the intersection of the $\Gamma$ pocket with itself (blue and orange contours at $\theta_k = n\pi/3$, $n = 1...6$) in Fig. 8 (a)]. The QFE theory is not applicable on these singular points since a vanishing curvature requires including next-to-leading order corrections in the saddle point approximation [40], but, in general, one can expect that they create an enhancement of the LDOS prefactor. It is then natural to ask if the emerging inflection points with vanishing curvature yield new peaks in the
prefactor curve that would appear as additional petals in the YSR LDOS.

Answering that question boils down to determining the set of real-space polar angles $\theta_r$ for which the corresponding critical points $\epsilon \tilde{k}_j(\theta_r)$ lay on the IPs of the Fermi contours. If the former differ from the $\theta_r$ with maximum prefactor at zero CDW potential, new petals will appear in the YSR LDOS. We recall that the Fermi velocity at the critical momenta is parallel to the polar vector along the observation direction in real space $\hat{r}$ [see Eq. (4)], therefore the Fermi velocity angle allows us to relate $\theta_r$ to $\theta_k$ for each band [Fig. 8 (b)]. The curves $\theta_{\hat{v}_F}(\theta_k)$ implicitly contain all the information about the curvature of the Fermi contour since the Fermi velocity vector is perpendicular to the latter by definition. In particular, a positive (negative) slope indicates a convex (concave) segment in the Fermi contour, and, consequently, the extrema of $\theta_{\hat{v}_F}(\theta_k)$ correspond to the IPs where the curvature vanishes. The sharp positive variations in the $\theta_{\hat{v}_F}(\theta_k)$ curve describe the highly curved convex segments of the Fermi contours that result from the CDW hybridization and, consequently, they always appear flanked by two IPs. Further, they coincide with sharp negative steps in the $\theta_{\hat{v}_F}(\theta_k)$ curve of the other contour involved in the contour-touching point (recall, e.g., the blue and orange contour touching at $\theta_k = n\pi/3, n = 1...6$). The sharp negative variations of $\theta_{\hat{v}_F}(\theta_k)$ also describe highly curved segments resulting from the CDW hybridization, but in this case, they do not involve any change in the curvature sign, and therefore, there are not any IPs associated with them.

The analysis in Fig. 8(b) demonstrates that, for a choice of $V_0$ relevant to 2H-NbSe$_2$, the $\theta_r$ corresponding to the IPs are distributed over a finite region around the high-symmetry points $\theta_r = (n + 1/2)\pi/3, n = 1...6$ where the theory of the QFE to lowest order of approximation fails. The full calculation of the LDOS prefactor in Fig. 8(c) shows that the prefactor curves with finite and vanishing CDW potential perfectly match within the region of validity. Crucially, the breakdown regions where an enhancement of the LDOS prefactor may be expected are aligned with the maxima of the prefactor curve calculated for $V_0 = 0$. Thus, we can conclude that the band structure modifications induced by the CDW ordering on 2H-NbSe$_2$ do not modify the spatial orientation of the Yu-Shiba-Rusinov states as it was argued above based on the ARPES observations [41].

Finally, it is worth noting that the emergence of inflection points illustrated here is not specific to 2H-NbSe$_2$, but applies to any band structure where CDW can be treated as a commensurate perturbation. In general, a larger value of the CDW potential amplitude $V_0$ will significantly modify the Fermi contours and the Fermi velocity direction at the emergent inflection points may substantially differ from its direction at the lowest-curvature points for $V_0 = 0$, thereby yielding additional focused petals in the YSR LDOS.

## Appendix E: Derivation of the Fourier transform of the LDOS

We derive here the Fourier transform of the LDOS at the YSR-state energy:

$$ 
\begin{aligned}
\tilde{\rho}(\mathbf{q}) &\equiv \int d\mathbf{r} e^{-i\mathbf{q}\cdot\mathbf{r}} \rho(\mathbf{r}) \\
&= \sum_{j,j'} \int_0^\infty dr \int_0^{2\pi} d\theta_r \Gamma_{j,j'}(\theta_r) \sum_{\epsilon, \epsilon'} \exp \left\{ r \left[ i \mathbf{p}_{\epsilon, \epsilon'}(\mathbf{q}, \theta_r) \cdot \hat{\mathbf{r}} - \frac{1}{\xi_{j,j'}(\theta_r)} \right] \right\} \quad (E1) \\
&= \sum_{j,j'} \int_0^{2\pi} d\theta_r \Gamma_{j,j'}(\theta_r) \sum_{\epsilon, \epsilon'} e^{i\varphi_{\epsilon, \epsilon'}} \frac{1}{i \mathbf{p}_{\epsilon, \epsilon'}(\mathbf{q}, \theta_r) \cdot \hat{\mathbf{r}} - 1/\xi_{j,j'}(\theta_r)},
\end{aligned}
$$

where $\mathbf{p}_{\epsilon, \epsilon'}(\mathbf{q}, \theta_r) = \epsilon \tilde{k}_j(\theta_r) - \epsilon' \tilde{k}_{j'}(\theta_r) - \mathbf{q}$, and Eq. (9) in the main text is recovered by setting

$$ \tilde{F}_{j,j'}(\mathbf{q}, \theta_r) \equiv \sum_{\epsilon, \epsilon'} e^{i\varphi_{\epsilon, \epsilon'}} \frac{1}{i \mathbf{p}_{\epsilon, \epsilon'}(\mathbf{q}, \theta_r) \cdot \hat{\mathbf{r}} - 1/\xi_{j,j'}(\theta_r)}. \quad (E2) $$

The maps in Figs. 3(b) and 3(c) are obtained by replacing the angular integral by a sum over a discrete set of $\{\theta_r\}$ values, for each of which the term $\sum_{j,j'} \Gamma_{j,j'}(\theta_r) \tilde{F}_{j,j'}(\mathbf{q}, \theta_r)$ is evaluated. Importantly, to ensure that $\tilde{\rho}(\mathbf{q})$ has the periodicity of the reciprocal lattice, the pairs of critical points $\epsilon \tilde{k}_j(\theta_r) - \epsilon' \tilde{k}_{j'}(\theta_r)$ have to be mapped onto the FBZ. Further, as it was shown in Appendix B, the $\varphi_{\epsilon, \epsilon'}$ phase depends on the nature of the scatterer. For simplicity, we set $\alpha = 1$ in Eq. (B13), which assumes the strong-coupling regime, i.e., $E_S = 0$.

<footnote>[1] A. I. Rusinov, Superconductivity near a Paramagnetic Impurity, Sov. J. Exp. Theor. Phys. Lett 9, 85 (1969).</footnote>

<footnote>[2] H. Shiba, Classical Spins in Superconductors, Progr.</footnote>
Theor. Phys. **40**, 435 (1968).

[3] L. Yu, BOUND STATE IN SUPERCONDUCTORS WITH PARAMAGNETIC IMPURITIES, Acta Phys. Sin. **21**, 75 (1965).

[4] B. W. Heinrich, J. I. Pascual, and K. J. Franke, Single magnetic adsorbates on s-wave superconductors, Progress in Surface Science **93**, 1 (2018).

[5] V. Kaladzhyan, C. Bena, and P. Simon, Characterizing p-wave superconductivity using the spin structure of shiba states, Phys. Rev. B **93**, 214514 (2016).

[6] V. Perrin, F. L. N. Santos, G. C. Ménard, C. Brun, T. Cren, M. Civelli, and P. Simon, Unveiling odd-frequency pairing around a magnetic impurity in a superconductor, Phys. Rev. Lett. **125**, 117003 (2020).

[7] P. O. Sukhachov, F. von Oppen, and L. I. Glazman, Tunneling spectra of impurity states in unconventional superconductors, Phys. Rev. B **108**, 024505 (2023).

[8] D.-J. Choi, C. Rubio-Verdú, J. De Bruijckere, M. M. Ugeda, N. Lorente, and J. I. Pascual, Mapping the orbital structure of impurity bound states in a superconductor, Nat. Commun. **8**, 15175 (2017).

[9] D.-J. Choi, C. G. Fernández, E. Herrera, C. Rubio-Verdú, M. M. Ugeda, I. Guillamón, H. Suderow, J. I. Pascual, and N. Lorente, Influence of Magnetic Ordering between Cr Adatoms on the Yu-Shiba-Rusinov States of the β-Bi₂Pd Superconductor, Phys. Rev. Lett. **120**, 167001 (2018).

[10] H. Huang, C. Padurariu, J. Senkpiel, R. Drost, A. L. Yeyati, J. C. Cuevas, B. Kubala, J. Ankerhold, K. Kern, and C. R. Ast, Tunnelling dynamics between superconducting bound states at the atomic limit, Nat. Phys. **16**, 1227 (2020).

[11] H. Huang, J. Senkpiel, C. Padurariu, R. Drost, A. Vilkas, R. L. Klees, A. L. Yeyati, J. C. Cuevas, B. Kubala, J. Ankerhold, K. Kern, and C. R. Ast, Spin-dependent tunneling between individual superconducting bound states, Phys. Rev. Research **3**, L032008 (2021).

[12] L. Schneider, P. Beck, J. Wiebe, and R. Wiesendanger, Atomic-scale spin-polarization maps using functionalized superconducting probes, Science Advances **7**, eabd7302 (2021).

[13] K. Vaxevani, J. Li, S. Trivini, J. Ortuzar, D. Longo, D. Wang, and J. I. Pascual, Extending the spin excitation lifetime of a magnetic molecule on a proximitized superconductor, Nano Lett. **22**, 6075-6082 (2022).

[14] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, and A. Yazdani, Proposal for realizing Majorana fermions in chains of magnetic atoms on a superconductor, Phys. Rev. B **88**, 020407 (2013).

[15] F. Pientka, L. I. Glazman, and F. von Oppen, Topological superconducting phase in helical Shiba chains, Phys. Rev. B **88**, 155420 (2013).

[16] B. Braunecker and P. Simon, Interplay between classical magnetic moments and superconductivity in quantum one-dimensional conductors: Toward a self-sustained topological Majorana phase, Phys. Rev. Lett. **111**, 147202 (2013).

[17] S. Nakosai, Y. Tanaka, and N. Nagaosa, Two-dimensional p-wave superconducting states with magnetic moments on a conventional s-wave superconductor, Phys. Rev. B **88**, 180503 (2013).

[18] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of Majorana fermions in ferromagnetic

atomic chains on a superconductor, Science **346**, 602-607 (2014).

[19] J. Röntynen and T. Ojanen, Topological superconductivity and high chern numbers in 2D ferromagnetic Shiba lattices, Phys. Rev. Lett. **114**, 236803 (2015).

[20] J. Li, T. Neupert, Z. Wang, A. MacDonald, A. Yazdani, and B. A. Bernevig, Two-dimensional chiral topological superconductivity in Shiba lattices, Nat. Commun. **7**, 12297 (2016).

[21] S. Jeon, Y. Xie, J. Li, Z. Wang, B. A. Bernevig, and A. Yazdani, Distinguishing a Majorana zero mode using spin-resolved measurements, Science **358**, 772-776 (2017).

[22] M. Ruby, B. W. Heinrich, Y. Peng, F. von Oppen, and K. J. Franke, Exploring a proximity-coupled Co chain on Pb (110) as a possible Majorana platform, Nano Lett. **17**, 4473 (2017).

[23] H. Kim, A. Palacio-Morales, T. Posske, L. Rózsa, K. Palotás, L. Szunyogh, M. Thorwart, and R. Wiesendanger, Toward tailoring Majorana bound states in artificially constructed magnetic atom chains on elemental superconductors, Science Advances **4**, eaar5251 (2018).

[24] A. Palacio-Morales, E. Mascot, S. Cocklin, H. Kim, S. Rachel, D. K. Morr, and R. Wiesendanger, Atomic-scale interface engineering of Majorana edge modes in a 2D magnet-superconductor hybrid system, Science Advances **5**, eaav6600 (2019).

[25] E. Liebhaber, L. M. Rütten, G. Reecht, J. F. Steiner, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Quantum spins and hybridization in artificially-constructed chains of magnetic adatoms on a superconductor, Nat. Commun. **13**, 2160 (2022).

[26] L. Farinacci, G. Reecht, F. von Oppen, and K. J. Franke, Yu-Shiba-Rusinov bands in a self-assembled kagome lattice of magnetic molecules, Nature Communications **15**, 6474 (2024).

[27] M. O. Soldini, F. Küster, G. Wagner, S. Das, A. Aldarawsheh, R. Thomale, S. Lounis, S. S. P. Parkin, P. Sessi, and T. Neupert, Two-dimensional shiba lattices as a possible platform for crystalline topological superconductivity, Nat. Phys. **19**, 1848-1854 (2023).

[28] A. Yazdani, F. Von Oppen, B. I. Halperin, and A. Yacoby, Hunting for Majoranas, Science **380**, eade0850 (2023).

[29] G. C. Ménard, S. Guissart, C. Brun, S. Pons, V. S. Stolyarov, F. Debontridder, M. V. Leclerc, E. Janod, L. Cario, D. Roditchev, P. Simon, and T. Cren, Coherent long-range magnetic bound states in a superconductor, Nat. Phys. **11**, 1013 (2015).

[30] H. Kim, L. Rózsa, D. Schreyer, E. Simon, and R. Wiesendanger, Long-range focusing of magnetic bound states in superconducting Lanthanum, Nat. Commun. **11**, 4573 (2020).

[31] U. Thupakula, V. Perrin, A. Palacio-Morales, L. Cario, M. Aprili, P. Simon, and F. Massee, Coherent and Incoherent Tunneling into Yu-Shiba-Rusinov States Revealed by Atomic Scale Shot-Noise Spectroscopy, Phys. Rev. Lett. **128**, 247001 (2022).

[32] S. Kezilebieke, M. Dvorak, T. Ojanen, and P. Liljeroth, Coupled yu-shiba-rusinov states in molecular dimers on nbse2, Nano Lett. **18**, 2311-2315 (2018).

[33] E. Liebhaber, S. Acero González, R. Baba, G. Reecht, B. W. Heinrich, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Yu-Shiba-Rusinov States in the Charge-Density Modulated Superconductor NbSe₂, Nano
Lett. 20, 339 (2020).

[34] A. Ptok, S. Glodzik, and T. Domański, Yu-Shiba-Rusinov states of impurities in a triangular lattice of NbSe$_2$ with spin-orbit coupling, Phys. Rev. B 96, 184425 (2017).

[35] D. Sticlet and C. Morari, Topological superconductivity from magnetic impurities on monolayer NbSe$_2$, Phys. Rev. B 100, 075420 (2019).

[36] L. M. Rütten, H. Schmid, E. Liebhaber, G. Franceschi, A. Yazdani, G. Reecht, K. Rossnagel, F. von Oppen, and K. J. Franke, Wave-function engineering on superconducting substrates: Chiral Yu-Shiba-Rusinov molecules, arXiv:2404.16513 10.48550/arXiv:2404.16513 (2024).

[37] J. Ortuzar, S. Trivini, M. Alvarado, M. Rouco, J. Zaldívar, A. L. Yeyati, J. I. Pascual, and F. S. Bergeret, Yu-Shiba-Rusinov states in 2D superconductors with arbitrary Fermi contours, Phys. Rev. B 105, 245403 (2022).

[38] A. Weismann, M. Wenderoth, S. Lounis, P. Zahn, N. Quaas, R. G. Ulbrich, P. H. Dederichs, and S. Blügel, Seeing the Fermi Surface in Real Space by Nanoscale Electron Focusing, Science 323, 1190 (2009).

[39] S. Lounis, P. Zahn, A. Weismann, M. Wenderoth, R. G. Ulbrich, I. Mertig, P. H. Dederichs, and S. Blügel, Theory of real space imaging of Fermi surface parts, Phys. Rev. B 83, 035427 (2011).

[40] M. Uldemolins, A. Mesaros, and P. Simon, Quasiparticle focusing of bound states in two-dimensional s-wave superconductors, Phys. Rev. B 105, 144503 (2022).

[41] D. J. Rahn, S. Hellmann, M. Kalläne, C. Sohrt, T. K. Kim, L. Kipp, and K. Rossnagel, Gaps and kinks in the electronic structure of the superconductor 2H-NbSe$_2$ from angle-resolved photoemission at 1 K, Phys. Rev. B 85, 224532 (2012).

[42] J. G. Rodrigo and S. Vieira, STM study of multiband superconductivity in NbSe$_2$ using a superconducting tip, Physica C 404, 306 (2004).

[43] M. Kuzmanović, T. Dvir, D. LeBoeuf, S. Ilić, M. Haim, D. Möckli, S. Kramer, M. Khodas, M. Houzet, J. S. Meyer, M. Aprili, H. Steinberg, and C. H. L. Quay, Tunneling spectroscopy of few-monolayer NbSe$_2$ in high magnetic fields: Triplet superconductivity and Ising protection, Phys. Rev. B 106, 184514 (2022).

[44] D. Sanchez, A. Junod, J. Muller, H. Berger, and F. Lévy, Specific heat of 2H-NbSe$_2$ in high magnetic fields, Physica B 204, 167 (1995).

[45] A. Nader and P. Monceau, Critical field of 2H-NbSe$_2$ down to 50 mK, SpringerPlus 3, 16 (2014).

[46] A. Mourachkine, Determination of the coherence length and the cooper-pair size in unconventional superconductors by tunneling spectroscopy, Journal of Superconductivity 17, 711 (2004).

[47] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Ising pairing in superconducting NbSe$_2$ atomic layers, Nat. Phys. 12, 139 (2016).

[48] M. Ugeda, A. Bradley, Y. Zhang, S. Onishi, Y. Chen, W. Ruan, C. Ojeda-Aristizabal, H. Ryu, M. Edmonds, H.-Z. Tsai, A. Riss, S.-K. Mo, D. Lee, A. Zettl, Z. Husain, Z.-X. Shen, and M. Crommie, Characterization of collective ground states in single-layer 2H-NbSe$_2$, Nat. Phys. 12 (2016).

[49] P. Dreher, W. Wan, A. Chikina, M. Bianchi, H. Guo, R. Harsh, S. Mañas-Valero, E. Coronado, A. J. Martínez-Galera, P. Hofmann, J. A. Miwa, and M. M. Ugeda, Proximity effects on the charge density wave order and superconductivity in single-layer NbSe$_2$, ACS Nano 15, 19430 (2021).

[50] F. Flicker and J. v. Wezel, Charge order from orbital-dependent coupling evidenced by NbSe$_2$, Nat. Commun. 6, 7034 (2015).

[51] S.-G. J. Ángel, O. Pablo, G. Francisco, and C. Eric, Electronic structure of 2H-NbSe$_2$ single-layers in the CDW state, 2D Materials 3, 035028 (2016).

[52] M. Ruby, Y. Peng, F. von Oppen, B. W. Heinrich, and K. J. Franke, Orbital Picture of Yu-Shiba-Rusinov Multiplets, Phys. Rev. Lett. 117, 186801 (2016).

[53] F. Massee, Q. Dong, A. Cavanna, Y. Jin, and M. Aprili, Atomic scale shot-noise using cryogenic MHz circuitry, Rev. Sci. Instrum. 89, 093708 (2018).

[54] M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B. W. Heinrich, and K. J. Franke, Tunneling Processes into Localized Subgap States in Superconductors, Phys. Rev. Lett. 115, 087001 (2015).