### Experimental procedures

Table of Contents

The experiments were performed in a commercially purchased SPECS STM system operated at *T* = 4.5 K, which is equipped with home-built ultrahigh-vacuum chambers for sample preparation^{46}. STM images were obtained by regulating the tunnelling current *I*_{stab} to a constant value with a feedback loop while applying a constant bias voltage *V*_{stab} across the tunnelling junction. For measurements of differential tunnelling conductance (d*I*/d*V*) spectra, the tip was stabilized at bias voltage *V*_{stab} and current *I*_{stab} as individually noted in the figure captions. In a next step, the feedback loop was switched off and the bias voltage was swept from −*V*_{stab} to +*V*_{stab}. The d*I*/d*V* signal was measured using standard lock-in techniques with a small modulation voltage *V*_{mod} (RMS) of frequency *f* = 1,097.1 Hz added to *V*_{stab}. The following measurement parameters have been used for the data presented in the main figures:

Fig. 1a: *V* = 50 mV, *I* = 1 nA, *V*_{mod} = 5 mV; Fig. 1c,d: *V* = 5 mV, *I* = 1 nA; Fig. 1e,f: *V*_{stab} = −100 mV, *I*_{stab} = 2 nA, *V*_{mod} = 2 mV; Fig. 2a,c: *V*_{stab} = −15 mV, *I*_{stab} = 4 nA, *V*_{mod} = 50 µV; Fig. 2b: *V*_{stab} = −15 mV, *I*_{stab} = 4 nA, *V*_{mod} = 100 µV.

d*I*/d*V* line profiles were acquired recording several d*I*/d*V* spectra along a one-dimensional line of lateral positions on the sample, respectively. Note that the tip was not stabilized again after each individual spectrum was acquired but the line profiles were measured in constant-height mode. This avoids artefacts stemming from a modulated stabilization height. At the chosen stabilization parameters, the contribution of multiple Andreev reflections and direct Cooper pair tunnelling to the superconducting tip can be neglected (see Supplementary Note 1). Throughout this work, we use Nb tips made from a mechanically cut and sharpened high-purity Nb wire. The tips were flashed in situ to about 1,500 K to remove residual contaminants or oxide layers. The use of superconducting tips increases the effective energy resolution of the experiment beyond the Fermi–Dirac limit^{47} but requires careful interpretation of the acquired d*I*/d*V* spectra. These are proportional to the convolution of the sample’s LDOS, the superconducting tip DOS and the difference in the Fermi–Dirac distributions of the tip and sample. Notably, the latter can play a large role when measuring at *T* = 4.5 K. We measure a value of *Δ*_{s} = 1.35 meV (Supplementary Note 1), which is similar to the gap of elemental Nb, *Δ*_{Nb} = 1.50 meV (refs. ^{3,48}), indicating a high interface quality between Nb and Ag. Details on the interpretation of SIS tunnelling spectra and on the determination of the tip’s superconducting gap *Δ*_{t} can be found in the next section, as well as in Supplementary Note 1.

### Sample preparation

A Nb(110) single crystal was used as a substrate and cleaned by high-temperature flashes to *T* ≈ 2,000 K. This preparation method yields the characteristic oxygen-reconstructed Nb surface observed in previous works^{49}, as can be seen in Extended Data Fig. 5a. Notably, a similar surface quality can be achieved by sputter-annealing cycles only, that is, without the need for challenging ultrahigh temperature flashes typically required for the preparation of clean c(1×1) Nb(110) surfaces^{48}. Ag was deposited from an e-beam evaporator using a high-purity rod at a deposition rate of about 0.1 monolayers (ML) min^{−1}. In agreement with previous studies, evaporation of Ag at elevated temperatures leads to the formation of two pseudomorphic monolayers of Ag followed by Stranski–Krastanov growth of large Ag(111) islands (Extended Data Fig. 5c). To obtain preferably thin islands, we grew Ag islands in a three-step process, starting with the deposition of 2ML at *T* ≈ 600 K, creating two closed wetting layers. In a second step, the temperature was reduced to *T* ≈ 400 K to limit the lateral diffusion of Ag on the surface and to create more nucleation centres for the Stranski–Krastanov islands. Under these conditions, another 2ML of Ag were deposited, followed by three further ML grown at *T*≈ 600 K, again to guarantee a well-annealed surface of the topmost layers.

