Yu-Shiba-Rusinov states in the charge-density modulated superconductor NbSe2

TL;DR

This study investigates how charge-density waves in NbSe2 influence Yu-Shiba-Rusinov states created by magnetic adatoms, revealing site-dependent effects and potential for topological nanostructure applications.

Contribution

It provides a detailed experimental and theoretical analysis of YSR states in NbSe2, highlighting the impact of incommensurate charge-density waves on these states.

Findings

Charge-density wave modulates YSR state energies.

YSR wavefunctions depend on adsorption site relative to CDW.

Model calculations explain site-dependent YSR properties.

Abstract

NbSe$_2$ is a remarkable superconductor in which charge-density order coexists with pairing correlations at low temperatures. Here, we study the interplay of magnetic adatoms and their Yu-Shiba-Rusinov (YSR) bound states with the charge density order. Exploiting the incommensurate nature of the charge-density wave (CDW), our measurements provide a thorough picture of how the CDW affects both the energies and the wavefunctions of the YSR states. Key features of the dependence of the YSR states on adsorption site relative to the CDW are explained by model calculations. Several properties make NbSe$_2$ a promising substrate for realizing topological nanostructures. Our results will be important in designing such systems.

Tables

Table 1. Table S1: Values of the fitting parameters t n subscript 𝑡 𝑛 t_{n} in Eq. ( S10 ) in meV.

	$t_0$	$t_1$	$t_2$	$t_3$	$t_4$	$t_5$
band 1	10.9	86.8	139.9	29.6	3.5	3.3
band 2	203.0	46.0	257.5	4.4	-15.0	6.0

Equations

$$E_{\text{YSR}}=\pm \Delta \frac{1-{\left(\pi JS{\nu }_0/2\right)}^2}{1+{\left(\pi JS{\nu }_0/2\right)}^2}.$$

$$H=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{c}_{\mathbf{k}}^{†}{c}_{\mathbf{k}}+\sum _{\mathbf{q}}{\omega }_{\mathbf{q}}{b}_{\mathbf{q}}^{†}{b}_{\mathbf{q}}+\sum _{\mathbf{k},\mathbf{q}}{g}_{\mathbf{q}}(b_{\mathbf{q}}+b_{-\mathbf{q}}^{†})c_{\mathbf{k}+\mathbf{q}}^{†}{c}_{\mathbf{k}}.$$

$$H=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{c}_{\mathbf{k}}^{†}{c}_{\mathbf{k}}+\sum _{\mathbf{k}}\sum _{\mathbf{Q}}{g}_{\mathbf{Q}}⟨b_{\mathbf{Q}}+b_{-\mathbf{Q}}^{†}⟩c_{\mathbf{k}+\mathbf{Q}}^{†}{c}_{\mathbf{k}}$$

$$H_{\mathrm{el}-\mathrm{ph}}=\int \mathrm{d}\mathbf{r}V(\mathbf{r}){\psi }^{†}(\mathbf{r})\psi (\mathbf{r}),$$

$$V(\mathbf{r})=\sum _{\mathbf{Q}}{e}^{-i\mathbf{Q}\cdot \mathbf{r}}{g}_{\mathbf{Q}}⟨b_{\mathbf{Q}}+b_{-\mathbf{Q}}^{†}⟩.$$

$${\mathbf{a}}_1$$

$${\mathbf{a}}_2$$

$${\mathbf{b}}_1$$

$${\mathbf{b}}_2$$

$${\mathbf{Q}}_1$$

$${\mathbf{Q}}_2$$

$${\mathbf{Q}}_3$$

$$\begin{array}{rl}E(\mathbf{k})& =t_0+t_1(2\cos \xi \cos \eta +\cos 2\xi )+t_2(2\cos 3\xi \cos \eta +\cos 2\eta )+t_3(2\cos 2\xi \cos 2\eta +\cos 4\xi )\\ & +t_4(\cos \xi \cos 3\eta +\cos 5\xi \cos \eta +\cos 4\xi \cos 2\eta )+t_5(2\cos 3\xi \cos 3\eta +\cos 6\xi )\end{array}$$

$$H_{\mathrm{exch}}=-\frac{JS/2}{3}\sum _{i,j}\sum _{s,s^{\prime }}{c}_{is}^{†}{\sigma }_{s{s}^{\prime }}^z{c}_{j{s}^{\prime }}$$

$$H_{\mathrm{pot}}=\frac{K}{3}\sum _{i,j}\sum _s{c}_{is}^{†}{c}_{js}.$$

$$H_{\mathrm{cdw}}=\sum _{\mathbf{r},\sigma }{c}_{\mathbf{r}\sigma }^{†}V(\mathbf{r})c_{\mathbf{r}\sigma },$$

$$V(\mathbf{r})=V_0\left[\cos ({\mathbf{Q}}_1\cdot \mathbf{r}+{\varphi }_1)+\cos ({\mathbf{Q}}_2\cdot \mathbf{r}+{\varphi }_2)+\cos ({\mathbf{Q}}_3\cdot \mathbf{r}+{\varphi }_3)\right],$$

$$V(\mathbf{r})=V_0\left[2\cos \left(\sqrt{3}\frac{q}{2}(x-x_0)\right)\cos \left(\frac{q}{2}(y-y_0)\right)+\cos \left(q(y-y_0)+{\varphi }_3\right)\right],$$

$$E=\pm \Delta \frac{1-{\alpha }^2+{\beta }^2}{\sqrt{4{\alpha }^2+(1-{\alpha }^2+{\beta }^2{)}^2}},$$

$$I(V)\propto {\int }_{-\infty }^{+\infty }\text{d}E{\rho }_{\text{t}}(E-eV,T){\rho }_{\text{s}}(E,T)[f(E-eV,T)-f(E,T)]|M_{\mu ,\nu }{|}^2.$$

$$\frac{\partial I}{\partial V}\propto {\int }_{-\infty }^{+\infty }\text{d}E\frac{\partial {\rho }_{\text{t}}(E-eV)}{\partial V}{\rho }_{\text{s}}(E,T)[f(E-eV,T)-f(E,T)]$$

$$+{\int }_{-\infty }^{+\infty }\text{d}E{\rho }_{\text{t}}(E-eV){\rho }_{\text{s}}(E,T)\frac{\partial f(E-eV,T)}{\partial V}.$$

$${\rho }_{\text{t}}(E)=\text{sgn}(E)\Re \left(\frac{E-i\Gamma }{\sqrt{(E-i\Gamma {)}^2-{\Delta }_{\text{t}}^2}}\right).$$

