Full Counting Statistics of Yu-Shiba-Rusinov Bound States

With the help of scanning tunneling microscopy (STM) it has become possible to address single magnetic impurities on superconducting surfaces and to investigate the peculiar properties of the in-gap states known as Yu-Shiba-Rusinov (YSR) states. However, until very recently YSR states were only investigated with conventional tunneling spectroscopy, missing the crucial information contained in other transport properties such as shot noise. Here, we adapt the concept of full counting statistics (FCS) to provide a very deep insight into the spin-dependent transport in these hybrid systems. We illustrate the power of FCS by analyzing different situations in which YSR states show up including single-impurity junctions, as well as double-impurity systems where one can probe the tunneling between individual YSR states. The FCS concept allows us to unambiguously identify every tunneling process that plays a role in these situations. Moreover, FCS provides all the relevant transport properties, including current, shot noise and all the cumulants of the current distribution. Our approach can reproduce the experimental results recently reported on the shot noise of a single-impurity junction with a normal STM tip. We also predict the signatures of resonant (and non-resonant) multiple Andreev reflections in the shot noise of single-impurity junctions with two superconducting electrodes. In the case of double-impurity junctions we show that the direct tunneling between YSR states is characterized by a strong reduction of the Fano factor that reaches a minimum value of 7/32, a new fundamental result in quantum transport. The FCS approach presented here can be naturally extended to investigate the spin-dependent superconducting transport in a variety of situations, and it is also suitable to analyze multi-terminal superconducting junctions, irradiated contacts, and many other basic situations.

With the help of scanning tunneling microscopy (STM) it has become possible to address single magnetic impurities on superconducting surfaces and to investigate the peculiar properties of the in-gap states known as Yu-Shiba-Rusinov (YSR) states. These systems are an ideal playground to investigate multiple aspects of superconducting bound states, such as the occurrence of quantum phase transitions or the interplay between Andreev transport physics and the spin degree of freedom, with profound implications for disparate topics like Majorana modes or Andreev spin qubits. However, until very recently YSR states were only investigated with conventional tunneling spectroscopy, missing the crucial information contained in other transport properties such as shot noise. In this work we adapt the concept of full counting statistics (FCS) to provide the deepest insight thus far into the spin-dependent transport in these hybrid atomic-scale systems. We illustrate the power of FCS by analyzing different situations in which YSR states show up including single-impurity junctions with a normal and a superconducting STM tip, as well as double-impurity systems where one can probe the tunneling between individual YSR states [Nat. Phys. 16, 1227 (2020)]. The FCS concept allows us to unambiguously identify every tunneling process that plays a role in these situations and to classify them according to the charge transferred in them. Moreover, FCS provides all the relevant transport properties, including current, shot noise and all the cumulants of the current distribution. In particular, our approach is able to reproduce the experimental results recently reported on the shot noise of a single-impurity junction with a normal STM tip [Phys. Rev. Lett. 128, 247001 (2022)]. We also predict the signatures of resonant (and non-resonant) multiple Andreev reflections in the shot noise and Fano factor of single-impurity junctions with two superconducting electrodes and show that the FCS approach allows us to understand conductance features that have been incorrectly interpreted in the literature. In the case of double-impurity junctions we show that the direct tunneling between YSR states is characterized by a strong reduction of the Fano factor that reaches a minimum value of 7/32, a new fundamental result in quantum transport. The FCS approach presented here can be naturally extended to investigate the spin-dependent superconducting transport in a variety of situations, such as atomic spin chains on surfaces or superconductor-semiconductor nanowire junctions, and it is also suitable to analyze multi-terminal superconducting junctions, irradiated contacts, and many other basic situations.

I. INTRODUCTION
Yu-Shiba-Rusinov (YSR) bound states are one of the most fundamental consequences of the interplay between magnetism and superconductivity at the atomic scale [1][2][3]. They appear when a magnetic impurity is coupled to one or several superconducting leads, and they have been extensively investigated in the context of STM experiments and impurities on superconducting surfaces , for recent reviews see Refs. [32,33]. The appeal of these hybrid junctions is manifold. On the one hand, they enable the investigation of a fundamental quantum phase transition [6,7,16,18,27]. They also provide the chance to study states that are the precursors of Majorana modes, which might potentially appear when magnetic impurities are combined to form magnetic chains [34][35][36][37][38][39][40]. On the other hand, these systems are also a new kind of Andreev spin qubit, a topic that is rapidly * david.ohnmacht@uni-konstanz.de growing in the community of superconducting qubits [41][42][43][44]. More important for the topic of this work is the fact that magnetic impurities on superconducting surfaces offer the possibility of exploring the role of the spin degree of freedom in situations where the electronic transport is dominated by Andreev reflections, which is a topic of great interest for the field of superconducting spintronics [45,46].
Until very recently, YSR states had been only investigated by the means of standard tunneling spectroscopy, i.e., via current or conductance measurements. This has already allowed us to elucidate many basic properties of these superconducting bound states. However, as it is well-known in mesoscopic physics, the analysis of other transport properties, such as shot noise, may provide very valuable information that is out of scope of conventional tunneling spectroscopy [47,48]. An experimental breakthrough illustrating this idea was recently reported in which the shot noise through a magnetic impurity featuring YSR states was measured [49]. These experiments revealed the power of shot noise measurements by accessing scales like the intrinsic YSR lifetimes, which are very arXiv:2305.04758v2 [cond-mat.mes-hall] 30 Sep 2023 difficult to obtain via conductance measurements. These experiments call for more comprehensive theories able to provide knowledge about the YSR states and their transport properties beyond the conventional calculations of the charge current. In this regard, the concept of full counting statistics (FCS), introduced some time ago in the context of quantum transport in mesoscopic systems [50,51], is clearly the most powerful tool that we have at our disposal to reach this goal. This concept has already provided unprecedented insight into the physics of superconducting junctions [52][53][54][55][56], which includes not only the current, shot noise and all possible cumulants of the current distribution, but also the possibility to unambiguously identify the contribution of every possible tunneling process and the charge transferred in it. Unfortunately, the concept of FCS has not been used in the context of the superconducting hybrid junctions featuring YSR states mainly due to technical difficulties implementing the spin-dependent scattering that occurs in these systems.
In this work we fill this void by providing a systematic study of the FCS in a variety of junctions featuring YSR bound states. Those cases include single-impurity junctions with one and two superconducting leads, and two-impurity systems in which the tunneling between individual YSR states has been recently reported for the first time [25,26]. In all these cases we show how the concept of FCS allows us to identify all the relevant tunneling processes and to classify them according to the charge transferred in them, providing so a very deep insight into the physics of YSR states. Moreover, from the knowledge of the FCS in these situations we obtain all the relevant transport properties: current, shot noise, and all the cumulants of the current distribution. Among the main result of this work we can highlight the full analysis of the shot noise through a single-magnetic impurity coupled to a normal and to a superconducting lead in excellent agreement with very recent experimental results [49]. This analysis reveals the possibility to access energy and time scales that are usually out of the scope of conductance measurements, and it provides a very fresh insight into the nature of a resonant Andreev reflection. We also present very concrete predictions for the shot noise and Fano factor in the case in which a magnetic impurity is coupled to two superconducting leads for arbitrary junction transparency. In particular, we elucidate the signatures of resonant and YSR-mediated multiple Andreev reflections (MARs) in both the current and shot noise and amend some misinterpretations related to these processes that have been reported in the literature. In the case of two-impurity junctions, we demonstrate that the direct quasiparticle tunneling between YSR states provides an unique signature in the noise and Fano factor. In particular, this tunneling can be identified by a strong reduction of the Fano factor at the resonant voltage at which the two states align and it can reach a minimum value of 7/32. We quantitatively explain this result in terms of the quasiparticle tunneling between two sharp levels, which constitutes an extreme example of quantum tunneling. It is also worth stressing that the FCS approach presented here can be readily adapted to analyze the charge transport properties in a plethora of superconducting nanoscale junctions including magnetic atomic chains, superconductor-semiconductor nanowire hybrid junctions, multi-terminal Josephson junctions, irradiated junctions, just to mention a few.
The rest of this paper is organized as follows. First, in Sec. II we remind the basics of FCS in the context of electronic quantum transport. Then, Sec. III presents the details on how the FCS can be obtained in practice in all the systems analyzed in this work with the help of a powerful Keldysh action. Section IV shows how this action can be used in combination with a mean-field model for the YSR states to describe the FCS in single-impurity junctions. The results of this combination are discussed in detail in Sec. V for the case in which an impurity is coupled to a superconducting lead and a normal reservoir. In particular, we focus on the analysis of the resonant Andreev reflection that can take place in the presence of YSR states and we present a detailed discussion on very recent shot noise experiments [49]. This analysis is extended to the case of single-magnetic impurities with two superconducting electrodes in Sec. VII. In this case, we focus on the prediction of the signatures of MARs in the noise and Fano factor to guide future experiments. Section VIII deals with the case of two-impurity junctions in close connection with recent experiments [25,26], and shows how the tunneling between individual YSR states is revealed in the current fluctuations. Finally, we present an outlook and our main conclusions in Sec. IX. Some of the technical details related to the Keldysh action and the corresponding Green's functions are discussed in Appendices A and B.

