Excitations in a superconducting Coulombic energy gap

Cooper pairing and Coulomb repulsion are antagonists, producing distinct energy gaps in superconductors and Mott insulators. When a superconductor exchanges unpaired electrons with a quantum dot, its gap is populated by a pair of electron–hole symmetric Yu-Shiba-Rusinov excitations between doublet and singlet many-body states. The fate of these excitations in the presence of a strong Coulomb repulsion in the superconductor is unknown, but of importance in applications such as topological superconducting qubits and multi-channel impurity models. Here we couple a quantum dot to a superconducting island with a tunable Coulomb repulsion. We show that a strong Coulomb repulsion changes the singlet many-body state into a two-body state. It also breaks the electron–hole energy symmetry of the excitations, which thereby lose their Yu-Shiba-Rusinov character.


INTRODUCTION
In a large superconductor, an adsorbed spin impurity binds to a quasiparticle screening cloud to form a state known as the Yu-Shiba-Rusinov (YSR) singlet, whose excitation energy with respect to the unbound doublet is below the superconducting energy gap, ∆ 1 .The miniaturization of the superconductor into an island reduces charge screening and introduces an energy gap for the addition of electrons, the Coulomb repulsion, E c (see Fig. 1a) 2,3 , with yet unexplored consequences on the ground state and the subgap spectrum.Such exploration is of relevance in the study of magnetic impurities adsorbed to superconducting droplets 4,5 , in quantum-dot (QD) readout of Majorana qubits based on superconducting islands [6][7][8] , and in realizations of superconducting variants of the multichannel Kondo model [9][10][11] .
In the absence of a spin impurity, the charging of a superconducting island (SI) depends on the ratio E c /∆, with E c /∆ < 1 leading to Cooper pair (2e) charging and E c /∆ > 1 to 1e charging 2,3 .In the latter case, even numbers of electrons condense as Cooper pairs, while a possible odd numbered extra electron must exist as an unpaired quasiparticle 3 .
Here we provide the first spectral evidence of the manybody excitations in a superconducting Coulombic gap.The spin impurity resides in a gate-defined QD in an InAs nanowire, and the SI is an Al crystal grown on the nanowire with gate-tunable Coulomb repulsion.Both QD-SI and SI-QD-SI devices are investigated in this work.We demonstrate that a strong Coulomb repulsion forces exactly one quasiparticle in the SI to bind with the spin of the QD in the singlet ground state (GS).The Coulomb repulsion also enforces a positive-negative bias asymmetry in the position of the excitation peaks which is uncharacteristic of YSR excitations.
Figures 1b-e summarize the energy dispersions which can arise when a QD is coupled to a superconductor.In Fig. 1b, the usual YSR case (E c = 0) with the QD gateinduced charge tuned to ν = 1 is depicted.The doublet GS and singlet excited state energies are independent of the gate-induced charge in the superconductor, n 0 , and excitations between these two states are electron-hole symmetric 12 .As shown in Fig. 1c, introducing E c > 0 in the superconductor produces a parabolic dispersion distorted by the hybridization (Γ) between the QD and the SI, which couples states of the same total charge.For odd n 0 , the energy of the doublet state is increased by ≈ E c (green dot), while for even n 0 it is the energy of the singlet state which is penalized by this amount.For odd n 0 and E c > ∆, the GS is a singlet even if Γ → 0. For E c < ∆, the singlet can be the GS if the YSR binding energy E B is large enough so that E c > ∆ − E B , which is achieved by increasing Γ 13 .
Due to U > 0, the dispersion against ν is approximately parabolic in both the E c = 0 (up to a constant) and E c > 0 cases, as shown in Figs.1d,e.For E c = 0 (Fig. 1d), the electron and hole excitations are symmetric due to the degeneracy of the even-parity parabolas.This ceases to be the case for E c > 0 (Fig. 1e).The asymmetry is maximal in the absence of additional QD levels.For ν > 1 (ν < 1), an extra electron (hole) must be stored in the SI with excitation energy ∆ + E c , but an extra hole (electron) can be added to either the QD or the SI, leading to a superposition of states with excitation energies −( d + U ) and −(∆ + E c ), where d is the energy level of the QD (details on Extended Data Fig. 1).The extra electron or hole either forms a quasiparticle or a Cooper pair, depending on the parity of the SI occupation of the initial state.Superconducting Coulombic  scheme shown in Fig. 2c 13 .The SI is conceived as several hundreds of electronic levels.Its charge is tuned by n 0 , equivalent to top gate voltage V S in the device.The corresponding Hamiltonian includes pairing between time-reversed states to produce the superconducting gap, ∆, and coupling to the QD, Γ, which is tuned by top gate voltage V 3 in the device.We consider constant Coulomb interactions U for the QD, E c for the SI, and V for the interdot charging due to the QD-SI inter-capacitance, C m (as in usual double QDs 14 ).The QD is itself modelled as an Anderson impurity, whose charge is tuned by ν, equivalent to top gate voltage V N .Other top gates (V 1 , V 5 ) control the couplings of the QD and SI to the source and drain, not included in the model.The output of the model is the energy spectrum of the system, consisting of a few low-lying many-body states and the edge of the continuum.These states are sketched in Fig. 2d between the source and drain.Table I shows device and model parameters.
To record the spectrum of excitations between the low-lying many-body states, the device is biased by a source-drain bias voltage, V sd , and the differential conductance, G, is measured at the grounded drain, as shown in Fig. 2d.Asymmetric source and drain couplings are needed for G to embody the energy asymmetry of the SCE.While we cannot account for the magnitude of G, we expect that the measured G(−V sd ) reflects the excitation energies at eV sd , where the negative sign in the argument is necessary to account for the voltage drops and polarity conventions.Symmetric barriers would instead result in a trivial bias symmetry.A peak at one polarity thus demonstrates the existence, at the corresponding excitation energy, of an excited many-body state with n GS + 1 electrons in the device, while a peak at the opposite polarity that of a many-body state with n GS − 1 electrons, where n GS is the total charge in the GS.
The zero-bias G signal exhibits a strong dependence on V S and V N , as shown in the diagram of Fig. 3a.Singlet↔doublet GS transitions are observed in the experiment when conductance lines are crossed, as at V sd = 0 these lines appear when the n GS and n GS +1 (or n GS −1) states in Fig. 2d are degenerate at zero energy.The repetition of the central hexagonal charge domain in the V S direction indicates filling of the SI.As a guide of the filling of the QD and the SI, we approximate their charge expectation values as integers n N , n S in each of the charge domains.This is an approximation as only n GS is integer with n GS = n N + n S 13 (see Extended Data Fig. 1).Small but resolvable 1,1 singlet domains (an example is enclosed in a dotted line) are seen between the 1,2 and 1,0 doublet domains.In contrast, the lines to the sides of the central hexagonal domains, which separate the 0,0 and 0,2 domains and the 2,0 and 2,2 domains, show no splitting at this resolution.The difference stems from finite Γ and V , which stabilize the 1,1 but not the 0,1 and 2,1 domains.The presence of the 1,1 Coulomb-aided YSR singlet is the key difference from a trivial double QD stability diagram 14 and from the ∆ = 0 case 9 .For instance, a raise of the interdot coupling in a double QD introduces molecular orbitals which show as avoided crossings at triple points (TPs), whereas in the QD-SI system the YSR singlet is a many-body state for these parameters.Finite Γ and V are also responsible for increasing the distance between the points of multiple degeneracy, for the acute angle between vertical and horizontal conductance lines and, in the case of Γ, for curving the conductance TABLE II.Effective and bare g-factors of the two components of the QD-SI device.Effective g-factors of the QD, gN, and of the SI, gS, extracted from the data in Fig. 3a (for methods, see Section I of Supplementary Information), and bare g-factors g0N, g0S used as input in the corresponding finite B calculation.
Our model of the system produces a diagram of GS transitions of the SCE that matches the gate position of the conductance lines, as shown in the comparison of the calculation to the experimental data in Fig. 3a.The quality of the match for model parameters approximately similar to the experimentally measured values (with ∆ as the only fit parameter) constitutes a first proof of the presence of SCE in our device.
We corroborate the spin (S z ) assignment done in Fig. 3a (right panel) at B = 0 from the variation of GS domain sizes with B = 0.3 T (in inset).Doublet domains are stabilized by B more than singlet domains, while triplet domains are stabilized further than doublet and singlet domains.The model fits the data using the g-factors as free parameters, and taking into account the GS transition from singlet to triplet in the 1,1 charge sector (charge parabolas are shown in Extended Data Fig. 2).The g-factors in the Hamiltonian are significantly larger than the measured effective g-factors (see Table II).These bare g factors produce Zeeman splittings E Z,QD = g 0S µ B B and E Z,SI = g 0N µ B B in the QD and the SI, where µ B is the Bohr magneton.The effective and bare g factors would be equal if the expectation values of the QD and the SI charges increased in steps of exactly 1e across the GS transition lines.For non-zero Γ, (non-integer) charge distributions occur between the QD and the SI on either side of the GS transition line (see Extended Data Fig. 1), hence the effective g-factors are some non-trivial function of the true (bare) g-factors which appear in the Hamiltonian 15 .
Following this comprehensive mapping, we show in Fig. 3b the G spectrum at finite V sd versus V N for fixed V S , at which the SI contains only Cooper pairs in the GS up to a good approximation.The SCE have a double-S shape, spanning V sd ≈ −0.37 → 0.37 mV.They are approximately inversion symmetric in position and in G intensity with respect to the electron-hole gate-symmetric filling point of the QD, which corresponds to the center of the 1,0 sector (indicated by a cross), from where removing/adding an electron from/to the QD are equally energetically unfavorable.G jumps in intensity when the SCE cross zero bias, as highlighted by the insert traces at gate points before (grey) and after (black) one of such changes.While the SCE are expected to appear as a pair at asymmetric positive and negative bias positions for a given gate voltage, in practice only one SCE is observed.A GS change brings discontinuously up to the continuum the other state, as charge is suddenly redistributed between the QD and the SI 13 .
Our model reproduces the position of the subgap resonances, as evidenced in the overlay of the calculated spectrum on the experimental data in Fig. 3b (see also Extended Data Fig. 3).Differences between the SCE spectrum and the spectrum in the Coulombic (∆ = 0) and Yu-Shiba-Rusinov (E c = 0) limits are shown in Extended Data Fig. 4. The Coulombic spectrum bears resemblance to that of an impurity in the paramagnetic Mott insulator described by the Hubbard model 16,17 , despite the differences in the Hamiltonian (local Hubbard interaction versus constant Coulomb repulsion in our model).In both cases the charge transfer from the impurity site to the bath costs energy corresponding to the total charge gap of the system in the absence of the impurity (≈ U/2 in the Hubbard model at half-filling, E c + ∆ in our device), and in both cases there is a (quasi)continuum of fermionic states extending above this gap (doublons/holons in a Mott insulator, and Coulomb quasiparticles with a mixed character of Bogoliubov quasiparticles due to ∆ in our device), leading to the same phenomenology.
To map these limits, we vary continuously the Coulomb repulsion in the superconductor in a second device.We first explore the role of the Coulomb repulsion on the stability of the YSR singlet as the GS.To this aim, we define two quantities, x = 1 − ∆/E c , and y, the YSR singlet GS size in units of e.In the Γ/U 1 regime, y = (E c − ∆ + E B )/E c .Figs. 4a,b explain how x and y are experimentally extracted.In the limits when E c → 0 and E c → ∞, x → −∞ and x → 1, respectively.When E c = ∆, then x = 0. Fig. 4c shows a measurement of y versus x in a device consisting of a QD coupled to two SIs with hybridization Γ L and Γ R (top inset in Fig. 4c).The SIs have charging energies E cL and E cR and superconducting gaps ∆ L and ∆ R , and their occupations are tuned with top gate voltages V SL and V SR .The advantage of this three-component device over the two-component one is that the presence of only one QD between the two SIs can be verified from stability diagrams similar to that in Fig. 3a against the pairs of gate voltages (V N ,V SL ) and (V N ,V SR ).In Fig. 4c, y characterizes the GS stability of the YSR state formed by the binding of the QD spin to the quasiparticle cloud in the right SI, and x = 1 − ∆ R /E cR .To employ the device as this two-component system, the left SI is kept either as a cotunnnelling probe at even occupation (for E cL > ∆ L ) or as a soft-gap superconducting probe (for E cL ∆ L and Γ L Γ R ).Three QD shells with different values of Γ R /U are analyzed.At the weakest Γ R /U , y shows a trivial linear dependence with a slope of 1 and with endpoints at (0,0) and (1,1), connected by a fitted solid line in the graph.In this regime, x only stabilizes the YSR singlet as the GS.At the other extreme, at the largest Γ R /U , x stabilizes the doublet more strongly than the YSR singlet, reducing y.In between these two extremes, at Γ R /U = 0.3, the behavior is intermediate.When x → 1, y converges to 1e independently of Γ R /U , as x stabilizes equally well the doublet and singlet states for even and odd gate-induced charges in the right SI.In the other limit, when x → −∞, y depends exclusively on Γ R /U , as in the usual E cR = 0 YSR regime.
Next, we describe how the Coulomb repulsion in the superconductor affects the dispersion of the excitations and how this is related to changes in the stability diagram.In Fig. 5, we show the evolution of the excitations produced by one QD shell on the right SI over a wide range of V SR , corresponding to a charge variation of ≈ 960 electrons.In this range, E cR /∆ R goes from 0 to 1.71, as measured from Coulomb-diamonds spectroscopy.The increase in E cR /∆ R is reflected on the stability diagram.In the usual E cR ≈ 0 YSR regime (Fig. 5a), the diagram shows two vertical dispersionless lines, and the spectrum consistently displays a YSR loop (for measurement details, see Methods).When the right SI enters into Coulomb blockade (Fig. 5b, E cR /∆ R = 0.36), the lines in the stability diagram wiggle as interdot charging and tunnelling effects enter into consideration.Consequently, the YSR loop in the spectrum gets skewed rightwards and increases its bias size as the energy gap includes now a Coulombic component.At E cR /∆ R = 0.75 (Fig. 5c), the entrance of the 1,1 YSR singlet GS breaks the stability diagram into several domains, and the excitations adopt a double-S shape.At this setting, the 1,0 doublet domain has a , where E BR is the YSR binding energy of the spin in the QD to the quasiparticle cloud in the right SI.This results in a maximum of the bias size of the double-S shape excitation.From then on, an increase in E cR /∆ R in Figs.5d-f reduces the energy of the doublet→singlet excitation and the double-S shaped feature shrinks in bias size, concomitantly with the stronger stabilization of the YSR 1,1 singlet GS in the stability diagram.
Throughout this article, we provided compelling evidence for the existence of superconducting Coulombic subgap excitations arising from states bound to a semiconductor-superconductor interface, and we showed how these are related to the usual electron-hole symmetric Yu-Shiba-Rusinov excitations.On one hand, we showed that a small Coulomb repulsion in the superconductor is enough to turn the excitations asymmetric in the polarity of the bias voltage.On the other hand, a strong Coulomb repulsion (E c → ∞) converts the YSR singlet many-body state into a two-body state formed by a spin in the QD and a single quasiparticle in the superconductor.
Though our model is successful at matching excitation energies, an extension which includes transport is needed to account for the magnitude of the conductance features and for their bias positions in devices with more symmetric source-drain barriers.The observation of current blockade in a regime of weaker Γ hints at elastic cotunnelling as the transport mechanism in our QD-SI de-vice (see Section II of Supplementary Information).The absence of zero-bias G in the ∆ > E c regime (e.g.Extended Data Fig. 5) indicates that Andreev reflection (2e charge transfer) is not a transport mechanism in our devices in this regime.Due to charge transfer between the SI and the QD, the model indicates that an upwards reconsideration of bare g-factor values extracted from experimentally-determined g-factors of subgap excitations is needed to match the experimental results 6,18 .Based on its success, the model can also inform on future developments, e.g.qubit and multi-channel devices which utilize the SI-QD-SI device, as outlined in Ref. 19 .
The system sketched in Fig. 1a may also be realized with magnetic adatoms on superconducting droplets (e.g.Pb on an InAs substrate), and probed with scanning tunnelling microscopy 4,5 .Several open questions could be answered with this technique: What is the spatial extension of the excitations in a superconducting Coulombic energy gap 20 ?Is there orbital structure in the excitations 21 ?How do the excitations behave in chains of magnetic adatoms when these chains are deposited on top of a E c > 0 SI 22 ?Do chains of magnetic adatoms deposited on finite E c SIs support Majorana excitations 23 ?

