Two-Impurity Yu-Shiba-Rusinov States in Coupled Quantum Dots

TL;DR

This paper demonstrates a superconducting analog of the two-impurity Kondo model using coupled quantum dots, revealing gate-tunable ground state transitions and quantum phase boundaries in YSR states.

Contribution

It introduces a novel superconducting device that models the two-impurity Kondo effect and explores the quantum phase transitions between different YSR ground states.

Findings

Gate-tunable transition from inter-dot singlet to YSR singlets.

Observation of quantum phase boundaries to doublet YSR phases.

Cusped current-voltage characteristic modeled by an overdamped circuit.

Abstract

Using double quantum dots as the weak link of a Josephson junction, we realize the superconducting analog of the celebrated two-impurity Kondo model. The device shows a cusped current-voltage characteristic, which can be modelled by an overdamped circuit relating the observed cusp current to the Josephson critical current. The gate dependence of the cusp current and of the subgap spectra are used as complementary ground-state indicators to demonstrate gate-tuned changes of the ground state from an inter-dot singlet to independently screened Yu-Shiba-Rusinov (YSR) singlets. In contrast to the two-impurity Kondo effect in normal-state systems, the crossover between these two singlets is heralded by quantum phase boundaries to nearby doublet YSR phases in which only a single spin is screened.

Figures

Equations

$$\frac{du}{d\tau }$$

$$\frac{d\varphi }{d\tau }$$

$$⟨\sin \varphi ⟩=\mathrm{Im}\frac{I_{1-iuA}(A)}{I_{-iuA}(A)}$$

$$\frac{V}{R_{\mathrm{s}}{I}_{\mathrm{c}}}-\frac{R}{R_{\mathrm{s}}}⟨u⟩=\mathrm{Im}\frac{I_{1-i⟨u⟩A}(A)}{I_{-i⟨u⟩A}(A)}.$$

$$⟨I⟩=I_c\mathrm{Im}\frac{I_{1-i⟨u⟩A}(A)}{I_{-i⟨u⟩A}(A)},$$

$$⟨I⟩\approx \frac{1}{2}\frac{I_c{A}^2⟨u⟩}{1+A^2⟨u{⟩}^2}$$

$$I_0=⟨I⟩(V_0)=I_{c}A/4$$

$$G=d⟨I⟩/d⟨V_{\mathrm{sd}}⟩(0)=A^2/2R$$

$${\omega }_{RC}\ll {\omega }_{RL}\lesssim k_{B}T/\hslash \ll {\omega }_{R{C}_{\mathrm{J}}},$$

$$Z_{eff}(\omega )={\left(1/R+i{C}_{\mathrm{J}}\omega +1/(i\text{L}\omega )\right)}^{-1},$$

$$S_{\varphi }={\int }_{-\infty }^{\infty }\frac{d\omega }{\omega }\frac{\mathrm{Re}{\mathrm{Z}}_{\mathrm{eff}}\left(\omega \right)}{2R_Q}\frac{1}{1-e^{-\hslash \omega /k_{B}T}}$$

$$I_{\mathrm{c}}^{\ast }\approx I_{\mathrm{c}}\exp (-S_{\varphi }).$$

$$H=\sum _i{H}_{\mathrm{imp},i}+H_{\mathrm{LR}}+\sum _i{H}_{\mathrm{BCS},i}+\sum _i{H}_{\mathrm{hyb},i},$$

$$H_{\mathrm{imp},i}=\frac{U_i}{2}{\left({\hat{n}}_i-n_i\right)}^2,$$

$$H_{LR}=\sum _{\sigma }{t}_{\mathrm{d}}{d}_{L,\sigma }^{†}{d}_{R,\sigma }+\text{H.c.},$$

$$H_{\mathrm{BCS},i}=\sum _{k,\sigma }{ϵ}_k{c}_{i,k\sigma }^{†}{c}_{i,k\sigma }+\sum _k\Delta c_{i,k\uparrow }^{†}{c}_{i,k\downarrow }^{†}+\text{H.c.},$$

$$H_{\mathrm{hyb},i}=t_i\sum _{k\sigma }{c}_{i,k\sigma }^{†}{d}_{i,\sigma }+\text{H.c.}$$

$$\begin{array}{r}I=\frac{4\pi e}{\hslash }|t_{\mathrm{d}}{|}^2{\int }_{-\infty }^{\infty }[f(\omega -{\mu }_{\mathrm{L}})-f(\omega -{\mu }_{\mathrm{R}})]\\ A_{\mathrm{L}}(\omega -{\mu }_{\mathrm{L}})A_{\mathrm{R}}(\omega -{\mu }_{\mathrm{R}})d\omega ,\end{array}$$

$$I=\frac{4\pi e}{\hslash }|t_{\mathrm{d}}{|}^2{\int }_{-e{V}_{\mathrm{sd}}/2}^{+e{V}_{\mathrm{sd}}/2}{A}_{\mathrm{L}}(\omega -e{V}_{\mathrm{sd}}/2)A_{\mathrm{R}}(\omega +e{V}_{\mathrm{sd}}/2)d\omega .$$

Full text

Magnetism relies on the presence of magnetic moments and their mutual exchange interactions. At low temperatures, local moments in metals may be screened by the Kondo effect and magnetism can be disrupted. This competition was first proposed by Mott Mott (1974) as a mechanism for the vanishing of magnetism at low temperatures in the $f$-electron metal CeAl3, and later explored by Doniach Doniach (1977) within a simple one-dimensional Kondo-lattice model, from which he established a phase diagram delineating the magnetic phase as a function of the ratio between the Kondo temperature, $T_{K}$, and the inter-impurity exchange. The essence of this competition was subsequently reduced to the vastly simpler two-impurity Kondo model, which exhibits an unstable fixed point separating a ground state (GS) of two local Kondo singlets from an inter-impurity exchange singletJones and Varma (1987). This competition remains a central ingredient in the current understanding of many heavy-fermion materials and their quantum critical properties Georges et al. (1996); Hewson (1997); Coleman (2007); Löhneysen et al. (2007); Bulla et al. (2008); Gull et al. (2011).

In a superconductor (S), the gap around the Fermi surface terminates the Kondo screening process before its completion, but local magnetic moments may still be screened by forming local singlets with BCS quasiparticles. As demonstrated by Yu, Shiba, and Rusinov (YSR) Yu (1965); Shiba (1968); Rusinov (1969), a local exchange coupling between a superconductor and a magnetic impurity leads to a subgap bound state. In a full quantum description, the corresponding bound state crosses zero energy and the singlet subgap state becomes the new GS at a coupling strength corresponding to $T_{K}\approx 0.3\Delta$, where $\Delta$ denotes the superconducting gap Satori et al. (1992); Bauer et al. (2007). This quantum phase transition reduces the spin by $\hbar/2$, quenching a spin-1/2 altogether.

Half-filled Coulomb-blockaded quantum dots (QDs) are model magnetic impurities, and can be readily coupled to superconductors in semiconductor nanowires Prada et al. (2019). In Josephson junctions (JJs) hosting an impurity, the above GS transition results in a phase change from $\pi$ to 0 Van Dam et al. (2006); Cleuziou et al. (2006); Grove-Rasmussen et al. (2007); Jørgensen et al. (2007); Maurand et al. (2012); Delagrange et al. (2015, 2016); Estrada Saldaña et al. (2018, 2019); Meden (2019), and constitutes a complementary experimental signature to bias spectroscopy of YSR states crossing zero energy Deacon et al. (2010); Kim et al. (2013); Lee et al. (2014); Jellinggaard et al. (2016); Lee et al. (2017); Gramich et al. (2017); Grove-Rasmussen et al. (2018); Gramich et al. (2017); Su et al. (2017); Estrada Saldaña et al. (2020).

In contrast to the normal-state two-impurity Kondo effect, the two-impurity YSR ground-state phase diagram depends strongly on the two different local exchange couplings and includes not only the two different singlets, but also a doublet GS in which only a single spin is screened. In a previous Letter, we used a nanowire device to demonstrate the Josephson effect in a serial double quantum dot (DQD) in the low-coupling regime, in which the inter-dot tunnel coupling, $t_{d}$, dominated over the dot-lead tunnelling rates, $\Gamma_{L},\Gamma_{R}$ Estrada Saldaña et al. (2018). Here, we investigate in the same device the possible GSs at stronger individual couplings to the leads. As the gate voltages controlling these couplings are tuned, the boundaries of the honeycomb charge stability diagram between GSs of different parity are erased and new two-impurity YSR GSs are accessed. This boundary deletion affects the dispersion in gate voltage of the extracted Josephson critical current, $I_{\mathrm{c}}$, and of a closely related and directly measurable cusp current Naaman et al. (2001); Rodrigo et al. (2004); Randeria et al. (2016); Cho et al. (2019); Peters et al. (2020); Liu et al. (2020), as well as the dispersion of the spectral YSR resonances, which render mutually consistent GS parity information.