II. FULL COUNTING STATISTICS: A REMINDER
The electronic transport in a quantum device can be viewed as a stochastic process which can be completely characterized by a probability distribution. Most experiments focus on the measurement of the electrical current, i.e. the average of that distribution. However, it has been shown in numerous systems that nonequilibrium current fluctuations (shot noise), i.e., the second cumulant of the current distribution, contain very valuable information out the scope of current measurements such as the charge of the carriers or the distribution of conduction channels [47,48]. Ideally, one would like to have access to all the cumulants of the current distribution to completely characterize the electronic transport. This idea was developed in the context of quantum optics and led to the introduction of the concept of photon counting statistics [57], which turned out to be key to characterize quantum states of the light such as those realized in a laser. In the 1990s, Levitov and Lesovik adapted this concept to mesoscopic electron transport [50,51], in which the electrons passing a certain conductor are counted. Later on, Nazarov and coworkers showed that these ideas could be combined with the powerful Keldysh-Green's function approach [48,58], which enabled the analysis of the electron counting statistics in numerous situations including those that involve superconducting electrodes [52][53][54]. In this section, we briefly remind the reader of the basics of full counting statistics in the context of quantum electronic devices, and in Sec. III we shall address how it can be computed in the situations in which we are interested in this work, namely superconducting hybrid junctions featuring YSR states.
Since the charge transfer in any junction is fundamentally discrete, one can aim at counting those charges. In a measuring time t 0 , there is a certain likelihood for N particles to cross the junction, which is described by a set of probabilities {P t0 (N )} N , the so-called full counting statistics (FCS), which contains all the information concerning the possible transport processes. In practice, the FCS can be obtained from the so-called cumulant generating function (CGF) A t0 (χ), which is given by where χ is the so-called counting field. Notice that because N P t0 (N ) = 1, the CGF is normalized: A t0 (0) = 0. The counting field χ is the conjugate variable to the charge number N and it is used as an auxiliary variable. By performing derivatives of the CGF with respect to the counting field, one obtains the cumulants C n as follows Thus, for instance, the first cumulant reads which is the expectation value of the charge number. Physically, it corresponds to the current I averaged over the time interval t 0 , I = eC 1 /t 0 . The second cumulant is given by which describes the variance. This second cumulant can be related to the zero-frequency shot noise S via S = 2e 2 C 2 /t 0 . An important measurable quantity is the socalled Fano factor, which is given by which is easily accessible in the framework of FCS. The Fano factor is a measure of the effective charge of the carriers in a system when the junction transmission is low and tunneling events are uncorrelated (Poissonian limit). In superconducting (SC) junctions, the Fano factor can be super-Poissonian (F * > 1), indicating that the transferred charge is larger than one due to the occurrence of Andreev reflections. In contrast, resonant tunneling might result in a so-called sub-Poissonian (F * < 1) Fano factor, where the interpretation of the Fano factor as an effective charge breaks down [47].
In the simple case of a two-terminal device featuring a single conduction channel characterized by an energyindependent transmission coefficient τ , the CGF at zero temperature reads (ignoring spin) where V is the bias voltage between the two terminals. The interpretation of this result is that the transport is dominated by single-electron tunneling with a probability P 1 = τ . This interpretation becomes apparent when considering the corresponding full counting statistics, which is given by This is a binomial distribution where we define the number of attempts M = ⌈eV t 0 /h⌉, where ⌈ ⌉ describes the next highest integer. This is justified by the fact that t 0 is chosen to be sufficiently large. Hence, it is evident that the transport can be described as N particles being transmitted in individual tunneling processes transferring a single electron charge with a probability τ . In the case of junctions containing SCs, there are additional multi-particle tunneling processes, namely MARs, whose probabilities can also be obtained from the knowledge of the CGF. As we shall show below, in the case where there is no spin-flip scattering the CGF of a superconducting contact can be expressed as [55,56] where P σ n (E, V ) is the spin-, energy-and bias-dependent probability of a transport process transferring n electron charges with spin σ. Thus, the evaluation of the CGF in the different physical situations that we shall address in this work will allow us to classify the different tunneling processes that can take place according to the charge they transfer, something that cannot be objectively done with any other theoretical method. Moreover, from the knowledge of those probabilities we can readily obtain the different transport quantities: current, noise, and higher-order cumulants. In the most general case of spin-flip, like in the example of Section VIII, the main difference will be the impossibility to spin-resolve the tunneling probabilities (see discussion below).

III. KELDYSH ACTION
Our goal is to obtain the probabilities of the different tunneling processes P n (E, V ) from the knowledge of the CGF, like in Eq. (8). Here, we shall make use of a remarkable result obtained by Snyman and Nazarov [60,61] in which the CGF of an arbitrary mesoscopic/nanoscopic device can be expressed in terms of two basic ingredientes: (i) the Green's functions of the electron reservoirs or leads, which can be superconducting, and (ii) the scattering matrix of the device in the normal state. Technically speaking, these researchers showed that the CGF of any type of junction (ignoring inelastic interactions) can be expressed as (in this section we drop the subindex t 0 and omit the prefactor t 0 /h) whereĜ(χ) contains the information of the reservoir Green's functions (GFs) andŜ is the normal-state scattering matrix whose structure will be explained in what follows. First,Ĝ(χ) is a GF in the lead-time-Keldyshspin-Nambu space, denoted by the symbol (ˆ). In this work, we focus on two-terminal settings for which this GF adopts the generic form which is a block-diagonal matrix with the two time-Keldysh-spin-Nambu GFs as entries on the diagonal. The counting field, similarly to the voltage, can be gauged away from the right onto the left lead. In what follows the right terminal will always be superconducting G R = G SC , while the left terminal can be either normal or superconducting. The GFs are in fact infinite matrices in time space, which means that the element (t, t ′ ) of the GF G L,R (χ) is the Keldysh GFǧ L,R (χ, t, t ′ ). In other words, the lead GF is an infinitely large matrix in time space whose entries are 16 × 16-matrices in lead-Keldysh-spin-Nambu space. For a detailed explanation of the structure of the GFs we refer to Appendix A. Equivalently, the normal state scattering matrixŜ is also an infinite matrix in time space whose (t, t ′ ) entry is a 16×16-matrix in lead-Keldysh-spin-Nambu space, which we denote as s(t, t ′ ), where the symbol (˜) indicates that the scattering matrix is expressed in lead-Keldysh-spin-Nambu space. However, the scattering matrix only depends on the relative time t rel = t − t ′ , meaning thats(t, t ′ ) =s(t rel ). To elaborate on the structure of the scattering matrix, in a two-terminal setting it can be written in lead space as followss whereť(t rel ) is the transmission matrix andř(t rel ) the reflection matrix in Keldysh-spin-Nambu space. Upon a Fourier transformation the scattering matrix reads The scattering matrix is unitary, meaning that s(E)s † (E) =1, where1 is the identity in lead-Keldyshspin-Nambu space. If the transport preserves timereversal symmetry, the scattering matrix is symmetric, meaning thats(E) =s T (E). The scattering matrix depends on the system under consideration and below we shall show how it can be obtained from models for the description of YSR states. Fourier transforming, we can write Eq. (9) in a Floquet representation as follows with the matricesG(χ, E) andS(E) both expressed in Floquet-lead-Keldysh-spin-Nambu space which are the Floquet representations of the GFĜ(χ) and the scattering matrixŜ respectively. In addition, we define the matrixQ(χ, E). The trace that initially went over time space goes now about the Floquet space and Eq. (13) can be rewritten in terms of the Floquet energy E ∈ [−eV, eV ] as follows (see Appendix B for more details) The trace now goes over the matrixQ(χ, E), which is infinite in the Floquet index with its entries being 16 × 16matrices in lead-Keldysh-spin-Nambu space. Note, that we can use that Tr lnQ = ln det(Q). In addition, in the case of single-impurity junctions the matrixQ(χ, E) is a block-diagonal matrix in spin space and thus the determinant factorizes into two contributions, a spin ⇑ and ⇓ contribution,Q ⇑ andQ ⇓ , which leads to the result where we have defined the charge-and spin-resolved counting polynomials P σ (χ, E). Note, that these counting polynomials follow from the determinant of infinitely large matrices, namelyQ σ . In addition, it holds that P σ (0, E) = 1. Therefore, it is evident that the counting polynomials take the following form where we encounter the spin-, charge-and energyresolved tunneling probabilities P σ n (E) with the respective counting factor e ınχ . Thus, we obtain the action in the form of Eq. (8) by adding the prefactor t 0 /h.

IV. SINGLE-IMPURITY YSR JUNCTIONS: MODEL AND SCATTERING MATRIX
The goal of this section is to use the concept of full counting statistics to describe the electronic transport properties of a junction featuring YSR states in which a single magnetic impurity is coupled to superconducting leads. As we showed in Sec. III, we need as an input a normal-state scattering matrix describing these types of junctions. To obtain it, we make use of a mean-field Anderson model with broken spin symmetry (see Ref. [62] and references therein for a discussion on its origin and range of applicability), and which has been very successful describing different transport characteristics [23,27]. This model, which is illustrated in Fig. 1, describes the experimentally relevant situation in which a magnetic impurity (an atom or a molecule) is coupled to a superconducting substrate (S) and to an STM tip (t), which can also be superconducting. The model used here is summarized in the following Hamiltonian Here, H j (with j = t, S) is the BCS Hamiltonian of the lead j given by where c † kjσ and c kjσ are the creation and annihilation operators, respectively, of an electron of momentum k, energy ξ kj , and spin σ =↑, ↓ in lead j, ∆ j is the superconducting gap, and φ j is the corresponding superconducting phase. On the other hand, H imp is the Hamiltonian of the magnetic impurity, which reads Here, n σ = d † σ d σ is the occupation number operator on the impurity, U is the on-site energy, and J is the exchange energy that breaks the spin degeneracy on the impurity. Finally, H hopping describes the coupling between the magnetic impurity and the leads and adopts the form where t j describes the tunneling coupling between the impurity and the lead j = t, S and it is chosen to be real. It is convenient to rewrite the previous Hamiltonian in terms of four-dimensional spinors that live in a space resulting from the direct product of the spin space and the Nambu (electron-hole) space. In the case of the leads, the relevant spinor is defined as while for the impurity states we definē Using the notation τ i and σ i (i = 1, 2, 3) for Pauli matrices in Nambu and spin space, respectively, and with τ 0 and σ 0 as the unit matrices in those spaces, one can show that the Hamiltonian in Eq. (17) can be cast into the form The k-dependent retarded and advanced GFs of the leads can be expressed in spin-Nambu space as a function of the energy as where the Dynes' parameters η j are introduced. The kdependence can be eliminated by summing over it and defining the spin-Nambu GF where ϕ j is the superconducting phase of the order parameter of electrode j. In the case of a normal metal, the previous expression reduces tō Equivalently for the impurity, its GF is given bȳ where η imp is the regularization constant of the impurity GF. It will later be evident that it can be set to zero. In the following, we focus on the electron space of Nambu space. Namely, on the first and third component of the spinors in Eq. (21) and Eq. (22). Of special interest are the electron self-energies of the two leads which are given byΣ r/a t,e =V imp,t,eḡ r/a t,eVt,imp,e where the index ( e ) refers to extracting the first and third component of the respective quantity in spin-Nambu (31) We define the electron coupling matrix in spin space with Γ t/S,e by taking the imaginary part of the self-energȳ Let us recall that in the usual regime in which the STM experiments are operated Γ t ≪ Γ S , this model predicts the appearance of a pair of fully spin-polarized YSR bound states in the limit J ≫ ∆ S , and they are inside the gap when also Γ S ≫ ∆ S . In this case, the energy of the YSR states (measured with respect to the Fermi energy) is given by [62] In the electron-hole symmetric case U = 0, the previous expression reduces to The entries of the electron scattering matrix, namely the electron reflection matrices r e , r ′ e and the transmission matrices t e , t ′ e , can be computed using the Fisher-Lee relations following Ref. [63] r e = σ 0 − 2ı Γ t,e 1/2Ḡ r imp,e (Γ t,e ) 1/2 (35) t The hole-components of the scattering matrix follow from the electron-components with [61] It is important to remark that in the action of Eq. (9), the scattering matrix is the normal-state one. Hence, to obtain the desired result, we have to evaluate the transmission matrix in the case in which both reservoirs are in the normal state, i.e., when the corresponding GFs are given byḡ r/a j =ḡ r/a N . Thus, for instance, it is straightforward to show that the transmission matrixt in the spin-Nambu basis (Ψ † ↑ , Ψ ↓ , Ψ † ↓ , −Ψ † ) is given by the following diagonal matrix where we have defined the tunneling rates Γ t/S = πN 0,t/S t 2 t/S , where N 0,t/S corresponds to the normal density of state of the corresponding electrode. The tunneling rates describe the strength of the coupling between the impurity and the corresponding lead (j = t, S). It is evident that the impurity regularization constant can be set to zero η imp ≈ 0 as the substrate tunneling rate is always orders of magnitude larger than it.
Due to the unitary of the scattering matrix, it can be written in terms of the transmission matrixt as follows which is a 8 × 8-matrix in lead-spin-Nambu space. Note that the scattering matrix is proportional to the identity in Keldysh space in this case, thus its Keldysh structure is not included in the above formulas.
There are some interesting limiting cases to be discussed. First, in the case in which the energy of the impurity level is U = 0, the electron-and hole-transmission function are related via and thus |t ⇑ e | 2 = |t ⇑ h | 2 , so the transmission is the same for electrons and holes. Another interesting case is the high-coupling regime where Γ L , Γ R ≫ U, J. In that case, the transmission at low energies adopts the form so the transmission matrix for electrons and holes are the same. In the case where Γ t = Γ S , it holds that r ⇑ e = r ⇑ h ≈ 0 and the junction behaves as a single-channel highly transmissive point contact.