METHODS
Devices fabrication and layout.QD-SI device (Fig. 2a).A 110-nm wide InAs nanowire picked with a micromanipulator was contacted by 5/200 nm Ti/Au (in yellow) source and drain leads.The ≈350-nm long, 7-nm thick epitaxially-grown Al SI covering three facets of the nanowire was defined by chemically etching the upper and lower sections of the nanowire before contacting.After insulating the nanowire and the leads with a 6-nm thick film of HfO 2 , five Ti/Au top gates were deposited along the nanowire.The QD was defined in the bare nanowire next to the SI by setting top gates 1 and 3 to negative voltage.A Si/SiO 2 substrate backgate was kept at zero voltage throughout the experiment.
SI-QD-SI device (Extended Data Fig. 5).The SI-QD-SI device was fabricated using a nanowire from the same growth batch.Two nominally identical 7-nm thick, ≈300nm long, epitaxially-grown Al SIs were defined by chemical etching.The nanowire was contacted by 5/200 nm Ti/Au leads, and then insulated by a 5-nm thick layer of HfO 2 from seven Ti/Au top gates deposited after.The QD was defined between the two SIs by setting top gates 3 and 5 to negative values.The substrate backgate was used to aid the top gates in depleting the device.E cR /∆ R was tuned by using an auxiliary QD (QD R aux ) defined between the right SI and the source lead.When QD R aux was put near resonance by sweeping V SR , E cR − ∆ R could be tuned to negative values, and when QD R aux was put in cotunnelling, E cR − ∆ R could be tuned to positive values.Similarly, E cL − ∆ L was tuned using an auxiliary QD defined between the left SI and the drain lead.
The critical B of the superconducting Al film was mea-sured to be B c = 2.1 T in nanowire devices made from the same batch of nanowires used in the fabrication of the present device 24,25 , which left ample room for B-resolved measurements in the superconducting state.In the QD-SI and SI-QD-SI devices, the presence of superconductivity at large B was determined from size differences of adjacent charge domains with odd and even occupation of the SI, observed up to B = 1.2 T and B = 1.5 T, respectively.Larger B was not explored.