The article is organized in sections as follows. In Section I, we show an overview of the DQD and its parameters. In Section II, we introduce the method used to extract $I_{\mathrm{c}}$. In Section III, we delineate the theoretical expectations for GS, $I_{\mathrm{c}}$ and subgap spectra in the S-DQD-S system. In Section IV, we show our main experimental results, making use of the concepts introduced in Section III. Finally, in Section V we present our conclusions and comment on the possible applications of our findings.

I Double quantum dot characterization

The measured device (Fig. 1a) consists of an Al-covered InAs nanowire Krogstrup et al. (2015) deposited on top of narrow gates, insulated by 20 nm of hafnium oxide from the wire. Al is etched away to form a Josephson junction with a bare segment of InAs nanowire as the weak link. The device has a Si/SiO2 substrate backgate which we operate at $V_{\mathrm{bg}}\sim 10$ V to observe a measurable supercurrent. Gates 1, 3 and 5 are set to negative voltages to define the DQD, while gates 2 (4) are used as the plunger $V_{\mathrm{gL}}$ ($V_{\mathrm{gR}}$) to load electrons into the left (right) QD. Slight changes in the voltage of gate 1 (5) modify $\Gamma_{\mathrm{L}}$ ($\Gamma_{\mathrm{R}}$) of the left (right) QD to the left (right) lead. Small changes in $V_{\mathrm{bg}}$ also affect these tunnelling rates. All measurements are done in an Oxford Triton® dilution refrigerator at $T_{\rm fridge}\approx 20$ mK, using standard lock-in techniques with a lock-in AC excitation of 2 $\mu$V at 116.81 Hz to obtain the differential conductance, $dI/dV_{\mathrm{sd}}$, superposed to a source-drain DC voltage, $V_{\mathrm{sd}}$, while simultaneously recording the current, $I$, with a digital multimeter. The data was corrected for an offset of 3.5 pA from the current amplifier.

Figure 1b shows a zero-bias differential conductance ($G$) colormap which represents a portion of the honeycomb stability diagram of the DQD in the superconducting state. Conductance lines appear at places in which the parity of the system changes, constituting a way to accurately map GS boundaries. We label charge sectors of the DQD shell selected for gate tuning (shell X) by $(N_{\mathrm{L}},N_{\mathrm{R}})$, where $N_{\mathrm{L}}$ and $N_{\mathrm{R}}$ are the (integer) charges in the highest level of the left and right QDs, respectively. Parity lines separating doublet and singlet regions according to the total number of electrons in both QDs alternate in spacing, consistent with shell filling Estrada Saldaña et al. (2018); Jørgensen et al. (2008).

From Coulomb-diamond spectroscopy of shell X, we obtain the charging energies of the left and right QDs, $U_{\mathrm{L}}\approx 1.9$ meV and $U_{\mathrm{R}}\approx 1.6$ meV. From this spectroscopy, we also find that the level spacing of both QDs is equivalent to their charging energies ($\Delta E_{\mathrm{L}}\approx 1.6$ meV and $\Delta E_{\mathrm{R}}\approx 1.8$ meV). $\Delta=0.27$ meV, the parent Al superconducting gap, is found from Coulomb-diamond spectroscopy in deep Coulomb blockade in an opaque regime. The fact that $U_{\mathrm{L}},U_{\mathrm{R}}>\Delta$ places the system firmly in the YSR regime Kiršanskas et al. (2015). From the (0,2)-(1,1) anticrossing in the charge stability diagram, we estimate the inter-dot coupling, $t_{\mathrm{d}}=0.03-0.05$ meV and the inter-dot charging energy, $U_{\mathrm{d}}=0.13-0.23$ meV. Measurements of DQD and superconductor parameters are provided in Appendix A (see Figs. 9, 10 and 11).

The DQD parameters $\Gamma_{\mathrm{L}}$ and $\Gamma_{\mathrm{R}}$ could not be independently measured due to the high critical field ($B_{\mathrm{c}}=2.1$ T) and high critical temperature ($T_{\mathrm{c}}=2.2$ K) of the superconducting device. We calculate in Section III for our approximate DQD parameters GS boundaries set at $\Gamma_{\mathrm{L}}\approx\Gamma_{\mathrm{R}}\approx 0.4<\Delta E_{\mathrm{L}},\Delta E_{\mathrm{R}}$, which provide $\Gamma_{\mathrm{L}},\Gamma_{\mathrm{R}}$ sufficiently small to satisfy the single-level approximation Van Dam et al. (2006). Multilevel QD JJs show phase changes at fixed dot charge, irrespective of charge occupation Van Dam et al. (2006); Delagrange et al. (2016), which is not observed in shells X, Y and Z.

Near charge depletion, QDs are typically strongly confined, but also show the smallest $\Gamma_{i}/U_{i}$ ratio, while a sizeable $\Gamma_{i}/U_{i}$ ratio for a given temperature is needed to observe a finite Josephson current De Franceschi et al. (2010). This sets a narrow ($\Delta E_{i}$, $\Gamma_{i}$, $U_{i}$) window which translates into a narrow gate range for single-level DQD JJ operation and $\Gamma_{i}$ tuning. Outside of this window either multilevel ($\Delta E_{i}\ll\Gamma_{i}$), resonant level ($\Gamma_{i}\gg U_{i}$), or too weakly-coupled regimes arise ($\Gamma_{i}\ll U_{i}$), which is why only a few shells near QD depletion could be studied here.

II Determination of critical and cusp currents

Next, we focus on the low-bias transport characteristic of the device in the superconducting state.

Initially intended for two-terminal bias spectroscopy of YSR excitations in a S-DQD-S system, for sufficiently strong $t_{\mathrm{d}}$ the device exhibits in $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ traces a bias-symmetric, 20 $\mu$V-wide zero-bias peak flanked by negative differential conductance (NDC) dips, exemplified in Fig. 2a. As shown in Fig. 2b, the corresponding $I-V_{\mathrm{sd}}$ traces have a finite slope between sharp cusps which occur at a current $I_{\mathrm{0}}$. The positive $V_{\mathrm{sd}}$ value at which the cusps occur, $V_{0}=9.1\pm 0.6$ $\mu$V, is gate-independent over 3 orders of magnitude of $G$ and over cusp currents, $I_{\mathrm{0}}$, ranging from 0.3 pA to 50 pA (see Fig. 12 in Appendix B). Additionally, $V_{\mathrm{0}}$ matches the $V_{\mathrm{sd}}$ value of NDC dips, and $I_{\mathrm{0}}$ shows a gate dependence similar to $G$ (for example, compare Figs. 14c,d in Appendix C), indicating common origins for both phenomena.

We relate this narrow, bias-symmetric, finite-sloped $I-V_{\mathrm{sd}}$ characteristic, which is cusped at $I_{\mathrm{0}}$ and locked at $V_{\mathrm{0}}$, to the supercurrent of the circuit-damped JJ formed by the S-leads and the DQD weak-link, schematized in Fig. 2c. This minimal model assumes the presence of a resistance $R$ at the level of the leads (indicated in green) which was not designed and which was not independently measured, and which provides Johnson-Nyquist white noise at a temperature $T$ which in principle can be larger than the temperature in the dilution refrigerator Pekola et al. (2010); Naaman et al. (2001). From its fitted value and position in the model, we speculate that $R$ comes from the contact resistance between the nanowire and the Ti/Au leads.

In Fig. 2c, the designed elements of the circuit at 20 mK are indicated in black. The JJ is assumed to have a current-phase relationship $I_{\mathrm{c}}\sin\phi$, where the critical current, $I_{\mathrm{c}}$, is unknown. This relationship is valid away from phase transitions Estrada Saldaña et al. (2018). Standard filtering and shielding elements include: 1) Two-stage low-pass RC filters with 30 kHz cut-off frequency, of total resistance $R_{\mathrm{s}}/2$ per line, where $R_{\mathrm{s}}=8.24$ k$\Omega$, and capacitance per stage $C_{\rm filter}=2.7$ nF. 2) 7-stage low-pass $\pi$ filters with 200 MHz cut-off frequency, positioned at 20 mK in series with the RC filters (not shown). 3) Lines to the sample are additionally attenuated at high frequency by Eccosorb® and encased by copper tape. 4) Radiation shielding is used at several stages in the dilution refrigerator, with the deepest encasing being at 20 mK at the sample level.