$$\overrightarrow{\frac{\partial I}{\partial V}}(V)\propto \mathbf{K}(E,V,T)\overrightarrow{{\rho }_{\text{s}}}(E).$$

$${\mathbf{K}}_{ij}(E_j,V_i,T)=\text{d}E\frac{\partial {\rho }_{\text{t}}(E_j-e{V}_i)}{\partial V}[f(E_j-e{V}_i,T)-f(E_j,T)]$$

$$+\text{d}E{\rho }_{\text{t}}(E_j-e{V}_i)\frac{\partial f(E_j-e{V}_i,T)}{\partial V}.$$

$$\text{LDOS}(E)=D_0+A_1{e}^{-(E-E_{\alpha }{)}^2/(2{\sigma }^2)}+A_2{e}^{-(E+E_{\alpha }{)}^2/(2{\sigma }^2)}$$

$$+B_1{e}^{-(E-E_{\beta }{)}^2/(2{\sigma }^2)}+B_2{e}^{-(E+E_{\beta }{)}^2/(2{\sigma }^2)}.$$

Full text

At the same time, recent experiments suggest that NbSe2 has attractive properties in the context of engineering topological superconducting phases Sticlet and Morari (2019); Głodzik and Ojanen (2019). Following the seminal work of Nadj-Perge et al. Nadj-Perge et al. (2014), a topological phase and Majorana bound states can appear in chains of magnetic adatoms on superconducting substrates. One way to understand these remarkable phases is based on coupling Yu-Shiba-Rusinov (YSR) bound states Nadj-Perge et al. (2013); Pientka et al. (2013) in the superconducting gap which are induced by individual adatoms Yu (1965); Shiba (1968); Rusinov (1969); Yazdani et al. (1997). Due to its layered structure, YSR states in NbSe2 fall off laterally much more slowly than in three-dimensional superconductors Ménard et al. (2015). This enhances the coupling between the YSR states of neighboring magnetic impurities Kezilebieke et al. (2018), potentially increasing the energy scale and thus the stability of possible topological phases. NbSe2 is also a relatively robust material which may enable one to build nontrivial adatom structures from the bottom up by manipulation using a scanning tunneling microscope tip.

The YSR states of different adatoms hybridize particularly strongly when the YSR energies of the isolated adatoms are identical and their wavefunctions extend sufficiently far along their connecting line Ruby et al. (2018); Choi et al. (2018); Kezilebieke et al. (2018); Kamlapure et al. (2018); Kim et al. (2018). For NbSe2, this is nontrivial to achieve due to the CDW ordering. Indeed, the CDW is incommensurate with the underlying atomic lattice so that adatoms which are nominally identical in terms of their atomic adsorption site, will in general be positioned differently with respect to the CDW. Here, we observe that both the energy and the wavefunction of YSR states depend sensitively on the adatom position relative to the CDW, and provide model calculations to understand these effects qualitatively. The dependence on the CDW imposes requirements which are crucial to take into account when designing adatom structures to realize topological superconducting phases.

To isolate the effect of the CDW on the YSR states, we focus on adatoms on clean $2H$-NbSe2 surfaces instead of buried impurities as used in previous experiments Ménard et al. (2015); Senkpiel et al. (2019). This gives us full control over the adsorption site and allows us to compare adatoms which differ only in their locations relative to the CDW. Indeed, numerous previous studies Yazdani et al. (1997); Ji et al. (2008); Ruby et al. (2015, 2016); Choi et al. (2017); Cornils et al. (2017); Kamlapure et al. (2018) on other substrates have shown that in the absence of CDW ordering, YSR spectra and wavefunctions are fully reproducible for a particular adsorption configuration of the adatom on the substrate.

A clean $2H$-NbSe2 surface is prepared by cleaving a bulk crystal (grown by iodine vapor transport Rahn et al. (2012)) in ultra-high vacuum. Scanning tunneling spectroscopy (performed by standard lock-in techniques at a temperature of $1.1\,$K) shows that the superconducting energy gap of the substrate is framed by characteristic coherence peaks in an energy range of $2.0-2.7\operatorname{\,\mathrm{meV}}$ (black $\mathrm{d}I/\mathrm{d}V$ spectrum in Fig. 1e). The origin of this peculiar quasiparticle density of states (DOS) has been discussed in terms of multiband or anisotropic superconductivity Hess et al. (1990); Yokoya et al. (2001); Noat et al. (2010, 2015); Rodrigo and Vieira (2004); Boaknin et al. (2003); Fletcher et al. (2007); Guillamón et al. (2008a, b). Notice that we use superconducting Pb tips to improve the energy resolution beyond the Fermi-Dirac limit. The convolution of the DOS of tip and substrate shifts all spectral features by the tip’s energy gap $\Delta_{\mathrm{t}}\approx 1.35\operatorname{\,\mathrm{meV}}$ (for details see Supporting Information (SI) SM ).

Scanning tunneling microscopy (STM) images of the clean surface clearly reveal the modulation of the local DOS induced by the CDW superimposed on the atomic corrugation (Fig. 1a). The CDW has a lattice constant $a_{\mathrm{cdw}}\gtrsim 3a$, making it incommensurate with the underlying atomic lattice. Thus, the phase of the CDW relative to the atomic lattice varies smoothly across the surface Gye et al. (2019); Guster et al. (2019). The close-ups show two extremal cases for which the maxima coincide either with a Se atom (chalcogen-centered, CC, red circle in Fig. 1a and Fig. 1b) or a hollow site (hollow-centered, HC, green circle in Fig. 1a and Fig. 1b SC (1)).

Fe atoms (deposited at temperatures below 12 K into the STM) appear as two stable species with different apparent heights (Fig. 1c,d). By atomic-resolution imaging (for details see SI SM ), we assign adatoms with small and large apparent height to adsorption in the two inequivalent hollow sites of the terminating Se layer, identified as hollow-centered (HC) and metal-centered (MC) sites, respectively (Fig. 1b) SC (3). $\mathrm{d}I/\mathrm{d}V$ spectra taken above the centers of HC and MC adatoms (blue and orange in Fig. 1e) show several YSR states inside the superconducting gap as well as in the energy range of the substrate’s coherence peaks as has been observed previously for buried impurities Senkpiel et al. (2019). The energy and intensity of the YSR states differ between the two species. Presumably, the splitting of the adatom $d$-levels is sensitive to the different local environments of the adsorption sites, which in turn affects the potential- and exchange-scattering strengths Ruby et al. (2016).

