Scanning Tunneling Spectroscopy of Subsurface Ag and Ge Impurities in Copper

We investigate single Ge and Ag impurities buried below a Cu(100) surface using low temperature scanning tunneling microscopy. The interference patterns in the local density of states are surface scattering signatures of the bulk impurities, which result from 3D Friedel oscillations and the electron focusing effect. Comparing the isoelectronic d scatterer Ag and the sp scatterer Ge allows to distinguish contributions from impurity scattering and the host. Energy-independent effective scattering phase shifts are extracted using a plane wave tight-binding model and reveal similar values for both species. A comparison with ab-initio calculations suggests incoherent sp scattering processes at the Ge impurity. As both scatterers are spectrally homogeneous, scanning tunneling spectroscopy of the interference patterns yields real-space signatures of the bulk electronic structure. We find a kink around zero bias for both species that we assign to a renormalization of the band structure due to many-body effects, which can be described with a Debye self-energy and a surprisingly high electron-phonon coupling parameter $\lambda$. We propose that this might originate from bulk propagation in the vicinity of the surface.


I. INTRODUCTION
Scattering at impurities in a solid crucially affects material properties like electrical or thermal conductivity. Almost a century ago, the empirical Linde rule already related the residual resistivity of a dilute alloy ρ with the valence difference ∆Z between host and impurity species by ρ ∝ (∆Z) 2 [1]. As compound in a dilute Cu alloy, e.g., Ge shows a 27 times higher residual resistivity than Ag [1][2][3]. Complementary to macroscopic transport measurements, the scanning tunneling microscope (STM) yields access to the atomic scale that allows to investigate the scattering properties of single impurity atoms. This was first performed in the early days of STM in 2D scattering experiments of noble metal surface states, e.g. on Cu(111) or Au(111) [4][5][6]. The resulting Friedel oscillations contain information about both the defect scattering properties as well as the host. For example, on the one hand, the STM characterized a surface magnetic impurity with its Kondo signature and its scattering phase shift in the scattering pattern of a metal surface state [7][8][9][10][11][12]. On the other hand, by measuring the energy-dependent wave length of the scattering pattern at step edges, the dispersions of the surface states of Cu(111) and Au(111) were observed in real space [4,6]. The Fourier analysis of standing wave patterns, termed as quasiparticle interference (QPI), characterizes a materials' surface 2D band structure as well as available scat- * martin.wenderoth@uni-goettingen.de tering channels [13][14][15][16]. In specific sample systems, 2D QPI can be used to study parts of the 3D band structure [17,18].
Of peculiar interest are properties exceeding a single electron dispersion that describe many-body effects, e.g. electronic interactions with phonons where self-energy corrections cause a characteristic kink around the Fermi level. Especially metal surface states as prototypical 2D system were investigated in detail [19][20][21], e.g. by angleresolved photoemission spectroscopy (ARPES) [22][23][24][25]. By means of STM, the life times [26,27] as well as the electron-phonon coupling parameter λ and the selfenergy Σ [28] of surface states were determined. QPI also gave access to other many-body effects including coupling to surface plasmons [29] and different phonon modes [30][31][32]. As another approach by STM, Landau level spectroscopy offered access to 2D sample systems and their renormalized bands due to electron-phonon coupling [33,34].
In 3D, first, subsurface properties in metals were explored through quantum well states [35][36][37]. Then electron focusing was found to be a powerful tool to access bulk properties with an STM in real space on the atomic scale. Due to the anisotropy of a metal's band structure, the electrons propagate directionally through the bulk and form interference patterns at the surface. This was found at single magnetic impurities characterizing the host's Fermi surface [38], at noble gas cavities [39][40][41][42] and point-like subsurface defects [43]. The electron focusing effect was theoretically described [44] and applied at impurities in different contexts for the spectral signature of the Kondo effect at single atoms [45][46][47], for the two-band superconducting gap of Pb [48], and for quantum well states within layered systems [49].
In this study, we map the surface signatures of nonmagnetic impurity atoms buried in Cu(100). The scattering patterns of the isoelectronic d scatterer Ag and the sp scatterer Ge are compared in order to distinguish contributions from the impurity and the host. We extract similar energy-independent effective scattering phase shifts for both impurity species. Scanning tunneling spectroscopy (STS) reveals the energy-dependence of the electron focusing patterns, which act as real-space signature of the isoenergy surfaces of the host's 3D band structure due to the spectral homogeneity of the scattering phase. Most of the detected patterns is reproduced by electronic structure simulations, except from spectral anomalies intriguingly located at the Fermi energy. These can be described as the signature of a Debye self-energy that is discussed as bulk electron-phonon coupling in the vicinity of the surface.

A. Topographic signature of subsurface impurities
The experiments were performed with home-built STMs operating at a temperature of T ≈ 6 K and a base pressure of p < 5 × 10 −11 mbar. We investigate dilute Cu surface alloys with < 1% of Ge or Ag impurities, respectively. In a first step of sample preparation, Cu(100) single crystals are cleaned by cycles of Ar + sputtering and annealing. Subsequently, Cu and the impurity material are simultaneously deposited from two electron beam evaporators. The sample preparation is described in detail in Supplementary Note S1.
Both species of non-magnetic impurities feature picometer-high Friedel oscillations as surface signature due to electron focusing. A large scale topography in Fig. 1a shows interference patterns of Ge impurities at the Cu(100) surface. Each pattern corresponds to a single subsurface impurity causing a standing wave pattern of scattered electrons that propagate to the surface. The surface signatures have four-fold symmetry as defined by the crystal surface orientation and the anisotropic Cu band structure. The latter determines a specific angle under which the electrons propagate from the impurity towards the surface. Thus, a laterally more extended surface pattern corresponds to an impurity located deeper below the surface. Strong bright contrast are surface adatoms while dark contrast is assigned to surface layer impurities (termed '1st layer impurities' hereafter). The topographic contrast in the local density of states (LDOS) for subsurface impurities is with a few picometers a small fraction of a monoatomic step. (For more information see Supplementary Note S1).
In Fig. 1b, we present the surface signatures of Ag and Ge impurities from 3rd to 7th monolayer (ML) in Cu(100). The topographic contrast of a single impurity can be identified and assigned a depth using atomic resolution data (c.f. Supplementary Note S2). With increasing depth, the picometer-high corrugations increase in lateral size. Each four-fold pattern has its characteristic shape of maxima and minima due to constructive and destructive interference. Comparing both species, we find the topographic signatures to be very similar in both shape and amplitude. We can reproduce the experimental data with a tight-binding plane wave model adapted from Ref. [38] which uses an effective scattering phase shift η to characterize the shape of the interference pattern (for detailed information, see Appendix A). The average best-fit effective scattering phase shifts of η Ag = (1.32 ± 0.04)π and η Ge = (1.23 ± 0.05)π match the similarity of both species as observed in topography data.