V. SINGLE-IMPURITY YSR JUNCTIONS: NORMAL CONDUCTING TIP
We now analyze the situation in which the STM tip is in the normal state, while the substrate is in the superconducting state. In this case the transport properties are a result of the competition between two tunneling processes, namely single-quasiparticle tunneling and an Andreev reflection, and the fact that both of them become resonant due to the presence of in-gap YSR states, see Fig. 2. In this case, the GFs only depend on the relative time and we just have to integrate over all energies E in the formula of the Keldysh action in Eq. (9). Moreover, in this single-impurity case the electrode GFs and the scattering matrix are block-diagonal in spin space. This allows us to carry out the calculation of the CGF analytically and the final result reads where P σ n (E, V ) corresponds to the probability of transferring n charges across the junction for spin σ. The single-quasiparticle tunneling probabilities from tip to substrate (n = 1) and from substrate to tip (n = −1) are given by Tunneling processes in the case of an impurity coupled to a normal and to a superconducting lead. In these diagrams the left density of states (DOS) corresponds to a normal STM tip and the right one to the impurity coupled to the superconducting substrate. The red lines correspond to electron-like quasiparticles and the blue ones to hole-like. In all cases, we indicate the threshold voltage at which the process starts to contribute to the transport. (a) Single-quasiparticle tunneling in which a quasiparticle tunnels either into the continuum DOS of the superconducting electrode (i) or resonantly into the excited YSR state inside the gap (ii). (b) Standard non-resonant Andreev reflection in which an electron is reflected as a hole. (c) A resonant Andreev reflection in which the electron that is retro-reflected impinges at the energy of a YSR state.
while the corresponding Andreev reflection probabilities from tip to substrate (n = 2) and from substrate to tip (n = −2) are given by In these expressions we have assumed that a bias voltage is applied to the normal metal, 2 are the GFs of the superconducting substrate, ∆ being the gap of the SC electrode and η the corresponding Dynes' parameter, ρ S = (g r − g a )/2 is the substrate density of states (DOS), f n (E) = f (E + neV ) with the Fermi function f (E) = 1/(e E/kBT + 1), and the normalization factor D σ (E) is given by Let us say that these probabilities reduce to the known result for a NS junction with energy-independent trans-mission and spin degeneracy [54]. Moreover, we have verified that they lead to the same results for the current as in Ref. [62]. These probabilities have the expected structure. Thus, for instance, the single-quasiparticle probabilities (n = ±1) are to a leading order proportional to the transmission coefficients (for electrons and holes) and to the DOS in the superconducting substrate. The Andreev reflection probabilities (n = ±2) are proportional to the product of the electron and hole transmission coefficients and to the Cooper pair density in the superconducting electrode (this will become more obvious in a moment). On the other hand, the denominators or normalization factors D σ are responsible for higher-order terms in the transmission and they contain the information of the energy of the bound states and their lifetimes. One can gain more insight into these expressions by considering the zero-temperature case where the above probabilities reduce to where −f a f r corresponds to the energy-dependent Cooper pair density. The other contributions are zero (P σ −1/−2 = 0) as there is no current flowing to the normal metal without thermal excitation. For negative voltages, we obtain where we see that now only the currents flowing from the substrate to the tip are nonzero and P σ 1/2 = 0. From the knowledge of the probabilities P σ n (E, V ) we can easily compute all the cumulants of the current distribution. Here, we shall focus on the analysis of the current and the noise, which can be obtained from the CGF of Eq. (44) using Eqs. (3) and (4), and are given by Let us now illustrate the results. We begin by analyzing the impact of the broadening or lifetime of the YSR states on the different transport properties. For this purpose, we choose the following system parameters: Γ S = 100∆, J = 80∆, U = 60∆, Γ t = ∆, and assume zero temperature. With these parameters the junction has a total normal-state conductance of G N = 0.025G 0 , where G 0 = 2e 2 /h. This corresponds to the tunnel regime in which the STM usually operates, and the corresponding YSR energy given by Eq. (33) is ϵ YSR = 0.41∆. In Fig. 3 we show the bias dependence of the differential conductance, shot noise and Fano factor for these parameters and for different values of the Dynes' parameter η of the substrate. The most salient feature in the conductance is the appearance of a peak (for both positive and negative voltages) exactly at the energy of the YSR states. The broadening of these peaks increases as η increases, as expected, and the height goes from being independent of the bias polarity for very small η to exhibit a very clear asymmetry when η is relatively large. In the case of the shot noise, see Fig. 3(b), the presence of the YSR states results in an abrupt increase of the noise at the energy of these states, while the value of η determines the rounding of the noise step. Finally, the impact of the YSR lifetime is most notable in the Fano factor, see Fig. 3(c). In this case, for very long lifetimes, F * exhibits a plateau for voltages below the YSR energy (e|V | < ϵ YSR ), while it adopts values very close to 1 for higher voltages. As the broadening of the YSR states increases (or their lifetime decreases), the Fano factor at low voltages progressively diminishes. Moreover, the dependence on the bias polarity also becomes more apparent. Notice also that values of F * < 1 become possible, in particular, inside the gap. Another thing that it worth mentioning is the absence of pronounced features at eV = ±∆ in all transport characteristics, contrary to what happens in the absence of YSR states. This is due to the fact that, due to the conservation of the number of states, the appearance of in-gap states is accompanied by the disappearance of the BCS gap edge singularities. This fact will become important in the next section when we compare with recent experimental results.
The previous results can be easily understood thanks to the unique insight that FCS offers us by identifying every individual tunneling process that contributes to the transport. In this regard, we show in Fig. 4 the results for the charge-resolved contributions to the differential conductance corresponding to the example of Fig. 3. In other words, we present in Fig. 4 the individual contributions to the conductance due to single-quasiparticle tunneling (G 1 ) and to the Andreev reflection (G 2 ). These processes are schematically represented in Fig. 2. The first thing to notice is that while the contribution of singlequasiparticle tunneling depends on the bias polarity, it is not the case for the Andreev reflection. Thus, any asymmetry in the transport characteristics must be produced by the contribution of the single-quasiparticle tunneling. The main message of this figure is that the Dynes parameter η determines the relative contribution of both tunneling processes: in the limit of large η the singlequasiparticle tunneling dominates the transport characteristics, while the Andreev reflection takes over in the opposite limit of very long-lived YSR states. This naturally explains the fact that the Fano factor is reduced upon increasing the broadening, which is simply due to the fact that in this case the transport is dominated by single-quasiparticle tunneling, see Fig. 2(a). This also explains the doubling of the Fano factor at low bias (e|V | < ϵ YSR ) in the limit of small η because in this case the Andreev reflection (transferring two electron charges) dominates the transport, see Fig. 2(b). A less trivial issue is to understand the origin of the abrupt jump of the Fano factor when eV = ±ϵ YSR and why in the limit of η → 0 the Fano factor gets much smaller than 2 in the voltage range ϵ YSR < e|V | < ∆ where the Andreev reflection completely dominates the transport. These in- With these parameters the junction has a normal-state conductance of 0.026G0 and the corresponding YSR energy is ϵYSR = 0.41∆.
teresting issues deserve a detailed explanation, as the one we are about to provide in the following paragraphs.
To clarify these issues we make use of the following analytical approximation for the probabilities of the two tunneling processes. Assuming zero temperature and focusing on energies close to the YSR energy, the probabilities for single-quasiparticle tunneling and the Andreev reflection for positive bias in Eq. (50) can be approximated as Lorentzians of the form where P max 1/2 are the energy-independent maxima of these two probabilities that is reached at E = ϵ YSR and W is the broadening of the YSR state in our model. If we assume electron-hole symmetry (U = 0) for simplicity, the factors P max 1/2 are given by where Notice that in the limiting case η → 0, it holds that Z = 0 and then, P max 2 = 1 and P max 1 = 0, i.e., there is no single-quasiparticle tunneling as expected. On the other hand, in the limit Γ L ≪ Γ R , J and η ≪ ∆, the broadening of the YSR states is given by With these approximate expressions, one can compute the current and shot noise in the voltage range ϵ YSR < e|V | < ∆. Considering first the ideal case of η = 0, where only the Andreev reflection contributes to the ingap transport, we obtain Thus, the corresponding Fano factor is F * = S/(2eI) = 1 in this voltage range, while it is easy to show that F * = 2 for e|V | < ϵ YSR (as long as η = 0). Thus, the abrupt reduction of the Fano factor at e|V | = ϵ YSR is a signature of the fact that the Andreev reflection becomes resonant because of the presence of the bound state, see Fig. 2(c). The occurrence of this resonant Andreev reflection can reduce the Fano factor all the way down to 1 inside the gap, see Fig. 3(c). This is due to the fact that the Andreev reflection can have a probability as high as 1 at the energy of the YSR states, which leads to a reduction of the Fano factor from 2 to 1 when crossing the bound state. In other words, as the Andreev reflection is resonant at the YSR energy, the transport is not longer in the Poissonian limit (with independent tunneling events) and the Fano factor cannot longer be interpreted as an effective charge. The fact that we observe in Fig. 3(c) that the Fano can be larger than 1 in the voltage range In this case, it can be shown that the probability of the Andreev reflection does not reach unity, see Fig. S2(d) in the Supplemental Material [71], and the Fano factor reduction upon crossing the YSR state is not complete, i.e., F * > 1. Actually, all of this is analogous to what happens in other resonant situations like in the case of a normal conductive double barrier structure. In that case, the Fano factor becomes 1/2 for a symmetric situation (two identical barriers), and it gets close to 1 for a very asymmetric case [47,64,65]. It is also worth mentioning that this crossover of the Fano factor, related to an Andreev reflection from 2 to 1 when crossing a bound state, has been reported theoretically in mesoscopic normal-superconducting structures [66]. To conclude this discussion, let us mention that one can also show using Eqs. (55) and (56) that, as long as U = 0, even in the case of a finite η, the zero-temperature Fano factor is equal to 1 in the voltage range ϵ YSR < e|V | < ∆.
Let us now analyze the evolution of the transport characteristics in the crossover between the tunneling regime and a highly transparent situation. For this purpose, we show in Fig. 5 the bias dependence of the differential conductance, shot noise and Fano factor for different values of the tunneling rate Γ t and Γ S = 100∆, J = 80∆, U = 60∆, η = 0.001∆, and T = 0. In this case the conductance evolves from exhibiting peaks at the YSR energies (see curve for Γ t = 0.1∆, which corresponds to G N = 0.0025G 0 ), to display a plateau inside the gap in which the conductance is close to 2G 0 (see curve for Γ t = 100∆, which corresponds to G N = 0.83G 0 ). This latter case essentially corresponds to a fully transparent standard (spin-degenerate) NS point contact. Notice that in this case the Fano factor has a complex evolution, namely it first increases at low bias upon increasing the transparency of the contact and then becomes sub-Poissonian (F * < 1) in the whole voltage range due to the reduction caused by the Pauli exclusion principle.
Again, these results can be easily rationalized making use of the charge-resolved contributions to the total differential conductance, which are displayed in Fig. 6 for this example. These results show that at the lowest transparency in this example the single-quasiparticle tunneling and the Andreev reflection give similar contributions. However, as the transmission increases the Andreev reflection completely dominates the transport inside the gap. For this reason, the evolution of the Fano factor can be solely understood from the evolution of the probability of the Andreev reflection. In particular, the sub-Poissonian Fano factor inside the gap for the highest transmission is simply due to the fact that the Andreev reflection reaches almost perfect transparency. Outside the gap, the Factor factors remains smaller than 1 due to the competition between the two tunneling processes that both give finite contributions. Finally, notice again that the contribution of the Andreev reflection is independent of the bias polarity and, thus, any asymmetry in the transport characteristics must be due to the contribution of single-quasiparticle tunneling.