A topographic image of a Ag island grown on NbO_{x}/Nb(110) is shown in Extended Data Fig. 5a. This sample features a Ag coverage of only 15%, enabling the identification of the substrate’s oxygen reconstruction (see refs. ^{48,49}) and of the islands’ apparent heights. Nearly all of such islands are found to have heights in the range 500–550 pm, consistent with a preferred double-layer growth. Low-energy d*I*/d*V* spectroscopy measurements (Extended Data Fig. 5b) reveal clean SIS tunnelling on both the NbO_{x}/Nb(110) substrate and the Ag double-layer island: sharp peaks of high differential tunnelling conductance appear at bias voltages *e**V* = ±(*Δ*_{t} + *Δ*_{s}), corresponding to quasiparticle tunnelling between the coherence peaks of the tip and sample, respectively. Also, weaker resonances are found at voltages *e**V* = ±(*Δ*_{t} − *Δ*_{s}). These are typically attributed to thermally activated tunnelling between the partially occupied and unoccupied coherence peaks of the tip and sample^{47}. From the positions of these peaks measured with different microtips, the tip and sample gaps can be unambiguously determined. Notably, there is no clear difference between the spectra measured on NbO_{x}/Nb(110) and on the Ag double layer, providing evidence that the interface quality between Nb and Ag is sufficient to open a full proximity gap in the Ag states.

As the Ag coverage is increased above 2ML, the Ag double layer is gradually closed and the formation of further large islands in the Stranski–Krastanov growth mode is observed (Extended Data Fig. 5c). For these samples, the closed double layer can be investigated in more detail. Characteristic defects of unknown chemical composition are found on the double layer, exhibiting a twofold symmetry (Extended Data Fig. 5d). This already suggests that the Ag film does not grow in a fcc(111) fashion for the first two layers. Instead, atomically resolved images of the double layer (Extended Data Fig. 5e) reveal a pseudomorphic growth of Ag on the bcc(110) surface of the underlying NbO_{x}/Nb(110). Previously, a similar growth mode has been reported for Ag layers on oxygen-reconstructed V(100) in ref. ^{50}: the pseudomorphic nature of the growth despite a reconstructed substrate surface has been explained by the enhanced mobility of oxygen in vanadium at elevated temperatures, leading to a substitution of oxygen atoms by Ag atoms during the growth and to a clean interface. We speculate that a similar growth mode is taking place for the first double-layer Ag/Nb(110). In contrast to this, atomically resolved images on the thicker islands (Extended Data Fig. 5f) reveal the hexagonal lattice expected for the energetically favoured fcc(111) growth of Ag and a very low number of impurities (typically 1–2 per 100 × 100 nm^{2}; see Extended Data Fig. 6a). These results are in agreement with earlier reports of the growth mode of Ag on Nb(110) by low-energy electron diffraction^{51} and STM^{15,52}.

On the islands, we find pronounced signatures of the well-known Shockley-type surface state of Ag(111), providing further evidence for a fcc(111) growth. In Extended Data Fig. 6b, an example of a constant-height d*I*/d*V* map measured in the area of Extended Data Fig. 6a is shown, measured at an energy above the onset of the surface state and outside the superconducting gap. It features a clear periodic modulation in agreement with pronounced quasiparticle interference of the surface-state electrons^{53,54}. A Fourier transformation of the map visualizes the circular Fermi surface of the surface state (inset of Extended Data Fig. 6b).

### Construction of QDs

As previously reported in ref. ^{55}, approaching the tip to a Ag(111) surface can lead to two processes: (1) single Ag atoms can be reproducibly pulled out of the surface by attractive tip–sample interactions, leaving a vacancy behind in the Ag lattice, and (2) single Ag atoms are dropped from a Ag-coated tip, which was coated previously by dipping the tip into the Ag(111) surface.

