Microwave excitation of atomic scale superconducting bound states

Magnetic impurities on superconductors lead to bound states within the superconducting gap, so called Yu-Shiba-Rusinov (YSR) states. They are parity protected, which enhances their lifetime, but makes it more difficult to excite them. Here, we realize the excitation of YSR states by microwaves facilitated by the tunnel coupling to another superconducting electrode in a scanning tunneling microscope (STM). We identify the excitation process through a family of anomalous microwave-assisted tunneling peaks originating from a second-order resonant Andreev process, in which the microwave excites the YSR state triggering a tunneling event transferring a total of two charges. We vary the amplitude and the frequency of the microwave to identify the energy threshold and the evolution of this excitation process. Our work sets an experimental basis and proof-of-principle for the manipulation of YSR states using microwaves with an outlook towards YSR qubits.

Magnetic impurities coupled to a superconductor give rise to Yu-Shiba-Rusinov (YSR) states, which are subgap states protected by parity (even/odd particle number conservation) [1][2][3].They exhibit a variety of interesting phenomena including (but not limited to) their resonant character, which enhances higher order processes in tunneling (Andreev processes) or their parity protection, which enhances their lifetime [4][5][6].Comparatively long coherence times can also be expected in YSR states, but work on coherent coupling of YSR states so far has been limited [4,7].The rst step towards coherent manipulation is the use of microwaves in a tunnel junction, which leads to microwave-assisted tunneling [8][9][10].However, parity conservation has to be considered when exciting a YSR state using microwaves.
Elementary excitations in a superconductor, i.e.Bogoliubov quasiparticles, come in pairs due to parity conservation, but only one quasiparticle is needed to excite the YSR state [11].The second quasiparticle can escape to the continuum, which requires excitation energies of at least the superconducting gap, or through a tunneling contact, where much lower excitation energies are sufcient.A scanning tunneling microscope (STM) provides such a tunneling contact o ering the ability to manipulate a YSR state with moderate excitation energies far below the superconducting gap.This makes the STM an ideal platform for the manipulation of YSR states as an extension of nondegenerate Andreev bound states to the atomic scale [12], which provides a starting point for YSR qubits [13][14][15][16].
Here, we demonstrate the excitation of YSR states using microwaves in the tunnel junction of an STM.We are able to separate di erent tunneling processes involving the YSR states, which allows us to identify a tunneling process that is only possible through the direct excitation of a YSR state by the microwave.We map out an amplitude threshold that has to be overcome to excite the YSR state.This threshold depends on the applied bias voltage, which allows for great exibility in di erent YSR excita-tion schemes.In this way, we provide a proof of principle for the excitation and manipulation of YSR states by microwaves in the presence of a tunnel junction, which is an important prerequisite for the preparation and control of complex YSR structures, for example, in the context of quantum simulations.
We use a scanning tunneling microscope (STM) with an external microwave antenna optimized for operation between 60 GHz and 90 GHz [18], which is schematically shown in Fig. 1(a).By controlled dipping of a vanadium tip in a V(100) surface, we create a YSR state at the apex of the tip [4,17], which is subsequently irradiated by microwaves.In Fig. 1(b), the di erential conductance (green line) through a YSR state in the absence of microwaves is shown.The salient features of the YSR state are two sharp peaks in the superconducting gap at  = ±( s +) corresponding to the electron and hole parts of the Bogoliubov quasiparticle ( is the bias voltage,  , is the superconducting gap parameter in tip and sample, and  is the YSR energy).
In recent experiments, microwaves have successfully been implemented in STMs with various applications, such as resolving the internal structure of complex tunneling processes.Initial experiments on clean superconductors [10] show good agreement with a theory for microwave-assisted tunneling [19,20], which we refer to in the following as Tien-Gordon (TG) theory.This theory predicts the formation of replicas of very sharp spectral features (e.g.coherence peaks, YSR states) at integer multiples of ℏ r / weighted by a squared Bessel function ( r is the microwave radiation frequency, ℏ is Planck's constant, and  is the elementary charge), which depends on the microwave amplitude.Further work has shown that this theory needs to be generalized beyond the tunneling regime for higher order processes such as the Josephson e ect or Andreev re ections [8].For a non-resonant transfer of  charges, replicas form at multiples of ℏ r / [21][22][23].Also, it has been demonstrated that replicas of YSR states can show asymmetries which are not contained arXiv:2303.13098v1[cond-mat.supr-con]23 Mar 2023  c) are ts to the data using the full Green's function model and two transport channels (one BCS and one YSR channel (cf.[4,17])).(d-g) Schematics illustrating ground state and excited state tunneling processes with and without microwaves.The schematics are drawn for the zero order processes, i.e. no net energy quanta transferred.Energy quanta may be absorbed/emitted in steps 1 / 3 leading to replicas at di erent bias voltages.
within the TG theory [9].This was corroborated by a simpli ed Green's functions approach [24].
The microwaves induce an alternating voltage  ac in the tunnel junction, which is on the order of 100 μV to 10 mV.The conductance spectrum with a YSR state irradiated by microwaves at a frequency of  r /2 = 60.05GHz and an amplitude of 570 V is shown in Fig. 1(c) (yellow green line).We note that the temperature of the junction only increases by a few mK, which we can safely assume to be constant in line with previous work [8].The interaction of the tunneling electrons with the microwave leads to both the absorption and emission of energy quanta by the tunneling electrons in integer multiples of ℏ r = 248.3μeV.In the simplest approximation, this interaction leads to the appearance of replicas of the spectral features in Fig. 1(b).In Fig. 1(c), the expected replicas of the YSR states are indicated by blue vertical lines at distances of 248.3 μV.However, we also observe a number of additional peaks marked by the red vertical lines, which appear at  = ±( s −) +ℏ r , where  is an integer.This might suggest a thermal origin, but the temperature of 560 mK is very low and no corresponding peak can be seen in the spectrum in the absence of microwaves (cf.Fig. 1

(b)).
To understand the origin of the di erent peaks seen in Fig. 1(c), we present schematics of the underlying tunnel-ing processes in Fig. 1(d)-(g).To induce tunneling through the YSR state without microwaves, we apply a bias voltage of  =  s +  as shown in Fig. 1(d).To illustrate this, we divide the tunneling process into three steps using the density of states picture.In the rst step (labelled 1 ), an electron is transferred across the tunnel junction.In the second step (labelled 2 ) a Cooper pair is split lling the hole, but leaving the YSR state excited.This excited quasiparticle then relaxes into the continuum (step 3' ) or tunnels across the junction as well (step 3 ).If the tunnel coupling is weak, quasiparticle relaxation in the YSR electrode dominates (step 3' ).As the tunnel coupling increases, step 3 becomes dominant transferring a total of two charges across the junction.This step is termed "resonant Andreev process" as its tunneling path involves a real state (the YSR state [6,9,25]) instead of a virtual state as in conventional Andreev re ections [26].We note that higher-order transfer processes appear in resonant tunneling processes at much lower conductances than for "conventional" tunneling, e.g.Andreev reections.Therefore, a theoretical description has to include these processes processes already at a conductance of 2.2 × 10 −3  0 , where Andreev re ections can still be neglected ( 0 = 2 2 /ℎ is the quantum of conductance).The event illustrated in Fig. 1  In the presence of microwaves, the tunneling process indicated by the blue arrow in the experimental spectrum in Fig. 1(c) is schematically shown in Fig. 1(e).We rst observe the conventional peak to appear at a bias voltage of  =  s + , which implies that a total of zero energy quanta are exchanged with the microwave during step 1 .However, energy quanta can be exchanged during step 3 , yet without shifting the position of the peak.
In fact, the peak position only changes if energy quanta are absorbed or emitted during step 1 such that they appear at di erent bias voltages  = ±( s + ) + ℏ r in the spectrum.Other than that, the process is analogous to the tunneling without microwaves (cf.Fig. 1(d)).In the following, we consider the two processes involving step 3 and 3' together and refer to this family of peaks as ground state tunneling.
The additional peaks seen as red lines in Fig. 1(c) cannot be explained by ground state tunneling (cf.Fig. 1(d) and (e)).They can be attributed to processes which originate from tunneling events in absence of microwaves at bias voltages of  =  s −  as depicted in Fig. 1(f), where we would expect them to occur via thermal activation.However, in our experiment, the Boltzmann factor exp (−   B  ) for a YSR state of energy  = 280 V (for Fig. 2-4) at a temperature of 0.56 K predicts a contribution of 0.03%, so that thermal excitations are strongly suppressed.Indeed, Fig. 1(b) shows no spectral feature, where the red arrow is pointing.When we turn on the microwaves, a strong and clear peak can be observed at the location of the red arrow, in contrast to a strong peak in presence of microwaves in Fig. 1(c).In this situation, the microwaves open new transfer channels as delineated in Fig. 1(g).The absorption of multiple energy quanta during step 1 induces an excited YSR state (step 2 ) and allows for subsequent relaxation into the continuum through step 3 .Multiple quanta being absorbed or emitted during process 3 then lead to a family of additional peaks marked by the red lines in Fig. 1(c) at bias voltages  = ±( s − ) + ℏ r .All the peaks of this family have in common that the excited state is aligned with the coherence peak through the bias voltage modulo an integer number of microwave quanta, which is why call these processes excited state tunneling.
In order to understand the evolution of the ground state and excited state tunneling more quantitatively, we measure di erential conductance spectra as function of the dimensionless microwave amplitude  =  ac /ℏ r .the quantum of conductance.We can clearly see many well de ned peaks, which we will assign to ground state or excited state tunneling in the following.In order to distinguish these peaks, we use the TG model to calculate the expected microwave amplitude dependence from the measured conductance spectrum without microwaves [10,19] where   () is the -th order Bessel function of the rst kind and  0 ( ) is the tunneling current without microwaves.The calculated image starting from the zero amplitude spectrum in Fig. 2(a) is shown in Fig. 2(b).We note that the TG model does not reproduce all of the experimentally observed peaks.The replicated peaks in Fig. 2(b) are entirely due to ground state tunneling, so that all additional peaks in Fig. 2(a) must be due to excited state tunneling.For comparison, we calculate the data set in Fig. 2(a) using the full Green's function theory taking into account microwaves, higher order tunneling processes (e.g.Andreev processes) as well as the interference between them [22] (for details see the Supplementary Information [27]).We found that due to the resonant tunneling through the YSR states the interplay between the microwave and the higher order tunneling processes become non-negligible such that approximative calculations fail and the full Green's function model has to be applied (for details see the Supplementary Information [27]).The calculation is shown in Fig. 2(c), which shows excellent agreement with the measured data in Fig. 2(a).Both ground state and excited state tunneling processes are reproduced with the full Green's function model.
To substantiate our claim that there are indeed two fam-ilies of processes, we present frequency dependent differential conductance spectra at a constant dimensionless amplitude of  = 3 in Fig. 3(a).The higher the order of the replica, the more tilted the spectral feature will appear in the map.An -th order replica of a feature at  0 moves as  =  0 + ℏ r .The replica and their dispersion are calculated from the full Green's function theory in Fig. 3(b) as well as presented schematically in Fig. 3(c).We can identify four vertical lines corresponding to zero order replicas, marked by the lines connecting panels (a), (b), and (c).The blue and red colors mark ground state tunneling ( 0 = ±( s + )) and excited state tunneling ( 0 = ±( s − )), respectively.We note that at  = 3, the microwave has enough power to excite the YSR state, such that excited state tunneling becomes possible.If the excited state replica actually did appear in the spectrum without microwaves, which is not the case (cf.Fig. 1(b)), the original spectrum would appear as in Fig. 3(d).In Fig. 3(d), the excited state tunneling peak is arti cially added (red line), where thermal tunneling would appear.However, microwaves could trigger these transfer processes to occur beyond a given threshold as discussed below.
In essence, the breakdown of the simple TG model (Fig. 2(b)) is expected because it leaves the ground state untouched and only considers the spectrum in the absence of microwaves without taking into account processes activated by the microwaves, such as excited state tunneling.The full model (Fig. 2(c)) agrees quantitatively with the experiment (cf.excellent t in Fig. 1(c)).An intuitive understanding of the mechanism behind excited state tunneling can be derived from a simpli ed model.Employing a perturbative approach including second order resonant Andreev processes, the excited state tunneling cur-rent  ex,e/h ( , ) appears as where e/h refers to the peak at negative/positive bias voltage  = ∓( s − ).The bare excited state tunneling current  0 ex,e/h ( ) is replicated by the microwave beyond an amplitude threshold (see below).This is described in the Supplementary Information along with details on the approximations being used [27].We further introduce a weight function which sums over all possible energy quanta that can be exchanged during step 1 , where  0 = 2 ℏ r is the minimum number of quanta needed to excite the YSR state (cf.Fig. 1(g)).The excited state tunneling current in Eq. ( 2) and the weight function in Eq. ( 3) show that step 1 in Fig. 1(g) only contributes to the magnitude of the current, but it does not generate any replica.This also explains why the replica are a distance ℏ r / apart despite two charges being transferred in the whole process.A very similar argument can be made for step 3 of the ground state tunneling in Fig. 1(e).However, in this case the sum condition in the weight function is  > − 0 , which does not introduce a new threshold, but just leads to a renormalization of the spectral weight.A number of di erent approximations between the full model and the simpli ed model in Eq. ( 2) can be made, e.g.[24], which are discussed in the Supplementary Information [27].
In Fig. 4(a), we separately calculate the ground state and excited state tunneling conductances using Eq. ( 2) and the corresponding formula for ground state tunneling [27], which are shown in blue and red, respectively.The stepped shaded area around zero bias voltage represents the threshold amplitude needed to activate excited state tunneling.The weight function for the  = −1, 0, 1 processes marked in Fig. 4(a) are plotted in Fig. 4(b) as function of dimensionless microwave amplitude .We see that the threshold is not a sharp cuto , but follows the leading edge of the lowest order Bessel function  2  2,3,4 (), respectively, enabling the process.For  = −1, 0, 1, the threshold is roughly at  ≥ 2, 3, 4, respectively, when the weight function becomes signi cant.To demonstrate the threshold e ect of the weight function, we plot the corresponding data from Fig. 2(a) at  = −1, 0, 1 in Fig. 4(c).The ts are shown with and without the weight function as solid and dashed grey line, respectively.We can directly see how the weight function imposes the threshold for small amplitudes and nicely follows the experimental data.
At this point, we emphasize that here YSR states are excited using energies much smaller than the minimal energy  >  s + , if the YSR state is connected to a tunnel junction.In fact, the excitation energy can be as low as 2 (cf.Fig. 1(g) and Eq. ( 2)), which we have demonstrated through the excited state tunneling process and the imposed activation threshold.Even though two electrons are transferred in the resonant tunneling process through the YSR state, the replica are spaced ℏ r / apart instead of ℏ r /2 as for conventional Andreev re ections.Hence, the spacing between replica cannot be used for inferring the number of charges being transferred.Our ability to excite YSR states with high precision can now be exploited for direct manipulation protocols.This opens up new possibilities for pump-probe schemes to address the nite lifetime of YSR states.
In summary, we have conducted a proof-of-princple experiment showing that the combination of a tunneling current and microwave radiation can excite YSR states without the need to cross the energy gap.In particular, microwave-assisted tunneling can be used as a tool not only for ground state tunneling, but also for excited state tunneling.The sub-gap excitation is attractive for future applications (such as information storage) as it does not introduce decoherence by coupling to the continuum to which the YSR state is coupled.Therefore, microwaves could pave the path towards coherent manipulation, similar to ESR-STM [28] or Andreev qubit architectures [14].Additionally, this work has shown replicas at multiples of ℏ r / as opposed to the ℏ r /2 that one would expect for a two-electron process.Pulse schemes or shot noise measurements [29,30] could shed further light on this process.

TIP AND SAMPLE PREPARATION
The V(100) sample was cleaned by repeated Argon ion bombardment and annealing to 700 • C. The typical appearance of the surface are square terraces with an oxygen reconstruction as shown in Fig. S1 The tip was made superconducting using eld emission (40 V bias voltage and 15 A current).By controlled dipping (4 nm dip at 100 mV), YSR states were created on the apex of the tip [1,2].For the present system, depending on the exact composition of the apex, the YSR states appear at di erent energies, allowing us to tune the YSR energy to the relevant frequency range between 60-100 GHz.

MICROWAVE TRANSMISSION
The microwave setup is similar to the one introduced in Ref. [3].The microwave source is a Keysight 8257D frequency generator (up to 20 GHz), whose frequency output is multiplied by a factor of six using a Virginia Diodes WR12SGX device.A millimeter wave 511E attenuator is used to tune the attenuation.In the vacuum chamber, we use semirigid Cu coaxial cables, which, starting at the 4 K stage, is replaced by a superconducting semirigid coaxial cable.Finally, the radiation is transmitted through vacuum to the tunnel junction using a custom made bow-tie antenna on a chip [3].To measure the transfer function, we use a feedback scheme as shown in Fig. S2.We broaden the peak by applying a lock-in amplitude and then reduce the attenuation until the peak drops below a threshold of 80 % of its maximum value.Then the ratio of the actual peak height   to the original peak height  0 is used to calculated the ac amplitude according to: ℏ , where  0 is the zeroth order Bessel function of the rst kind.

General theory
In this section we show how the general theory of photon-assisted tunneling in superconducting junctions developed in Ref. 22 can be adapted to the description of the microwave-assisted tunneling through a magnetic impurity coupled to superconducting leads.Our goal is to calculate the current through a voltage biased superconducting tunnel junction in the presence of a monochromatic radiation of frequency  r .For simplicity, we focus on the case of a single channel contact.We assume that the external radiation produces an e ective timedependent voltage  () =  +  ac sin  r .The task now is to extend the theory for multiple Andreev re ections (MARs) to the case of such a time-dependent voltage, for which the so-called Hamiltonian approach is a convenient starting point [5].The irradiated single channel superconducting tunnel junction can be described by means of the following tight-binding-like Hamiltonian [5] where  L,R are the Hamiltonians following Bardeen-Cooper-Schrie er (BCS) theory for the isolated elec-trodes.In the coupling term, L and R stand for the outermost sites of each electrode, and  is a hopping parameter describing the coupling between these sites.This parameter determines the normal state transmission coe cient  of this model in a way that depends on the nature of the tunnel junction.For instance, in a tunnel junction formed by two conventional BCS superconductors, the transmission adopts the form with   being the electrodes' density of states at the Fermi level [5].
In this model the current evaluated at the interface between the two electrodes adopts the form The nonequilibrium expectation values in Eq. ( S3) can be expressed in terms of the Keldysh-Green functions Ĝ+− Thus, the current can be now written as ) where τ3 is the corresponding Pauli matrix in Nambu space, Tr denotes the trace in Nambu space and the t's are given by Here,  () =  0 +  0  + 2 cos  r  is the time-dependent superconducting phase di erence.In this expression,  0 is the dc part of the superconducting phase di erence,  0 = 2 /ℏ is the Josephson frequency, and the constant  =  ac /(ℏ r ) measures the strength of the coupling to the electromagnetic eld, and is proportional to the square root of the radiation power.
Using the relation where   () is the Bessel function of order , one can write the time dependence of the hopping as follows In order to determine the Green functions we follow a perturbative scheme and treat the coupling term in the Hamiltonian (cf.Eq. (S1)) as a perturbation.The unperturbed Green functions ĝ, correspond to the uncoupled electrodes in equilibrium, where the superscript  ,  denotes the retarded and advanced components, respectively.Following Ref. [5], one can express the current in terms of a  -matrix, rather than in terms of the Green functions.The  -matrix associated to the time-dependent perturbation of Eq. (S6) is de ned as where the • product is a shorthand for integration over intermediate time arguments.With this de nition, it is easy to show that As shown in Ref. [5], the current in terms of the T-matrix components reads In order to solve the  -matrix integral equation, it is convenient to Fourier transform with respect to the temporal arguments: (S14) Thus, one can nally write down the current as where the current amplitudes    can be expressed in terms of the T-matrix Fourier components, T   () ≡ T ( +  + ℏ r ,  +  + ℏ r ), in the following way We are interested in the dc current  dc .In general, this current is the sum of two contributions  dc =  B +  Shapiro , where  B ≡  0 0 (cf.Eq. (S15)) is a background current and  Shapiro = ,      0  ( −    ) is the current from Shapiro steps contribution at discrete voltages    = (/)ℏ r /2.We shall ignore the contribution from the Shapiro steps in the following and focus only on the background current.
Using the relations which can be demonstrated using the corresponding  -matrix equations for these components, we can write the dc current exclusively in terms of T   ≡ T  0 LR,0 as follows Finally, the T   ful ll the following set of linear algebraic equations where the di erent matrix coe cients adopt the following form in terms of the unperturbed Green functions , where we have used the shorthand notation (  L, ) , =   L,, ( +  + ℏ r ), where ,  = 1, 2 are indexes in Nambu space.