VI. SINGLE-IMPURITY YSR JUNCTIONS: COMPARISON WITH SHOT NOISE MEASUREMENTS
Very recently, Thupakula and coworkers reported the first measurements of the shot noise through a magnetic impurity on a superconducting substrate featuring YSR states [49]. To be precise, these researchers employed shot-noise scanning tunneling microscopy to measure the nonequilibrium current fluctuations through a magnetic impurity deposited on 2H-NbSe 2 at temperatures around 0.7 K, way below the critical temperature of this superconductor. The STM tip was made of a normal metal and the main observation was the appearance of a Fano factor above 1. This fact was interpreted as the evidence of the contribution of an Andreev reflection to the charge transport. Moreover, those authors presented a theoretical analysis that suggested that the broadening of the YSR state was of the order of 1 µeV, which was clearly below the thermal energy (k B T ) in this experiment demonstrating that shot noise can probe energy scales that are not accessible in conventional tunneling spectroscopy measurements. The goal of this section is to provide a thorough analysis of those experimental results in the light of the theory presented in Sec. V.
In Fig. 7 we show an example of the experimental re- sults reported in Ref. [49], which corresponds to the conductance and Fano factor measured in the tail of the YSR states of a magnetic impurity. The most notable feature in the conductance is the appearance of two peaks inside the gap, which correspond to the YSR states. Notice also the asymmetry in the height of the peaks for positive and negative bias. With respect to the Fano factor, the main observation is the appearance of values above 1 for negative bias, but inside the gap region. This signature was taken as the evidence of the occurrence of Andreev reflections. Notice also that the Fano factor is asymmetric and for positive bias it remains below 1 for most voltages inside the gap. The Fano factor at very low bias was not reported simply due to the difficulty of measuring the very low currents in this voltage range. These experimental results were analyzed in Ref. [49] in the light of a model for the YSR state similar in spirit to ours, but just focusing on voltages inside the gap. To be precise, these authors used an approximation similar to that summarized in our Eqs. (55) and (56). Although such a simplified model qualitatively captures the physics of the YSR states, it is obvious that it cannot provide an overall correct picture of the experimental results. Actually, the major problem is that no model based on a single magnetic channel can describe the simultaneous appearance in the conductance of YSR peaks and pronounced coherent peaks as those in Fig. 7(a). As we explained above, these models do not predict the appearance of the coherent peaks simply due to the conservation of the number of states. So, as a way out to provide an overall consistent fit of the experimental results, we propose that there are at least two channels or pathways for the current: one due to the magnetic impurity and another non-magnetic channel, which probably results from the direct tunneling of the STM tip to the superconducting substrate. This is actually a solution that has been proposed before to explain the reported conductance spectra, see for instance Refs. [23,27]. There is another subtlety that we need to take into account, namely 2H-NbSe 2 is not a regular BCS superconductor. In fact, it has been shown by several groups that 2H-NbSe2 can be described as a two-band superconductor [20,67,68]. So, we propose here to explain the experimental results assuming that the transport takes place via two independent channels, one magnetic that is described by our YSR model of the previous sections and a non-magnetic channel that proceeds from the normal tip to the 2H-NbSe 2 substrate that we describe as a two-band superconductor. Moreover, we shall assume that this non-magnetic channel can be described in the tunneling regime, i.e., taking only into account the quasiparticle tunneling at the lowest order in transmission. We proceed now to describe the details of such a two-channel model.
First, for the non-magnetic channel we describe the superconductivity in 2H-NbSe 2 using the two-band model described in Ref. [68]. In this model one considers that the SC order parameter has two energy-dependent components ∆ 1 (E) and ∆ 2 (E) that can be computed by solving the following self-consistent set of equations where ∆ BCS 1/2 describes the bare SC gap of the separate bands. The density of states for the two bands can be computed as follows where the DOS at the Fermi energy is adjusted to fit the experimental data and α is a measure of the band anisotropy. The total DOS then follows Using the procedure explained in Ref. [20], we solved the algebraic system of Eqs. (61) and (62) and found that the best set of parameters is given by ∆ BCS 1 = 1.23 meV, Γ 12 = 0.27 meV, ρ 1 (E F ) = 1, ∆ BCS 2 = 0.29 meV, Γ 21 = 1.25 meV, ρ 2 (E F ) = 0.18 and α = 0.2. As mentioned above, we assume that the non-magnetic channel operates in the tunnel regime such that its transport properties are solely determined by single-quasiparticle tunneling, whose probabilities can be computed adapting Eqs. (45) and (46) to a non-magnetic situation in the tunneling regime, i.e., where ρ S is given by Eq. (64). We have made use of the fact that all the transmission coefficients are the same (for electrons and holes and for spin up and spin down) and equal to |t| 2 . This transmission coefficient was adjusted to fit the experimental results for the conductance and we obtained a value of |t| 2 = 0.0044. For the second, magnetic channel we simply used the theory of Sec. V and adjusted the different parameters to describe as well as possible the in-gap conductance due to the YSR states.
In particular, our best fit was produced with the following set of parameters: ∆ = 1.23 meV, Γ S = 123 meV, Γ t = 0.014 meV, J = 99.32 meV, U = 49.2 meV, η = 0.0011 meV, and we used the temperature of the experiment T = 0.7 K. The results of the best fit with our two-channel theory are shown in Fig. 7 alongside with the experimental results for both the conductance and Fano factor. As a reference, we also include the corresponding theory results obtained taking only into account the contribution of the magnetic channel. As one can see in panel (a), the theory captures very well the salient features of the differential conductance, namely the YSR peaks inside the gap and the coherent peaks related to the doublegap structure. Notice that the magnetic channel by itself reproduces well the YSR peaks, but it is unable to properly describe the conductance close to the gap edges and outside the largest gap.
Concerning the Fano factor, which was simply calculated from the parameters extracted from the fit of the conductance, there is an excellent agreement for positive voltages, see Fig. 7(b). The Fano Factor is super-Poissonian for voltages below the YSR energy and becomes sub-Poissonian for voltages higher than the bound state energy. For negative voltages there seems to be fine structure that we are not able to perfectly reproduce, but overall the agreement is quite satisfactory. In particular, we are able to reproduce the fact that the Fano factor is always bigger than 1 inside the gap and it tends to 1 for higher voltages. Here, it becomes even more apparent that a two-channel model is necessary. The contribution of the magnetic channel alone does not reproduce very well the structure of the Fano factor and it supports our hypothesis on the need of an additional contribution. Let us also say that at very low bias (not shown here) the Fano factor becomes very large simply because the current fluctuations are dominated by a finite thermal noise, while the current tends to zero.
Something that is very important to emphasize is the fact that, as we showed in Sec. V, the Fano factor is very sensitive to the broadening or lifetime of the YSR states. In our fit we obtained a value of η = 1.1 µeV for the Dynes' parameter in the superconducting substrate. Using Eq. (58) and the rest of the parameter values extracted from the fit, we obtain that W ≈ η = 1.1 µeV, which is much smaller than the thermal energy in this case (k B T ≈ 60 µeV). Thus, as stated in Ref. [49], the analysis of the noise and the corresponding Fano factor allows us to have access to energy and time scales that are usually out of the scope of conventional conductance measurements due to thermal broadening.