An example of this process of adatom gathering is shown in Extended Data Fig. 7a–e. Once the Ag atoms are moved to a region without surface-state contributions—for example, inside a QD structure (Extended Data Fig. 7f)—the d*I*/d*V* spectra on top of Ag atoms (Extended Data Fig. 7g) show clean SIS tunnelling without signs of Yu–Shiba–Rusinov states^{31}. This provides further evidence that the adatoms used for our QDs are indeed non-magnetic, as expected for Ag/Ag(111). Subsequently, after gathering a sufficient number of Ag atoms, the Ag QDs were constructed by lateral atom manipulation^{56,57} at low tunnelling resistances of *R* ≈ 100 kΩ. Because the Ag walls of the QDs have a finite transparency for the surface-state electrons, a second wall of Ag atoms is constructed around the central QD wall to screen the interior from surface-state modes located outside the structure.

### Modelling of QD eigenmodes

The wavefunctions of the eigenmodes of the QDs can be well modelled by hard-wall particle-in-a-box simulations to a first approximation, as it was done already in the pioneering work by Crommie et al.^{11} in 1993. The eigenmodes of an infinitely high rectangular potential wall with dimensions *L*_{x} and *L*_{y} are the well-known analytical solutions:

$$\varPsi ({n}_{x},{n}_{y})={\varPsi }_{x}({n}_{x})\times {\varPsi }_{y}({n}_{y})$$

(4)

with

$${\varPsi }_{j}({n}_{j})=\sqrt{2/({L}_{j}-\delta )}\times \{\begin{array}{c}\sin ({\rm{\pi }}{n}_{j}/({L}_{j}-\delta )\times j)\,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{e}}{\rm{v}}{\rm{e}}{\rm{n}}\,{n}_{j}\\ \cos ({\rm{\pi }}{n}_{j}/({L}_{j}-\delta )\times j)\,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{o}}{\rm{d}}{\rm{d}}\,{n}_{j}\end{array}$$

(5)

*j* = *x*, *y* and the quantum numbers *n*_{x} and *n*_{y}, corresponding to the number of antinodes of a certain eigenfunction. These correspond to eigenenergies of

$$E({n}_{x},{n}_{y})=\frac{{\hbar }^{2}}{2{m}_{{\rm{eff}}}}[{({\rm{\pi }}{n}_{x}/({L}_{x}-\delta ))}^{2}+{({\rm{\pi }}{n}_{y}/({L}_{y}-\delta ))}^{2}]+{E}_{0}\,.$$

(6)

Note that the parameter *δ* is introduced to renormalize the effective dimensions of the QDs, because the distances seen by the scattered quasiparticles are not necessarily given by the distances of the adatoms of the walls. Here, *m*_{eff} = 0.58*m*_{e}, the surface-state band edge *E*_{0} = −26.4 meV and *δ* = −0.28 nm are used, as motivated in Supplementary Note 2. The LDOS patterns presented in Figs. 1e and 2b have been calculated as a sum of the individual eigenfunctions with a finite Lorentzian broadening of *Γ* = 3 meV acting on their eigenenergies:

$${\rm{L}}{\rm{D}}{\rm{O}}{\rm{S}}(E)=\sum _{{n}_{x},{n}_{y}}\frac{{|\varPsi ({n}_{x},{n}_{y})|}^{2}}{1+{(E-E({n}_{x},{n}_{y}))}^{2}/{\varGamma }^{2}}.$$

(7)

### MSS model

We start by considering a system of a single spatially extended spin-degenerate level coupled to an *s*-wave superconducting three-dimensional bath, being a generalization of the model introduced in ref. ^{12}:

$$\begin{array}{l}{\mathscr{H}}=\sum _{{\bf{k}},\sigma }{{\epsilon }}_{{\bf{k}}}{c}_{{\bf{k}},\sigma }^{\dagger }{c}_{{\bf{k}},\sigma }+\sum _{{\bf{k}},\sigma }(\mathop{V}\limits^{ \sim }({\bf{k}}){c}_{{\bf{k}},\sigma }^{\dagger }{{d}_{\sigma }+{\mathop{V}\limits^{ \sim }{({\bf{k}})}^{\ast }d}_{\sigma }^{\dagger }c}_{{\bf{k}},\sigma })+\sum _{\sigma }{E}_{{\rm{r}}}{d}_{\sigma }^{\dagger }{d}_{\sigma }\\ \,\,-\,{\varDelta }_{{\rm{s}}}\sum _{{\bf{k}}}({c}_{{\bf{k}},\uparrow }^{\dagger }{c}_{-{\bf{k}},\downarrow }^{\dagger }+{c}_{-{\bf{k}},\downarrow }{c}_{{\bf{k}},\uparrow }).\end{array}$$