To isolate the influence of the CDW, we focus on Fe atoms which are all sitting in the same atomic adsorption site (HC) in the following. Six different atoms (labeled by I-VI) are shown in Fig. 2a. The corresponding $\mathrm{d}I/\mathrm{d}V$ spectra (Fig. 2b) reveal that the energy and intensity of their YSR states differ strongly even though their atomic adsorption sites are identical with respect to the unperturbed lattice.

YSR states in the energy range of the superconducting coherence peaks are difficult to disentangle from the background. To avoid this complication, we focus on deep-lying YSR states, specifically the two lowest YSR pairs labeled as $\pm\alpha$ and $\pm\beta$ in the close-up view presented in Fig. 2c. $\mathrm{d}I/\mathrm{d}V$ maps recorded at the corresponding bias voltages $V_{\alpha}$ and $V_{\beta}$ are shown in Fig. 3a for adatoms I-VI. The main panels show the extended patterns while the insets focus on the immediate vicinity of the adatom. All extended maps show patterns with oscillating intensity SC (2). The overall symmetries of the patterns clearly differ between adatoms I-VI.

The symmetry of the patterns associated with YSR states is expected to originate from the anisotropy of the Fermi surface Salkola et al. (1997); Flatté and Byers (1997) and the local crystal field Ruby et al. (2016). The threefold-symmetric atomic adsorption site together with the sixfold symmetry of the Fermi surface should therefore lead to $D_{3}$ symmetry (threefold rotation as well as three mirror axes, cf. dashed lines in Fig. 1b). Interestingly, this $D_{3}$ symmetry is observed for both YSR states $\alpha$ and $\beta$ of adatoms I and III, which reside at a maximum and a minimum of the CDW, respectively (Figs. 2d, 3a). However, this symmetry is lost in both, the long-range and the immediate vicinity of the adatoms for the other atoms shown in Fig. 3a. For adatoms II, IV, and V, the symmetry is reduced to $D_{1}$ symmetry (single mirror axis) while no symmetry axis can be discerned for adatom VI.

The symmetry reductions of the YSR patterns coincide with the reductions of the local symmetry of the adsorption sites by the CDW (Fig. 2d). Adatoms I and III are positioned at a maximum or a minimum of the CDW, so that the CDW respects the atomic $D_{3}$ symmetry (HC structure of the CDW). Atoms II, IV, and V are located on one of the three equivalent symmetry axes connecting the high-symmetry positions, but unlike adatoms I and III do not directly fall on an extremum (CDW in the CC domain). Thus, the CDW breaks the atomic $D_{3}$ symmetry, leaving only $D_{1}$ symmetry consistent with the observed YSR patterns. Finally, the position of adatom VI is totally asymmetric with respect to the CDW, which is reflected in the absence of any symmetry in the corresponding $\mathrm{d}I/\mathrm{d}V$ maps.

The microscopic physics behind these symmetry reductions can be understood theoretically within a phenomenological mean-field description of the CDW combined with a tight-binding model of the $2H$-NbSe2 band structure. Foregoing a realistic description of the Fe $d$-orbital physics, we model the adatom as a classical impurity with isotropic potential and exchange couplings to the substrate (see SI SM for details). Due to the electron-phonon interaction, the CDW acts on the electrons as a weakly incommensurate, static periodic potential. This potential reduces the symmetry of the YSR wavefunctions in a manner consistent with the experimental results. Corresponding numerical results are shown in Fig. 3b for a subset of the adsorption sites investigated experimentally (see SI SM for the remaining sites and further details). We find that the wavefunctions (and the number) of the YSR states are quite sensitive to details of the band structure and the neglected $d$-orbital physics, so that one expects agreement only for qualitative aspects of symmetry.

So far, the assignments $\alpha$ and $\beta$ for the YSR states simply refer to the two deepest resonances. Analyzing the observed YSR patterns allows us to track these two YSR states as a function of position relative to the CDW. For the YSR state $+\beta$ of adatom I, the short-range pattern exhibits a triangular shape (Fig. 3c) with high intensity on its sides and slightly smaller intensity at its vertices (orange circles). In contrast, the YSR state $+\alpha$ for adatom I lacks intensity at the vertices (blue circles) and thus exhibits pronounced nodes. Signatures of this distinction persist at the lower-symmetry positions II, IV, and V. While the intensity distributions change for both states, we still observe distinct intensities at one of the vertices which distinguish between the $+\alpha$ and $+\beta$ resonances (while present, the distinction is somewhat less pronounced for adatom V, see also SI SM ). A similar distinction cannot be made for adatom III. We attribute this to the fact that in this case, both the $\pm\alpha$- and the $\pm\beta$-resonances lie very close to the Fermi level. This prevents us from individually resolving them within our energy resolution (Fig. 2c, green trace).

With these characteristics, we can track the resonances as a function of position relative to the CDW and probe its effect on the YSR energies. As the CDW transforms smoothly between HC and CC domains, we can study numerous positions along its symmetry axis (such as positions II, IV, and V in Fig. 2d, but also in between). Altogether, we analyze approximately $90$ adatoms adsorbed close to one of the three equivalent CDW symmetry axes, combining data from several samples and Pb tips. For all adatoms, we extract the energies of the YSR states by deconvolving the $\mathrm{d}I/\mathrm{d}V$ spectra to remove the influence of the Pb tip (see SI SM ). Figure 4a plots the energies of both the $\alpha$ and $\beta$ states as a function of adatom position along a high-symmetry axis of the CDW. Here, we assume that the energy dependences of both resonances follow the same trend which implies that the energy of the $\alpha$ resonance changes sign and crosses the quantum phase transition from a screened-spin to a free-spin ground state Choi et al. (2017); Franke et al. (2011); Farinacci et al. (2018); Malavolti et al. (2018) as a function of adatom position. Note that some regions (shaded in Fig. 4a) of adatom positions remain inaccessible for Fe atoms in HC sites, because the CDW maxima avoid MC domains SC (1).

A qualitatively similar dependence of the energy of the YSR states is found in our model calculations (Fig. 4c and SI SM ). The YSR energies are (anti)correlated with the CDW potential and the local DOS, with the type of correlation depending on whether the impurity is in the weak- or strong-coupling regime and on the strength of potential scattering by the impurity SM . (We find that the CDW potential is anticorrelated with the local DOS.) Except for very strong potential scattering (see SI SM ), our results are consistent with a simple model for the energy of YSR states Yu (1965); Shiba (1968); Rusinov (1969),

[TABLE]

