Yu–Shiba–Rusinov screening of spins in double quantum dots

A magnetic impurity coupled to a superconductor gives rise to a Yu–Shiba–Rusinov (YSR) state inside the superconducting energy gap. With increasing exchange coupling the excitation energy of this state eventually crosses zero and the system switches to a YSR ground state with bound quasiparticles screening the impurity spin by ħ/2. Here we explore indium arsenide (InAs) nanowire double quantum dots tunnel coupled to a superconductor and demonstrate YSR screening of spin-1/2 and spin-1 states. Gating the double dot through nine different charge states, we show that the honeycomb pattern of zero-bias conductance peaks, archetypal of double dots coupled to normal leads, is replaced by lines of zero-energy YSR states. These enclose regions of YSR-screened dot spins displaying distinctive spectral features, and their characteristic shape and topology change markedly with tunnel coupling strengths. We find excellent agreement with a simple zero-bandwidth approximation, and with numerical renormalization group calculations for the two-orbital Anderson model.


SUPPLEMENTARY NOTE 1: THEORETICAL MODEL
The double-dot system is described by a two-orbital Anderson model, where L refers to the leads, D to the dots, and T to the tunnel coupling between the different parts. Since our system has the two dots in series, we shall refer to them by either QD N or QD S , or sometimes even just N or S, i.e. with the label of the closest N (normal) or S (superconducting) lead (cf. Supplementary Fig. 1). The three parts of the Hamiltonian are where ∆ N = 0 and ∆ S = ∆ denotes the (BCS) superconducting gap in lead S, and ξ αk denotes the energy relative to the chemical potential of an electron with momentum k in lead α. Operators c finite dot spin [1]. Since also t S t N , we expect the Kondo temperature due to the normal lead to be very small, and throughout our analysis, we shall merely treat the normal lead as a weakly coupled tunnel probe.

SUPPLEMENTARY NOTE 2: ADDITIONAL DATA ON N-DQD-S SYSTEMS
This section presents additional data on normal metal-double quantum dot-superconductor (N-DQD-S) systems for two devices (A and B). We first elaborate on the data set in Fig. 3 in the partly screened regime, which was measured in device A. In particular, we clarify how the parameters used in the fit of the zero-band-width (ZBW) model are estimated by analyzing stability diagrams and bias spectroscopy plots. Furthermore, we comment on the behavior of other double-dot shells observed in the stability diagram. Finally, we present stability diagrams and sub-gap spectroscopy plots for device B in the honeycomb regime, where the tunnel couplings between the dots and to the superconductor are small.   are observed. The two-orbital double-dot shell A in the top left corner shows a honeycomb diagram representing the familiar honeycomb regime (see ZBW model calculations in Fig.   2f in the article). Interestingly, the stability diagram also shows other two-orbital shells, in particular D-F, which all resemble the behaviour of the partly screened regime. The double-dot shell D was analyzed in Fig. 3 in the article. The double-dot shells B-C display less regular behaviour, which we interpret as the behaviour expected close to the transition between the honeycomb and partly screened regimes (see for instance the ZBW and NRG generated stability diagrams versus tuning of t S in Supplementary Figs. 7c,f and 7g,j).  Fig. 3b), the features are smeared, which is attributed to the strong coupling to the superconducting electrode. We note that in studies of the two-impurity Kondo effect with double dots strongly coupled to leads, the honeycomb pattern appears similarly smeared [2]. When the electrode turns superconducting, the features sharpen and a conductance pattern resembling the partly screened regime is observed ( Supplementary Fig. 3a). Similar characteristic mirrored arc patterns have also been observed in SQUID nanotube samples in the S-DQD parallel configuration (where these patterns are also expected for appropriate parameters) [3].
The parameters used in the ZBW modelling are shown in Supplementary Table 1. We now explain how they were obtained. For easy comparison, columns 1 and 3 (c-e & i-k) in Supplementary Fig. 3 show sub-gap spectroscopy plots, also displayed in the main text yields an upper limit of around 1 meV. Since the important parameter for the model is t 2 S /U S (not U S ), an error in the estimate of the charging energy may be compensated by adjusting the tunnel coupling t S , and therefore we do not have a fully independent estimate for this parameter. We use t S and t d as fitting parameters to fit the 6 sub-gap spectroscopy plots in     by a single quasiparticle, corresponding to a perturbative solution based on Yosida's variational wavefunction [5,6]. Here we extend this ZBW approach to deal with the two-orbital Anderson model (1). As we shall see below, the resulting sub-gap spectrum compares very well with the exact results obtained from numerical renormalization group (NRG) calculations, even in the strong coupling regime with a fully YSR screened ground state, and with a number of experimental details observed throughout the 9 different charge states of the stability diagram.
The ZBW approximation corresponds to replacing the Hamiltonian for the S-lead by where c † σ may be thought of as the creation of an electron in the local frontier orbital of the superconducting lead (cf. Appendix C in Ref. 7), which is most strongly coupled to the dot. For now, we focus on the sub-gap spectrum, and the normal lead will be neglected altogether until we return to discuss the actual conductance.