Non-designed on-chip elements of known origin in the experiment are indicated in blue. Using a parallel-plate capacitor model, and from the device geometry and materials, we evaluate the junction capacitance at $C_{\mathrm{J}}\approx 3$ aF, and the capacitance of the large-area bonding pads through the backgate at $C\approx 9$ pF. The latter estimation assumes that the resistance R is somewhere in the leads between the bonding pads and the device. These capacitances place the JJ in the overdamped regime, but knowledge of their values is not needed to extract $I_{\mathrm{c}}$.

In the overdamped and voltage-biased regime this circuit gives rise to a three-parameter ($I_{\mathrm{c}}$, $R$, $T$) voltage-current relation, as demonstrated by Ivanchenko and Zil’berman Ivanchenko and Zil’berman (1969). When $T>\hbar I_{\mathrm{c}}/2ek_{B}$, corresponding to the limit of large noise amplitude Naaman et al. (2001), $V_{\mathrm{0}}$ is $I_{\mathrm{c}}$-independent as in the experiment and its value $V_{\mathrm{0}}=2eRk_{B}T/\hbar=10$ $\mu$V is close to the measurement. Additionally, $I_{\mathrm{c}}$ can be rescaled to $I_{\mathrm{0}}$ and $G$, as shown in Fig. 13 in Appendix B, explaining their connection in the experiment.

We fit $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ and $I-V_{\mathrm{sd}}$ curves to obtain $I_{\mathrm{c}}$. Since the device is measured in a two-terminal configuration, $R_{\mathrm{s}}=8.24$ k$\Omega$ is subtracted from the raw bias voltage, $V_{\mathrm{sd}}=V-IR_{\mathrm{s}}$, and raw differential conductance data, $dI/dV_{\mathrm{sd}}$=$(dI/dV)/(1-R_{\mathrm{s}}dI/dV$), before inputting them into the formula. Despite the circuit assumptions and the crudity of the model employed, the fits to the data are good for $I_{\mathrm{0}}$ between 0.2 pA and 50 pA, and $I_{\mathrm{c}}$ ranging from 0.03 nA to 0.8 nA. Examples of the fit are shown as blue curves in Figs. 2a,b.

The parameters $R=3$ k$\Omega$ and $T=80$ mK are extracted through the fitting procedure by initially keeping the parameters $I_{\mathrm{c}}$, $R$ and $T$ free. As long as $RT=240$ k$\Omega\cdot$mK, $R$ and $T$ can adopt any value and still produce a good fit to the data. A rescaling by a dimensionless number $\gamma$ to $R/\gamma$ and $T\gamma$ globally will merely change $I_{\mathrm{c}}$ to $\sqrt{\gamma}I_{\mathrm{c}}$ without any qualitative change in its gate dependence. This is crucial for a robust interpretation of the latter, but it prevents absolute quantitative $I_{\mathrm{c}}$ assessments. Nevertheless, $T$ is chosen from an upper-bound estimate of the electron temperature from the thermal broadening of a Coulomb peak, which also produces a reasonable $I_{\mathrm{c}}$ value, given the maximum imposed by $\Delta$ and the DQD couplings. The fitted value of $R$, in turn, is compatible with previous measurements of contact resistance in highly transparent nanowire junctions van Weperen et al. (2013). Once the noise $RT$ was set, $I_{\mathrm{c}}$ was kept as the only free parameter in the fit to $V_{\mathrm{sd}}-I$ curves at other gate voltages. Extraction of an unambiguous $I_{\mathrm{c}}$ value, impeded here by lack of measurement of $R$ Naaman et al. (2001), should benefit from a fully designed circuit Steinbach et al. (2001). A derivation of the model and additional details are shown in Appendix B.

III Theoretical insights

We introduce below the phase diagram of the S-DQD-S system, which distills the main consequence of the YSR spin-screening mechanism: changes in tunnelling rates $\Gamma_{L}$ and $\Gamma_{R}$ drive GS parity transitions. We also provide calculations of the parity stability diagram, of $I_{c}$, and of the subgap $dI/dV_{\mathrm{sd}}$ YSR spectra, of importance to undertake the experimental results in the next section.

In our GS naming convention, the parity is specified by $D$ for doublet and $S$ for singlet. Indexes depend on QD occupation, as follows. If the QD levels contain each one electron, we use two indexes; the first (second) denotes if the left (right) spin is screened (1) or not (0) (Fig. 3a). If only one QD contains a single electron and the other is empty/full, we use a single index, L or R, to indicate if the left or right QD level is occupied (Fig. 4a). Finally, if both QD levels are empty/full, we use no index. For example, in Fig. 3a $D_{\mathrm{01}}$ denotes screening of the right-QD spin for the (1,1) charge state.

In Fig. 3b, we show the two-impurity phase diagram of the S-DQD-S system vs. $\Gamma_{\mathrm{L}}$, $\Gamma_{\mathrm{R}}$ with dot occupations fixed at $n_{\mathrm{L}}=n_{\mathrm{R}}=1$. This diagram was calculated by the Numerical Renormalization Group (NRG) technique Wilson (1975); Lee et al. (2014) for $t_{\mathrm{d}}=0.1$ meV, charging energies of the QDs $U_{\mathrm{L}}=U_{\mathrm{R}}=2$ meV, and $\Delta=0.25$ meV, parameters which are used in all NRG calculations shown in this work. For simplicity, $U_{\mathrm{d}}$ was kept at zero. Black lines denote GS parity boundaries between $S_{\mathrm{00}}$, $S_{\mathrm{11}}$ and $D_{\mathrm{10}}$, $D_{\mathrm{01}}$. These boundaries avoid each other due to $t_{\mathrm{d}}$. An NRG calculation of the one-impurity phase diagram is in turn shown in Fig. 4b. The one-impurity case corresponds to $n_{\mathrm{L}}=1$, $n_{\mathrm{R}}=0,2$ (red) and to $n_{\mathrm{L}}=0,2$, $n_{\mathrm{R}}=1$ (green). In contrast to two-impurity GS parity boundaries (black dashed lines overlaid from the calculation in Fig. 3b), one-impurity parity boundaries are straight lines.

GS parity boundaries can also be plotted as function of $n_{\mathrm{L}}$, $n_{\mathrm{R}}$ for fixed $\Gamma_{\mathrm{L}}$, $\Gamma_{\mathrm{R}}$, which encourages qualitative comparison to experimental charge stability diagrams such as the one in Fig. 1b. Three examples calculated by NRG are shown in Fig. 5. Similarly to the charge stability diagram of shell X in Fig. 1b, the GS diagram of Fig. 5a exhibits all possible parity boundaries; its $\Gamma_{\mathrm{L}}$, $\Gamma_{\mathrm{R}}$ coordinates situate it at the bottom-left corner of the two and one-impurity phase diagrams of Figs. 3b and 4b. In Fig. 5b, at stronger $\Gamma_{\mathrm{L}}$, the left spin is screened and the GS changes in the (1,0) and (1,2) charge sectors from $D_{\mathrm{L}}$ to $S_{\mathrm{L}}$, and in the (1,1) sector from $S_{\mathrm{00}}$ to $D_{\mathrm{10}}$, deleting the horizontal parity boundaries which previously existed in Fig. 5a between $S_{\mathrm{00}}$ and $D_{\mathrm{R}}$, and between $D_{\mathrm{L}}$ and $S$. In Fig. 5c, at stronger $\Gamma_{\mathrm{R}}$ than in Fig. 5b, the remaining doublet sector is screened away into a singlet, resulting in the disappearance of vertical parity boundaries.