Differential conductance measurements.
A standard lock-in technique was used to measure the differential conductance, G = dI/dV sd , of the QD-SI device by biasing the source with an AC excitation of 5 µV at a frequency of 223 Hz on top of a DC source-drain bias voltage, V sd , and recording the resulting AC and DC currents on the grounded drain lead.In the case of the SI-QD-SI device, G was measured at the grounded drain with a 5 µV lock-in excitation applied at the source at 84.29 Hz.The measurements were performed in an Oxford Triton dilution refrigerator at 30 mK for the QD-SI device and 35 mK for the SI-QD-SI device, such that k B T E c , where k B is the Boltzmann constant and T is the refrigerator temperature.
Calculation of subgap and continuum excitations.The calculations were done using the densitymatrix renormalization group approach (details in Section III of Supplementary Information).The quantum numbers are the total number of electrons in the system, n, the z-component of the total spin, S z (see Extended Data Fig. 2), and the index for states in a given (n, S z ) sector, i = 0, 1, . ... The superconducting Coulombic excitation energies are given by E = E(n = n GS ± 1, S z = S z,GS ± 1/2, i = 0) − E(n = n GS , S z = S z,GS , i = 0).The edges of the continuum excitations are given by E edge = E(n = n GS ± 1, S z = S z,GS ± 1/2, i = 1) − E(n = n GS , S z = S z,GS , i = 0).Due to the finite size of the SI, the continuum is in truth only a quasi-continuum of states.The nature of these states and the excitation energies depend on the values of ∆ and E c .For ∆ = 0, the quasiparticles are free-electron states.For ∆ = 0, these are Bogoliubov quasiparticles with pronounced inter-level pairing correlations c † i↑ c † i↓ c j↓ c j↑ .If E c = 0, the excita-tion spectrum is not affected by the number of preexisting particles in the superconductor (up to finite-size effects).
If E c = 0, the particle-addition and particle-removal energies are affected by the charge repulsion (parabolas).The calculations do not provide direct results for the differential conductance of the system, only information about the energies of the GS and the low-lying excitations.
Spectral measurements in Fig. 5.To obtain sharp spectral features visible over the continuum background, we tuned the left SI into a superconducting probe (E cL = 0, Γ L ≈ 0).The strong hybridization of the left SI with the drain needed to achieve E cL = 0 resulted in an unintended soft gap in this probe, which produced faint replica of the main excitations.For example, in Fig. 5a, black dotted lines correspond to the YSR loop coming from the QD-right SI being probed by the coherence peaks of the probe.The loop is thus followed by negative differential conductance (NDC) and appears at ∓eV sd = ±∆ L + E(n = n GS ± 1, S z = S z,GS ± 1/2, i = 0) − E(n = n GS , S z = S z,GS , i = 0), reaching ∆ L at GS transitions.The gray dotted lines highlight a YSR replica probed by the soft gap of the probe, thus an order of magnitude weaker in conductance and without associated NDC.This replica appears at ,GS , i = 0), and therefore crosses zero bias at GS transitions.When E cR /∆ R increases, the relationships between the excitations and the bias positions of the conductance features become approximations due to nonideal Γ source , Γ drain asymmetry in the device, which also results in negative-slope features included in the dotted lines but not accounted by the model, as the model does not consider transport.The interpretation of features in the continuum at higher bias is outside the scope of this work.The asymmetry comes from having only one level in the QD which can be filled with up to 2 electrons.As a consequence, in the nGS + 1 state an additional electron must go into the SI.In contrast, in the nGS − 1 state an electron can be removed from the SI or from the QD.In general, the charge excitations are electron-hole asymmetric (Coulombic-like) for ν = 1.G(e 2 /h) 0 0.07