B. Spectroscopic signature of subsurface impurities
We resolve the energy-dependence of the electron focusing patterns by using scanning tunneling spectroscopy. For each voltage, STS probes the interference pattern of the subset of states with wave vectors of the corresponding energy. Data was processed with topography normalization to obtain spectroscopic information for a constant height contour. For analyzing the change in LDOS due to the impurity, we subtract averaged spectra of the pristine Cu surface far away from the buried atom from the differential conductance around the impurity. As the topographic patterns only show a height contrast of a few picometers, the remaining spectroscopic signal shown in the following corresponds to approximately 1% of LDOS modulation obtained in the raw data.
We find the spectroscopic signatures of buried impurities to be very similar for both species, Ag and Ge, which agrees with the findings for the topography data. In Fig. 2 we show a 7th ML Ag impurity along with its spectroscopic signature. In the topography, the four-fold symmetry is well-resolved. A spectroscopic line section through the pattern along the crystal's [010] axis is presented in Fig. 2c. The color-code depicts the difference in differential conductance ∆dI/dU with respect to the pristine Cu surface. The main characteristic is an inwards movement of the positions of maxima and minima due to energy dispersion. A higher electron energy leads to a shorter wave length which fulfills the same interference condition for shorter distances. Dispersion is also the dominating effect in simulated spectroscopy data, calculated using the tight-binding model (Fig. 2b). The general shape of the simulation, which is given by the positions of maxima and minima lines, matches well with the experiment (Fig. 2c).
A spectroscopy data set of a 5th ML Ge impurity is shown in Fig. 3 with topography, spectral section and tight-binding simulation. As for Ag, we find the spectroscopic signature of a Ge impurity to be governed by electron dispersion with good agreement between experimental and calculated patterns. A single effective scattering phase can describe the data in the whole energy range without any strong additional resonances at element-specific energies. As for the topographies, the spectroscopy data for Ag and Ge is very similar for the same impurity depths.
An additional spectroscopic feature is found for both species at the Fermi level. In Fig. 2c and Fig. 3c, a bending of the interference pattern can be seen around 0 eV bias voltage where the pattern shifts towards the center in an immediate jump as if it experiencing a small phase shift. Considering a larger energy interval, the dispersion is unaffected by this anomaly. We observed these characteristics in multiple data sets for both species and different depths, including a Ge impurity buried in the 18th ML (c.f. Supplementary Note S4). The additional feature is not found in simulations of the tight-binding model, so it is beyond the scope of this single-particle band structure calculation.

III. AB-INITIO CALCULATIONS
Complementary to the experiment and the simplified tight-binding model, we performed ab-initio calculations using density functional theory (DFT) as implemented in the full potential Korringa-Kohn-Rostoker Green function (FP-KKR-GF) method [50] with the local density approximation (LDA) as parameterized by Vosko, Wilk and Nusair [51]. The theoretical and numerical approach is similar to that used in Refs. [38,44,47,49] and is described in more detail in Supplementary Note S3.
In Fig. 4, we show surface spectral sections along the crystal symmetry axis as calculated from the KKR-DFT potentials. Figure 4a and Fig. 4b show a 7th ML Ag impurity and a 5th ML Ge impurity, corresponding to the experimental sections presented in Fig. 2c and Fig. 3c.
The calculated spectra show richer patterns than the tight-binding model because the full crystal potential is considered. The pattern shape shows more detailed features than an almost linear dispersion and the LDOS amplitude modulates with energy. The total amplitude of the signal for Ge is multiple times higher than for Ag.
Comparing ab-initio calculation and experiment, the Ag data shows good qualitative agreement reproducing the overall pattern shape. With only small corrections to the ab-initio calculations' parameters (e. g. a shift of the energy scale by 150 mV or slightly corrected scattering phase shifts), the positions of the pattern's minima and maxima for the Fermi level match also quantitatively (c.f. Supplementary Note S3). For Ge, the patterns differ from each other. The calculation's rich features in the center including the contrast inversion from positive to negative with increasing energy are not found in the experiment, which presents a pattern that does not change qualitatively with energy and only shows a small dispersive shift. The pattern shape of the ab-initio data in the outer region (r > 8Å) resembles better the overall experimental pattern shape and is also similar to the Ag data. This different behavior for Ag and Ge data holds also for spectroscopic sections of other depths of both species.