VII. SINGLE-IMPURITY YSR JUNCTIONS: SUPERCONDUCTING TIP
We now analyze the case in which a magnetic impurity is coupled to two superconducting leads, which for simplicity we assume to have the same gap ∆. In this case, the novelty with respect to the previous case is the possibility of having MARs. As a reference for our discussions below, we illustrate the relevant tunneling processes in this case in Fig. 8.
Technically speaking, this case is considerably more complicated than the NS one because the GFs of the terminals are not longer diagonal in energy space, which is due to the ac Josephson effect. This means in practice that we have to treat the problem in the Floquet language, as explained in Appendix A. For this reason, it is not possible to provide a complete analytical solution for the CGF and we have to resort to numerics. The technical details are presented in Appendix A,B and here we shall focus on the discussion of the physical results.
As discussed in Sec. III, in this case the CGF reads where P σ n (E, V ) corresponds to the probability of transferring n charges across the junction for spin σ and E is the Floquet energy. Notice that the main difference with respect to Eq. (44) is the fact that n now runs from −∞ to ∞ due to the occurrence of MARs and we integrate over a finite energy range. Again, from the knowledge of the probabilities P σ n (E, V ) we can easily compute the current and the noise, which are given by the standard formulas of a multinomial distribution We begin by analyzing the impact of the broadening or lifetime of the YSR states on the different transport properties. For this purpose, we choose the following system parameters: Γ S = 100∆, J = 80∆, U = 60∆, Γ t = 5∆, η t = 0.001∆ and assume zero temperature. With these parameters the junction has a normal-state conductance of G N = 0.12G 0 and the corresponding YSR energy given by Eq. (33) is ϵ YSR = 0.41∆. Figure 9 displays the conductance, shot noise, and Fano factor for these parameters and different values of the Dynes parameter of the substrate η S . Notice that the conductance exhibits a very rich subgap structure due to the occurrence of MARs. In any case, the most visible conductance peaks appear at eV = ±(∆ + ϵ YSR ), which as we show below are due to both single-quasiparticle tunneling and the lowest-order resonant Andreev reflection. Notice also that the conductance depends on the bias polarity because we are dealing with a situation with electronhole asymmetry (U = 60∆). The additional conductance peaks appear at eV = ±2∆/n, which are due to conventional (non-resonant) Andreev reflections, and at eV = ±(∆ + ϵ YSR )/n, whose origin will be discussed below. Overall, the role of the Dynes' parameter is to determine the width and the height of all these conductance peaks, as expected. On the other hand, the shot noise exhibits steps at the voltages at which the conductance peaks appear and these steps are progressively more pronounced as η S decreases, see Fig. 9(b). Again, this parameter has a strong impact in the Fano factor, see Fig. 9(c), which now exhibits super-Poissonian values well above 2 as η S decreases. This is obviously a signature of the occurrence of MARs. It is also important to notice that the Fano factor does not exhibit integer values inside the gap, as it occurs in the absence of in-gap states [69,70], which suggests that no tunneling process completely dominates the transport at any subgap voltage. Let us also say that we do not report the results at very low bias because the current becomes exceedingly small and it is not measurable in practice.
Again, we can rationalize these results with the analysis of the individual contributions of the different tunneling processes classified according to the charge they transfer. In Fig. 10, the charge-resolved currents are shown for different values of η S (we plot the absolute value of those currents for clarity). The corresponding charge-resolved conductances are shown in the Supplemental Material (Fig. S7) [71]. The light grey line in-  dicates the total current, whereas the others correspond to the contributions I n due to the processes transferring n electron charges. Notice that for the smallest value of the Dynes' parameter, η S = 0.1∆ in panel (a), the singlequasiparticle tunneling and Andreev reflection dominate the transport and the MAR contributions with n > 2 are negligible for all voltages. However, the resonant Andreev reflection illustrated in the upper graph of Fig. 8(c) gives a sizeable contribution to the conductance peak at eV = ±(∆ + ϵ YSR ) = ±1.41∆. Decreasing η S results in the suppression of the quasiparticle current (I 1 ) because the SC DOS in the gap region is lowered. For η S = 0.01∆ in panel (b), we see that the lowest-order Andreev reflection takes over at voltages ∆ < e|V | < 2∆ and it is responsible for the conductance peaks at eV = ±∆ (non-resonant Andreev reflection, see Fig. 8(b)) and at eV = ±(∆ + ϵ YSR ) = ±1.41∆ (resonant Andreev reflection, see Fig. 8(c)). Notice, on the other hand, that the third-order MAR starts to play a fundamental role in the subgap transport and it is responsible for the conductance peak at eV = ±(∆ + ϵ YSR )/2 = ±0.705∆ and eV = ±(∆ + ϵ YSR )/3 = ±0.47∆, this latter one exhibiting moreover a pronounced negative differential conductance. The conductance peak at eV = ±(∆ + ϵ YSR )/2 = ±0.705∆ is extremely interesting because it has been observed in a number of experiments and its origin has been attributed to a second-order Andreev reflection that starts or ends in a YSR state (transferring 2 charges) [12,16,26,30]. This process, which we term YSRmediated Andreev reflection, is illustrated in the upper diagram of Fig. 8(d) and, it has in fact a threshold voltage equal to eV = ±(∆ + ϵ YSR )/2. However, this interpretation is incorrect in our case, as it is evidenced by the charge-resolved currents and by the fact that the Fano factor is clearly above 2 in this voltage region. Such a conductance peak originates indeed from a MAR of order 3 that becomes resonant precisely at eV = ±(∆ + ϵ YSR )/2, as we illustrate in the lower diagram of Fig. 8(c). It is easy to show that such a resonant MAR requires the YSR energy to fulfill ϵ YSR ≥ ∆/3, which is satisfied in this example. Another fact that confirms our interpretation is the absence of negative differential conductance at that bias, which would be expected from a YSR-mediated Andreev reflection, but not from a resonant MAR. This discussion illustrates again the unique insight of the FCS approach, without which it would be hard to draw the right conclusion in this case.
Coming back to the features at eV = ±(∆+ϵ YSR )/3 = ±0.47∆, they also originate from a third-order MAR, but in this case it is a MAR process that involves the tunneling into an YSR state (i.e., a YSR-mediated MAR), which explains the negative differential conductance. This process is illustrated in the lower graph of Fig. 8(d). For η S = 0.001∆ in panel (c), the current steps are very abrupt leading to very high conductance peaks. In this case, the conductance peaks at eV = ±(∆ + ϵ YSR )/3 can now be mainly attributed to the fifth-order YSRmediated MAR. For the smallest η S = 0.0001∆ in panel (d), the total current exhibits pronounced negative differential conductance whenever the tunneling into a YSR state is involved. An important observation in Fig. 9 is that, contrary to the case of a normal tip, the contribu- tion of the different Andreev reflections does depend on the bias polarity. Thus, asymmetries in the the currentvoltage characteristics cannot be solely attributed to the contribution of single-quasiparticle tunneling.
Let us explore now the impact of the junction transmission by changing the tip tunneling rate Γ t , while maintaining constant the other parameters: Γ S = 100∆, J = 80∆, U = 60∆, η S = 0.001∆, η t = 0.001∆, and T = 0. In Fig. 11, the conductance, shot noise, and Fano factor are shown as a function of the voltage and for different tip tunneling rates Γ t . Considering the conductance, for the smallest value of Γ t = 0.1∆ (G N = 0.0025G 0 ), the first YSR resonance at eV = ±(∆ + ϵ YSR ) dominates the subgap structure. Increasing Γ t firstly broadens the conductance peaks, and secondly shifts the position of the conductance peaks slightly. This is caused by the renormalization of the YSR energies due to the finite coupling to the tip. For Γ t = ∆ (G N = 0.025G 0 ), the first YSR resonance at eV = ±(∆ + ϵ YSR ) becomes increasingly pronounced and additional subgap structure starts to appear, namely at the Andreev reflection onsets at eV = ±∆ but also at the onset of the YSR-mediated and resonant Andreev reflections eV = ±(∆ + ϵ YSR )/n (with n > 1). For Γ t = 5∆ (G N = 0.12G 0 ), the conductance peaks are increasingly broadened and, in particular, the peaks at eV = ±(∆ + ϵ YSR )/3 become clearly visible. For Γ L = 10∆ (G N = 0.22G 0 ), the YSR peaks at eV = ±(∆ + ϵ YSR ) are almost completely washed out. The peaks at eV = ±∆ are clearly visible and corre-spond to the onset of the regular Andreev reflection. At eV = ±(∆ + ϵ YSR )/3, we see high conductance peaks that correspond to YSR-mediated Andreev reflections. If one continued increasing the coupling (not shown here), the system would resemble a highly transparent superconducting point contact [72], and the YSR resonances are no longer resolvable. On the other hand, the shot noise exhibits steps at voltages corresponding to onsets of the different types of Andreev processes. Concerning the Fano factor, for large values of Γ t = 5, 10∆ it exhibits very high super-Poissonian values due to the occurrence of MARs. In addition, we see the characteristic drop-off of the Fano factor from roughly 2 to almost 1 at the first YSR resonance at eV = ±(∆ + ϵ YSR ), which has the same origin as in the NS case just shifted by the gap energy ∆. Decreasing the tip coupling decreases the Fano factor because the MAR contributions are suppressed. Notice that for small Γ t = 0.1∆, the Fano Factor never reaches over 2 indicating that MARs do not contribute to the transport. In addition to the characteristic drop-off at eV = ±(∆ + ϵ YSR ), there are other signatures in the Fano factor which correspond to the MAR onsets at eV = ±(2∆/n) and to the onset of both resonant MARs and YSR-mediated Andreev reflections at eV = ±(∆ + ϵ YSR )/n. Again, we can pinpoint the origin of every feature in the transport characteristics by considering the chargeresolved currents, as we illustrate in Fig. 12 (see Fig. S8 for the corresponding charge-resolved conductances [71]). Notice that for the smallest coupling Γ t = 0.1∆ in panel (a) the single-quasiparticle tunneling and the lowestorder Andreev reflection dominate the transport and give similar contributions at eV = ±(∆ + ϵ YSR ). However, as the transmission increases, the quasiparticle tunneling becomes progressively more irrelevant in the subgap transport. Thus, for instance, for Γ t = ∆, the structure at eV = ±(∆ + ϵ YSR ) is solely due to the resonant Andreev reflection. For Γ t = 5∆ in panel (c), the subgap structure becomes much more pronounced with current steps not only at the standard MAR onsets eV = ±2∆/n, but also at voltages eV = ±(∆ + ϵ YSR )/2 and eV = ±(∆ + ϵ YSR )/3 which, as explained above, are due to resonant MARs and to YSR-mediated Andreev reflections, respectively. Thus, for instance, the structure at eV = ±(∆ + ϵ YSR )/2 is mainly due to the third-order MAR, whereas several YSR-mediated MARs contribute to the structure at eV = ±(∆ + ϵ YSR )/3. The steps in the charge-resolved currents are smoothed out when Γ t is further increased towards Γ t = 10∆ in panel (d).
In this case the first YSR resonance becomes so broad that it can barely be identified as a resonance, whereas the higher-order YSR resonances are still very sharp and become increasingly easy to detect.
To conclude this section, let us say that one can provide some analytical insight in the tunnel regime when Γ t ≪ Γ S . In this limit, and with the help of Ref. [62], one can derive perturbative expressions for the current contributions of single-quasiparticle tunneling (I 1 ) and the lowest-order Andreev reflection (I 2 ). These expressions are given by Here, ρ t (E) is the BCS DOS of the superconducting tip, ρ imp (E) is the total DOS of the impurity coupled to the superconducting substrate, and F (E) is the corresponding anomalous Green's function of the impurity coupled to the substrate. The impurity DOS is given by where and the impurity anomalous Green's function is given by By construction, these perturbative current formulas are only valid in the tunnel regime, but they also require the broadening of the YSR states to be large in comparison with the tip tunneling rate Γ t . To test these approximate formulas, we present in Fig. 13 a comparison with the exact numerical results for the case in which Γ S = 100∆, J = 80∆, U = 60∆, η t,S = 0.01∆, Γ t = 0.01∆, and T = 0. Notice that the analytical formulas nicely reproduce the numerical results, except for the Andreev reflection contribution at very low bias when higher-order terms in the transmission are expected to play a role.