FIG. 1 .FIG. 2 .
FIG. 1.From Yu-Shiba-Rusinov to Coulombinfluenced excitations.a Idealized system displaying Coulomb-influenced subgap excitations.An impurity with Coulomb repulsion U , carrying a spin degree of freedom when occupied by one electron, is coupled with hybridization Γ to a superconducting island with Coulomb repulsion Ec and energy gap ∆.A quasiparticle is plucked away from the Cooper pair condensate to form a YSR singlet bound state with the spin, with the competing doublet state destabilized by Ec.In a device, the impurity can be a quantum dot, and the QD and SI gate-induced charges ν and n0 can be tuned with gate voltages.b,c Calculated charge parabolas versus δn0, which is n0 referenced to an even integer value (with ν = 1).d,e Same as (b,c), but now sweeping ν (with δn0 = 0).Ec = 0 in (b,d) and Ec = 1.3∆ in (c,e).Γ/U = 0.05 in (b-e).The parabolas are tagged by their total charge nGS referenced to an even integer value.A green dot indicates the destabilization of the doublet state by Ec for n0 = 1 from (b) to (c).Red (blue) arrows indicate addition (removal) excitations.For simplicity, continuum parabolas are not included.

FIG. 3 .
FIG. 3. Superconducting Coulombic subgap excitations and Coulomb-aided Yu-Shiba-Rusinov singlet.a Left.Zero-bias G versus VS and VN at B = 0 measured in the QD-SI device.Other gates are set at V1 = −350 mV, V3 = −52 mV, V5 = −169 mV.An unwanted gate glitch is indicated by an asterisk.The Coulomb-aided YSR singlet domain is encircled.Right.Calculation of GS transitions (blue lines) versus charges induced in the QD, ν, and in the SI, n0, overlaid on a duplicate of the experimental data for N = 800 and parameters indicated in Table I.The graph is a collage of five identical plots with n0 ranging from 799.5 to 801.5.Inset.Zero-bias G versus VS and VN at B = 0.3 T. The same colorscale and gate settings as the B = 0 diagram are used.A calculation of GS transitions (blue lines) versus ν and n0 is overlaid on the experimental data for B = 0.3 T, N = 200 and parameters indicated in Tables I and II.b Left.Colormap of G versus VN and V sd , with VN swept along the blue arrow in (a).The color scale is saturated to highlight SCE.The overlaid grey and black traces, set to the same G scale and shifted for clarity, are taken at the VN positions indicated by arrows of the same color.Right.Calculated low-energy SCE spectrum (addition: red, removal: blue dots) for n0 = 800 overlaid on the data.The continuum edge is indicated by dashed lines.Approximate QD, SI charges are given in red, while their (a) GS and (b) excitation spin Sz is given in blue.