Here, $\Delta$ is the superconducting gap, $S$ the (classical) spin, $J$ the exchange coupling, and $\nu_{0}$ the local DOS of the superconductor. This predicts that YSR energies are (anti)correlated with the local DOS for impurities in the strong (weak) coupling limit. Experimentally, we record a constant-height $\mathrm{d}I/\mathrm{d}V$ map at zero bias above $T_{\mathrm{c}}$ and remove the atomic corrugation by an FFT filter (cf. inset Fig. 4b). The remaining signal is proportional to the local DOS due to the CDW at the Fermi level. We find that the YSR energies correlate with the local DOS as seen in Fig. 4a,b. This suggests that the YSR states are in most of the adsorption sites in the strong-coupling regime. Only when the Fe atoms are located in the minima of the CDW, resonance $\alpha$ undergoes the quantum phase transition to the free-spin ground state.

In conclusion, we find that the properties of YSR states are strongly affected by CDW ordering. The CDW shifts the energy of YSR states and reduces the symmetry of YSR wavefunctions both in the immediate vicinity of and further from the adatoms. These observations impose stringent constraints on the design of strongly coupled magnetic adatom structures on NbSe2 substrates. Neverthless, NbSe2 remains a highly attractive substrate in this context, owing to the long-range nature of the YSR states and its suitability for adatom manipulation. Our results extend to truly two-dimensional NbSe2 substrates with their promise of device applications. Indeed, single-layer NbSe2 also exhibits an incommensurate CDW Guster et al. (2019); Ugeda et al. (2016). Importantly, pairing correlations in single-layer NbSe2 are remarkably robust against in-plane fields as a consequence of spin-valley locking (Ising superconductivity) Xi et al. (2016). This may enable measurements of YSR states at magnetic fields which polarize the adatom spin without suppressing superconductivity. The interaction of YSR states with superconductivity, charge-density order, and external magnetic fields entails substantial flexibility to tune into or control a topologically non-trivial superconducting phase. Finally, our study suggests that YSR states provide a sensitive probe of these competing interactions which should extend to other systems.

Supporting Information

S1 Theoretical considerations

S1.1 Mean-field description of the charge-density wave

We resort to a phenomenological mean-field description. We start with a generic Fröhlich-type Hamiltonian describing electrons coupled to phonons Grüner (1994),

[TABLE]

Here, $c_{\bf k}$ annihilates electrons with momentum ${\bf k}$ and energy $\epsilon_{\bf k}$, $b_{\bf q}$ annihilates a phonon with wavevector ${\bf q}$ and frequency $\omega_{\bf q}$, and $g_{\bf q}$ denotes the strength of the electron-phonon coupling. Within mean-field theory, we assume that certain phonon modes (denoted by ${\bf Q}$) go soft due to their coupling to the electronic system and develop a finite expectation value. Restricting to the lowest Fourier components and neglecting phonon dynamics beyond the static charge-density wave (CDW) distortions, we find the mean-field Hamiltonian

[TABLE]

for the electronic degrees of freedom. Within this approach, the CDW acts on the electrons as a periodic potential $V({\bf r})$. Indeed, the electron-phonon interaction term can also be expressed in real space as

[TABLE]

where the mean-field CDW potential takes the form

[TABLE]

To find the vectors Q, we consider a triangular lattice of Nb atoms with lattice vectors

[TABLE]

and bond length $a=3.445\,$Å. The corresponding reciprocal lattice vectors ${\bf b}_{i}$ satisfying ${\bf a}_{i}\cdot{\bf b}_{j}=2\pi\delta_{ij}$ are

[TABLE]

The unit cell of the CDW has a linear dimension which is approximately three times larger than the unit cell of the atomic lattice, ${\bf a}^{\rm cdw}\approx 3{\bf a}_{i}$ (for $i=1,2$). Correspondingly, the reciprocal lattice vectors of the CDW are approximately a factor of three smaller, ${\bf b}^{\rm cdw}_{i}\approx{\bf b}_{i}/3$. Then, the first harmonics of the CDW have wavevectors

[TABLE]

where $q=\frac{4\pi}{3\sqrt{3}a}$, and $\delta\ll 1$ accounts for the fact that the CDW is not exactly commensurate with the lattice.

S1.2 Tight-binding model

To illustrate the effect of the CDW on the YSR state wavefunction, we perform model calculations within an effective tight-binding description. Specifically, we implement a model for the NbSe2 band structure first introduced in Ref. Smith et al. (1985) and later used in Refs. Rossnagel et al. (2005); Inosov et al. (2008, 2009); Rahn et al. (2012); Ménard et al. (2015). The model uses that the relevant states near the Fermi energy have mostly Nb character and therefore focuses on one atomic $d$-orbital per niobium atom. (The model neglects an additional band centered at the $\Gamma$ point which predominantly derives from Se orbitals and has strongly three-dimensional character Bawden et al. (2016).) Reducing the band-structure problem to the triangular Nb lattice makes it necessary to include hopping up to fifth nearest neighbors to reproduce the NbSe2 band structure. The resulting band structure

[TABLE]

is sixfold symmetric about the $\Gamma$ point. Here, we defined $\xi=k_{x}/2$ and $\eta=\sqrt{3}k_{y}/2$. The hopping strengths $t_{n}$ are used as fitting parameters to reproduce the NbSe2 band structure. The values of the parameters $t_{n}$ are given in Table S1 Rahn et al. (2012) for the two bands. The parameters for band 1 reproduce the inner cylindrical bands of NbSe2, while the parameters for band 2 reproduce the outer bands. Figure S1a shows the Fermi surfaces of the two bands, with the Brillouin zone outlined in gray.

While this model accounts for the symmetries of the band structure and the approximate shapes of the Fermi surfaces, it is limited in other ways. The model neglects interlayer coupling which, however, is expected to be weak compared to the intralayer couplings. In this approximation, it suffices to consider a single (tri)layer of $2H$-NbSe2, consisting of a Nb layer sandwiched between two Se layers with equal orientations of the triangles. The model also neglects spin-orbit coupling Xi et al. (2016). This precludes a one-to-one identification of the two sets of bands in Fig. S1a with specific spin directions.

The elementary layer of bulk $2H$-NbSe2 consists of two trilayers with the orientations of the Se lattices rotated by $180^{\circ}$ with respect to one another. Neglecting interlayer coupling, the band structure of the second trilayer is identical to the first (but differs in the spin assignments). Thus, the bulk crystal has twofold degenerate bands, consistent with the fact that it has an inversion center located between the trilayers. The adatom predominantly couples to one of the trilayers and the band structure in Eq. (S10) should be a good starting point. In the following, we follow Ref.  Ménard et al. (2015) in neglecting the splitting between the two sets of bands and model the system either by band 1 or band 2, taking these bands as spin degenerate.