(8)

The Hamiltonian given in equation (1) and in ref. ^{12} is a special case of equation (8) for a perfectly localized impurity level (\(\widetilde{V}({\bf{k}})=V={\rm{constant}}\)). Here and in the following, we set *ħ* = 1.

We aim to calculate the LDOS at the local level. For that, we use the Green’s function equations of motion in energy space^{58}

$$E{G}_{{a}_{i,}{a}_{j}^{\dagger }}(E)={\delta }_{ij}+\langle \langle [{a}_{i},{\mathscr{H}}]\,;{a}_{j}^{\dagger }\rangle \rangle ,$$

(9)

in which \({G}_{{a}_{i,}{a}_{j}^{\dagger }}(E)=\left\langle \left\langle {a}_{i}\,;{a}_{j}^{\dagger }\right\rangle \right\rangle \) is the shorthand notation for the usual retarded Green’s function^{58}, for which *a*_{i} is one of the operators *d*, *c*_{k} or their adjoint. The LDOS at the local level is

$${\rm{LDOS}}(E)=-\frac{1}{{\rm{\pi }}}{\rm{Im}}\left[{G}_{{d}_{\uparrow },{d}_{\uparrow }^{\dagger }}(\omega )+{G}_{{d}_{\downarrow },{d}_{\downarrow }^{\dagger }}(\omega )\right],$$

(10)

in which *ω* = *E* + i*δE* and *δE* is a small and positive real number approximating the experimentally observed energy broadening. We obtain the Green’s function by solving the system of equations of motion in equation (9) for the Hamiltonian in equation (8) after linearizing the dispersion around *E*_{F}, that is, *ϵ*_{k} = *v*_{F}(*k* − *k*_{F}), with *v*_{F} and *k*_{F} being the Fermi velocity and momentum, respectively:

$${G}_{{d}_{\sigma },{d}_{\sigma }^{\dagger }}(\omega )=\frac{\omega +{E}_{{\rm{r}}}-\sum _{{\bf{k}}}|{\mathop{V}\limits^{ \sim }({\bf{k}})|}^{2}\frac{(\omega -{{\epsilon }}_{{\bf{k}}})}{({\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2})}}{{G}_{1}-{G}_{2}}$$

(11)

with

$${G}_{1}=(\omega +{E}_{{\rm{r}}}-\sum _{{\bf{k}}}\frac{|{\mathop{V}\limits^{ \sim }({\bf{k}})|}^{2}(\omega -{{\epsilon }}_{{\bf{k}}})}{({\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2})})\,(\omega -{E}_{{\rm{r}}}-\sum _{{\bf{k}}}\frac{|{\mathop{V}\limits^{ \sim }({\bf{k}})|}^{2}(\omega +{{\epsilon }}_{{\bf{k}}})}{({\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2})})$$

and

$${G}_{2}=(\sum _{{\bf{k}}}\frac{{\varDelta }_{{\rm{s}}}\mathop{V}\limits^{ \sim }({\bf{k}})\mathop{V}\limits^{ \sim }(\,-\,{\bf{k}})}{({\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2})})\,(\sum _{{\bf{k}}}\frac{{\varDelta }_{{\rm{s}}}{(\mathop{V}\limits^{ \sim }({\bf{k}})\mathop{V}\limits^{ \sim }(-{\bf{k}}))}^{\ast }}{({\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2})})$$