FIG. 4 .
FIG. 4. From a many-body Yu-Shiba-Rusinov state to a two-body singlet.a,b Parameter extraction from charge parabolas for (a) an empty QD (ν = 0) and (b) a half-filled QD (ν = 1).For simplicity, the Γ/U 1 case is illustrated.Dashed (solid) parabolas indicate doublet (singlet) states.a For x > 0, scaling the horizontal solid bar by the dashed one provides x = 1 − ∆/Ec.The short vertical bar corresponds to Ec − ∆, while the long one equates to Ec + ∆, from which x can be determined when x < 0. b Scaling the horizontal solid bar by the dashed one provides y, a measurement of the YSR singlet ground state stability in units of e. Γ increases the solid bar size (red arrows).c y versus x for different ΓR/U values, measured in a SI-QD-SI device (top inset, scale bar=100 nm.) with full tunability of the Coulomb repulsion in the superconductors.The left SI (of even occupation) is used as a cotunnelling probe for the QD-right SI by setting ΓL/U = 0.02, EcL/∆L = 1.5 (bottom curve), ΓL/U = 0.3, EcL/∆L = 1.6 (middle curve), or as a soft-gap superconducting probe by setting ΓL/U ≈ 0, EcL/∆L ≈ 0 (top curve).The YSR singlet is stabilized by x for weak ΓR/U , but it is hindered for large ΓR/U .When x → 1, y converges to 1e (dashed line) independently of ΓR/U , at which point exactly one quasiparticle is bound to the spin of the QD, as sketched in the inset.y(x = 0) provides the YSR binding energy in units of 1e.The bottom right inset shows the zero-bias conductance stability diagram for ΓR/U = 0.01 at the lowest x, to exemplify parameter extraction for x > 0. x and y are measured from the thin and thick bars in the diagram, drawn as horizontal bars in (a, b).For x < 0, x is extracted from bias spectroscopy by measuring the quantities indicated by vertical bars in (a).Vertical and horizontal error bars correspond to the sum of full widths at half maximum (FWHM) of the conductance lines delimiting measurement bars for y and x.Bias spectroscopy data and full stability diagrams are shown in Extended Data Fig. 5.