Approximations
In general, one has to solve Eq. (S19) numerically to then evaluate the current via Eq.(S18).However, in lowtransmission junctions there are a number of approximations that one can make.In the deep tunnel regime (when the tunnel coupling is the smallest energy scale), one can use the following approximation for the solution of Eq. (S19): This leads to the standard Tien-Gordon result for the tunneling of single quasiparticles (see below).
If we want to consider at least the lowest order Andreev re ection, the next approximation is Using this approximation in Eq. (S18), we get the lowestorder approximation for the contributions of both the quasiparticle current (| | 2 ) and the Andreev re ection (| | 4 ).Additionally, one gets terms like a higher order contribution for the quasiparticle current.
The previous two approximations are perturbative in nature and may lead to divergencies, if they are not properly regularized.This is what happens, for instance, when there is a bound state inside the gap with a very long lifetime.In those cases, one can x that problem by solving the following closed system for T  1 and T  −1 : whose solution is Note that in Eq. (S23) the di erent matrices have to be understood as big matrices in microwave space.It is worth remarking that this approximation exactly reproduces the results for the YSR problem for the typical transmissions of the experiments.Actually, there are intermediate approximations that seem to work very well.For instance, to regularize the quasiparticle term the following approximation su ces Here, one can ignore the o -diagonal elements (in Nambu space) of Ê, .
The minimal approximation to regularize the Andreev term is given by where again one can ignore the o -diagonal elements (in Nambu space) of Ê, .