Modifications of the GS diagram result in changes of the gate dependence of $I_{\mathrm{c}}$ and subgap YSR spectra which are also experimentally resolvable. An NRG calculation of $\absolutevalue{I_{\mathrm{c}}}$ is shown in Fig. 5d for an $n_{\mathrm{L}}$ trajectory which crosses two GS parity boundaries in the case of the blue trace, and no GS parity boundaries in the case of the red one. The blue trace shows asymmetric peaks located at the two $D_{\mathrm{L}}-S$ GS parity changes, whereas the red one evolves into a smooth and broad peak centered at $n_{\mathrm{L}}=1$. These are the well-known signatures of a $\pi$ and a [math] QD JJ, respectively. In Figs. 5e-h, we compare the gate dependence of these traces (Figs. 5e,f) to a calculation of subgap $dI/dV_{\mathrm{sd}}$ (Figs. 5g,h), which is detailed in Appendix E. In Figs. 5g,h, the calculation shows $dI/dV_{\mathrm{sd}}$ colormaps in which the innermost red lines appear at $eV_{\mathrm{sd}}=\pm(\Delta+E_{\mathrm{L}})$, corresponding to quasiparticle tunnelling from the gap edge of the right superconductor via the spinless right QD, to the YSR state at energy $E_{\mathrm{L}}$ in the spinful left QD, with tunnelling rate $|t_{\mathrm{d}}|^{2}$. Since BCS peaks probe sub-gap bound states, these peaks are followed by dips of NDC. This transport channel assumes a quasiparticle relaxation rate in the left superconductor, which is faster than the inter-dot tunneling rate Ruby et al. (2015a). In Fig. 5g, the YSR $dI/dV_{\mathrm{sd}}$ peaks form the well-known split-loop which kinks downwards at parity changes coinciding with peaks in $\absolutevalue{I_{\mathrm{c}}}$ in Fig. 5e. In the absence of parity changes as in Fig. 5h, the YSR peaks evolve smoothly with minimum splitting at $n_{\mathrm{L}}=1$, which coincides with a broad peak in $\absolutevalue{I_{\mathrm{c}}}$ in Fig. 5f.

Since the high magnetic field and/or temperature needed to drive the device into the normal state impede independent measurement of $\Gamma_{\mathrm{L}}$ and $\Gamma_{\mathrm{R}}$, these calculations serve only as a qualitative guide for the gate-tuned behavior of experimental charge stability diagrams and of the gate dependence of $I_{\mathrm{c}}$, $I_{\mathrm{0}}$ and YSR excitations in the section below.

IV Results

By following the deletion of parity transition lines in the charge stability diagram for DQD-shell X, while knowing the departure ground states of each charge configuration $(N_{\mathrm{L}},N_{\mathrm{R}})$ of the original diagram at weak-coupling, it is possible to track GS changes of the system at stronger coupling to the leads, to which we access by changing gate voltages other than $V_{\mathrm{gL}}$ and $V_{\mathrm{gR}}$. Before describing the data, we discuss three points not included in the NRG model above:

While the model is based on a single DQD shell, the DQD device has a staircase of shells. GS changes demonstrated here for shell X are not necessarily concurrent in shells Y and Z, as it is typically the case in nanowire-defined QD devices Jellinggaard et al. (2016); Estrada Saldaña et al. (2019). 2. 2.

Parity transition lines in the experimental charge stability diagram of shell X have an acute angle with respect to each other due to trivial $V_{\mathrm{gL}}$, $V_{\mathrm{gR}}$ cross-talk van der Wiel et al. (2002), instead of the nearly perpendicular angle which they exhibit in the NRG calculation. While knowledge of the exact angle is not relevant, it is important to keep track of the slope of the parity transition lines of both QDs during the experiment, to avoid confusing gate-controlled parity changes with the formation of a single QD or with the introduction of a third QD, both of which should have parity transition lines of different slope. We therefore keep the same total gate variations $\Delta V_{\mathrm{gL}}$ and $\Delta V_{\mathrm{gR}}$ and aspect ratio in the stability diagrams of Figs. 6, 7 and 8, and use red dashed lines to mark the slope of the left and right QD parity lines. These lines have the same slope in all plots, showing that the device stays as a DQD through the whole gating procedure. 3. 3.

As discussed below, electron-hole asymmetric GS parity changes in shell X are observed for YSR screening of the right-QD spin, which can be modelled by a dependence of $\Gamma_{\mathrm{R}}$ on $V_{\mathrm{gL}}$ (see Appendix F). This cross-coupling is not unreasonable in view of the slight horizontal rightwards shift of the bottom gates in Fig. 1a with respect to the nanowire channel Estrada Saldaña et al. (2019).

In Fig. 6a, we show a zero-bias conductance colormap which represents the charge stability diagram of DQD shell X at slightly different gate settings than in the overview in Fig. 1b. In these settings, the DQD has been brought to the verge of a GS transition due to YSR screening of the left-QD spin. We tune the leftmost gate of the device, $V_{\mathrm{g1}}$, to gradually merge (Fig. 6b) and then, in combination with slight changes in $V_{\mathrm{g3}}$, $V_{\mathrm{g5}}$ and $V_{\mathrm{bg}}$, to erase (Fig. 6c) the parity transition lines of the stability diagram corresponding to the loading of a spin-1/2 in the left QD, consistent with an increase in $\Gamma_{\mathrm{L}}$. The end result is that the singlet-doublet-singlet GS sequence $S-D_{\mathrm{L}}-S$ along the solid blue line in Fig. 6a changes into an all-singlet sequence $S-S_{\mathrm{L}}-S$ along the solid red line in Fig. 6c.

To support our interpretation, we show in Figs. 6d-h $I_{\mathrm{c}}$, $I_{0}$ and subgap $I-V_{\mathrm{sd}}$ and $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ data taken with $V_{\mathrm{gL}}$ swept along these solid lines of corresponding color at different instances of the GS transition; $V_{\mathrm{gR}}$ is also swept to compensate for cross-capacitance. $I_{\mathrm{c}}$ and $I_{\mathrm{0}}$ traces (in inset) in Fig. 6d behave similarly to NRG calculations of $\absolutevalue{I_{\mathrm{c}}}$ in Fig. 5d. Asymmetric peaks in the blue trace come together in the intermediate green trace and merge into a single broad resonance in the red trace as the coupling $\Gamma_{\mathrm{L}}$ is increased by tuning the $V_{\mathrm{g1}}$. The asymmetric peaks in $I_{\mathrm{c}}$ are consistent with a phase-shift of $0-\pi-0$ added to the current-phase relationship when the GS parity changes as even-odd-even. In turn, the culminating single-broad peak is compatible with an enhanced supercurrent from screening of the spin of the left QD, as seen in earlier S-QD-S devices Kim et al. (2013); Cleuziou et al. (2006); Maurand et al. (2012). The absence of two asymmetric peaks in this trace is related to the even parity of the GS for all the relevant charge states.

Figures 6e,f compare the evolution of low-bias $I-V_{\mathrm{sd}}$ traces in gate voltage taken along solid lines of corresponding color in Figs. 6a,c, to Figs. 6g,h, which show $dI/dV_{\mathrm{sd}}$ colormaps vs. $V_{\mathrm{gL}}$ taken along the same respective gate trajectories, but with $V_{\mathrm{sd}}$ swept along a larger bias window within $\pm 2\Delta$, so as to capture the behaviour of YSR subgap $dI/dV_{\mathrm{sd}}$ peaks. In Fig. 6g, the lowest-lying, gapped peaks exhibit a small split-loop structure with kinks at -3.05 V and -3.1 V, which align in gate voltage to peaks in $I_{\mathrm{0}}$ in the $I-V_{\mathrm{sd}}$ traces from Fig. 6e, as predicted by our NRG calculations above, and as previously observed in S-QD-S Josephson junctions Kim et al. (2013). The kinks vanish in Fig. 6h, indicating absence of parity changes. Instead, the peaks exhibit a smooth point of inflexion at $V_{\mathrm{gL}}$=-3.43 V which aligns with a broad resonance in the $V_{\mathrm{gL}}$ dependence of $I_{\mathrm{0}}$ in the $I-V_{\mathrm{sd}}$ traces from Fig. 6f.