The dot part of the Hamiltonian can be rewritten as
with n i = n i↑ + n i↓ andñ i = 1/2 − i /U i . The parametersñ i are hence rescaled gate voltages defined so that in the limit of weakly coupled dots one has n i ≈ñ i . One can also introduce the related parametrisation δ i = 1 + 2 i /U i which quantifies the departure from the particlehole symmetric point at δ i = 0. The two sets are connected through δ i /2 = 1 −ñ i . The tunneling Hamiltonian reduces to This simple ZBW Hamiltonian is readily diagonalized numerically. Focusing on the sub-gap spectrum, we may now compare to the experimental data. In Supplementary Fig. 6 we show our best ZBW fit (b, e) to the experiment (a, d). To further verify that the ZBW model indeed yields an acceptable qualitative picture, we also performed NRG calculations and found excellent agreement, see Supplementary Fig. 6(c, f). In addition, we systematically compared ZBW and NRG in different regimes. In Supplementary Fig. 7 the stability diagrams calculated by both methods for different couplings to the superconducting lead show the evolution from the honeycomb regime via the partly screened to the screened regime.
The contour lines indicate an excitation energy of zero (red) and 0.01 meV (black), respectively, while the colours indicate the ground state (blue for singlet and brown for doublet).
It is a priori unclear how to relate the hopping t S in the ZBW to the hybridisation strength Γ in the NRG, since in ZBW the continuum of states is represented by a single level. In our NRG calculations we used a flat density of states with a Γ that does not depend on the energy. Empirically we found that the best agreement between both approaches as regards the positions of the YSR states, which is our main focus here, is obtained for Γ = c t 2 S with c ≈ 2π. not observed in the experiments presented in Fig. 4 in the main text, although the triplet state is theoretically predicted to exist and the transition is allowed by the selection rules.
Supplementary Fig. 8a shows the stability diagram in the screened regime, where the ground state changes from singlet (blue) to doublet (brown) to singlet as the electrochemical potential on QD N is swept, almost independently of the filling of QD S . To examine the sub-gap excitation spectrum, we make a cut along the n N = 1 line shown in Supplementary Fig. 8a by an arrow. The resulting spectral density plot along this line is shown in Supplementary   Fig. 8b. The spectral density plot reveals two excitations corresponding to the singlet and triplet states at low and high excitation energy, respectively. It is clearly seen that the triplet excitation has significantly lower spectral weight than the singlet excitation, which may give an indication why the triplet state was not experimentally observed. For finite coupling to N, the sub-gap states broaden and this effect is further pronounced.
with perturbation theory using the full quantum spin [1], so with S = 1/2, we end up with the following tunneling Hamiltonian Neglecting J NS and W NS , this Hamiltonian leads to sub-gap YSR bound states with energies with α = 3πν FS J SS /2 and β = πν FS W SS , where ν FS denotes the normal-state density of states in the superconductor. Within the polarized-spin approximation, the S-electron self-energy is exact and the Dyson equation is readily solved to find the exact Green's function which has poles at ω = ± 0 . Following Refs. 8 and 12, we now restrict the Green's function to |ω| < ∆ and expand around ω ≈ ± 0 to arrive at the retarded 4 × 4 Nambu Green's function expressed in terms of 2 × 2 blocks written in the Nambu-spinor basis and with u and v defined as Notice that, unlike 0 , u and v are sensitive to the sign of β, i.e. to the sign of δ S , which is opposite for gate-voltages on opposite sides of the particle-hole symmetric point, S = −U S /2.
Switching sign of δ S effectively interchanges u 2 and v 2 . Turning back on the coupling to the normal lead, this effective resonant-level Green's function is readily dressed by cotunneling amplitudes J NS and W NS , and analytically continued to arrive at where the total broadening, Γ = Γ r + Γ e + Γ h , is given in terms of Real parts of the tunneling induced self-energy are neglected here, assuming that 0 is much smaller than the normal lead bandwidth [13]. The additional term, Γ r , appearing in the broadening is put in by hand to capture quasiparticle relaxation in the superconductor, as discussed in Refs. 8 and 12. Quasiparticle relaxation could be due to electromagnetic noise or interactions among quasiparticles, say, but here we shall simply model Γ r as the tunneling rate to an additional fermionic bath with an unspecified bias-independent chemical potential µ r . Since the problem has been reduced to an effective resonant-level problem, the non-equilibrium Green's functions and the resulting current are readily calculated, and one finds that giving the zero-temperature differential conductance which no longer contains the unspecified µ r . When Γ r vanishes this expression becomes symmetric in bias-voltage, V , since transport is then carried exclusively by Andreev reflections proportional to 2Γ e Γ h in which incoming electrons are reflected as holes, thereby using both the electron and hole components of the YSR state and transporting two charges (cf. Supplementary Fig. 9b). Notice that aligning the normal-lead Fermi level instead with the negative energy YSR state (labeled v 2 in the figure), it is a Cooper pair at zero energy which splits into the two YSR states and leaves into the empty states of the normal lead. For finite Γ r , on the other hand, transport may also take place via the bias-asymmetrical relaxation process proportional to Γ e/h Γ r in which an electron tunnels from the normal lead to a YSR state, which then relaxes into the quasi-particle continuum (cf. Supplementary Fig. 9a).