Superconductivity is included by incorporating the tight-binding Hamiltonian into a Bogoliubov-deGennes Hamiltonian with conventional (isotropic) $s$-wave pairing. We choose the pairing strength to be $\Delta=1$ meV. In experiment, one finds multiple coherence peaks for $2H$-NbSe2. Since we focus on YSR states far from the gap edge, we choose a representative value for $\Delta$ which falls into the range of the peak distribution ($0.7-1.4$ meV) found in experiment.

The multiple coherence peaks originate from multiband or anisotropic superconductivity, or both. Our modeling does not account for these effects. In principle, an anisotropic gap function affects the wavefunction patterns of the YSR states. However, the (observable) spatial extent of YSR states is smaller than the coherence length through which effects of the pairing function would enter. We thus expect that the effect of an anisotropic pairing function on the YSR wavefunctions is weak and therefore account for superconductivity via an isotropic pairing function.

S1.3 Coupling to the magnetic impurity

While the origin of the magnetic moment and the resulting exchange coupling can be understood starting with an Anderson model and its mapping to a Kondo Hamiltonian in the absence of valence fluctuations (see Ref. Hewson (1993) for an extensive pedagogical discussion), YSR states are typically well described within a model which treats the impurity spin in the Kondo model as classical Yu (1965); Shiba (1968a); Rusinov (1969). For YSR states in NbSe2, this model was previously used in Ref. Ménard et al. (2015). Following this reference, we model a magnetic impurity as a classical spin $S$ which interacts with the substrate electrons through the exchange interaction $JS\sigma_{z}/2$ and accompanying potential scattering $K$. In our experiment, we focus on magnetic impurities located above the center of a Nb triangle (hollow-centered, see Fig. 1b of the main text). The impurity induces not only separate exchange and potential couplings to its three neighbors, but also nonlocal exchange couplings which scatter electrons between the three sites via the impurity. It is essential to retain the latter to ensure that the impurity induces only one pair of YSR resonances. If we assume identical hopping amplitudes between the three sites and the impurity and assume that the impurity spin is aligned along the $z$-direction, we obtain the exchange coupling

[TABLE]

and potential scattering terms

[TABLE]

Here, the sums over $i$ and $j$ run over the three nearest-neighbor Nb sites of the impurity, $\sigma^{z}$ denotes a Pauli matrix, and $c_{is}$ annihilates an electron on site $i$ with spin $s$.

Figures S1b and c show the probability densities of the YSR states for (a spin-degenerate version of) bands 1 and 2, respectively, for a specific choice of impurity parameters. The spatial distribution of the probability density is threefold symmetric in both cases. Due to the more isotropic character of band 1, the corresponding YSR state is also more spatially isotropic than its band 2 counterpart, which exhibits a more pronounced anisotropy. In a more realistic modeling including the effects of spin-orbit coupling on the band structure (but within the approximation that the impurity is coupled only to a single trilayer), the YSR states would be simultaneously coupled to spin-textured, but polarized versions of both bands. While this implies that our results for the YSR wavefunctions are only qualitative, the symmetry of the YSR wavefunctions should not be affected by this approximation. Similarly, we expect that the qualitative variation of the YSR energies with the location of the adsorption site relative to the CDW would also be unaffected.

We note that in our approximation, there is only a single pair of YSR states for a classical magnetic impurity. In fact, this conclusion would remain valid even when spin-orbit splitting of the bands was taken into account. On the one hand, it is known that even in multiband situations, a classical magnetic impurity induces only a single pair of YSR states. In fact, if we write the matrix of exchange couplings $J_{\alpha\alpha^{\prime}}$ between bands $\alpha$ and $\alpha^{\prime}$, it rather robustly satisfies the relation $J_{11}J_{22}=J_{12}J_{21}$. With this relation, a classical magnetic impurity generates only a single pair of YSR states. For instance, this was shown explicitly for a generic model motivated by the multiband superconductor MgB2 Moca et al. (2008) and for a quantum dot coupled to multiple superconducting leads Kirs̆anskas et al. (2015). In the presence of spin-orbit coupling, the same conclusions can be shown to remain valid. This can be seen by expressing the Hamiltonian in terms of creation and annihilation operators for Kramers pairs of electron states rather than conventional spin states.

S1.4 Coupling to the charge-density wave

As explained above, the effect of the CDW on the electrons of the substrate can be described through a periodic potential. Within the tight-binding model, we include the CDW as a modulation of the on-site potential,

[TABLE]

where the sum runs over the lattice sites of the triangular Nb lattice. The CDW potential $V(\mathbf{r})$ is constructed from the main Fourier components of the CDW $[$see Eq. (S9)$]$,

[TABLE]

where $V_{0}$ denotes the amplitude of the CDW potential. Except for $\phi_{1}=\phi_{2}=\phi_{3}$ and similar fine-tuned cases, typical choices of the phases $\phi_{i}$ yield CDW potentials with the desired symmetry and shape, exhibiting an absolute maximum, an absolute minimum, and a local minimum/saddle point. For definiteness, we choose $\phi_{1}=\phi_{2}=0$.

Due to the small deviations from commensurability, the CDW shifts slowly as a function of position relative to the underlying atomic lattice. On the scale of the YSR states, these shifts can be considered constant so that the CDW potential can be written as

[TABLE]

where $q=4\pi/3\sqrt{3}$. The offsets $x_{0}$ and $y_{0}$ originate from the deviation $\delta$ from commensurability, and describe the local shift of the CDW relative to the atomic lattice. In the experiment, the measured maxima of the CDW correspond to regions where the lattice deformation compresses the ions, creating an attractive potential for the electrons. For this reason we choose the potential $V(\mathbf{r})$ to have minima where the measured tunneling density of states has maxima, i.e., position I of the adatom corresponds to a minimum of $V(\mathbf{r})$ and position III corresponds to a maximum. In principle, the CDW may also modulate the hopping parameters of the tight-binding model, but the effective long-range hopping processes of the model make this cumbersome to include.

We note that this model also includes configurations of the CDW relative to the atomic lattice which are not observed experimentally (see Fig. S18). Theoretical results for these unobserved configurations are only included for completeness. The results for observed configurations are not affected since the Moire-like variations occur on relatively large scales.

S1.5 Numerical results

Within this model, we can study how the wavefunction and the energy of YSR states depend on the adatom position relative to the CDW. As discussed above, the model cannot reproduce the experimental results quantitatively as it neglects various relevant effects. Nevertheless, the model is expected to capture the symmetries of the YSR wavefunctions and to provide insights into the qualitative dependence of the YSR energy on the adsorption site relative to the CDW.