As in our simple model in Figs. 5e,f, the lowest-lying peaks appear split in bias voltage by a gapped region Lee et al. (2014); Jellinggaard et al. (2016); Deacon et al. (2010); Pillet et al. (2013); Lee et al. (2017); Li et al. (2017), are followed by NDC Kim et al. (2013); Pillet et al. (2010); Grove-Rasmussen et al. (2009); Eichler et al. (2007); Lee et al. (2012); He et al. (2020), and ascribe to the expected shape in gate voltage Lee et al. (2014); Deacon et al. (2010); Gramich et al. (2015); Estrada Saldaña et al. (2020). However, closer inspection reveals additional complexity in the form of reduced gap values and peak replicas, which we will discuss in some detail. If the Bardeen model in Section III is used to interpret the bias position of the lowest-lying gapped peak as $\Delta^{*}_{\mathrm{R}}+E_{\mathrm{L}}$, the measurements in Figs. 9c,d in Appendix A show that the superconducting gap is effectively reduced with respect to that of the parent Al superconductor ($\Delta$) and that it is different for the left ($\Delta^{*}_{\mathrm{L}}=0.110$ meV) and right ($\Delta^{*}_{\mathrm{R}}=0.140$ meV) leads. By an effectively reduced gap, we mean that the gapped region in the weak tunnelling regime at $V_{\mathrm{sd}}<2\Delta$ observed at $V_{\mathrm{bg}}=-15$ V in Fig. 9c is replaced in Fig. 9d, at the backgate voltage of operation of the device in shell X, $V_{\mathrm{bg}}\approx 10$ V, by conductance resonances above $V_{\mathrm{sd}}>\Delta^{*}_{\mathrm{L}}+\Delta^{*}_{\mathrm{R}}$. In electrostatic simulations of gating in InAs nanowire/Al hybrids, reduction of the effective induced gap in the nanowire occurred at more positive gate voltages due to weaker semiconductor/superconductor hybridization, which could explain our measurement Antipov et al. (2018).

Additionally, the colormaps display significantly less conducting, bias-symmetric replicas at $V_{\mathrm{sd}}>\Delta^{*}_{\mathrm{R}}+E_{\mathrm{L}}$, as well as one or more faint replicas at $V_{\mathrm{sd}}<\Delta^{*}_{\mathrm{R}}+E_{\mathrm{L}}$, the lowest of which can cross zero bias at parity crossings (see Appendix C). These replicas are often followed by NDC. In addition to replicas, on occasion a few straight horizontal lines cross the colormaps, such as in Fig. 6h, which shows a pair of bias-symmetric lines at $V_{\mathrm{sd}}=\Delta^{*}_{\mathrm{R}}$. The colormaps also display additional lines with opposite curvature to that of the lowest-lying gapped peaks; a clear example of this is seen in Fig. 6h around $V_{\mathrm{gL}}=-3.43$ V as split-loops terminating at the main YSR peaks. These additional lines can also exhibit replicas. The intricate replica behavior has been previously related to multiple peaks inside the superconducting gap of the hybrid nanowire-superconductor leads Su et al. (2018); Deng et al. (2016), to multiple Andreev reflection Nilsson et al. (2011); Villas et al. (2020), quasiparticle relaxation, thermally excited transport Ruby et al. (2018); Huang et al. (2020); Kumar et al. (2014) and/or to inelastic Andreev tunnelling Gramich et al. (2015). Understanding the origin of these replicas, some of which may be visible in this work due to the unprecedented resolution of our data (the full width at half maximum of YSR peaks can be as low as 10 $\mu$V, as shown in Fig. 14b), is outside of the scope of this work. In what follows, we will focus only on the gate dependence of the curvature of the lowest-lying pair of gapped peaks, and on its relation to $I_{\mathrm{0}}$.

We now bring shell X back to the honeycomb regime, this time with gate settings which put it close to a GS transition due to YSR screening of the right-QD spin. Fig. 7a shows a colormap of zero-bias conductance which represents the charge stability diagram of shell X in this regime. The sequence of colormaps in Figs. 7a-f shows tuning of the stability diagram from the honeycomb pattern into a pattern lacking parity lines from the right QD. We interpret the change as stemming from an effective increase in $\Gamma_{\mathrm{R}}$ due to our tuning of the rightmost gate voltage, $V_{\mathrm{g5}}$. The change is more subtle as $V_{\mathrm{gL}}$ also affects $\Gamma_{\mathrm{R}}$, as mentioned above. Due to this, the upper section of the stability diagram transitions into a new GS faster than the lower one. We are able to model this qualitatively by introducing a linear dependence of $\Gamma_{\mathrm{R}}$ on $n_{\mathrm{L}}$ within a zero-bandwidth approximation of the superconducting continuum Grove-Rasmussen et al. (2018), as shown in Appendix F, and find that this is a faithful demonstration of a $S_{\mathrm{00}}$ to $D_{\mathrm{01}}$ GS transition in the (1,1) charge sector, and of a $D_{\mathrm{R}}$ to $S_{\mathrm{1}}$ GS transition in the (0,1) and (2,1) charge sectors. In Appendix F, we also show that GS transitions obtained by increasing $\Gamma_{\mathrm{L}}$, as those shown Fig. 6, are not significantly affected by introducing a dependence of $\Gamma_{\mathrm{R}}$ on $n_{\mathrm{L}}$.

Finally, departing from the pattern in Fig. 7f, in which the right-QD spin of DQD-shell X is screened, we demonstrate simultaneous screening of the left-QD spin, and therefore, full YSR screening of the spins of shell X. Figure 8a-c shows the step-by-step modification of this stability diagram, which contains only parity lines of the right QD in shell X in Fig. 8a, to a diagram lacking parity lines of this shell in Fig. 8c, which marks a $D_{\mathrm{01}}$ to $S_{\mathrm{11}}$ GS transition. In Fig. 8c, after the GS transition occurs, the charge stability diagram of shell X shows a broad peak of conductance as the remains of the larger availability of conduction channels for Cooper pair transport in the (1,1) charge sector. The GS transition is also reflected in the gate dependence of $I_{\mathrm{c}}$, $I_{\mathrm{0}}$ and YSR peaks of largest conductance in Figs. 8d-h, obtained along gate paths given by solid lines of corresponding color in Figs. 8a-c, which cross the (1,1) charge sector. As in Fig. 6d, Fig. 8d shows that $I_{\mathrm{c}}$ and $I_{\mathrm{0}}$ (in inset) gradually change from a split-peak structure in the blue trace to a broad peak in the red trace. As in Figs. 6e-h, parity changes manifested in Fig. 8g in two kinks in the split-loop YSR structure which align to $I_{\mathrm{0}}$ peaks in $I-V_{\mathrm{sd}}$ traces in Fig. 8e are replaced by an all-singlet parity sequence, which expresses itself in the absence of kinks in Fig. 8h, and in the smooth inflection point at $V_{\mathrm{gL}}=-1.9$ V which coincides in Fig. 8f with a broad $I_{\mathrm{0}}$ peak in $I-V_{\mathrm{sd}}$ traces.

From the succession of Figs. 7a-f and Figs. 8a-c, we can see that the GS of the (1,1) charge sector, which started in Fig. 7a as an inter-dot singlet, $S_{\mathrm{00}}$, has been gate-tuned in Fig. 8c into independently screened YSR singlets, $S_{\mathrm{11}}$. This GS transition occurred through an intermediate YSR doublet phase, $D_{\mathrm{{01}}}$, as shown in Figs. 7f and 8a.

V Conclusion

In summary, we have demonstrated one- and two-impurity YSR physics in a shell of a DQD hybrid nanowire. We obtain the GS from the gate dependence of $I_{\mathrm{c}}$, of $I_{\mathrm{0}}$ and of the gate dispersion of YSR subgap conductance peaks, supported by the step-by-step tuning of the stability diagram to different endpoints of the phase diagrams. We find a reasonable qualitative agreement between experiment and theory. However, technical difficulties which prevent the measurement of $\Gamma_{L}$ and $\Gamma_{R}$, together with cross-couplings in the device, have as a result that the S-DQD-S model cannot be quantitatively established in relation to the experiment without free fitting parameters. Quantitative comparison between experiment and theory is also complicated by lack of independent measurement of $R$ within our model for the characteristic of the JJ circuit, which leads to extraction of $I_{\mathrm{c}}$ only up to a factor $\sqrt{\gamma}I_{\mathrm{c}}$, and by the complexity of the reduced hybrid nanowire-superconductor gap, which is reflected in multiple replicas of the YSR state filling the gap between the YSR state and the parent Al superconducting gap edge.

As a spectroscopic probe of parity changes, the narrow zero-bias conductance peak of the effectively voltage-biased Josephson junction Naaman et al. (2001); Rodrigo et al. (2004); Randeria et al. (2016); Cho et al. (2019); Peters et al. (2020); Liu et al. (2020) maintains the sharpness of the relevant features in all charge, or more accurately, parity stability diagrams independently of the tunnelling rates. This is in stark contrast to the case of magnetic impurities (spinful QDs) coupled to normal metals, which broaden the conductance features at strong hybridization Jeong et al. (2001); Bork et al. (2011); Chorley et al. (2012); Spinelli et al. (2015); Chang and Chen (2009); Georges and Meir (1999). The zero-bias conductance peak also provides very direct access to $I_{\mathrm{c}}$, just as it provides access to the local superfluid density in scanned Josephson-tunneling microscopy, which has been used recently to detect Cooper-pair density waves on surfaces of Bi2Sr2CaCu2O8+x and NbSe2 Hamidian et al. (2016); Liu et al. (2020). Unlike scanning tunnelling spectroscopies of dimers of magnetic adatoms on superconducting surfaces Ji et al. (2008); Ruby et al. (2018); Kamlapure et al. (2018), our DQD realization comprises two spin-1/2 states which are completely screened by individual superconducting channels. The Kondo-YSR analogy breaks down towards zero temperature, as confirmed by the existence of doublet domains in the phase diagram at $k_{B}T\ll\Delta$.