Bias asymmetry and quasiparticle relaxation
In this section we fit experimental differential conductance with our classical spin model by using the parameters obtained with ZBW. To estimate the parameters we first analyze a b Supplementary Figure 9. Relaxation mechanism. Sketch of the two conducting processes in Supplementary Eq. 22. In a we see the relaxation process leading to bias asymmetrical conductance as one can probe the two YSR peaks independently and in b the Andreev process which is always bias symmetric as both peaks are always probed together.
the particle-hole symmetric point at which Γ e = Γ h . From Supplementary Eq. 21 we obtain the following finite-temperature conductance for the positive-bias peak, This equation is a convolution of a Lorentzian with a thermal distribution. Experimentally a temperature of T el ∼ 80 mK ∼ 7 µV is obtained (we here use a conservative estimate for the electron temperature based on measurements of another sample in the Coulomb blockade regime), while the full width of the YSR state is measured to be approximately 40 µV.
For the relevant range of parameters, this integral is very well approximated by adding the temperature to the width of the zero-temperature Lorentzian, that is with conductance peak height This peak height is measured to be approximately 0.11 e 2 /h, which together with the width 2(Γ + T ) ∼ 40 µV implies that Γ e ∼ 0.55 µV. Writing out Γ e as we find the following estimate for the dimensionless normal-lead exchange-cotunnel ampli-  which indicates that we are well away from the Kondo regime, since the Kondo temperature [7] are exponentially suppressed by a factor of ∼ exp(−1/0.005) ∼ 10 −87 . Since Γ r = Γ − 2Γ e , we may finally estimate the quasiparticle relaxation rate to be Γ r ∼ 12 µeV.
Bias asymmetry at the particle-hole symmetric point.
With this model we are able to explain the observed bias asymmetry by inclusion of a relaxational transport channel, but an additional asymmetry appears in the experimental data. For both (n S , n N ) = (1, 0) and (1, 2), the positive-bias peak conductance is different from the negative-bias peak conductance at the particle-hole symmetrical points, δ S = 0 (cf. Supplementary Fig. 10, panel e). This asymmetry is almost exactly opposite for charge states (1, 0) and (1,2), and therefore the explanation is likely to depend on the inter-dot Coulomb interaction, U d , by which the charge of the N-dot would effectively gate the S-dot.
Carrying out the Schrieffer-Wolff transformation to arrive at the effective Kondo model (13), we have indeed excluded all corrections of order U d /U N/S , and redoing the Schrieffer-Wolff transformation without this omission we find a number of corrections to the cotunnelling amplitudes (12). The second order amplitudes, J SS and W SS , are still given by Supplementary Eq. 12 when the N-dot is empty, but with a doubly occupied N-dot, the energy denominators are shifted, and the amplitudes may be obtained from Supplementary Eq. 12 by simply replacing δ S by δ (2) S = δ S − 4U d /U S . This corresponds to a shift in the apparent deviation from the particle-hole symmetric point. The third order inter-lead amplitudes, J NS and W NS , undergo similar shifts, and expanding to leading order in U d /U N , we obtain J (0) W where the superscript refers to the N-dot occupancy. These correction terms are indeed down by a factor of U d /U N , but whereas the potential scattering terms are W NS = 0 for δ S = 0, their new correction terms take a finite value. Therefore including U d /U N to lowest order in our effective coupling between metal and superconductor renders conductance asymmetric at the particle-hole symmetric point, since the positive peak depends on Γ e ∝ (3/2J NS + W NS ) 2 and the negative on Γ h ∝ (3/2J NS − W NS ) 2 . Finding peak conductance with this more complete Schrieffer-Wolff transformation yields the same asymmetry as seen in experiment, with positive-bias peaks enhanced and negative-bias peaks suppressed for charge state (1,0), and opposite for (1,2)(cf. Supplementary Fig. 10, panels g,h).