Figure S2 shows the probability density of the YSR wavefunction for the six different positions of the adatom relative to the CDW, as discussed in the main text. The results are obtained for the parameters of band 2. The adatoms in positions I and III are at the maximum and the minimum of the CDW, respectively, i.e., at points with a threefold symmetric environment. Consequently, the corresponding density plots (Fig. S2a and c) present threefold symmetry. Positions II, IV, and V exhibit reflection symmetry about an axis that passes through the adatom. For the corresponding plots shown in Fig. S2b, d, and e, this axis is aligned along the vertical direction. Finally, position VI has no symmetries with respect to the CDW and consequently, the state shown in Fig. S2f also exhibits no symmetries. Interestingly, the YSR state seems to retain the original threefold symmetry (without CDW) to some degree (compare with Fig. S1c).

Figure S3 shows equivalent plots for band 1. Due to the more isotropic character of band 1, the threefold symmetry is less pronounced in these plots. Remarkably, position VI (Fig. S3f) does not preserve any of the original threefold symmetry in the vicinity of the adatom. More generally, a comparison between Figs. S2 and S3 shows that the detailed YSR wavefunctions depend quite sensitively on the band structure parameters.

Finally, we discuss the correlation of the energy of the YSR state with the modulation of the local density of states by the CDW, as shown in Fig. 4 in the main text. Within a simple model of a classical magnetic impurity, the energy of YSR states depends on the strength of the exchange coupling as measured by the dimensionless parameter $\alpha=\pi\nu_{0}JS/2$, where $\nu_{0}$ is the density of states, $J$ the exchange coupling, and $S$ the impurity spin. Thus, one indeed expects the energy of the YSR state to depend on the density of states $\nu_{0}$. The direction of the energy shift depends on whether the impurity is weakly or strongly coupled, corresponding to unscreened or (partially) screened impurity spins, respectively. The energy of the positive-energy YSR state decreases (increases) with the density of states, depending on whether the impurity spin is unscreened (screened). Variations in the local density of states due to the CDW should thus be reflected in the energy of the YSR states. Within this picture, the observed correlations suggest that in our experiment the impurity is in the strong-coupling limit.

This interpretation is consistent with our theoretical results. However, we find that in general the CDW affects the energy of the YSR state through two mechanisms. In addition to the density of states effect, the CDW potential has a second effect which can be roughly rationalized as the CDW affecting the strength of potential scattering from the impurity and thereby shifting the YSR energy. Including potential scattering within the simple model, the energy of a YSR state is given by Yu (1965); Shiba (1968b); Rusinov (1969)

[TABLE]

where $\beta=\pi\nu_{0}K$ is a dimensionless measure of the strength of the potential scattering $K$. The transition between weak and strong coupling takes place at $\alpha^{2}=1+\beta^{2}$ where the energy of the YSR state passes through zero. Roughly, the CDW potential $V$ seems to shift $K$ to $K+V(\mathbf{r}_{0})$, where $\mathbf{r}_{0}$ denotes the adsorption site of the impurity. For comparison with experiment, it is important to note that an increase in the CDW potential decreases the density of states. This can be seen in e.g. Fig. S4g, which shows numerical results for the average local density of states of the three sites coupled to the impurity (orange) in the presence of the CDW potential (blue) as obtained from a straight-forward band structure calculation. Thus, the experiment implies an anticorrelation between the CDW potential $V$ and the energy of the YSR states.

To illustrate these two effects of the CDW potential, we calculate the YRS energy for different adatom positions and various values of $K$. Consider first the results for weak coupling. We find that depending on the strength $K$ of potential scattering, the CDW potential is correlated or anticorrelated with the energy of the YSR state. For $K=180$ meV and $K=0$, the energy is correlated with the CDW potential, see Fig. S4c and d. In contrast, we find anticorrelations for $K=-180$ meV, see Fig. S4b. The correlations observed at $K=0$ are consistent with the effect of the CDW potential on the density of states. At this $K$, the energy shift of the YSR state due to the change in potential scattering, $K$ to $K+V(\mathbf{r}_{0})$, is weak since the energy of the YSR states is only weakly dependent on $K$ near the minimum at $K=0$. This is illustrated in Fig. S4a which shows the dependence of the energy of the YSR state on $K$ in the weak-coupling regime. Specifically, the black curve averages data for the energy of the YSR state over different adatom positions. The red curve corresponds to Eq. (S16) with the same values for $JS$ and $K$ and an appropriately chosen density of states $\nu_{0}$. Both curves exhibit similar ‘parabolic’ behaviors.

For the strongly coupled regime, we find anticorrelations between the CDW potential and the energy of the YSR state – and thus correlations between the local density of states and the YSR energy – for all values of $K$, see Fig. S4f,g and h. These correlations between the density of states and the YSR energy are consistent with experiment. To start with, the anticorrelations at $K=0$ are consistent with the expected density of states effect according to Eq. (S16). For strong coupling, the dependence of the YSR energy on $K$ is generally weaker compared to the weak coupling case, see Fig. S4e. For this reason, the shift of $K$ by the CDW potential no longer overcomes the density-of-states effect. Thus, unless the impurity is a surprisingly strong potential scatterer, we conclude that in experiment the impurity is in the strong coupling limit.

S2 Experimental details

S2.1 Superconducting tip

In order to increase the energy resolution beyond the Fermi-Dirac limit, all measurements presented in the manuscript were performed using a superconducting Pb tip fabricated by indenting a W tip into a clean Pb(111) surface until the tip exhibits a bulk-like superconducting gap. Details of the preparation procedure for tip and sample are described in Ref. Ruby et al. (2015a). The tunneling current is proportional to the convolution of the density of states of the sample $\rho_{\text{S}}$ and tip $\rho_{\text{t}}$ and can be expressed as Tersoff and Hamann (1985):

[TABLE]

Here, $f(E,T)$ is the Fermi-Dirac-distribution function accounting for the thermal occupation of the states in the tip and the sample, and $|M_{\mu,\nu}|$ is the tunneling matrix element between the initial state $\mu$ and final state $\nu$.

This leads to the differential conductance:

[TABLE]

Before investigating the NbSe2 sample, the tip was characterized on a clean Pb(111) sample. A spectrum taken on bare Pb(111) after preparation of the Pb tip is shown in Fig. S5. The two superconducting energy gaps of the substrate are well resolved Ruby et al. (2015b). The Pb tip was characterized by the fit procedure described in detail in Ref. Ruby et al. (2015b) where a BCS-like density of states is assumed for the tip:

[TABLE]