YSR states + microwaves
To describe the tunneling through an YSR impurity we use the mean-eld Anderson impurity model put forward in Refs.[1,6].Within this model the Green's functions of the electrodes are given as follows.For the left electrode, which is superconducting, we use the standard BCS Green's functions: where  0,L is the density of states at the Fermi energy of the left electrode in the normal conducting state.On the other hand, the Green functions for the right electrode features a superconducting electrode with the impurity, adopt the form [6] ĝR where Here, we have de ned the tunneling rate  R =  0,R  2 R (a similar rate  L =   0,L  2 L describes the strength of the tip-impurity coupling).
Let us recall that the condition for the appearance of superconducting bound states is  () = 0.In particular, the spin-induced YSR states appear in the limit | |  R (and they are inside the gap when also  R  R ).In this case, there is a pair of fully spin-polarized YSR bound states at energies ±, where which in the electron-hole symmetric case  = 0 reduces to In this case, using the approximation of Eq. (S20) we arrive at the following expression for the quasiparticle current to the lowest order in the tunnel coupling: where   () is the density of states of electrode  and  () is the Fermi function.This is simply the standard Tien-Gordon result.Thus, because of the presence of YSR states inside the gap (with energy  > 0), one expects the microwaves to give rise to a series of conductance peaks at  =  S +  + ℏ r with a height that should evolve with the microwave power as  2  ().Using the approximation of Eq. (S21) and selecting the contribution to the resonant Andreev re ection, we arrive at the following expression for the current due to the resonant Andreev re ection (to lowest order in the tunnel coupling): where ( R ) 12 () is the anomalous Green function at the impurity site and it is given in Eq. (S27).Equations (S31) and (S32) nicely explain the physics of the experimental observation.The last remaining thing is to establish what are the simplest expressions that regularize these equations when the YSR states are very long lived.After some careful analysis, we have arrived at the following regularized expressions: where  R, () = (1/)Im ( and Notice that the only di erence with respect to the perturbative results above is the presence of a denominator that regularizes the eventual divergencies.For a comparison of these approximations with the full Green's function model, see Section .In Fig. S3, we show an example of the di erential conductance as a function of the bias voltage  and  computed with Eqs.(S33) and (S34), i.e. we computed the total current as the sum of those two contributions.The parameters are given in Table SI.As one can see, these equations reproduce all the salient features of the experiment.
Finally, to make contact with recent results [7], we dene the energy-dependent tunneling rates for electrons  e and holes  h as where the coherent factors ũ2 and ṽ2 in our model are .
With these de nitions, the anomalous Green's function in the impurity can be approximated by (for energies close to the YSR energy) Thus, the perturbative expression for the current contribution of the resonant AR of Eq. (S32) becomes (S38) Figure S3: Calculated di erential conductance map as function of microwave amplitude Calculation of the YSR state replica as function of bias voltage and microwave amplitude using the regularized quasiparticle and Andreev currents (Eq.( S33) and (S34)).The excited state tunneling is clearly visible for higher amplitudes.
This expression has to be compared with Eq. (40) in Ref. [7].A list of parameters that were used to obtain the calculated spectra in the di erent gures is given in Table SI.

VALIDITY OF THE MODELS
Resonant tunneling easily involves higher order tunneling long before higher orders become signi cant in nonresonant tunneling.This is particularly the case when resonant tunneling is combined with the interaction with microwaves.We, therefore, evaluate up to which transparencies such approximations are valid in di erent models.Here, we compare three models: 1. Full Green's functions model (Eq.(S18) with Eq.
(S19), black lines) 2. Green's function model to rst order in Andreev reections (Eq.(S18) with Eq. (S23), blue lines) 3. Regularized Andreev model (Eq.(S33) and (S34), red lines) To assess the agreement of model A referenced to model B with functions  A () and  B (), we evaluate the mean squared deviation referenced to the mean squared deviation from zero: We plot the evolution of the spectrum without microwaves as a function of conductance for each model in Fig. S4(a).Every spectrum is normalized by the normal state tunneling conductance.For each spectrum, we calculated the deviation referenced to the full Green's function model and plot this in Fig. S4(b).The regularized Andreev model deviates by more than 5% from the full calculation at a transparency of 4 × 10 −2 , whereas the deviations of the rst order model only become relevant at a transparency of 1 × 10 −1 .
In contrast to that, when the microwaves are included, the three models become inconsistent much faster.Figure S5(a) shows four spectra for each model calculated at different conductances.At the highest conductance, the regularized Andreev model shows deviations in peak height, amplitude and even peak position.The rst order approximation still performs much better.This means that the interference of higher order processes is crucial to properly describe the spectrum under microwave irradiation.In Fig. S5(b), the deviations start about two orders of magnitude sooner for the regularized Andreev model, which crosses the 5% mark at a transparency of 5 × 10 −4 .The rst order approximation crosses the 5% mark at a transparency of 8 × 10 −2 .
The comparison of the behaviour with and without microwaves leads to an important conclusion for the experimental data.Firstly, for measurements without microwaves, the three models remain consistent up to about  = 2 × 10 −2 , which corresponds to roughly 12 nA for a set point bias voltage of 4 mV.This means that for the lifetime broadening of 0.6 μeV (cf.Table SI) used here and for typical setpoint currents of  (100 pA), higher order contributions are not relevant.If the lifetime broadening is smaller, higher order contributions will become relevant at lower transparencies (i.e.smaller setpoint currents) [1,8].In contrast to that, the regularized Andreev model with microwaves shows signi cant disagreement already at a transparency of 5 × 10 −4 , corresponding to about 150 pA.This means that for typical measurements with microwaves and YSR states a full Green's function approach is necessary.The resonances due to the interactions with the microwave lead to a failure of the lowest order approximation [7].A list of parameters that were used to obtain the calculated spectra in the di erent gures (both main text and supplementary information) is given in Table SI.