A. Contributions of orbital components
Surface interference patterns are detected by STM even for atoms buried many layers below the surface despite the fact that non-magnetic impurities are weak scatterers. By the specific design of this experiment, we probe only the contributions of coherently scattered electrons for which the incoming and the scattered waves interfere constructively or destructively at the surface. The shape and size of a pattern are determined by two main ingredients: firstly, the scattering event, i.e., the scattering properties of the impurity atom and the electronic states it interacts with and, secondly, the electron propagation between surface and impurity.
In the experiment, we find the topographic signature of both impurity species to be very similar. From the literature, one would expect to find severely different scattering characteristics as the species differ strongly in both residual resistivity as well as the predicted scattering channels. In scattering theory, the scattering process at an impurity is described by its influence on host states being scattered. These electron wave functions are expanded in spherical harmonics, usually limited to s, p, and d orbitals. The influence on their asymptotic behavior by the impurity potential is parametrized by orbital-dependent scattering phases which are the main contribution to the scattering matrix T . These considerations are essential to the DFT-KKR ab-initio calculations, where the difference of element-specific scattering matrices determines by its amplitude how strong the scattering is in the respective channel for a specific impurity in a specific host. In other words, defects scatter especially in those channels, where the impurity differs most from a host atom. In this framework, strong scattering means that the host's states of a specific orbital character are especially disturbed and it does not describe scattering during bulk state propagation in the pure host or how many states of the respective orbital character are present. For the transition metal Ag in the isoelectric Cu, scattering at a 3d virtual bound state gives rise to the low residual resistivity [52]. As the resonance state is located far from the Fermi energy [53], leading to only little interaction with the electrons at the Fermi level, a small scattering amplitude of d scattering is predicted. For the sp impurity Ge in Cu, the perturbations of the crystal potential due to the significant additional local valence charge [53] of three additional electrons in s and p orbitals constitutes the main scattering process, leading to higher scattering rates especially for electrons of s and p character.
These predicted scattering intensities agree with the scattering amplitudes calculated in our ab-initio calculations (see Supplementary Note S3). While the scattering amplitude for Ge in Cu is three to four times higher for s and p orbital characters than for the d contribution, all Ag scattering amplitudes are half or less of the smallest value for Ge. This is reflected in the calculated LDOS amplitudes (see Fig. 4) where for Ge the signal is up to a factor of 7 higher than for Ag in the same depth. Even more strikingly, also the calculated pattern shape strongly differs between both species (c.f. Supplementary Figure S6).
While ab-initio calculations and experiment match for Ag impurities, there are large deviations for Ge. In order to explain these differences, the importance of the orbitals for the scattering process leads us to investigate if LDOS features of the simulations can be assigned to specific orbital components. For this, it has to be noted that the ab-initio calculations contain two orbital-dependent components, namely the impurity scattering as well as the propagation of Cu bulk states within the crystal lattice. For the latter, we can apply the tight-binding model (c.f. Appendix A). We determine the orbital character of the host metal's Bloch states in order to analyze how different orbital contributions are distributed over the Fermi surface and how this affects the electron focusing patterns at the crystal surface.
In Fig. 5a the contributions of Cu sp and d states in an exemplary interference pattern are shown as line sections as calculated with the tight-binding model. For each orbital, both propagation paths to and from the impurity are described by the orbital-specific propagator. For equal scattering amplitudes for sp and d states, the surface signal of d electrons is around 60% of the whole signal and, therefore, larger in amplitude than the combined sp contributions. At the outmost angles, the maximum at 12Å is only found in the d signal. This is representative for the signatures of impurities in all depths. In the center of the pattern, contributions from all electronic orbitals superimpose, while in the most outward region, the surface pattern is strongly dominated by d contributions.
In more detail, we can link this observation on a simulated topography to the host's Fermi surface: the interference patterns are constituted by bulk electrons that propagate from the impurity to the crystal surface in the direction of their group velocity. Electron focusing describes the fact that flat regions of the Fermi surface accumulate states with a certain propagation angle and, thus, cause a strong surface LDOS signal [38,44]. Figure 5b shows that for Cu(100) the focusing regions are especially the flat facets near the 110 directions (between the necks) with a focusing angle of 30 • to 40 • . The regions in 100 direction and around the necks ( 111 ) show smaller focusing angles. From the tight-binding model, we determine the orbital character of states at different regions of the Fermi surface as shown in Fig. 5c in agreement with previous calculations [52,54]. For all areas contributions from all orbitals are present with varying intensity. We find that the regions dominated by d character (bright) coincide with the flat focusing regions, whereas the neck and belly regions show stronger contributions from s and p (dark) whose smaller focusing angle points more towards the center of the interference pattern. Consequently, the orbital-resolved characteristics of the Fermi surface support what has already been visualized in the simulated orbital-resolved interference pattern.
With this knowledge, we can interpret the differences in the data of buried Ge atoms. In the ab-initio calculations, the outer, weak-contrast dispersion-like features are assigned to electrons with d character. In the center area, a very dominant feature from sp (because of strong sp scattering at the Ge impurity) masks a weak d electron background. As in the experiment such a sp feature is not visible, the signal of the experimental Ge topography is attributed to d electrons, which is supported by DFT in the outer regions. As our experiment is only sensitive to coherently scattered electrons interfering at the surface, we conclude that Ge shows such coherent scattering only for Cu electrons with d character and attribute the missing sp signal to additional incoherent processes in the sample system, which are not being considered in the ab-initio calculation.
We propose that these incoherent scattering processes occur at the impurity while we assume incoherence introduced during electron propagation to play only a minor role. Previous ab-initio studies have found strongly kdependent relaxation times for bulk electrons regarding electron-phonon interactions for different regions of the Fermi surface [54] with shorter times at necks and in 100 directions and longer times at flat regions which could explain stronger damping for the low angle regions than for high angles. However, as the relaxation times for room temperature range from 18 fs to 54 fs for different parts of the Fermi surface, the corresponding mean free paths are up to an order of magnitude larger than the propagation distance between surface and impurity in our low-temperature experiment. At low temperatures and including electron-electron interaction, bulk electron coherence has been assumed to be several dozens of nanometer even for energies up to 1 V from the Fermi energy [55,56]. Therefore, the coherence lengths are too large for the sp contributions to be damped by scattering within the propagation distance of only few nanome- ters. Instead, we suggest that the sp scattering process at the Ge impurity is partially not coherent which would strongly reduce the signal in our experiment. Within these considerations, the substantially increased residual resistivity of Ge impurities in Cu with respect to Ag could be attributed to decoherent scattering of sp electrons at the impurity.
Having discussed the origin of the differences between ab-initio calculations and experimental data for Ge, it remains the question why Ag and Ge show almost identical surface signatures in the experiments. For Ag, our DFT calculations and literature describe sp scattering to be less dominant than d scattering and a similar effective phase shift for all orbital contributions. This agrees with the simulated LDOS data that does not exhibit such a strong feature as Ge, but instead shows dispersive lines over the whole lateral range. These consist of all orbital characters and especially d electrons as these are the strongest orbital contribution in the electron focusing propagation. For Ag, the DFT calculation matches with the experimental data. When comparing both species, the KKR scattering phases for Ag are similar to the Ge scattering phase for the d orbital which is, in this case, the only imaged orbital channel. Hence, this can explain why the experimentally acquired topographies resemble each other for both species.