By virtue of $t_{\mathrm{d}}$, non-magnetic, gap-protected superpositions of the two singlet states found, the exchange singlet $S_{\mathrm{00}}$ and the independently screened YSR singlet $S_{\mathrm{11}}$, could be prepared in future works for parameters close to the anti-crossing in the phase diagram of Fig. 3 Wineland (2013), with the purpose of using them as qubits DiVincenzo (2000). In addition, the demonstrated gate control of a two-site quantum dot chain in superconducting proximity is a crucial step towards the implementation in our hybrid wires of the YSR analog of Doniach’s Kondo necklace Doniach (1977) and of the Kitaev chain Kitaev (2001); Sau and Sarma (2012); Fulga et al. (2013), complementing ongoing research of emergent manifestations of topology Mourik et al. (2012); Nadj-Perge et al. (2014); Ruby et al. (2015b); Pawlak et al. (2016); Deng et al. (2016); Grivnin et al. (2019); Manna et al. (2020).

Data availability. All data needed to evaluate the conclusions in the paper are present in the paper. Raw data used to produce the experimental figures in the paper can be found at the repository ERDA of the University of Copenhagen at https://doi.org/10.17894/ucph.c83b3677-34ff-4d10-8f5a-36785d1ed59e.

Acknowledgements.

We thank A. Jellinggaard, M.C. Hels, and J. Ovesen for experimental assistance. The project received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 832645. We additionally acknowledge financial support from the Carlsberg Foundation, the Independent Research Fund Denmark, QuantERA ’SuperTop’ (NN 127900), the Danish National Research Foundation, Villum Foundation project No. 25310, and the Sino-Danish Center. P. K. acknowledges support from Microsoft and the ERC starting Grant No. 716655 under the Horizon 2020 program. R. Ž. acknowledges support from the Slovenian Research Agency (ARRS) under Grants No. P1-0044 and J1-7259.

Appendix A: Extraction of DQD parameters

To measure $\Delta$, we deplete the device by setting $V_{\mathrm{bg}}=-15$ V (which is $\approx 25$ V more negative than $V_{\mathrm{{bg}}}$ used in the measurements on shell X in the main text) and perform Coulomb-diamond spectroscopy as shown in Fig. 9a. The repetitive Coulomb-blockade pattern is indicative of a multilevel QD regime. The Coulomb diamonds appear split by a forbidden bias window. Zooming in on this window, as shown in Fig. 9b, a pair of bias-symmetric lines which cross all diamonds is identified at $eV_{\mathrm{sd}}=\pm 2\Delta$. This feature stems from the BCS coherence peak from the left lead at energy $-\Delta$ probing the coherence peak in the right lead at energy $+\Delta$ via co-tunnelling through the QDs. At parity crossings, which correspond to the apexes of split Coulomb diamonds in this opaque regime, a few pairs of bias-symmetric subgap states come down, the lowest of which reaches exactly $\Delta$. This occurs when the coherence peak at $\Delta$ in the right lead probes a YSR excitation in the left lead which crosses zero energy at parity crossings. At $e|V_{\mathrm{sd}}|<\Delta$, the conductance is blocked at all $V_{\mathrm{g1}}$ values. In Fig. 9c, we plot a linecut at the center of one Coulomb diamond to display the position of the quasiparticle tunnelling peaks at $2\Delta$ and their extremely narrow width Huang et al. (2020). We measure $2\Delta=0.54$ meV, which provides a value of $\Delta$ identical to the one obtained from S-QD-N devices made from the same batch of nanowires Estrada Saldaña et al. (2020). The position of horizontal lines outside $2\Delta$ depends on the Coulomb diamond chosen, as can be seen in Fig. 9b; we ascribe these lines to the probing of QD states in co-tunnelling by the sharp superconducting singularities. The trace in Fig. 9c is compared to two traces in Fig. 9d displaying YSR excitations in shell X which are extracted from Fig. 6g to show graphically that $\Delta^{*}_{\mathrm{R}}<\Delta$ and that $\Delta^{*}_{\mathrm{L}}+\Delta^{*}_{\mathrm{R}}<2\Delta$, where $\Delta^{*}_{\mathrm{L}}$ and $\Delta^{*}_{\mathrm{R}}$ are the effectively reduced gaps in the left and right leads at the backgate voltage at which the device is operated, as in Section IV. If we assign a bias position $E_{\mathrm{L}}+\Delta^{*}_{\mathrm{R}}$ to the innermost gapped YSR peaks, we measure $\Delta^{*}_{\mathrm{R}}=0.140$ meV from the black trace, which is taken at a parity crossing where $E_{\mathrm{L}}=0$. In turn, we find $\Delta^{*}_{\mathrm{L}}+\Delta^{*}_{\mathrm{R}}=0.250$ meV from the green trace, which is taken in the (0,2) charge state, where $E_{\mathrm{L}}$ approaches asymptotically to $\Delta^{*}_{\mathrm{L}}$. From this, we deduce $\Delta^{*}_{\mathrm{L}}=0.110$ meV.

To measure $U_{\mathrm{L}}$, $U_{\mathrm{R}}$, we perform Coulomb-diamond spectroscopy on each QD in shell X as shown in Fig. 10. In spite of the larger couplings of the DQD in this regime, we observe faded Coulomb diamonds from which we can trace the actual diamonds as denoted by dashed lines. $U_{\mathrm{L}}$, $U_{\mathrm{R}}$ correspond to the bias difference between the apex of the central Coulomb diamond and the edge of the superconducting gap at $2\Delta$. An indication of the order of magnitude of the tunnelling rates of the DQD can be obtained from the full-width-at-half-maximum of Coulomb lines outside the gap between the pairs of arrows in Figs. 10b,c, which is of 0.6 meV in both cases.

To estimate $t_{\mathrm{d}}$ and $U_{\mathrm{d}}$, we use the curvature of the parity lines at the (1,1), (0,2) charge transition in the charge stability diagram of shell X. We first prolong the parity lines of the left and right QDs, as indicated by the yellow dashed lines in Fig. 11. We join their two intersections by a green line which separates the (1,1) and (0,2) charge sectors. Afterwards, by projecting the green line into the two white lines which indicate the charging energies of the QDs, we convert the gate scale of the green line into energy. The red lines along the green line indicate the distance from the curved parity lines in the colormap to each of the two intersections of the yellow dashed lines and are each equal to $2t_{\mathrm{d}}$; the green line summed to the two red lines is $\sqrt{2}U_{\mathrm{d}}+4t_{\mathrm{d}}$. This procedure assumes that the parity stability diagram in the superconducting state is similar to that of the normal state. In this way, we measure $t_{\mathrm{d}}\approx 0.03-0.05$ meV and $U_{\mathrm{d}}\approx 0.13-0.23$ meV.

Appendix B: Model of the circuit

In this appendix we elaborate on the extended RCSJ model necessary to explain the $I-V_{\mathrm{sd}}$ and $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ curves observed in the experiment. The extension, compared to standard RCSJ, consists of an additional series resistance $R$ and shunt capacitance $C$. This model has been applied to ultra-small Josephson junctions Vion et al. (1996) and S-QD-S setups Jørgensen et al. (2007). Disregarding $C_{\mathrm{filter}}$ in the circuit depicted in Fig. 2c, Kirchhoff’s laws and the Josephson relation yield two coupled Langevin equations:

[TABLE]