DERIVATION OF THE SIMPLIFIED GROUND STATE AND EXCITED STATE TUNNELING MODEL
In the following, we will derive a simpli ed model to highlight the roles of the replicas in step 1 and step 3 of ground state and excited state tunneling (cf.schematic in Fig. 1(e) and (g) of the main text).We simplify the tunnel- Table SI: Fitting Parameters Table of t parameters that were used to calculate the spectra in the corresponding gures.The parameters are given in meV, except for , which is dimensionless, and  L,R , which is measured in μeV.The temperature was set to be 0.56 K, the coupling of the impurity to the substrate  R = 100 meV, and an overall Gaussian broadening was chosen to be 12.5 μeV.The right electrode (R) carries the YSR state, while the left electrode (L) features an empty gap.For the t in Fig.  ing and focus on the interplay of the Bessel functions and exchange of energy quanta.We start with Eq. (S38).In the case of long-lived YSR states, i.e. very small   , we can approximate the Lorentzian of the YSR state by a Dirac delta-function This step solves the integral in Eq. (S38) and the Andreev current becomes Each tunneling rate  e,h has two peaks at ±, such that we have a total of four peaks in the spectrum.To separate out these peaks, we use the Heaviside step function  () to de ne  ± e,h () =  (±) e,h () so that we can split  e,h () into  e,h () =  + e,h () +  − e,h ().