FIG. 5 .
FIG. 5. From Yu-Shiba-Rusinov to superconducting Coulombic excitations.a-f Left.Differential conductance, G, in the SI-QD-SI device versus V sd and VN.Right.Zero-bias G, versus VSR and VN.The measured EcR/∆R ratio is indicated for each pair of panels.For each bias spectrum (left panels), the VSR and VN are swept through the black line on the corresponding stability diagram (right panels), but only VN is indicated.For ease of comparison, the spectra and diagrams span equal bias and gate relative ranges, and share G scales.Other device parameters are: U ≈ 0.6 − 1.2 meV, ΓR/U ≈ 0.2, ΓL ≈ 0, EcL = 0, ∆L = 0.17 meV.With the left SI tuned to be a soft-gap superconducting probe, black dotted lines highlight the negative-bias side of excitations probed by coherence peaks at ∆L, while grey dotted lines correspond to replica probed by residual density of states at zero energy.The excitations progressively transform from the characteristic YSR (split) loop in (a) at EcR/∆R ≈ 0 to the double-S shape of SCE at EcR > ∆R in (b-f ).The change is concomitant to the apparition of YSR singlet domains for odd QD occupation in the corresponding stability diagrams.Other gate settings are: VSL = −50 mV, V3 = −732.5mV, V5 = −300 mV, V bg = −1.49V. Right SI parameters are: (a) EcR ≈ 0, ∆R = 0.15 meV, (b) EcR = 0.05 meV, ∆R = 0.14 meV, (c) EcR = 0.12 meV, ∆R = 0.16 meV, (d) EcR = 0.255 meV, ∆R = 0.215 meV, (e) EcR = 0.255 meV, ∆R = 0.175 meV, (f ) EcR = 0.24 meV, ∆R = 0.14 meV.In the spectra, the true zero bias (indicated by a horizontal dashed line) is offset by V sd = −18 µV, and the stability diagrams are measured at V sd = −18 µV.