VIII. DOUBLE-IMPURITY JUNCTIONS: TUNNELING BETWEEN YSR STATES
As mentioned in Sec. I, recently it has been experimentally demonstrated that a superconducting STM tip can be functionalized with a magnetic impurity that then features YSR states [25]. More importantly, it was shown that this YSR-STM can be used to probe other magnetic impurities deposited on a superconducting substrate that also features YSR states. In this way, these experiments realized for the first time the tunneling between individual superconducting bound states at the atomic scale, which is the ultimate limit for quantum transport. In particular, the current-voltage characteristics in these double-impurity junctions were shown to exhibit huge current peaks inside the gap (with an extremely pronounced negative differential conductance). These current peaks have been interpreted as the evidence of tunneling between YSR states (both direct at low temperatures and thermally excited at high temperatures) [25,26,73]. In this section we want to analyze this unique situation from the FCS point of view and, in particular, provide very concrete predictions for the shot noise and Fano factor in these junctions. In turn, this problem gives us the opportunity to show how the FCS approach can be used in a case in which the scattering matrix is not diagonal in spin space.

A. Model and scattering matrix
In the spirit of the Keldysh action of Eq. (9) we now need a scattering matrix describing these double impurity junctions. For this purpose, we follow the model put forward in Ref. [73], which has been very successful FIG. 14. Schematic representation of our two-impurity model. Two magnetic impurities are respectively coupled to a superconducting substrate and to an STM tip that is also superconducting. The tunneling rates Γt and ΓS measure the strength of the coupling of the impurity to the tip and substrate, respectively, and v is the hopping matrix element describing the tunnel coupling between the impurities. These impurities have magnetizations Jt and JS whose relative orientation is described by the angle θ.
describing the experimental observations for the currentvoltage characteristics [25,26]. We briefly describe this model and then proceed to the determination of the corresponding normal-state scattering matrix.
The double-impurity model of Ref. [73] is schematically illustrated in Fig. 14. In this model, and inspired by the experiments of Ref. [25], every impurity is strongly coupled to a given superconducting lead (tip and substrate) and both impurities are coupled via a tunnel coupling. The Hamiltonian describing this model is given by [73] H = H lead,t + H imp,t + H hopping,t + H lead,S + H imp,S + H hopping,S + V, where H lead,j describes the Hamiltonian of lead j = t, S, H imp,j describes the respective impurity which is coupled to lead j via the Hamiltonian H hopping,j and V describes the coupling between the two impurities. We choose a global z-quantization axis in which to describe the creation and annihilation operators wherec † k,j corresponds to lead j andd † j to impurity j. We can express the different terms of the Hamiltonian in their respective basis as follows where the respective matrices in spin-Nambu space take the formH with the on-site energies U j and the exchange energies J j . The tunneling term between the two magnetic impurities reads with the coupling matrix where v is the tunneling coupling between the impurities. To simplify things, it is convenient to transfer the dependence on θ j and φ j to the coupling term V in Eq. (75) and work with Hamiltonians describing the subsystems in which the corresponding spin points along its quantization axis. For this purpose, we introduce the combined unitary transformationR j = e ıθj σ2/2 ⊗ e −ıφj τ3/2 , where θ j is the angle formed between the exchange energy J j and the global z-quantization axis. Upon performing this unitary transformation, the Hamiltonian matrices of Eqs. (81)-(83) become now while the coupling matrices in Eq. (85) adopt now the formV Notice that these coupling matrices effectively describe a spin-active interface in which there are spin-flip processes whose probabilities depend on the relative orientation of the impurity spins described by θ.
To obtain the scattering matrix, we shall make use again of the Fischer-Lee relations. For this purpose, we first notice that the central region is now twice as big as in the single-impurity case, and the corresponding Hamiltonian describing this region reads in the rotated basis (d † t , d † S ) ≡ (R td † t ,R Sd † S ). As in previous cases, only the electron-components of the matrix representation are needed for the computation of the scattering matrix, namely the first and third component of the spinor basis in Eq. (77). Explicitly, that part of the Hamiltonian of the central system reads . Notice that we have included the bias voltage V in the tip subsystem as a shift of the corresponding impurity level. This way we describe a situation in which the voltage entirely drops between the two impurities, as it happens in the STM experiments. On the other hand, the couplings can be expressed as where the matrices V j,e are given by With these couplings, the self-energies of the tip and substrate can be readily computed as we are only interested in the normal state scattering matrix, and thus the lead GFs are normal metal ones. In particular, the (normalstate) electron part of the GFs are given by g so they are 4×4 matrices in impurity-spin space and they are not of full rank. The dressed central GFs can then be computed from the Dyson equation with the regularization factor η imp which can again be chosen as zero, as in the single impurity case. To compute the components, we also need the scattering rate matrices Γ t ≡ ℑ(Σ a t,e ) and Γ S ≡ ℑ(Σ a S,e ). Finally, the scattering matrix then follows from the Fisher-Lee relation using the 4 × 4 dressed central GFs, namely where ρ 3 is the third Pauli matrix in lead space. For the hole part of the scattering matrix, it holds that and the total scattering matrix reads