(S41)
Without microwaves, the four principal peaks correspond to two ground state and two excited state tunneling peaks at  = ±( + ) and  = ±( − ), respectively.The derivations for all of these four peaks are very similar, so that in the following we derive the behavior for one peak, which can be easily extended to the other peaks.
Derivation for excited state tunneling electron peak These peaks are located at bias voltages of  = −( − ) − ℏ r , i.e. at the bias voltage where  + e () is resonant.We assume that    , such that the Fermi function can be approximated by a step function.In order to observe this peak, the following conditions have to be fullled: 1 where we have de ned the weight function  (, ) as and where  0 = 2 ℏ , where is the ceiling function (  is de ned as  rounded to the next larger integer).The weight function does not change the position nor the number of the replicas.It only modi es the amplitude of the peak.This nicely explains the appearance of replica at integer multiples of ℏ/ instead of ℏ/2.The weight function also introduces a threshold through the condition  >  0 − , which means that  0 quanta of ℏ have to be absorbed from the microwave in order to excite the YSR state.The leading edge of the weight function determining the onset of the peak as function of microwave intensity is given by the lowest order Bessel function  2  0 − ().This means in particular that the bare tunneling current  0 ex,e ( ) as de ned above cannot be observed when the microwave is turned o .
Simpli ed tunneling equations for ground state and excited state tunneling We can straightforwardly extend the above derivation for all four peaks.We nd for the bare tunneling currents where the rst index (gr,ex) refers to ground state and excited state tunneling and the second index (e,h) refers to electron and hole tunneling, respectively.From these bare tunneling currents, which have one peak each, we nd the following equations to calculate the spectra with microwaves where  0 = 2 ℏ r is the minimum number of quanta needed to excite the YSR state.Interestingly, we nd that for ground state tunneling the weight function w (, ) does not impose a threshold for the activation of the tunneling process, because the condition  ≥ − 0 (for  = 0) always includes the zeroth order Bessel function, such that resonant Andreev processes are always possible without microwaves as has been discussed before [8].