From this fitting procedure we determine the values for the superconducting energy gap $\Delta_{\text{t}}$ and the depairing factor $\Gamma$ of the tip. The size of the energy gap $\Delta_{\text{t}}\approx 1.35\operatorname{\,\mathrm{meV}}$ does not vary for bulk-like tips, whereas the depairing factor was found to vary between $\Gamma\approx 10$-$20\operatorname{\,\mathrm{\mu eV}}$ for different tips.

S2.2 Numerical deconvolution of the $\mathrm{d}I/\mathrm{d}V$ spectra

In order find the exact position of the YSR energies on the NbSe2 sample, all data was numerically deconvolved as described in the following (similar to the procedure described in Ref. Choi et al. (2017)). Equation (S18) can be discretized into matrix form and then reads:

[TABLE]

The vector on the left hand side contains the differential conductance data, $\overrightarrow{\rho_{\text{s}}}(E)$ is the local density of states of the sample. The matrix $\mathbf{K}$ can be determined by comparing Eq. (S18) and Eq. (S20) to

[TABLE]

Thus, with the knowledge of the depairing factor and the energy gap of the tip, $\rho_{\text{s}}\left(E\right)$ can be calculated by finding the pseudoinverse of $\mathbf{K}$. In order to assure that the tip properties did not change during sample exchange from Pb to NbSe2, we checked the accurateness of the tip’s energy gap by deconvolution and recalculation of spectra on several Fe adatoms. If an inaccurate value of the energy gap is assumed for the tip, the energies of the thermally excited YSR resonances are found in a wrong position in the recalculated trace. Before deconvolving all spectra are normalized and smoothed using a Savitzky-Golay filter. Deconvolved spectra of atoms I-VI (cf. Fig. 2 in the main text) are shown in Fig. S6.

S2.3 Identification of adsorption sites

A dilute coverage of Fe atoms ($\approx 50\,\mathrm{atoms}/(100\times 100\operatorname{\,\mathrm{nm}}^{2})$) was deposited by evaporation on the NbSe2 sample in the STM at temperatures below $12$ K and a base pressure of $1\cdot 10^{-10}\operatorname{\,\mathrm{mbar}}$. Two distinct apparent heights of the Fe atoms are observed (Fig. S7 and Fig. 1c,d in the main text). Furthermore, the two adatom types differ in the $d$-level resonances (see blue and orange data in Fig. S8a). Here, only minor differences arise within different atoms of one species.

To identify the precise adsorption site of the two species, we recorded atomic-resolution images as shown in Fig. S8c,d. With the atomic lattice superposed on the STM images we find that all Fe atoms reside in hollow sites of the terminating Se layer. The Fe atoms with large apparent height solely occupy hollow sites, where the neighboring Se atoms form a triangle pointing in one direction, while the atoms with low apparent height are enclosed by a Se triangle rotated by 180∘. Inspection of the atomic structure of NbSe2 shows that there are two distinct hollow sites which differ by the presence or absence of a Nb atom beneath (metal-centered, MC, or hollow-centered, HC).

As discussed in the main text, the CDW is incommensurate with the lattice ($a_{\text{cdw}}\simeq 3.05a$ Soumyanarayanan et al. (2013)) which implies that the CDW smoothly shifts along the atomic lattice of the surface (see Fig. 1a of the main text and Fig. S11). Hence, the maxima of the CDW can be found at different sites with respect to the atomic lattice. Some striking motifs occur when the CDW maximum coincides with a hollow site (HC) or with a chalcogen site (CC) as illustrated in Fig. 1a of the main manuscript. Interestingly, the configuration with a CDW maximum on a Nb atom (MC) is not observed Gye et al. (2019); Guster et al. (2019). This local configuration is energetically unfavored Gye et al. (2019); Zheng et al. (2018); Guster et al. (2019); Lian et al. (2018); Cossu et al. (2018). To avoid the MC structure smooth grain boundaries are built into the CDW modulation. As illustrated in Fig. S9a,b Fe atoms with small apparent height can be found exactly in the minimum and maximum of the CDW. Hence, this species is assigned to adsorption in the HC position. Consequently, atoms with large apparent height are adsorbed in the MC sites.

S2.4 Determination of Fe-atom adsorption sites relative to the charge-density wave

The adsorption site of Fe atoms relative to the CDW was determined for a set of approximately 90 atoms (used for the data compilation in Fig. 4a of the main manuscript). The center of the Fe atoms as well as the location of the CDW can be determined from the STM images (see superposed CDW grid on the STM topographies in Fig. S10). The error in placing the CDW grids is estimated as $\delta x=\pm\,0.05\,\text{a}_{\text{cdw}}$, which accounts for the inaccuracy in placing the CDW grid, possible drift effects and for the inaccuracy in finding the center of the atom.

S2.5 Pinning of the charge-density wave

It has been shown that the CDW can be pinned by crystal defects Cossu et al. (2018); Langer et al. (2014); Flicker and van Wezel (2015); Chatterjee et al. (2015); Arguello et al. (2014). Therefore, the phase of the CDW is dictated by the complex interplay of lattice defects and long-range order of the CDW. One such example of a local modification of the CDW in the presence of a defect can be seen in Fig. S11a in the bottom left, where the CC structure of the CDW is pinned around a defect (visible as depression). Hence, one may argue that the Fe atoms can affect the CDW. However, a pinning of the CDW by the adatoms is not observed in the experiment. The adatoms seem not to favor a specific CDW domain as can be seen in Fig. S11a and b. Here, pure HC and CC domains are marked. The adatoms are located in the smooth transition region.

S2.6 Spatial variation of $\mathrm{d}I/\mathrm{d}V$ spectra

Figure S12 shows false-color plots of the $\mathrm{d}I/\mathrm{d}V$ spectra taken along high-symmetry directions of the YSR patterns across atoms I and IV of the main text. The YSR states exhibit long-range intensity oscillations with no detectable variations in energy within our experimental resolution.

We note that in close vicinity to the atom (Figure S13), the intensity variations are particularly strong and some YSR states are accompanied by a negative value of the differential conductance (NDC) at their high-bias flank. The NDC results from the convolution of YSR states with the tip density of states and can lead to a strong suppression of conductance. In an extreme case, an additional (weak) YSR resonance may be hidden in this part of the spectrum. Examples, where YSR resonances are suppressed below zero are shown by arrows for atom IV, V and VI in Fig. S13. The local variation of the NDC also affects the spatial appearance of the YSR states in the $\mathrm{d}I/\mathrm{d}V$ maps as is discussed in the following section.