B. Results
The calculation of the Keldysh action in this doubleimpurity case is very similar to that of a single impurity discussed in Sec. VII. Again, we have to treat the problem in the Floquet language and resort to numerics to compute the charge-resolved probabilities. On a conceptual level, the main differences are that the scattering matrix has a different energy dependence and it is non-diagonal in spin space. Therefore, it can be shown that the CGF in this case reads where P n (E, V ) corresponds to the probability of transferring n charges across the junction and E is the Floquet energy. Notice the absence of the spin index, when compared to Eq. (67). We can then compute the current and the noise, which are given by the standard formulas of a multinomial distribution To make a connection with the experiments of Refs. [25,26], we shall mainly focus here on the analysis of the results in the regime of weak coupling between the impurities in which the transport properties are mainly determined by the tunneling of quasiparticles. Moreover, we shall assume that the two superconducting electrodes are identical and the corresponding gap is given by ∆.
To illustrate the results in the weak-coupling regime, we consider an example in which the different model parameters are given by: Γ S,t = 100∆, J t = 60∆, U t = 0, J S = 60∆, U S = 60∆. With these parameters, the YSR energies are ϵ t = 0.48∆ for the SC tip and ϵ S = 0.64∆ for the substrate. The Dynes' parameter is chosen to be the same for both SCs η t,S = 0.001∆. Moreover, we choose the impurity coupling v = ∆ such that the normal state conductance is relatively low (≈ 4 × 10 −4 G 0 ) and the quasiparticle tunneling dominates the transport. Finally, the spin-mixing angle is chosen as θ = π/2, which has been shown to reproduce the experimental results [26]. Actually, this value was interpreted as a way to describe the average transport properties in a situation in which the two spins are freely rotating, see Ref. [26] for details. The results for the current, shot noise, and Fano factor as a function of the bias voltage are shown in Fig. 15 for two cases: zero temperature (k B T = 0) and k B T = 0.2∆. Let us first discuss the zero-temperature results. The main salient feature in the current is the appearance of very pronounced peaks (with a huge negative differential conductance) inside the gap region at voltages given by eV = ±|ϵ S + ϵ t | ≈ ±1.12∆, see Fig. 15(a), which reproduces the main observation in Refs. [25,26]. As explained in Refs. [25,26,73], these current peaks are due to the resonant quasiparticle tunneling between the lower and upper YSR states of the impurities in the tip and substrate, as we illustrate in Fig. 15(g). Such a tunneling process, which we shall refer to as direct Shiba-Shiba tunneling, can occur at any temperature and it only requires the impurity spins to be nonparallel (i.e., θ ̸ = 0), otherwise it would be forbidden due to the full spin polarization of the YSR states. Notice that in Fig. 15(a) we compare the exact result taking into account any possible contribution (also from Andreev reflections) and the result only including single-quasiparticle tunneling (n = ±1). Such a comparison shows that the current is dominated in this case by quasiparticle tunneling. With respect to the shot noise, see Fig. 15(b), it exhibits the same voltage dependence as the current and it shows two pronounced peaks at eV = ±|ϵ S + ϵ t | as a result of the direct Shiba-Shiba tunneling. Again, the noise is dominated by the contribution of single-quasiparticle tunneling. Turning now to the Fano factor, see Fig. 15(c), the main feature is the strong reduction at the voltages of the direct Shiba-Shiba tunneling. In this particular case, the Fano factor reaches a value around 0.5 that actually depends on the bias polarity. Notice that the Fano factor reduction is also dominated by single-quasiparticle tunneling. Obviously, the origin of this pronounced reduction of the Fano factor must be due to the resonant character of this tunneling process between the two individual YSR states. This will be explained below. Notice, on the other hand, that the Fano factor is close to 1 away from the resonant voltage and it reaches values above 1 inside the gap due to the contribution of Andreev reflections.
Turning now to the high temperature case (k B T = 0.2∆), the main novelty in the current is the presence of two additional current peaks inside the gap that appear at eV = ±|ϵ t − ϵ S | ≈ 0.16∆, see Fig. 15(d). These current peaks, which were in fact experimentally observed [25,26], are due to the thermally activated tunneling between the two upper (or excited) YSR states, as we illustrate in Fig. 15(h). As explained in Refs. [25,26,73], this thermally activated tunneling process, which we shall refer to as thermal Shiba-Shiba tunneling, can occur when the temperature is high enough to have a partial occupation of the excited YSR states and it requires that the impurity spins are not antiparallel (i.e., θ ̸ = π), otherwise it would be forbidden due to the full spin polarization of the YSR states. Again, we note that the current in this case is also dominated by single-quasiparticle tunneling, see Fig. 15(d). The thermal Shiba-Shiba tunneling is also visible in the shot noise in the form of two peaks at eV = ±|ϵ t − ϵ S |. In the case of the Fano factor, the values at the voltages eV = ±|ϵ t + ϵ S | associated to direct Shiba-Shiba tunneling have increased as compared to the zero-temperature case, while there is no visible feature at the biases of the thermal Shiba-Shiba tunneling because in this example it is masked by the rapid increase of the thermal noise that dominates the low-bias regime in the Fano factor at finite temperatures.
From the discussion of the example of Fig. 15 it is obvious that the most important feature related to the tunneling between two YSR states is the Fano factor reduction at the resonant bias of the direct Shiba-Shiba tunneling. To understand its origin and magnitude we first need to analyze more systematically how this reduction depends on the different parameters of the model. In particular, the Fano factor related to the direct Shiba-Shiba tunneling is expected to depend on the relative value between the effective tunneling rate (related to the hopping element v) and the broadening (or lifetime) of the bound states. In fact, this issue is very much related to another interesting observation reported in Ref. [25], namely the fact that the height of the direct Shiba-Shiba current peaks (and their area) undergoes a crossover be-  (d-f) The same as in panels (a-c), but for a temperature kBT = 0.2∆ with additional dotted vertical lines indicating the position of the voltages eV = ±|ϵS − ϵt| of the thermal Shiba-Shiba tunneling. (g,h) Schematic representation of the tunneling processes between the two YSR states. Here, the left electrode is a magnetic impurity coupled to a SC tip and the right one is another impurity coupled to a superconducting substrate both featuring YSR states. The green shaded areas represent the corresponding Fermi functions of both electrodes. The diagrams follow the same convention as in Fig. 2. Panel (g) corresponds to the direct Shiba-Shiba tunneling enabled by a finite spin-mixing angle θ and panel (h) to the thermal Shiba-Shiba tunneling enabled by a finite temperature and a spin-mixing angle θ ̸ = π. tween a linear regime at very low transmission (or normal state conductance) and a sublinear regime at higher transmission when the STM tip with its impurity was brought closer to the impurity on the substrate. To understand the impact of the junction transmission, we considered the zero-temperature example of Fig. 15 and computed the value of the current, noise and Fano factor at the direct Shiba-Shiba voltage eV = |ϵ S + ϵ t | as a function of the hopping matrix element v describing the coupling between the impurities, while the rest of the parameters are kept constant. The results are displayed in Fig. 16, where we show both the exact results and those computed taking only into account the contribution of single-quasiparticle processes. As it can be seen in panel (a), the peak height indeed exhibits the type of crossover observed experimentally. For small values of v, i.e., deep into the tunnel regime, the current peak height increases linearly with the transmission (which in this regime is proportional to v 2 ) and it crosses over to a sublinear regime at higher couplings. Moreover, this crossover is mainly about the quasiparticle current and only at high enough couplings (v ≳ ∆) we start to see additional contributions due to Andreev reflections. Interestingly, the shot noise exhibits a similar crossover, but with an important difference, namely the fact that it enters the sublinear regime more quickly than the current and even exhibits a local minimum at a certain value of v. This peculiar behavior is clearly reflected in the evolution of Fano factor, see Fig. 16(b). The Fano factor for very small couplings is equal to 1, as expected from standard (non-resonant) tunneling, it is progressively reduced as the coupling increases reaching a minimum value of around (∼ 0.22), and it increases again monotonically for even higher couplings. Notice also that the quasiparticle contribution dominates the Fano factor up to coupling values clearly beyond that at which the minimum takes place.
To understand the origin of this crossover, the existence and magnitude of the Fano factor minimum, and whether this behavior is universal, we put forward the following toy model. We consider two energy levels ϵ t and ϵ S that are coupled to their respective electron reservoirs. Due to this coupling, the quantum levels acquire a broadening, which we assume to be equal for both levels and given by η S = η t = η. Additionally, these two levels are coupled via a tunnel coupling v eff . This model is illustrated in Fig. 17(a). Notice that we do not even need to invoke the existence of superconductivity. The energy and bias dependent electron transmission in this toy model is given by [63] T where g t,S (E) are advanced GFs related to both energy levels and given by and ρ t/S (E) = (1/π)ℑ(g L/R ) are the DOS associated to those two levels. Without loss of generality, we choose ϵ t < 0, ϵ S > 0, and a positive bias. Since we want to emulate the transport properties at the direct Shiba-Shiba energy, we shall set from now on eV = |ϵ t | + |ϵ S |. At this bias, the transmission adopts the form where γ = v eff /η. Within this model, the corresponding zero-temperature current and noise are given by Assuming that η ≪ |ϵ t,S | and v eff < |ϵ t | + |ϵ S |, one can compute analytically the previous integrals to obtain The corresponding Fano factor is then given by These expressions nicely summarize the type of crossover that we discussed above and show that it is simply controlled by the parameter γ, i.e., ratio between the tunnel coupling and the level width. Before comparing with the actual results of the two impurity system, let us analyze the Fano factor in different limiting cases. First, in the weak coupling regime γ ≪ 1, the Fano factor becomes 1, which is the expected result for non-resonant tunneling situations. Second, in the high coupling regime γ ≫ 1, the Fano factor tends to 1/2, which is the generic expected result for two very narrow energy levels. We shall see that this regime cannot be easily achieved with the YSR states because for very large v the Fano factor is altered by the Andreev reflections. More interestingly, the Fano factor can be much smaller than 1/2. Thus, for instance, for γ = 1 the Fano factor is exactly 1/4. In particular, the shot noise exhibits a local minimum at this value of γ, which we propose as a new method to extract the intrinsic lifetime of the states. However, the Fano factor of 1/4 does not mark the lowest Fano factor. It is easy to show that the minimum Fano factor value is achieved when γ = 5/3 > 1, which results in a Fano factor equal to 7/32. This is a unique value that, to our knowledge, is not found in any other generic situation.
To establish a quantitative comparison between the toy model and the direct tunneling between YSR states, we still need to identify the effective tunneling coupling in the latter situation. In the presence of YSR states the effective coupling must take into account both the spin-mixing angle θ and the coherent factors that determine the height of the DOS associated to the YSR states. These coherent factors are given in our model by (115) with j = t, S. Thus, the effective tunnel coupling v eff describing the direct Shiba-Shiba tunneling can be expressed in terms of the bare v hopping matrix elements as follows Finally, we are in position to establish the desired comparison, which is shown in Fig. 17(b,c). There we present both the analytical results obtained with the toy model for the current, noise, and Fano factor at the resonant bias, as well as the exact results shown in Fig. 16(a,b) where the tunnel coupling has been renormalized according to Eq. (116). The results are presented as a function of the ratio v eff /η. Notice that the toy model is able to quantitative reproduce the results of our two-impurity model over a broad range of values of that ratio. In particular, the toy model nicely reproduces the crossover in all transport properties up to coupling values in which Andreev reflections start playing a role. Interestingly, we now can see that the local minimum in the shot noise exactly corresponds to the case v eff = η when the effective coupling is equal to the broadening of the YSR states. Moreover, the minimum of the Fano factor is shown to reach the value 7/32 predicted by the toy model, and the Fano factor becomes 1/4 exactly when the shot noise is locally minimal. On the other hand, notice that the prediction of the toy model of a Fano factor equal to 1/2 when v eff ≫ η is not reached in the direct Shiba-Shiba tunneling case due to the onset of the Andreev reflections, as one can see in Fig. 17(c). Finally, the most important conclusion of the impressive agreement with the toy model is the universality of the results when considered as a function of the the ratio v eff /η, as long as the transport is dominated by the tunneling of single quasiparticles. This universality is illustrated in Ref. [71] where we present the results for this crossover for different sets of values of the model parameters, see Fig. S15.
It is important to stress that our FCS approach is valid for arbitrary junction transparency and can also describe situations in which the transport is dominated by MARs. Actually, this is a very interesting subject since, as it has been theoretically shown in Ref. [73] for the current, the transport properties are expected to exhibit an extremely rich subgap structure due to the occurrence of a variety of different types of MARs, some of which have no analog in the single-impurity case. However, due to the richness of that physics, and to keep the length of this manuscript under control, we shall postpone the analysis of that regime to a forthcoming publication. In any case, we have included some examples of the high transparency regime in Ref. [71] to illustrate once more the power of our FCS approach, see Figs. S10-S14.