1 Extended Data Figure 1 .
Charge redistribution in the quantum dot and in the superconducting island.Calculated occupation expectation values versus gate-induced charge in the QD, ν, for (a) Ec = 0 and V = 0 (other parameters are the same as those in TableI) and (b) device parameters in TableI, and corresponding QD and SI filling schematics for specific ν values.The calculation is done for N = 800 and n0 = 800.The occupation expectation values are those of the QD, nN , and of the SI, δ nS , where δ indicates that the number shown is the variation in occupation over 800 electrons.The values shown correspond to the GS and the first addition and removal excitations.The addition and removal expectation values overlap in (a), but do not do so in (b), leading in the latter case to the double-S shape of the SCE.In the QD and SI charge filling schematics, black (white) arrows represent electron (hole) spins.These schematics show that, whereas at ν = 1 both (a) and (b) are Yu-Shiba-Rusinov (YSR)-like with electron-hole symmetric excitation expectation values (with double-headed arrows in (b) reflecting the fractional expectation values for the QD and SI occupations), away from ν = 1 (e.g., right after the charge fluctuation point at ν = 1.55) there is a strong asymmetry in addition and removal occupation values in (b) not present in (a).

Extended Data Figure 2 .
Charge parabolas at finite B and exchange interaction.(a,b) Calculated charge parabolas versus δn0 at ν = 1 for Zeeman energies of the QD and the SI indicated in each panel, where δn0 corresponds to the gateinduced charge in the superconductor, n0, shifted by N = 200.Each parabola is tagged by its Sz number.A red bar in (a) indicates the singlet-triplet exchange splitting.Black bars in (b) indicate doublet and triplet Zeeman splittings.The Sz = 0 triplet state is not included in the calculation.The sketches on the right show spin states color-coded as the states shown in (b) for δn0 = 1.