Here $u=V_{\mathrm{sd}}/RI_{\mathrm{c}}$ denotes the dimensionless voltage across the shunt capacitor $C$ and we have introduced the dimensionless time $\tau=\omega_{RL}t$ with $\omega_{RL}=R/\textup{L}$, where $\textup{L}=\Phi/I_{\mathrm{c}}$ is the self-inductance of the Josephson junction set by the critical current, $I_{\mathrm{c}}$, and $\Phi=\hbar/2e$. The two different $RC$ frequencies, $\omega_{RC_{\mathrm{J}}}=(RC_{\mathrm{J}})^{-1}$ and $\omega_{RC}=(RC)^{-1}$, enter these equations via the dimensionless parameters $\alpha=\omega_{RL}/\omega_{RC}=R^{2}CI_{c}/\Phi$ and $\alpha_{0}=\omega_{RC_{\mathrm{J}}}/\omega_{RL}=\Phi/R^{2}C_{\mathrm{J}}I_{\mathrm{c}}$. Finally, $L(\tau)$ is a dimensionless stochastic parameter describing white-noise voltage fluctuation of $V_{\mathrm{sd}}$. $L(\tau)$ is composed of Nyquist noise from the series resistor ($R$) at temperature $T^{*}$ and a stray voltage noise characterized by a variance, $K$, to account for imperfect filtering, such that $\langle L(\tau)L(\tau+\Delta\tau)\rangle=\left(\frac{2k_{B}T^{*}}{\Phi I_{c}}+\frac{K}{RI_{c}\Phi}\right)\delta(\Delta\tau)$. Noise from the source resistance $R_{\mathrm{s}}$ is assumed to be excluded by the RC filters. As $K$ and $T^{*}$ can not be independently measured in this setup we ascribe all noise to the series resistor ($R$) with an effective temperature $T=T^{*}+K/(2Rk_{B})$ which can be different from the fridge temperature.

For $\alpha\gg 1$, $u$ changes slowly with time, which allows us to solve for $\phi$’s equilibrium distribution keeping $u$ constant. If furthermore $\alpha_{0}\gg 1$, solving the stochastic equation for the $\phi$-distribution is equivalent to the original RCSJ problem Ivanchenko and Zil’berman (1969) with the solution,

[TABLE]

where $A=\Phi I_{\mathrm{c}}/k_{B}T$ and $I_{n}(z)$ is a Bessel function. Inserting this solution into Eq. (1) a constant average voltage, $\dot{\langle u\rangle}=0$, may be enforced by satisfying the equation

[TABLE]

Two limits of $R/R_{\mathrm{s}}$ are of interest. For small $R/R_{\mathrm{s}}$, depending on $I_{\mathrm{c}}$ and $T$, this equation has three solutions for $\langle u\rangle$ arising from the non-monotonic behavior of the right-hand side, and the setup is effectively current biased, exhibiting hysteretic $\langle I\rangle-V$ curves. For large $R/R_{\mathrm{s}}$, the solution for $\langle u\rangle$ is unique and the system is effectively voltage biased and can exhibit NDC. The latter case is consistent with the experiment.

Experimental $I-\langle V_{\mathrm{sd}}\rangle$ traces can now be fitted to

[TABLE]

which is identical to an earlier formula derived by Ivanchenko and Zil’berman Ivanchenko and Zil’berman (1969), while the derivative of this formula with respect to $\langle V_{\mathrm{sd}}\rangle$ allows fitting of $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ traces. In a typical Josephson junction $A\gg 1$ and the cusp voltage scale is set by $V_{\mathrm{0}}\approx RI_{\mathrm{c}}$. In our measurements $V_{\mathrm{0}}$ is constant over a wide range of cusp current, $I_{\mathrm{0}}$, as we show in Fig. 12 where extracted $I_{\mathrm{0}}$, $V_{\mathrm{0}}$ and fitted $I_{\mathrm{c}}$ is shown for two linecuts in the experiment.

This is consistent with the large-noise limit, $A<1$, of Eq. (5) where one can perform a Taylor expansion of Eq. (5) around $A\approx 0$ to obtain Ivanchenko and Zil’berman (1969),

[TABLE]

from which it is clear that the cusp voltage $V_{\mathrm{0}}=\operatorname{arg}\max\langle I\rangle(\langle V_{sd}\rangle)=Rk_{B}T/\Phi$ is solely determined by noise and is independent of $I_{\mathrm{c}}$. Cusp current,

[TABLE]

and zero-bias conductance,

[TABLE]

can also easily be identified in this limit. We also identify a scaling $\gamma$ with $\bar{R}=R/\gamma$, $\bar{T}=T\gamma$ and $\bar{I}_{\mathrm{c}}=\sqrt{\gamma}I_{\mathrm{c}}$ which keeps $V_{\mathrm{0}}$ and Eq. (6) unchanged.

If the limits above are satisfied we indeed have a model where: 1) The junction is overdamped and non-hysteretic, 2) the junction is effectively voltage biased allowing measurements of NDC and 3) the $I-V_{\mathrm{sd}}$ curves are cusped with an $I_{\mathrm{0}}$-independent cusp voltage $V_{\mathrm{0}}$. We will now discuss the validity of the above limits. By fixing $T=80$ mK, consistent fits of Eq. (5) can be done for all gate ranges yielding a constant series resistance $R=3$ k$\Omega$ and variable $I_{\mathrm{c}}$ in the range $I_{\mathrm{c}}\sim 0.02-3$ nA which roughly corresponds to $\alpha\sim 10-100$ and $\alpha_{0}\sim 10^{4}-10^{5}$ using the geometric estimates of $C$ and $C_{\mathrm{J}}$, consistent with our assumption that $\alpha,\alpha_{0}\gg 1$. For these parameters Eq. (4) has only one solution, consistent with the junction being voltage biased. Lastly, for these parameters $A\leq 1$ for the fitted range of $I_{\mathrm{c}}$ with $A\approx 1$ for the largest $I_{\mathrm{c}}$. In Fig. 13, we compare the simple large-noise-limit expressions for $I_{\mathrm{0}}$ and $G$ using $I_{\mathrm{c}}$ from the fitting with measured $I_{\mathrm{0}}$ and $G$ and find excellent agreement with the largest discrepancies occurring for large $I_{\mathrm{c}}$ values consistent with our Taylor expansion in $A\propto I_{\mathrm{c}}$.

As the above analysis is a fully classical treatment we will now discuss the impact of quantum fluctuations. The characteristic frequencies of the circuit are estimated to be $\omega_{RC}\sim 0.01$ GHz, $\omega_{RL}\sim 0.1-10$ GHz, and $\omega_{RC_{\mathrm{J}}}\sim 10^{5}$ GHz. With $T=80$ mK, i.e. $k_{B}T/\hbar\sim 10$ GHz, and $I_{\mathrm{c}}=0.02-3$ nA, we arrive at the following hierarchy of frequencies:

[TABLE]

For these values, one may safely neglect the large capacitance, $C$, in the effective circuit impedance experienced by the junction

[TABLE]

and which determines the mean-squared phase fluctuations, $S_{\phi}=\langle(\phi-\langle\phi\rangle)^{2}\rangle$, which is given by Refs. Grabert et al., 1998; Grabert and Ingold, 1999; Joyez, 2013,

[TABLE]

with $R_{\mathrm{Q}}=h/4e^{2}\approx 6.5$ k$\Omega$. Within a fully quantum mechanical treatment, this correlation function in turn leads to an approximately exponential reduction of the critical current Grabert et al. (1998); Joyez (2013):

[TABLE]

As a function of temperature, $S_{\phi}$ is roughly constant and due to quantum fluctuations for $k_{B}T<\hbar\omega_{RL}$. For $\hbar\omega_{RL}<k_{B}T$, corresponding to the parameters found above, the fluctuations are largely classical and $S_{\phi}$ increases linearly with $T$, consistent with the classical treatment employed above. Since $k_{B}T$ is just barely larger than $\hbar\omega_{RL}$, and $R$ just barely smaller than $R_{Q}$, a further reduction of $I_{\mathrm{c}}$ due to quantum fluctuations must be expected to reduce the actual values for $I_{\mathrm{c}}$, which we deduce from the strictly classical analysis above and report in the figures throughout the main text. Nevertheless, the marked gate dependence of $I_{\mathrm{c}}$ remains and this is what provides the real value of the zero-bias conductance peak as a probe of YSR screening and the concomitant quantum phase transitions reported in this paper.

Summarizing the circuit analysis, the magnitudes of $I_{\mathrm{c}}$ obtained by fitting $I-V_{\mathrm{sd}}$ and $dI/dV_{\mathrm{sd}}-V_{\mathrm{sd}}$ traces are expected to be slightly overestimated by leaving out quantum, and retaining only thermal fluctuations of the electromagnetic environment. Additionally, $T$ has been fixed by an independent measurement which is only an estimate of the effective temperature experienced by the circuit. Relaxing this constraint on $T$, we can still obtain consistent fits utilizing a global scaling of $R$ and $I_{\mathrm{c}}$. The value of these measurements therefore lies in the gate dependence of $I_{\mathrm{c}}$, not its amplitude which can include a global scaling, in good qualitative agreement with NRG calculations.