Quasiparticle tunneling from the ground state
For completeness, we note that quasiparticle tunneling has to be considered, when modeling ground state tunneling, since the lifetime of the YSR state is not in nite in practice.In the deep tunneling regime, quasiparticle tunneling can be calculated from the Tien-Gordon model  qp ( , ) = ∑︁   2  ()  ( + ℏ r /, 0) .(S56) As Andreev processes become more dominant with increasing tunneling conductance, the quasiparticle current reduces (cf.also regularized quasiparticle current in Eq. (S33)) [1,8].Excited state tunneling is a two-electron tunneling process, so that quasiparticle tunneling does not apply for that process.The full Green's function model naturally includes all current contributions.

RESONANT VS. NONRESONANT ANDREEV PROCESSES
In the previous section, we have derived a simple model that nds a spacing of ℏ/ between the replica of the resonant Andreev processes despite two charges being transferred.By contrast, the replica of nonresonant Andreev re ections are spaced by ℏ/2.We can explain this difference in behavior by deriving a simpli ed equation for the regular Andreev re ection starting from the same Eq.(S32) as for the resonant Andreev processes.The main difference is that the anomalous Green's function ( R ) 12 () is no longer given by a resonance, but by the standard result which in the following we will approximate by a constant, such that |( R ) 12 ()| 2 =  2 R .This is justi ed because the relevant part of the anomalous Green's function that is probed here is very close to zero energy.We further de ne (S60) which is again a Lorentzian that nicely shows how the replica are spaced by ℏ/2.This demonstrates that depending on the presence of a resonance inside the superconducting gap, the spacing between replica changes accordingly.In this case, the spacing between replica cannot be used for inferring the number of charges being transferred.