IX. DISCUSSION AND CONCLUSIONS
As explained in Sec. I, the concept of FCS has already been used to understand the transport properties of several basic superconducting junctions. However, the Keldysh action described in Sec. III [60,61], which has remained fairly unnoticed, puts this concept at a whole new level because it allows us to describe the coherent transport in situations where the scattering matrix of the system may depend on energy, spin, and even involve an arbitrary number of superconducting terminals. Here, we have shown how this action can be used in practice by combining it with model Hamiltonians, and we have illustrated its power in the several examples concerning the spin-dependent transport in hybrid superconducting systems that exhibit YSR states. Such a combination opens up the possibility to shed new light on the coherent transport in numerous superconducting nanostructures. Thus, for instance, among the natural extensions of this work we can mention the study of Majorana physics in magnetic atomic chains from the FCS point of view [34][35][36][37][38][39][40]. In fact, it has been recently suggested that shot-noise tomography in hybrid magnetic-superconducting wire systems could be used to distinguish Majorana modes from other trivial fermionic states [74]. Another context in which one could straightforwardly apply our FCS approach is that of superconductor-semiconductor nanowire junctions in which one can realize single and double quantum dot systems that exhibit YSR states, see Refs. [75][76][77] and references therein. An interesting situation in which it would be very interesting to apply the FCS approach is that of microwave-irradiated junctions [78][79][80]. It has been recently shown that one can gain some new insight into the interplay between YSR states and Andreev transport by studying the single-impurity systems discussed in this work with the assistance of a microwave field [81][82][83]. Probably the most exciting possibility of the approach put forward in this work is the analysis of multiterminal superconducting systems. As mentioned above, the Keldysh action described here is valid for arbitrary number of terminals. In this sense, we are very much interested in extending this work to study different aspects of superconducting multiterminal systems that are currently attracting a lot of attention such as the generation of Cooper quartets [84][85][86][87][88], the study of the Josephson effects [89][90][91][92], or the possibility to engineer Andreev bound states with interesting topological properties [93][94][95][96][97][98][99][100][101][102][103][104]. Last but not least, nothing prevents us from using this FCS approach in the case of junctions involving unconventional superconductors.
It is worth remarking that our whole analysis of the FCS in these impurity systems has been done making use of mean-field models. In that regard, it would be interesting, albeit very challenging, to investigate the role of electronic correlations or spin fluctuations in the transport properties discussed in this work, most notably in the noise. This clearly goes beyond the scope of this work and it cannot be done with the action of Eq. (9) used here as starting point. The success of these mean-field models is indeed quite remarkable, as we have shown in Sec. VI and has been demonstrated in numerous publications [9, 23, 25-27, 73, 81-83]. So, it will be interesting to see if future shot noise experiments reveal the presence of significant correlation effects.
To summarize the present work, let us say that we have shown here that the FCS approach can be extended to describe the spin-dependent transport in systems involving magnetic impurities coupled to superconducting leads. In particular, with the analysis of different situations we have illustrated how the concept of FCS provides an unprecedented insight into the interplay between YSR states and electronic transport. Among the lessons and predictions put forward here, we can highlight the following ones. First, in the case of single-impurity junc-tions with only one superconducting electrode, we have shown that the whole subgap transport can be understood as a competition between single-quasiparticle tunneling and a resonant Andreev reflection. Such a competition is especially reflected in the shot noise and Fano factor, which allow us to access energy scales that are out of the scope of conventional conductance measurements, most notably the YSR lifetimes. In particular, we have illustrated this fact with the analysis of very recent experiments that have reported for the first time shot noise measurements in these hybrid atomic-scale systems [49]. Second, in the case of single-impurity junctions with two superconducting leads, we have shown how the FCS concept allows us to unambiguously identify the contribution of all the tunneling processes, including multiple Andreev reflections. In particular, this has helped us to correct common misinterpretations on the origin of certain features in the subgap structure of the conductance. In particular, we have discussed the unique signatures of the so-called resonant multiple Andreev reflections (multiple versions of the known lowest-order resonant Andreev reflection), which should enable their experimental identification, something that to our knowledge has not been reported thus far. Moreover, we have presented extensive predictions on how the occurrence of these resonant Andreev reflections, and YSR-mediated Andreev reflections, are reflected in the shot noise and Fano factor. Finally, in connection to recent experiments on the tunneling between two individual YSR states in two-impurity systems [25,26], we have shown that such a tunneling leads to an unambiguous signature in the Fano factor in the form of a strong reduction at the resonant bias voltage at which the two YSR states are aligned. In particular, we have demonstrated that the direct Shiba-Shiba tunneling at low temperatures can exhibit a Fano factor as small as 7/32, which results from the resonant quasiparticle tunneling between two YSR states. This constitutes a novel result in quantum transport that has not been realized in any other system. So, in short, we think that these examples and predictions will motivate other theoretical groups to use the concept of FCS to provide a new point of view on the superconducting transport in numerous situations. Moreover, our work clearly shows the importance of going beyond conductance measurements to truly understand the Andreev transport in the presence of superconducting bound states. In this sense, we are convinced that this work will trigger off the realization of new shot noise measurements in these atomic-scale superconducting systems exhibiting YSR states.
in Nambu space result in a dependence of the GF on the average time t av = (t + t ′ )/2. We can thus express the GF as a function of relative time t rel and average time t av = (t + t ′ )/2, namelyǧ(t rel , t av ). For a normal metal, the off-diagonal entries in Nambu space are zero and the GFs only depend on the relative time t rel . Then, the GFs can be trivially Fourier transformed back with respect to t rel to the energy space. However, because of the occurrence of t av , this is not possible in the SC case. Hence, we introduce the Wigner representation of the GFs following Ref. [106]. The nth moment of the Keldysh GF is defined by dt av e ıEt rel +ınΩtavǧ (t rel , t av ), where Ω = 2eV is the fundamental energy (periodicity of the average time argument) and τ = 2π/Ω is the corresponding period. This way, we restrict the Floquet energy to the first Floquet Brioullin zone, namely E ∈ [−Ω/2, Ω/2] = [−eV, eV ]. From the Wigner representation of the GF, one can compute the components of the so-called Floquet matrix G which are given by (G) mn (E) ≡ǧ mn (E) ≡ǧ m−n (E + (m + n)eV ) , (A9) which explicitly written out for the case of the voltagebiased SC GF in Keldysh space reads (we omit the spin structure as the GF is identity in spin space) −e ıχ g −+ 2m−1 δ m,n e ıχ g +− 2m+1 δ m,n f +− 2m+1 δ m+1,n g ++ 2m+1 δ m,n e ıχ f ++ 2m+1 δ m+1,n f +− 2m−1 δ m,n+1 −e −ıχ g +− 2m−1 δ m,n e −ıχ f ++ 2m−1 δ m,n+1 −g ++ 2m−1 δ m,n     . (A10) where we have used the notation g n ≡ g(E + neV ). Hence, the GF is an infinitely large matrix in Floquet space with each entry being an 8×8 Keldysh-spin-Nambu GF. In the case of a constant voltage, there are off-diagonal contributions in the Floquet space, characterized by factors of δ m+1,n and δ m,n+1 . In the case when the SC is not voltage-biased, there are no off-diagonal entries in Floquet space and we obtaiň so the Floquet GF is diagonal in Floquet space. In the case of a normal metal, g r/a N = ±1 and f r/a = f r/a = 0 and the voltage-biased normal metal Floquet GF follows easily by taking these limits in Eq. (A10). In particular, the voltage-biased normal metal Floquet GF is diagonal in Floquet space. Let us remark that in the case of a normal metal and one SC, the voltage can be gauged onto the normal metal. Hence, both of the Floquet GFs are fully diagonal in Floquet space and the problem can be treated by just integrating over the whole energy range.
The other ingredient for the action in Eq. (9) is the normal state scattering matrix. The structure of the scattering matrix is explained in detail in Sec. IV for the single-impurity and in Sec. VIII for the double-impurity case. The scattering matrix is a function of the energy E. Thus, upon a Fourier transformation, the scattering matrix only depends on relative times =s(t rel ) where (˜) signifies that the quantity is expressed in lead-Keldysh-spin-Nambu space. Thus, its Wigner representation is trivial and its Floquet matrixS is diagonal in Floquet space. In particular, the element (n, m) of the Floquet matrix is given by (S) mn ≡s mn (E) =s(E + 2meV )δ m,n , with the energy-dependent scattering matrix in lead-Keldysh-spin-Nambu space. detail. We start by the Keldysh action Tr lnQ(0), (B1) whereŜ is the normal-state scattering matrix and G(χ) = diag(G L (χ), G R ) is a block-diagonal matrix containing the lead GFs G L/R . In this formalism, the reservoir GFs are expressed as infinitely large matrices in time-space. Namely, the element (t, t ′ ) of the GF G L/R (χ) is the Keldysh GFǧ L/R (χ, t, t ′ ) from Eq. (A7). In addition, the scattering matrix is expressed in the same way as a matrixŜ with its element (t, t ′ ) being the scattering matrixs(t, t ′ ). Note that by expressing the GFs and scattering matrix as infinitely large matrices in time space, the trace and logarithm in Eq. (B1) imply convolutions over intermediate arguments. In particular, it holds where G L G R describes a matrix multiplication in time-Keldysh-spin-Nambu space between two GFs G L and G R in time-Keldysh-spin-Nambu space. The subindex ( t,t ′ ) refers to extracting the element (t, t ′ ). As a reminder, the Floquet matrix representation of the GFs G L is G L (E) with its element (n, m) being the Keldysh matrixǧ nm L (E) from Eq. (A10) and equivalently for the GF G R and its Floquet matrix G R . An important mathematical relation of the Floquet matrix representation in Eq. (A9) is that upon mapping the GFs in time space to the Floquet space, the algebraic structure of a generalized convolution is preserved [106] (ǧ L ⊗ǧ R )(t, t ′ ) = (G L G R ) t,t ′ ⇔ ∞ l=−∞ǧ ml L (E)ǧ ln R (E) = (G L (E)G R (E)) mn , whereǧ ml L andǧ ln R are the components of the Floquet matrices of the Keldysh functionsǧ L andǧ R respectively. It is seen that a multiplication of two Keldysh functions in time-Keldysh-spin-Nambu-space can just be translated into a multiplication of Floquet matrices in Floquet-Keldysh-spin-Nambu-space. Thus we can rewrite the action in Eq. (B1) using Floquet matrices Practically, the Floquet matrices cannot be expressed as infinitely large matrices. Thus, when computing the determinant in Eq. (B5) numerically, one introduces a cutoff in the Floquet index large enough for the cumulants to converge.
In the case of single-impurity junctions, the matrixQ is block-diagonal in spin space, thus detQ = detQ ⇑ detQ ⇓ and the charge-and spin-resolved tunneling probabilities are computed using the inverse Fourier transformation for σ =⇑, ⇓. In the two-impurity case, the matrixQ is not block-diagonal in spin space and the charge-resolved tunneling probabilities can be obtained as follows P n (E, V ) = 1 2π