Appendix C: Asymmetries in cusp current

While the model displays $I_{\mathrm{0}}$ symmetric in $V_{\mathrm{sd}}$, in some instances the experiment fails to do so, with the discrepancy attributed to conduction through low-lying YSR subgap peaks which approach $V_{\mathrm{sd}}=0$. In Fig. 14 we show two examples in which subgap peaks approach zero bias, and the effect which this has on the symmetry of $I_{\mathrm{0}}$ in $V_{\mathrm{sd}}$. The extreme case when a pair of subgap peaks reaches all the way to $V_{\mathrm{sd}}=0$ at two singlet-doublet parity crossings is shown in Fig. 14b, and its effect on $I-V_{\mathrm{sd}}$ traces is shown in Fig. 14a. Blue arrows indicate the subgap peak of largest conductance among the ones at positive and negative bias for a given gate voltage. To the left of state $D_{\mathrm{L}}$, this peak is at negative bias, whereas to the right it is at positive bias; i.e., the peak of largest conductance is antisymmetric with respect to the center of the state $D_{\mathrm{L}}$. Unsurprisingly, the same antisymmetry is observed for $I_{\mathrm{0}}$, with the negative-$I$ (positive-$I$) cusp occurring at larger $\absolutevalue{I}$ than its positive-$I$ (negative-$I$) counterpart to the left (right) of state $D_{\mathrm{L}}$. This occurs between the pairs of gate points marked by blue pins. This behavior contrasts with the hysteresis of underdamped Josephson junctions, where a sweep of $V_{\mathrm{sd}}$ from negative to positive values, as it is always the case in our device, would result in larger switching current for positive than negative $I$. Along the gate voltage delimited by the pair of red pins, the subgap states remain close to $V_{\mathrm{sd}}=0$, affecting the positive-negative symmetry of the $I$ cusps, which are just barely well-defined in this region. To the left and right of the parity crossings, in the $S$ states, YSR subgap peaks are far enough away from the supercurrent zero-bias conductance peak, and $I-V_{\mathrm{sd}}$ traces are consequently fully symmetric in $V_{\mathrm{sd}}$.

Figs. 14c,d, show a different example with subgap peaks appearing near $V_{\mathrm{sd}}=0$ in Fig. 14d. Nevertheless, these are far enough away to not affect the cusp symmetry in the respective $I-V_{\mathrm{sd}}$ traces in Fig. 14c.

It is natural to expect that the appearance of such additional sub-gap conductance channels close to zero bias affects the resistance of the junction, and therefore the effective relation between $\langle I\rangle$ and $\langle V_{\mathrm{sd}}\rangle$. The modelling of this additional complication, however, is beyond the scope of this work.

Appendix D: Two-impurity model

We model the double quantum dot system using the following superconducting generalization of the two-impurity Anderson impurity model:

[TABLE]

where the impurity Hamiltonians for quantum dots $i=L$ (left) and $i=R$ (right) with on-site energies $\epsilon_{i}$ and electron-electron repulsion (Hubbard) parameters $U_{i}$ are

[TABLE]

where the impurity occupancy (charge) operators are $\hat{n}_{i,\sigma}=d^{\dagger}_{i,\sigma}d_{i,\sigma}$, $\hat{n}_{i}=\hat{n}_{i,\uparrow}+\hat{n}_{i,\downarrow}$, while $n_{i}=1/2+\epsilon_{i}/U_{i}$ is the impurity energy level $\epsilon_{i}$ (controlled by the corresponding gate voltage) expressed in the dimensionless units of electron number. The inter-dot tunnelling amplitude, $t_{\mathrm{d}}$, enters through a tunnelling term,

[TABLE]

the leads are modelled using BCS Hamiltonians

[TABLE]

$U_{\mathrm{d}}$ is set to zero, and the QD-lead tunneling amplitudes, $t_{\mathrm{{L/R}}}$, enter through the hybridisation terms

[TABLE]

In these expressions, $d_{i\sigma}$ are QD operators, while $c_{i,k\sigma}$ are lead operators. Finally, the tunnelling rates are defined as $\Gamma_{i}=\pi\rho_{i}t_{i}^{2}$, where $\rho_{i}$ is the density of states in lead $i$ in the normal state (i.e., in the absence of superconductivity).

Appendix E:

Finite-bias conductance at weak inter-dot coupling

For weak $t_{\mathrm{d}}$, $\mathrm{d}I/\mathrm{d}V_{sd}$ can be estimated using the tunneling (Bardeen) approach, similar to that used in the analysis of scanning tunneling microscopy (STM) experiments. We assume that the voltage drops entirely across the weak link between the two dots, and that each quantum dot is in thermodynamic equilibrium with the neighboring lead. The chemical potential shifts in the leads are $\mu_{\mathrm{{L,R}}}=\pm eV_{\mathrm{sd}}/2$, and the local spectral functions on the dots are equally shifted. The superconductors are described in the semiconductor model, neglecting all coherence effects. The current is then given by

[TABLE]

where $f(x)=[1+\exp(x/k_{B}T)]^{-1}$ is the Fermi function. At zero temperature, this becomes

[TABLE]

We computed the spectral functions $A_{\mathrm{L,R}}(\omega)$ for the two subsystems with $t_{\mathrm{d}}=0$ using the density-matrix numerical renormalization group algorithm with discretization parameter $\Lambda=2$, and broadened raw spectra using the broadening parameter $\alpha_{\mathrm{b}}=0.2$. After convolution of the spectral functions to obtain the current, $dI/dV_{\mathrm{sd}}$ was computed by taking numerical differences. This single-particle tunneling approach, valid only to second order order in $t_{\mathrm{d}}$, neglects Andreev reflections altogether and presumes that relaxation rates of quasiparticles in the sub-gap states are larger than the inter-dot tunneling rate.

Appendix F: Effects of finite inter-dot charging energy and gate-dependent tunnel coupling

In this section we investigate the effects of a possible influence of $V_{\mathrm{gL}}$ on $t_{\mathrm{R}}$. For efficient parameter exploration, we do this using a zero-bandwidth (ZBW) approximation Grove-Rasmussen et al. (2018) of the superconducting leads in Eq. (13), taking $H_{\mathrm{BCS_{i}}}\approx\sum_{i}\Delta c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}+\operatorname{h.c.}$ and numerically diagonalizing the Hamiltonian. We model the tunnel couplings, $t_{i}$, as $t_{\mathrm{L}}=\sqrt{2\Delta\Gamma_{\mathrm{L}}/2}$ and $t_{\mathrm{R}}=\sqrt{2\Delta\Gamma_{\mathrm{R}}/2}(1+n_{\mathrm{L}}/4)$, where $t_{\mathrm{R}}$ now depends on $V_{\mathrm{gL}}$ through $n_{\mathrm{L}}$. To more realistically mimic the experimental scenario we also include an inter-dot Coulomb repulsion, $U_{\mathrm{d}}$, such that $H_{\mathrm{imp}}=\sum_{i=\mathrm{L,R}}U_{i}/2(\hat{n}_{i}-n_{i})^{2}+U_{\mathrm{d}}(\hat{n}_{\mathrm{L}}-n_{\mathrm{L}})(\hat{n}_{\mathrm{R}}-n_{\mathrm{R}})$ with $U_{\mathrm{d}}=0.5$ meV and $\hat{n}_{i}=\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}$ being the electron counting operator. Other parameters are identical to parameters used for the NRG calculations presented above. Comparing Fig. 16a to Fig. 5a, the main effect of $U_{\mathrm{d}}$ is to provide an angle on the $n_{\mathrm{R}}$, $n_{\mathrm{L}}$ parity lines.

In Figs. 15, 16, and 17 we show GS parity transition lines (black lines) of the S-DQD-S charge stability diagrams for various values of $\Gamma_{\mathrm{L}}$ and $\Gamma_{\mathrm{R}}$, corresponding qualitatively to data in Figs. 6, 7 and 8. The gate-dependent coupling is seen to account in Fig. 16 for the asymmetric screening around the electron-hole symmetry point occurring when $\Gamma_{\mathrm{L}}$ is increased, which is necessary to explain the experimental findings in Fig. 7, that is, the observation that the doublet to singlet transition reflected by the vanishing of the two parity transition lines happens faster on the top ($V_{\mathrm{gL}}\approx-2.1$V) than at the bottom ($V_{\mathrm{gL}}\approx-2.3$V).