Figure 1 :
Figure 1: Tunneling mechanisms of YSR states under microwave irradiation.(a) Schematic drawing of the experimental setup.(b) Di erential conductance measured without microwaves.Ground state tunneling is indicated by a blue arrow.No excited state tunneling is observed (red arrow).(c) Di erential conductance measured with microwaves at 61 GHz.The energy exchange with the microwave induces replicas.The zero order ground state tunneling is indicated by a blue arrow.Excited state tunneling induces additional peaks, with the zero order peak indicated by a red arrow.The dashed lines in (b) and (c) are ts to the data using the full Green's function model and two transport channels (one BCS and one YSR channel (cf.[4,17])).(d-g) Schematics illustrating ground state and excited state tunneling processes with and without microwaves.The schematics are drawn for the zero order processes, i.e. no net energy quanta transferred.Energy quanta may be absorbed/emitted in steps 1 / 3 leading to replicas at di erent bias voltages.
(d) leads to a spectral peak indicated by the blue arrow in Fig. 1(b).

Figure 2 :
Figure 2: Di erential Conductance as function of bias of bias voltage and microwave amplitude.(a) Experimental data measured at a setpoint of 500 pA at 3 mV with a microwave frequency 61 GHz.(b) Calculation based on the spectrum at zero amplitude in panel (a) using the Tien-Gordon (TG) model.The features connected to excited state tunneling are missing.(c) Full Green's function model (FM) calculation showing all details as in the experimental data.