Extended Data Figure 3 .Extended Data Figure 4 .
Changes in the VS dependence of superconducting Coulombic states across a nN = 2 → 1 transition.(a-e) Colormaps of G versus source-drain bias, V sd , and SI gate voltage, VS, at various settings of the QD gate voltage, VN, indicated by (a) red, (b) green, (c) purple, (d) yellow and (e) cyan vertical arrows in the stability diagram on the left, which is a duplicate of Fig.3ain the main text.The tails of the black arrows are attached to a subgap state, while their heads point to the direction of the evolution of said state with varying VN.The color scale is saturated to highlight subgap excitations.Charge numbers of the QD and the SI, nN and nS, are indicated in (a) and (e).An unwanted gate glitch is indicated by an asterisk in (b).Calculated SCE spectra using the same model parameters as in Fig.3ain the main text are overlaid as blue dots on each panel for comparison.In the calculation, ν is fixed to the values indicated on top of each plot, and n0 is swept.The calculation matches the position of positive-slope SCE, but does not account for negative-slope features, negative G features, and continuum features.Spectral comparisons with Coulombic and usual Yu-Shiba-Rusinov (YSR) limits.(a-c) Calculated diagram of GS transitions (collage of two identical plots) versus the gate-induced charge in the QD, ν, and the gate-induced charge in the SI, n0, for (a) our QD-SI device parameters, (b) Coulombic limit (∆ = 0), and (c) YSR limit (Ec = 0).Bars indicate the vertical size of charge domains in electron charge units.Whereas in (a) ∆ modulates this size depending on ν and favors 2e charging away from ν ≈ 1, in (b) the size is strictly 1e independently of ν and a honeycomb pattern is obtained.Arrows in (a) are anchored to GS transition lines which split when the ratio Ec/∆ is increased.In (c), absence of charging results in no GS dependence on n0.(d-f ) Calculated first-excitation spectrum versus ν for n0 = 800 for each of the cases shown in (a-c).Dashed lines of different color highlight the position of the different energy gaps in each case: (d) Ec + ∆, (e) Ec, and (f ) ∆. Aside from this difference, since Ec is sizeable in our device, the spectrum in (d) resembles qualitatively the Coulombic limit in (e).However, due to the presence of an equally sizeable ∆ in our device, the differences are very marked near n0 = 801, as shown for the GS in (a, b).The energy-asymmetric and discontinuous spectra in (d, e) are to be contrasted with the usual YSR limit in (f ).Parameters for (a, d) are noted in Tab.I.For (b, e): U/Ec = 4, Γ/U = 0.05, V = 0.For (c, f ): U/∆ = 4, Γ/U = 0.05, V = 0. N = 800 in all cases.

5 .
Tuning the Coulomb repulsion in a two-island device.a Scanning electron micrograph of the SI-QD-SI device, comprising an InAs nanowire with Al SIs (below top gates SL and SR).Top gates 3 and 5 confine a QD below top gate N. Top gates 1 and 7 are shorted to top gates SL and SR, respectively.Scale bar is 100 nm.b Electrostatics of the device.c, d, e Zero-bias G versus gate voltage of the right SI, VSR, and gate voltage of the QD, VN, from which the data plotted in Fig. 4 (indicated here by vertical bars) was extracted.Other gates are set to (c) VSL = −761 mV, V3 = −505 mV, V5 = −330 mV, V bg = −0.05V, (d) VSL = −536 mV, V3 = −443 mV, V5 = −170 mV, V bg = −2.69V, and (e) VSL = −260 mV, V3 = −830 mV, V5 = 150 mV, V bg = −0.48V, with V5 controlling ΓR by accumulating electrons in the nanowire.VSL corresponds to nSL =even, where the left SI does not play a role in defining the GS.Approximate charges in the parts of the device (nSL,nN,nSR) are indicated.Tuned by VSR, EcR − ∆R increases progressively from bottom to top in (c) and from top to bottom in (d, e), expanding nSR = 1 sectors.Charge sectors are identified by the approximately integer occupations in the three Coulomb-blockaded parts of the device (nSL,nN,nSR) indicated in red.The singlet assignment in (e) is verified by the B dependence of the stability diagram.f G versus VSR and source-drain bias voltage, V sd , with VSR swept along the dashed line in (d), where nQD = 2.The true zero bias (indicated by a dashed line) is offset by V sd = −18 µV.When EcR − ∆R < 0, the vertices of the Coulomb diamonds do not touch zero bias, while for EcR − ∆R > 0 small diamonds emerge.Bars indicate how EcR − ∆R relates to these features.Asterisks denote gate glitches.

TABLE I .
Parameters of the QD-SI device and model.Top-row parameters are estimates obtained from measurements (for methods, see Section I of Supplementary Information), while the bottom-row ones represent best fits of the model output to the experimental data based on the measured parameters as the initial input for subsequent fine tuning.