S2.7 Additional $\mathrm{d}I/\mathrm{d}V$ maps of YSR resonances (both bias polarities and thermal excitations)

Figure S14 shows the complete set of the $\mathrm{d}I/\mathrm{d}V$ maps of atoms I-VI (the positive bias voltage was presented in the main manuscript). $\mathrm{d}I/\mathrm{d}V$ maps at negative bias polarity are depicted in the lower two rows. Opposite bias polarities image the $\left|u\right|^{2}$ and $\left|v\right|^{2}$ components of the YSR states, respectively Rusinov (1969).

As expected, the $\left|u\right|^{2}$ and $\left|v\right|^{2}$ components lead to distinct scattering patterns. As discussed above, the convolution of tip and sample density of states may lead to distortions of the intensity of individual YSR states, culminating in their suppression in the NDC region of a nearby state. Consequently, also the $\mathrm{d}I/\mathrm{d}V$ maps at the energy of a specific YSR state may be affected by a close-lying YSR state.

In particular, the $\pm\beta$-states of atom II, IV and V are affected by the NDC of the $\alpha$-resonances, which are close in energy (Fig. S14). One trick to determine and eliminate the effect of NDC is to investigate the $\mathrm{d}I/\mathrm{d}V$ maps of thermally excited YSR states. Thermal excitation of quasiparticles leads to additional resonances within the energy gap of the tip (cf. grey shaded area in Fig. 2b,c in the main text). Whereas the original resonances are found at a bias of $eV_{\pm\alpha,\pm\beta}=\pm|\Delta_{\text{t}}+E_{\alpha,\beta}|$, the thermally excited states are found at $eV_{\pm\alpha^{\star},\pm\beta^{\star}}=\mp|\Delta_{\text{t}}-E_{\alpha,\beta}|$ Ruby et al. (2015a). The thermally excited YSR states are of much less intensity and, therefore, exhibit a small or negligible region of NDC. Figure S15 shows the $\mathrm{d}I/\mathrm{d}V$ maps recorded at the corresponding voltages of the thermally excited YSR states for atoms I-VI. In these maps the similarity of the $\pm\beta^{\star}$ states of atoms I, II and IV is more striking than in the original maps of $\pm\beta$. To illustrate the impact of the NDC from the $\pm\alpha$-YSR resonances on the $\pm\beta$ state, we simulate the map of $\pm\beta$ by subtracting a fraction of the $\pm\alpha$ maps (mimicking the NDC) from the $\pm\beta^{\star}$ maps. The result is compared to the maps at the energies corresponding to $\pm\beta$ (Fig. S16) and shows remarkable similarity, illustrating the impact of NDC.

S2.8 Determination of YSR energies

To determine the YSR energies, we fit the deconvolved spectra (cf. Fig. S6) in the low-energy region with a set of 4 Gaussians (2 pairs of peaks):

[TABLE]

Here, $D_{0}$ is an offset, $E_{\alpha}$ and $E_{\beta}$ are the YSR energies, $A_{1,2}$ and $B_{1,2}$ are four amplitudes and $\sigma$ is the width of the Gaussian peaks.

The top panel of Fig. S17c shows the extracted energies of the YSR states of approximately 90 atoms with respect to their location relative to the CDW. The error bar in x-direction results from the inaccuracy in finding the exact position relative to the CDW as discussed along with Fig. S10. The error margins of the YSR energy include the uncertainty in the exact value of the tip’s energy gap $\Delta_{\text{t}}$, the lock-in modulation of the bias voltage and the standard deviation of the fit routine. Additionally, an error of $\pm\sigma/2$ was added to the energy uncertainty if the peaks are not well separated (i.e. $E_{\alpha}-E_{\beta}\leq 2\sigma$). This is mainly the case for atoms close to the CDW minimum, where we cannot distinguish the $\alpha$- and $\beta$-resonances as they are both very close to zero energy (e.g. atom III in the main text). Some of these atoms can be fitted with only one pair of Gauss peaks as the YSR states overlap. In this case both energies have the same value ($E_{\alpha}=E_{\beta}$).

As discussed in the main manuscript, we find a correlation of the YSR energy with position along the CDW symmetry axis. $E_{\alpha}$ and $E_{\beta}$ are largest at the maximum of the CDW and decrease with distance from the maximum. By default, the fit outputs for $E_{\alpha}$ and $E_{\beta}$ are positive. Assuming that $E_{\alpha}$ and $E_{\beta}$ should follow a similar trend, we suggest that the $\alpha$-resonance crosses zero energy, i.e. undergoes the quantum phase transition Franke et al. (2011). Therefore, $E_{\alpha}<0$ was assumed for the atoms in close vicinity to the CDW minimum (cf. atom III, the set of atoms with $E_{\alpha}<0$ is marked in orange and grey). Having passed the minimum (cf. atom IV) the $\alpha$- and $\beta$-states can be separated and identified again in the $\mathrm{d}I/\mathrm{d}V$ maps (see Fig. 3 in the main manuscript and Figs. S14, S15, S16), i.e. $E_{\alpha}>0$.

The top panel of Fig. 4a in the main manuscript is obtained from the data in Fig. S17c by averaging over all atoms that are found within an interval of $\pm 0.05$ a${}_{\text{cdw}}$. This corresponds to the uncertainty in the determination of the atoms’ position with respect to the CDW (see above). The error margins in energy are the standard deviations of the averaging including the average of the errors shown in Fig. S17c.

As seen in Fig. S17 and Fig. 4 of the main text, some regions along the CDW cannot be probed by Fe adatoms (grey areas). To understand the absence of data in these regions, we illustrate the phase change of the CDW (black rhombi) with respect to the atomic lattice (grey lattice) along the main symmetry axis in detail in Fig. S18. A HC adsorption site of an Fe atom is marked by a colored triangle on the atomic lattice. The green triangles indicate observed lattice sites, whereas the red triangles indicate that we did not find these adsorption sites relative to the CDW. The latter cases correspond to configurations, where the maximum of the CDW is located on top of a Nb atom. These so-called MC configurations are energetically disfavored as shown by DFT calculations and do not occur on the surface Gye et al. (2019); Zheng et al. (2018); Guster et al. (2019); Lian et al. (2018); Cossu et al. (2018). Hence, the corresponding Fe adsorption sites cannot be found either. This explains the absence of data points in the grey areas with only one exception. The third high symmetry position, where an atomic HC adsorption site would sit exactly in the local minimum of the CDW (marked with the red arrow in Fig. S18b), was also not found in experiment. We may speculate that the energy landscape of the CDW can be affected by the Fe adatoms for this specific adsorption site.