Figure 3 :
Figure 3: Frequency dependence of the spectra at constant microwave amplitude .(a) Di erential conductance spectra measured as function of frequency at constant microwave amplitude  =  ac ℏ r = 3.(b) Calculated spectra in the same range as (a) (full model).(c) Theoretical location of normal state replicas (blue) and excited states replicas (red).(d) Base spectrum without microwaves (blue) and excited states added manually (red).The zeroth order replicas (vertical lines) connect the panels.
Figure 2(a) shows the di erential conductance measured at a microwave frequency of 61 GHz and a normal state conductance of  N = 2.2 × 10 −3  0 , where  0 = 2 2 /ℎ is

Figure 4 :
Figure 4: Excitation threshold for the YSR state.(a) Di erential conductance calculation (simpli ed model) to illustrate the origin of the tunneling processes.Ground state tunneling is shown in blue and excited state tunneling in red.(b) The weight function for excited state tunneling for di erent orders  = −1, 0, 1 as indicated by the vertical lines in (a).The initial threshold is clearly visible.(c) Slices of di erential conductance data as function of microwave amplitude for excited state tunneling at di erent orders  = −1, 0, 1 as indicated by the vertical lines in (a).The simpli ed model (SM) (solid gray line) ts well to the experimental data and nicely demonstrates the cuto due to the threshold at lower amplitudes compared to when the weight function is not considered (gray dashed line).The vertical blue lines indicate the threshold  = 2, 3, 4, respectively.

Figure S1 :
Figure S1: Topography of the V(100) surface.The data was obtained at a set point of 100 pA with a bias voltage of 3 mV.

Figure S2 :
Figure S2: Illustration of the algorithm for the transfer function determination.The plot shows di erential conductance spectra of a coherence peak at di erent values of microwave attenuation.The spectra are o set in voltage for clarity.Starting from a previously determined attenuation, the attenuation is reduced in steps of 1 dB until the peak height is below the threshold of 80% of the original peak.The data was measured at a frequency of 61.1 GHz, and a setpoint current of 100 pA at a bias voltage of 3 mV.

Figure S4 :
Figure S4: Comparison of the di erent models without microwaves.(a) Di erential conductance spectra calculated using the di erent models as function of selected junction transparencies  with the microwaves turned o .The full Green's function model is labeled "Full model", the Green's function model is labeled "1st order FM", and the regularized Andreev model is labeled "Reg.AM".The regularized Andreev model is calculated from the sum of the quasiparticle current and the Andreev current.(b) Deviations of the approximations from the full Green's function model as function of junction transparency .Both approximations only fail for high transparencies when higher order processes become relevant also for nonresonant tunneling processes.

Figure S5 :
Figure S5: Comparison of the di erent models with microwaves.(a) Di erential conductance spectra calculated using the di erent models as function of selected junction transparencies  with the microwaves turned on.The full Green's function model is labeled "Full model", the Green's function model is labeled "1st order FM", and the regularized Andreev model is labeled "Reg.AM".The regularized Andreev model is calculated from the sum of the quasiparticle current and the Andreev current.(b) Deviations of the approximations from the full Green's function model as function of junction transparency .The rst order model fails at somewhat lower transparencies than without microwaves, but the regularized model fails for two orders of magnitude lower transparencies than before.This indicates that interactions/interference between the resonant processes and the microwaves cannot be neglected for a quantitative and even a qualitative agreement.
Owing to the second condition, we approximate  h () by a constant, i.e.  h =  h ( ).