B. Signatures of Many-Body Effects
So far the analysis has shown that the surface signatures of buried impurities are also determined to a large extent by the bulk electron properties. This allows to compare scattering at two weak scatterers with different properties, here Ag and Ge impurities, and to separate the contributions from the specific impurity and general properties of the host. Thus, we use the real-space spectroscopic data to characterize the electronic structure of the Cu host around the Fermi level.
Similarly to the analogous STS studies for 2D surface states, here, we perform an experiment for standing waves of 3D Cu bulk electrons reflected at sub-surface point defects. In previous studies for single bulk impurities in Cu [45,46], the Kondo signatures of buried magnetic atoms were investigated. As magnetic impurities induce a phase shift around zero bias due to the Kondo resonance, an analysis of features from the host with weaker impact was not possible.
Analysing the spectroscopic data of the non-magnetic atoms Ag and Ge in Fig. 2 and and Fig. 3, respectively, the main observation in the range from −300 mV to 300 mV is a reduction of the interference pattern's lateral size for increasing bias voltage. This can be directly attributed to the dispersion relation of Cu bulk electrons as visualized in the tight-binding model (c.f. Figs. 2b, 3b). The reduction of the wave length for higher energies, which leads to a smaller ring structure, is the real space mapping of increasing electron momentum for higher energy.
In addition to this expected observation, the experimental data shows significant deviations from the simple, almost-linear bulk dispersion around the Fermi level. A bending of the pattern shape around zero bias is detected for both species in an energy range of approximately ±30 mV. This energy scale is common for lattice vibrations and it is known for the dispersion in momentum space that electron-phonon interaction leads to a renormalization of electron band structure, often manifested as 'kink' in the corresponding energy interval around the Fermi level [57]. We propose that in our experiment we map the real space signature of electronic coupling to phonon modes of the host.
In many-body theory, electron-phonon coupling is described by electron-like quasiparticles that include the interaction with phonons in their dispersion relation by self-energy corrections. These are determined by the Eliashberg function α 2 F (ω) or the averaged electron-phonon coupling parameter λ [57] whose cou-pling strength is given on the one hand by the available electron and phonon states, i.e. the respective densities of states, and on the other hand by the matrix element linking the two. While the real part of the self-energy shifts the center energy for a given k-value, the imaginary part causes a broadening in the spectral function due to finite quasiparticle lifetimes. In order to estimate the effect of such a band structure renormalization on the real space signatures, we implement the real part of the self-energy into the tight-binding model. For the Cu phonons, we use the Debye model with an Eliashberg function α 2 F (ω) = λ(ω/ω D ) 2 for ω < ω D and zero elsewhere, and the Debye frequencyhω D =30 meV [58,59]. The influence of the imaginary part of the self-energy is investigated in a simple 1D toy model with linear dispersion. We find broadening especially for states with energies E ≥ ω D , but no indications for strong changes in λ or the overall pattern shape (see Appendix B). Figure 6 shows the results and a comparison to the experimental data of 7th ML Ag and 5th ML Ge impurities. The tight-binding simulations including electronphonon coupling reveal distinctive features, especially for |E| <hω D (Fig. 6b and 6e). The patterns exhibit strong deviations from the simulations without many-body effects ( Fig. 2b and 3b) around the Fermi level, which are the real space signature of the well-known kink in momentum space. Comparing calculated and experimental data, we find good agreement for both species. For example, for the Ge pattern (Fig. 6d), the maxima almost join at around 25 meV before bending outward again, which is reproduced in the calculation (Fig. 6e). In the simulations, discontinuities are present at ω D , which are linked to the used Debye model that introduces an abrupt cut in the phonon density of states. For a real Cu phonon spectrum, the self-energy shows a smoother course that leads to weaker features at ω D (c.f. Supplementary Note S6). Furthermore, the simulation assumes a temperature of T = 0 K, while the experimental data is smoothened by thermal broadening on the scale of a few millivolts, blurring possible features at the Debye energy. The comparison in Fig. 6c and Fig. 6f shows the shift of the spatial position of the pattern's maximum with energy, as labelled in the 2D data. Good agreement is obtained with almost constant values for |E| >0.1 eV and a steep step within [−ω D , ω D ].
To describe the experimental data, we use an electronphonon coupling parameter of λ = 5.0, which is very high compared to the Cu bulk value λ Cu = 0.15 or even materials with strong coupling like PbBi alloys (λ ≈ 3.0) [57,60]. A calculation similar to Fig. 6 with the bulk value λ Cu does not reproduce the experimental data (c.f. Supplementary Note S5). Although the tight-binding model does not include all interactions and therefore is not expected to quantify exactly the experimentally found phenomena, this large deviation needs to be discussed.
Instead of coupling with delocalized bulk phonons, the spectroscopic features could be linked to local bulk vi-brational modes induced by the impurity atom at the impurity site or extended over several nanometers around the defect. However, as Ag is significantly heavier than Ge, different local vibrational modes would be expected, but both species' signatures show a similar strength of deviation from the single-particle dispersion. Furthermore, also electronic resonances, e.g. originating from a hybridization of bulk states with the impurity, which generally cause a phase shift, would very unlikely reveal a resonance exactly at the Fermi level for both species.
In the literature, the Cu bulk coupling parameter λ shows a directional variance of a factor of 2-3, leading to values as λ(k) = 0.08...0.23 [54,61,62]. Although our model uses the simplifications of no k-dependence and particle-hole symmetry in the self-energy, this alone cannot account for a factor of 30 between λ = 5.0 and the bulk value. According to ab-initio calculations, the bulk Eliashberg function α 2 F (ω) generally follows the phonon density of states without unusually large coupling of particular parts of the phonon spectra to electronic states [63]. Yet, our experiment is very specific about electron k values, as every point (x, y) probes only a specific region of the isoenergy contour with the corresponding angle of group velocity, so the general statement might not hold for an electron focusing signal.
We propose that the vicinity of the surface enhances the electron-phonon coupling parameter for our experiment. Missing translational symmetry breaks the momentum conservation perpendicular to the surface which opens up additional scattering processes with more available (bulk and surface) phonon modes as well as more available final electronic states.
Electron-phonon coupling of metal surface states has extensively been studied by surface-sensitive techniques [20,21]. For various systems, Mo(110) [23], Cu(111) [22,24,64,65], Cu(110) [66] as well as the Pb(110) bulk state [24], values of λ similar to the respective λ bulk were obtained by ARPES [57]. In STM measurements, many-body effects due to electron-phonon coupling of the Ag(111) surface state and other 2D systems were quantified by the coupling parameter λ and the self-energy which were found to match the literature values from ARPES and the bulk [28,33,34].
Further studies have shown that electron-phonon coupling can be enhanced compared to an average parameter λ, because the coupling of specific electronic and phonon modes can show strong variations as well as an energy-dependence [67,68]. Different interactions within the electronic surface band, with bulk electrons as well as with bulk phonon states and surface modes contribute unequally to the lifetime of surface states [26,27,69], e.g. for Cu(111) with strong coupling to the surface Rayleigh mode and interaction with specific bulk phonons [65,70]. By extracting the Eliashberg function an enhanced λ was obtained for the Be(1010) surface due to coupling to low energy surface modes [25]. These findings were crucial to understand the very high value of λ ≈ 1 for the Be(0001) surface [13,[71][72][73][74][75], four times higher than the bulk value [57]. Strong variations of the electron-phonon coupling parameter are also discussed for quantum well systems, e.g. thin Ag films [76][77][78][79][80].
It becomes clear that for a full picture of electronphonon interaction, all available final electronic states, available phonon modes, and the matrix elements for the scattering processes have to be taken into account at a state-specific level. This can play a key role in our experiment because we investigate a complex 3D system, which is dominated by the bulk band structure and additionally influenced by the vicinity of the surface and the atomic defects. Although the origin of the high electronphonon coupling value of λ ≈ 5 remains unclear, such an enhanced phase space of possible interactions might be a reason why the simulations including a self-energy from many-body effects can describe the experimentally observed features.

V. CONCLUSIONS
In summary, we perform a 3D scattering experiment at single non-magnetic impurities in a metal, detecting surface signatures of impurities buried up to 17 monolayers below the Cu(100) surface. The interference patterns at the surface characterize both impurity scattering properties of Ag and Ge as well as the host electronic structure of Cu. Analysing the topographic patterns with a plane wave tight-binding model, we extract similar effective scattering phase shifts for both impurity species. Following a comparison with ab-initio calculations, we propose that incoherent sp scattering at Ge impurities is crucial for the main scattering channel of this dilute alloy. As a result, the interference patterns of Ge impurities resemble those of Ag.
Due to the spectral homogeneity of impurity scattering, the surface signatures act as real space probe of the host's electronic dispersion, which we are able to map in good agreement with electronic structure calculations. We resolve small corrections in the form of a distinct bending of the interference pattern around the Fermi level that can be described with a Debye self-energy. It visualizes in real space the quasiparticle band structure of renormalized Cu bulk bands including many-body interactions. The high value of the electron-phonon coupling parameter λ could result from state-specific interactions and the vicinity of the surface in the investigated bulk system.
The used experimental approach can be versatilely applied to other systems to investigate scattering properties of subsurface impurity atoms as well as electronic properties of the host material. While for magnetic impurities the Kondo signature masks the host's features at the Fermi level, non-magnetic impurities or inherent defects can be functionalized to probe the electronic structure of both bulk and surface nanostructures even close to zero bias. For a simple material such as a noble metal, we unveil non-trivial features that we propose to be emerging from phonons coupling to electrons. As the miniaturization of nanoelectronic devices continues, the presented approach can be a useful tool to understand fundamental processes on the atomic scale in order to identify electronic interactions.
Appendix A: Description of surface signatures with tight-binding model The surface signatures are analyzed and reproduced with a tight-binding plane wave (TB) model adapted from Ref. [38]. It describes the interference of an outgoing electron wave and an incoming wave scattered at the impurity (see Fig. 7). The Bloch states are assumed to be plane waves without considering the lattice periodic part, but including the dispersion of the Cu bulk band structure, as calculated in a tight-binding calculation [81,82]. For simplicity, the scattering event is reduced to an energy-independent, orbital-independent effective scattering phase shift η imp that the wavefunction collects at the impurity and a corresponding unitary scattering matrix T = exp(i · η imp ). The change in LDOS with respect to the pure Cu surface is given with the Green's function G 0 as propagator and x imp as the position by The LDOS(x, ) at the surface exponentially decays into the vacuum towards the STM tip. We use the effective scattering phase shift η imp of the TB model to describe and quantify the shape of the experimental surface signatures. It can take values between 0 and 2π and is assumed to be constant for different bulk impurity depths. This is reasonable as for bulk conditions the surroundings of an impurity are always identical, independent of location within the single crystal. The other free parameters of the model are the height of the tip h above the sample, the distance a surf between surface layer and vacuum, and the depth of the impurity d which only takes discrete integer values (Fig. 7).
The typical distances in scanning tunneling microscopy are about 5Å to 8Å where shorter distances lead to more laterally high-frequency features because of larger transmission of k || [83]. Taking the bias voltage and current set points into consideration, we choose h = 5Å for Ag and h = 7Å for Ge data sets which proves to reproduce well the general appearance of interference patterns (cf. Figs. 8 and 9). The tip-sample distance h does not have strong influence on the simulated topographies as a large change of h only leads to a change of roughly 0.1π in phase shift which is about the size of the error bar for the extracted phase shift.  7. Sketch of the tight-binding plane wave model. Surface interference patterns result from interference between incoming (blue) and outgoing (red) electron waves, which oscillate within the crystal and experience a scattering phase shift ηimp at the impurity (yellow). The tip probes the surface signature with a distance-dependent exponential decay. This distance is labeled with h. The depth of the impurity is given by d.
The surface layer is referred to as the 'first layer'.
The surface layer spacing is set to a surf = a [100] = 181 pm in a rough approximation of the crystal potential [84]. Small changes in depth can be compensated by a shift in effective phase shift η imp because the phase contribution of additional propagation distance (Fermi wave length λ Cu = 4.6Å) cancels with the new η imp . A shift of 0.1 ML in depth corresponds to approximately −0.14π in phase shift. Therefore, as the parameters of the TB model are not independent of each other, there is not exactly one resulting parameter set. However, the lateral size of the surface signature which has to match with the focusing directions confines this flexibility to very few monolayers. The depth of an impurity d is deduced from images with atomic resolution (for more information see Supplementary Note S2). The depth is restricted to discrete integer values as impurities only sit substitutionally in Cu atomic lattice points.
The effective phase shift in the TB model includes all phase contributions within the solid, even if they do not occur at the impurity. Firstly, this is a phase shift from the transmission from the vacuum to the metal crystal [85]. Secondly, unlike the phase shift in the KKR framework, the TB model assumes propagation within the impurity cell which leads to a phase contribution. These contributions are identical for both impurity species, so the phase shift difference between them is independent of the contributions' absolute values. In Fig. 8 the surface signatures of Ag impurities from 3rd to 7th ML in Cu(100) are shown. The left column depicts the experimental topographies already presented in Fig. 1b in the main text. With increasing depth, the picometer-high corrugations increase in lateral size. Each four-fold pattern has its characteristic shape of maxima and minima due to constructive and destructive interference. The center column (Fig. 8b) depicts simulated topographies calculated with the TB model with a best-fit scatter- ing phase, reproducing their experimental counterparts in size and shape. The sections (Fig. 8c) along the [010] direction also indicate accurate agreement between experimental and simulated data. The deviations for the 7th ML originate from the atomic resolution superposed with the bulk focusing signal. We find an average effective scattering phase shift of η Ag = (1.32 ± 0.04)π. Figure 9 shows the surface signatures of another nonmagnetic impurity, Ge. Figure 9a shows experimental topographies in monolayers 3-7. The simulated 2D TB topographies (Fig. 9b) and sections of experiment and calculation (Fig. 9c) show good agreement. We obtain an average effective scattering phase shift of η Ge = (1.23 ± 0.05)π. The tight-binding model proves to be a reliable tool to describe the surface patterns of subsurface, nonmagnetic impurities. This can be especially useful when ab-initio calculations are not available or for deeper layers where computational cost increases strongly for DFT.

Appendix B: Electron-phonon coupling toy model
In the implementation of electron-phonon coupling in the tight-binding model, as described in the main text, we have only included the real part of the self-energy, which is responsible for the shift in center energy in the spectral function with respect to the single-particle system. Here, we analyze the role of the imaginary part for the spectral function in real space using a simple 1D toy model. We approximate the quadratic dispersion of a free electron to be linear in the considered energy interval around the Fermi energy (E F =7 eV). The interaction with Debye phonons is introduced via a self-energy calculated from an Eliashberg function α 2 F similar as described in the main text.
The resulting renormalized spectral function in momentum space is shown in Fig. 10a for λ = 5.0. The typical kink at the Debye energyhω D =30 meV as well as the broadening due to quasiparticle interactions for energies E > 0 eV are found as known from the literature [57].
Via Fourier transform, we obtain the real space representation of the Green's function propagator and calculate a simulated STS pattern for an impurity in 7ML with effective scattering phase η = 1.0π. Because the Friedel oscillations' amplitude in the 1D model is constant even for large distances, we include an additional damping term ∼ r −3 (with distance r from impurity to surface position) into the LDOS. This accounts for a 3D propagation from impurity to the surface like in the experimental patterns.
For a propagator with only the real part of the selfenergy, as implemented in the tight-binding model, we obtain a pattern of constant intensity with energy with a characteristic kink around the Fermi level (Fig. 10b). As this toy model is spatially isotropic and, thus, does not favor a certain electron focusing angle, there are multiple oscillations in the surface pattern. The smaller amplitude with larger distance only results from the additional 3D damping term.
Including the full self-energy (Fig. 10c), the imaginary part reduces the signal strength due to the finite lifetime. For E > ω D , the amplitude is attenuated by approximately 50% while close to the Fermi level, the patterns remains almost unchanged. The whole pattern shape (i.e. the position of minima and maxima) is almost identical to the case with only e(Σ), as depicted in Fig. 10d where for each line the spectral function is normalized to 1.
From the study of the toy model, we conclude that our implementation of many-body effects into the Cu tight-binding model can describe the basic effects of electron-phonon coupling for the STS spectra. We do not claim exact quantitative accuracy for the tight-binding model. After all, despite describing well the overall topographic and spectroscopic signatures for different impurity species, it is a simple model that does not include all effects of the solid.

S1. SURFACES OF DILUTE CU ALLOYS
The samples, Cu(100) with dilute Cu alloys in the topmost layers, were prepared in a home-built UHV preparation chamber with a base pressure of p = 5 × 10 −11 mbar. The Cu(100) single crystals were cleaned by cycles of sputtering by Ar + ions with an energy of 700 eV and annealing to 680 K. In the last cycles the sputtering and annealing times were reduced and the temperature was lowered to 630 K in order to flatten the surface without causing segregation of bulk defects towards the surface. The surface quality was checked by low-energy electron diffraction, Auger electron spectroscopy and STM which presents the sample as clean with terraces of a width of up to a few 100 nm. Dilute alloys in the 20 topmost layers are fabricated by co-deposition of 20 ML copper and < 1% silver and germanium, respectively, by two electron beam evaporators. Cu is evaporated at 1 ML/min.
In order to achieve very small ratios of composition, the impurity evaporator is equipped with a stepper motor which allows shutter opening times of 100 ms. During evaporation, the sample is held at a temperature of T≈100 K. In order to restore a flat surface that can be investigated by STM, the sample is flashed to T = 520 K for 3 seconds after growth.
Subsequently, it is immediately transferred into the low temperature STM.
The surfaces of dilute Cu alloys show different patterns at the surface corresponding to different defects. An example of a topography with Ag impurities is shown in Fig. S1 which is similar to the large-scale topography for Ge presented in the main text (cf. Fig. 1a). Many second layer Ag impurities are visible as strong contrast white dots with surrounding black rings. In the left center two impurities are buried, one rather close to the surface and one deeper impurity with clear four-fold pattern. Additionally, signatures of clusters presenting richer patterns are found (e.g. center, bright white contrast) which can either be linked to subsurface impurity clusters or nanocavities of residual argon gas from the sputtering process.
The four-fold signature for single, buried impurities are the consequence of electron focusing. Electrons within the crystal propagate corresponding to their group velocity which is given by the gradient of band structure [1]. Hence, for flat regions on the isoenergy con- directions, leading to a weak signal at the surface. This effect is depicted for two sections of the Cu Fermi surface in Fig. S2. The electrons dominantly propagate in directions where group velocity vectors accumulate, leading to surface interference patterns above buried impurities. In semiconductors or other materials with comparably low charge carrier density, subsurface dopant atoms can be mapped by STM in a different mechanism. The donor's extended Coulomb potential around the impurity site can be modified by switching the dopant's charge state by tip-induced band bending [2,3]. This mechanism is excluded for the impurities investigated here, buried in the metal Cu. The charge carrier density is significantly higher and the Thomas-Fermi screening length as short as 0.55Å, so that any charged region is screened on short length scales. In Fig. S3 three different types of single impurities in Cu(100) are shown in topography and line sections. The signature of first monolayer impurities is a depression for Ge (and also previously, e.g. Co [4]) with an amplitude of typically 30 pm. Second ML defects, only covered by one monolayer of Cu, reveal a bright spot surrounded by a dark ring.
The electronic contrast is typically 10-15 pm. As a reference also a 7th ML impurity with picometer-high signature is included in the data set. Due to very good imaging conditions in the presented data set, the typical heights for the features are exceeded.

A. Determination of impurity depth
The depth of an impurity is deduced from STM topographies with atomic resolution.
The center position of the surface interference pattern with respect to the surface lattice determines if the impurity is located in an odd or even crystal monolayer. The surface pattern is always centered around the buried atom due to symmetry. If the surface signature coincides with a surface atom, then the impurity sits perpendicular to the surface below this surface atom. Hence, with Cu being an fcc single crystal, the impurity is located in an odd monolayer. Here, the surface layer is counted as first layer. If the surface signature's center is located between surface atoms, then the impurity is positioned in an even crystal monolayer. An example for such an analysis is shown in Fig. S4a. The surface interference pattern of this impurity is centered around a surface atom (marked with a cross), so the impurity is located in an odd crystal monolayer.
Additional information for the determination of an impurity atom's depth is taken from the lateral pattern size which is restricted by the focusing angle, i.e., the main directions of electron focusing. Because of this, for a full data set in which all depths are available, the patterns can be ordered by size and then assigned a crystal monolayer. Alternatively, tight-binding calculations (see Appendix A in the main text) can be performed for different depths, for each fitting the best effective scattering phase shift. As the impurity positions are limited to a discrete number of crystal layers and by atomic resolution data further reduced to only even or odd layers, the comparison of simulation and experimental data reveals the impurity's depth. For the exemplary depth analysis in Fig. S4, in (b) the experimental data is compared with a tight-binding simulation for a third monolayer ('3ML') impurity. The patterns match far better than the candidates of 5th ML and 7th ML in subfigure (c). For deeper layers, the main weight of the interference pattern is shifted towards larger distances, whereas for the 3rd ML the main signal is located in the center peak.
For data sets where no atomic resolution was obtained, the depth is determined by comparison with tight-binding simulations of various depth with the average effective scattering phase shift, assuming a constant scattering phase shift for all bulk layers. Once the layer is identified (e.g. comparison of lateral size, ratios in contrast pattern), in a next step an  [5][6][7]. The interference patterns for residual sputter gas (Ar, Ne) can be distinguished from the patterns used in our analysis as they show richer surface interference patterns with larger lateral extent (e.g. in Fig. S1).
For the dilute alloys, the concentration of impurities is usually < 1%, so that impurity atoms can hardly form clusters. In some topographies we find clusters of Ag or Ge with rich surface signatures which are easily distinguished from point-like scattering centers which show clearly defined four-fold interference patterns confined to approximately 1.5 Friedel oscillations [1].
Nevertheless, the tip quality is crucial for resolving the surface signatures of buried impurities. They show apparent heights of few picometers in topography which correspond to LDOS modulations of around only 1% in spectra. Data sets with reproducible, well-defined interference patterns as shown in the main text can be recorded with a suitable sharp tip.
For dull tips, we obtain broader features which make it more difficult to distinguish between subsurface structures and other defect contrasts, which is why we exclude these from the analysis.

S3. AB-INITIO CALCULATIONS AND SCATTERING PHASE SHIFTS
The ab-initio calculations are performed using DFT as implemented in the full potential Korringa-Kohn-Rostoker Green function (FP-KKR-GF) method [8] with the local density approximation LDA [9]. First, the electronic structure of 18 layers fcc Cu slab with two additional vacuum regions (3 vacuum layers on each side of the Cu slab) are taken along the (001) direction. The decimation technique [10,11] is adopted to simulate the semi-infinite substrate to avoid size artifacts in the charge density. The experimental lattice parameter a = 3.61Å was considered without surface relaxations, which are negligible. Then each impurity is embedded underneath the surface together with its neighboring atoms, defining a cluster of atoms, where the charge is allowed to be updated during self-consistency. The induced charge density is then computed in the vacuum at h = 3.61Å above the surface, proportional within the Tersoff-Hamman approach [12] to the tunneling signal measured with STS. While the self-consistent calculations required 40 × 40 k-points in the two-dimensional Brillouin zone, the theoretical STS spectra were obtained with a set of 200 × 200 k-points.
The calculations are performed to obtain simulations of the experimental spectroscopy data and energy-dependent scattering phase shifts at the different impurity species. As the former results are discussed in the main text, here we mainly focus on the phase shifts at the Fermi energy.
We obtain scattering phases δ for host and impurity species consistent with previous studies [13]. The values of phase shifts at the Fermi energy are listed in Tab. I. A comparison between 5th ML and 3rd ML shows only tiny changes for the different depth. This indicates a depth-independent scattering phase shift which supports the assumption in our analysis for extracting the phase shift from experimental data.
With the tight-binding model (c.f. Appendix A in the main text), we parametrize the experimental surface patterns with a single effective scattering phase shift η imp which describes the phase relation of incoming and scattered wave. This parameter is determined by the scattering potential of the impurity. The corresponding quantity used in quantum scattering theory is the scattering phase shift δ which is linked to the scattering matrix t by t = −(1/ √ E) sin(δ) exp(iδ). In a Cu crystal with Bloch electron states, the scattering amplitude and phase due to an impurity are given by the complex number ∆t = t imp − t Cu which is the difference of the impurity scattering matrix with respect to the host [14].  The experimentally determined effective scattering phase shift η imp corresponds to arg(∆t). The scattering amplitude, i.e., the amplitude of the complex scattering matrix difference, is contained in the height of electronic contrast measured by STM. However, electronic contrast is very difficult to quantify as transmission strongly depends on absolute tip-sample distance and tip shape. In the calculations, Ge shows as an impurity in Cu a three to four times higher scattering amplitude for s and p orbitals than the d orbital contribution. For Ag, the orbital components only vary by a factor of two, with all amplitudes being half of the Ge's value for d orbitals or less. The energy-dependence of scattering phases for s, p, and d orbital characters are shown in Fig. S5 for Cu, Ag, and Ge. The phase around the Fermi level is approximately constant.
The strongest energy-dependence from −1 eV to +1 eV show the s and p orbitals of Ge. The scattering amplitude, i.e., the absolute value of difference of the complex scattering matrices of Ge and Cu, of these orbitals is 3-4 times higher than for the d orbitals.
Resulting from the ab-initio calculations, we obtain simulated spectroscopic sections as shown in Fig. S6 which also comprises the data shown in the main text. While the lateral sizes of the surface signatures are similar for both species for a specific depth, Ag and Ge differ strongly in LDOS amplitude and pattern shape.
As discussed in the main text, for Ag we find good agreement between ab-initio calculations and the experimental data. In order to obtain a match for the overall shape and lateral size of the interference pattern, in the ab-initio calculations either the scattering phase shifts have to be slightly corrected or the energy scale has to be shifted by 150 mV. Such an energy shift is known from the literature to be possible for metallic electronic states as obtained from LDA with respect to those measured experimentally [15]. In Fig

S4. DEEPLY BURIED IMPURITIES
Due to the electron focusing effect, the signature of buried impurities is detected by an STM despite the screening effects of the host's high charge density. While in previous data we have shown interference patterns corresponding to impurities of up to the 7th ML, in Fig. S8 we show a deeply buried Ge impurity. In Fig. S8a, the topographic surface pattern of a 18th ML impurity is presented. The four-fold ring structure of more than 3 nm diameter is clearly discriminable despite the adatoms and other buried impurities in the surrounding surface region. Because the focusing pattern is limited to distinct 1.5 wave lengths and clusters would lead to richer surface structures, we assign the pattern to a single impurity.
We identify the impurity's depth by comparison with tight-binding simulations (Fig. S8b) where we obtain good agreement.
In Fig. S8c and Fig. S8d, the spectroscopic signature is shown for two energy ranges.
For ±300 mV (Fig. S8c), the dominant feature is the shrinking of the ring-like structures for higher energies. This is understood by electronic dispersion. For high absolute values of Spectroscopic section similar to (c) with smaller energetic range. Around the Fermi level, an additional bending of the pattern is found that is assigned to electron-phonon coupling.
energies, the signal is more difficult to detect as the tunnel current noise is also increasing.
For these long propagation distances, there might be additional effects by electron-electron scattering for energies further from the Fermi level. As the impurity is located 30Å below the crystal's surface, the electrons in the interference pattern have travelled coherently 7 nm through the Cu crystal. In Fig. S8d, we show a spectroscopic section for an energy range from −100 mV to +100 mV. As in the data before, we find an additional bending at the Fermi level due to electron-phonon coupling. Even for this depth, the renormalized band structure of interacting electrons is resolved in spectroscopy.

S6. SPECTROSCOPIC SIGNATURES AND PHONON SPECTRUM
In Fig. S10, a comparison of the spectroscopic signatures of a 7th ML Ag impurity and a 5th Ge impurity are shown for different phonon spectra. The tight-binding simulations including a self-energy based on the Debye model (Fig. S10a,d) are also shown in Fig. 6 in the manuscript. Here, they are compared with calculations using a self-energy based on the Cu phonon spectrum. One finds that, firstly, the Debye model approximation of an abrupt cut-off at the Debye energy leads to discontinuous line shape. Therefore, the self-energy and the simulated STS pattern for the real Cu phonon spectrum reveal a smoother transition around the Debye energy. Secondly, the simulation assumes a temperature of T = 0K. The experimental data is recorded at T = 6K, so the experimental data is thermally broadened, which is included in the spectroscopic sections shown in Fig. S10(c,f).