The impurity, which is described by its coupling to the substrate *V*(**r**), has a localization length *L*_{imp} corresponding to the size of the QD, and drops to zero for |**r**| ≫ *L*_{imp}. We can therefore reasonably set the corresponding Fourier transform \(\widetilde{V}({\bf{k}})\) constant for momenta *k* = |**k**| in the interval [*k*_{F} − *β*/*L*_{imp}, *k*_{F} + *β*/*L*_{imp}], in which *β* is on the order of one, whereas its concrete value depends on the spatial details of the impurity coupling *V*(**r**). In the following order-of-magnitude approximation, we set *β* = 1. From equation (11), we find that the physics of a spatially extended impurity does not differ from that of a localized impurity if \(\frac{1}{{\omega }^{2}-{{\epsilon }}_{{\bf{k}}}^{2}-{\varDelta }_{{\rm{s}}}^{2}}\ll 1\) at *k* = *k*_{F} ± 1/*L*_{imp}. Combining the last two formulas, we find that an extended impurity can be considered as localized if *ω* is within a few *Δ*_{s} from *E*_{F}, which is the case for the experiment in the main text, and if \(\frac{{v}_{{\rm{F}}}}{{\varDelta }_{{\rm{s}}}}=\xi \gg {L}_{{\rm{imp}}}\), in which *ξ* = *v*_{F}/*Δ*_{s} is the proximitized superconducting coherence length in the Ag islands. For Ag, the Fermi velocities range from 0.518 to 1.618 × 10^{6} m s^{−1} (ref. ^{59}), resulting in *ξ* = 253 to 789 nm, which is considerably larger than the maximal extent of our QDs reaching *L*_{imp} = 24 nm. The QD level can therefore be treated as a localized impurity. In this limit, the Green’s function can be written as

$${G}_{{d}_{\sigma },{d}_{\sigma }^{\dagger }}(\omega )=\frac{\omega +{E}_{{\rm{r}}}+\frac{\varGamma \omega }{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{\omega }^{2}}}}{{\omega }^{2}\,\left(1+\frac{2\varGamma }{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{\omega }^{2}}}\right)-{E}_{{\rm{r}}}^{2}-{\varGamma }^{2}},$$

(12)

in which *Γ* = π*V*^{2}*D*, with *D* = \({k}_{{\rm{F}}}^{2}W\)/(2π^{2}*v*_{F}) being the density of states per spin species of the substrate above the critical temperature at *E*_{F} and *W* is its volume. The LDOS is given by

$${\rm{L}}{\rm{D}}{\rm{O}}{\rm{S}}(E)=-\frac{2}{{\rm{\pi }}}{\rm{I}}{\rm{m}}\,\left[\frac{\omega +{E}_{{\rm{r}}}+\frac{\varGamma \omega }{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{\omega }^{2}}}}{{\omega }^{2}\left(1+\frac{2\varGamma }{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{\omega }^{2}}}\right){-{E}_{{\rm{r}}}}^{2}-{\varGamma }^{2}}\right].$$

(13)

We note the emergence of in-gap states as found in ref. ^{12}. The energy *ε*_{+} of this in-gap state for a range of values *E*_{r} and *Γ* is plotted in Extended Data Figs. 2 and 3. Recently, an LDOS of a localized impurity including further magnetic scattering has been derived^{60}. In contrast to a metallic bath, in which the scattering results in a spectral broadening of the local level, the superconducting bath induces superconductivity by proximity to the local level. Hence, when *E*_{r} lies within the gap of the superconductor, the state at *E*_{r} splits into two particle–hole-symmetric ones around *E*_{F}. Notably, for energy scales *E*_{r} sufficiently larger than *Δ*_{s}, equation (13) reduces to a typical Lorentzian LDOS of width *Γ* at position *E*_{r}, as observed in the experiment.

The obtained spin-degenerate single-level Hamiltonian with proximity-induced pairing (equation (2)) is equivalent to the Green’s function approach above to the second order in the coupling constant \(V\propto \sqrt{\varGamma }\), as we show in the following.

### Derivation of the effective Hamiltonian

In this section, we derive an effective low-energy model for the electronic level valid when the bare energy of the spin-degenerate electronic level is close to the Fermi energy and the coupling to the superconducting bulk is smaller than the superconducting gap. We find that the level obtains proximity pairing and a correction in its chemical potential.

The Hamiltonian of a spin-degenerate electronic level locally coupled to a Bardeen–Cooper–Schrieffer *s*-wave superconductor is given in equation (1), which we repeat here for convenience

$$\begin{array}{l}{\mathscr{H}}=\sum _{{\bf{k}},\sigma }{{\epsilon }}_{{\bf{k}}}{c}_{{\bf{k}},\sigma }^{\dagger }{c}_{{\bf{k}},\sigma }+\sum _{{\bf{k}},\sigma }V({c}_{{\bf{k}},\sigma }^{\dagger }{{d}_{\sigma }+{d}_{\sigma }^{\dagger }c}_{{\bf{k}},\sigma })+\sum _{\sigma }{E}_{{\rm{r}}}{d}_{\sigma }^{\dagger }{d}_{\sigma }\\ \,\,-\,{\varDelta }_{{\rm{s}}}\sum _{{\bf{k}}}({c}_{{\bf{k}},\uparrow }^{\dagger }{c}_{-{\bf{k}},\downarrow }^{\dagger }+{c}_{-{\bf{k}},\downarrow }{c}_{{\bf{k}},\uparrow }),\end{array}$$

(14)

in which *c*_{k,σ} are the annihilation operators in the superconducting bulk with momentum **k** and spin *σ*, *d*_{σ} the annihilation operator of the electronic level with spin *σ*, *ϵ*_{k} the dispersion relation in the bulk, *E*_{r} the electric potential of the electronic levels and *V* quantifies the local coupling between the electronic levels and the superconducting bulk.

To derive the low-energy model, we use the Schrieffer–Wolff transformation

$$S=\sum _{{\bf{k}},\sigma }{{\rm{s}}{\rm{g}}{\rm{n}}(\sigma )A}_{{\bf{k}}}{d}_{\sigma }{c}_{{\bf{k}},-\sigma }+{B}_{{\bf{k}}}{d}_{\sigma }{c}_{{\bf{k}},\sigma }^{\dagger }-{\rm{h.}}\,{\rm{c.}}$$

(15)

where h.c. is the Hermitian conjugate, with sgn(↑) = 1, sgn(↓) = −1, and

$${A}_{{\bf{k}}}=\frac{-V{\varDelta }_{{\rm{s}}}}{{{\epsilon }}_{{\bf{k}}}{{\epsilon }}_{-{\bf{k}}}+{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}},$$

(16)

$${B}_{{\bf{k}}}=\frac{V({{\epsilon }}_{-{\bf{k}}}-{E}_{{\rm{r}}})}{{{\epsilon }}_{{\bf{k}}}{{\epsilon }}_{-{\bf{k}}}+{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}},$$

(17)

to obtain the effective Hamiltonian

$${{\mathscr{H}}}^{{\prime} }={e}^{S}{\mathscr{H}}{e}^{-S}={{\mathscr{H}}}_{{\rm{D}}}^{{\prime} }+{{\mathscr{H}}}_{{\rm{S}}{\rm{C}}}^{{\prime} }+{\mathscr{O}}({V}^{3}).$$

(18)

The physics inside the superconducting gap is contained in the effective Hamiltonian \({{\mathscr{H}}}_{{\rm{D}}}^{{\prime} }\), which is that of a spin-degenerate electronic level with proximity-induced superconductivity

$${{\mathscr{H}}}_{{\rm{D}}}^{{\prime} }=\sum _{\sigma }({E}_{{\rm{r}}}+{E}_{{\rm{s}}{\rm{h}}{\rm{i}}{\rm{f}}{\rm{t}}}){d}_{\sigma }^{\dagger }{d}_{\sigma }-{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}({d}_{\uparrow }^{\dagger }{d}_{\downarrow }^{\dagger }+{d}_{\downarrow }{d}_{\uparrow }),$$

(19)

with the induced gap *Δ*_{ind} and the shift *E*_{shift} in the chemical potential

$${\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}=-\,\sum _{{\bf{k}}}V{A}_{{\bf{k}}},$$

(20)

$${E}_{{\rm{s}}{\rm{h}}{\rm{i}}{\rm{f}}{\rm{t}}}=\sum _{{\bf{k}}}V{B}_{{\bf{k}}}.$$

(21)

We approximate equations (20) and (21) by linearizing the dispersion relation *ϵ*_{k} close to the Fermi momentum *k*_{F} by

$${{\epsilon }}_{{\bf{k}}}={v}_{{\rm{F}}}(k-{k}_{{\rm{F}}}),$$

(22)

in which *v*_{F} is the Fermi velocity of the superconductor and we only consider momenta within the range [*k*_{F} − *Λ*, *k*_{F} + *Λ*]. For a three-dimensional host superconductor, we find

$${\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}=2{V}^{2}D\,\arctan \,\left(\frac{\varLambda {v}_{{\rm{F}}}}{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}}}\right)\frac{{\varDelta }_{{\rm{s}}}}{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}}}\lesssim {V}^{2}D\,{\rm{\pi }}\frac{{\varDelta }_{{\rm{s}}}}{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}}}=\varGamma \frac{{\varDelta }_{{\rm{s}}}}{\sqrt{{\varDelta }_{{\rm{s}}}^{2}-{E}_{{\rm{r}}}^{2}}},$$

(23)

$${E}_{{\rm{s}}{\rm{h}}{\rm{i}}{\rm{f}}{\rm{t}}}=-\,{E}_{{\rm{r}}}\frac{{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}}{{\varDelta }_{{\rm{s}}}},$$

(24)

in which *Γ* and *D* are defined as in the main text. We infer that the effective Hamiltonian of the spin-degenerate electronic level close to *E*_{F} obtains a proximity-induced superconducting pairing. From equation (19), we calculate the energy *ε* of the level and the hole weight |*v*|^{2} of the negative-energy eigenvalue to be

$${\varepsilon }=\pm \sqrt{{E}_{{\rm{r}}}^{2}{(1-{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}/{\varDelta }_{{\rm{s}}})}^{2}+{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}^{2}},$$

(25)

$${|v|}^{2}=\frac{1}{2}-\frac{{E}_{{\rm{r}}}\left(1-\frac{{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}}{{\varDelta }_{{\rm{s}}}}\right)}{2{\varepsilon }}=\frac{1}{2}-\frac{\sqrt{{{\varepsilon }}^{2}-{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}^{2}}}{2{\varepsilon }},$$

(26)

in which we have neglected orders of \({\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}^{3}\) and higher in the last step. If the spin-degenerate electronic level originally lies at the Fermi energy, that is, *E*_{r} = 0, its effective Hamiltonian only contains induced superconductivity and its Bogoliubov quasiparticles have 50% particle and 50% hole content. Moreover, the resonances are located at ±*ε*_{min} = ±*Δ*_{ind} for *E*_{r} = 0. Thus, the proximity-induced pairing strength can be readily inferred from measuring the value of *ε*_{min}.

Using |*u*|^{2} = 1 − |*v*|^{2} for the particle weight, the Bogoliubov angle, which conveniently measures the amount of particle–hole mixing, takes the form

$${\theta }_{{\rm{B}}}({\varepsilon })=\arctan (\sqrt{{|u|}^{2}/{|v|}^{2}})=\arctan \,\left(\sqrt{\frac{1+\sqrt{{{\varepsilon }}^{2}-{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}^{2}}/{\varepsilon }}{1-\sqrt{{{\varepsilon }}^{2}-{\varDelta }_{{\rm{i}}{\rm{n}}{\rm{d}}}^{2}}/{\varepsilon }}}\right).$$

(27)

Notably, equation (27) is independent of *Δ*_{ind} if the energies *ε* are normalized by *Δ*_{ind}. This is the reason why equation (27) is used to plot the theoretical curve in Fig. 4. For the Bogoliubov angle of the MSSs based on the LDOS given in equation (13), the energy-dependent *θ*_{B}(*ε*) varies with *Γ* and is thus different for each eigenmode. In Extended Data Fig. 4, we compare the Bogoliubov angle for a single superconducting level (equation (27)) with the expected Bogoliubov angle of MSSs using the expression for the LDOS calculated in equation (13). In the low-energy limit, both theories agree well, verifying that the anticrossing of the MSSs is evidence for superconducting pairing in the spin-degenerate level. For higher energies, the MSSs approach the coherence peak of the bulk gap and their asymmetries decrease again and finally converge to zero (equivalent to *θ*_{B} approaching π/4 at *Δ*_{s}, marked by the dashed blue lines in Extended Data Fig. 4). This leads to an even better agreement with the experimental data and demonstrates that the observed resonances indeed behave like MSSs.