Multiterminal transport spectroscopy of subgap states in Coulomb-blockaded superconductors

Subgap states are responsible for the low-bias transport features of hybrid superconducting--semiconducting devices. Here, we analyze the local and nonlocal differential conductance of Coulomb-blockaded multiterminal superconducting islands that host subgap states with different spatial structures. The emerging patterns of their transport spectroscopy are used to characterize the possible topological nature of these devices and offer the possibility of controlling their transport properties. We develop a next-to-leading order master equation to describe the multiterminal transport in superconductors with both strong Coulomb interactions and multiple subgap states, coupled with metallic leads. We show that the nonlocal differential conductance characterizes the spatial extension of the subgap states and signals the presence of degenerate bound states with a finite support on different parts of the device. Additionally, it displays sharp sign changes as a function of the induced charge of the superconductor, signaling energy crossings among its lowest excited states.


I. INTRODUCTION
Hybrid systems fabricated with superconducting and semiconducting materials are considered a key ingredient for the development of quantum technologies and have demonstrated a potential for quantum information storage and processing [1][2][3][4].In particular, the advent of topological superconductors has raised the hope of realizing protected quantum devices, robust to decoherence induced by local perturbations [5].
In one dimension, topological superconductors host Majorana subgap modes [6], predicted to be nonabelian anyons and to give rise to several nonlocal properties [7][8][9][10].These Majorana modes are exponentially localized at the system edges but, importantly, each pair of them defines a nonlocal quasiparticle state, with vanishing energy in the ideal case.In this context, distinguishing between local and nonlocal subgap states becomes relevant for identifying devices that are potentially in the topological regime.
When the superconductor is not grounded but exhibits a strong charging energy, the intricate interplay between superconductivity and many-body interactions yields a rich phenomenology.In the topological regime, Coulomb blockade effects provide effective tools to initialize, read out, and measure coherence times in topological qubits, thus constituting a key element to design platforms for quantum information devices [43][44][45][46][47][48][49][50][51][52][53].
Transport properties of floating superconducting islands have been extensively studied experimentally in semiconducting nanowires, with the goal of extracting quasiparticle poisoning times [54,55], showing the doubling of Coulomb peaks in the high-field regime [56][57][58][59], and investigating the exponential protection [56] and spin polarization [60] of subgap states.These systems have been explored using interferometry [61] and reflectometry techniques [62,63].However, most of the studies so far, both experimental and theoretical, focused on local conductance measurements of two-terminal devices; nonlocal transport properties of multiterminal Coulombblockaded devices remain largely unexplored.
A more complex typology of devices is based on two nanowires coupled via a floating superconductor, sharing a common charging energy.In the case where the device couples to more than two leads, the topological Kondo effect is predicted to appear, and its realization would be an unequivocal signature of the Majorana nonlocality [64][65][66].As such, systems with Coulomb-blockaded double nanowires are at the focus of recent experiments [67,68] and theoretical investigations [69].
In this work, we propose multiterminal quantum transport measurements as a way to distinguish between local and nonlocal subgap states in Coulomb-blockaded superconducting devices.We focus on two geometries, namely a single and a double superconducting nanowire device.We employ a second-order master equation approach that allows for studying both sequential and cotunneling transport in superconducting devices with strong Coulomb interaction, multiple subgap states, and the coupling to two or more metallic leads.In particular, we show that sequential transport at low bias voltages is blocked whenever the lowest-energy subgap state does not couple two of the leads.Indeed, even though the Coulomb interaction is highly nonlocal, transport is dominated by the local projection of subgap states, providing indications about their spatial structure.Moreover, a three-terminal Coulomb-blockaded device displays abrupt changes in the current direction at fixed bias potential as a function of the charge n g induced on the superconducting island, due to energy crossings in the excitation spectrum, similarly to nonsuperconducting Coulomb-blockaded quantum dots [70].Concerning the double nanowire geometry, we investigate the local and nonlocal conductance associated to different spatial structures of the subgap states, aiming to facilitate the interpretation of future experiments.
The rest of the paper is organized as follows.In Sec.II we present the effective zero bandwidth model we use to describe the low-energy behavior of the devices and the formalism used to compute the sequential and cotunneling current signals.In Sec.III we investigate both the two-terminal and three-terminal conductance of a single nanowire with one or multiple subgap states.In Sec.IV we analyze the transport signatures associated to different structures of the subgap states in the double nanowire geometry.We draw our conclusions and discuss future developments in Sec.V.The appendices contain more supporting results and a detailed derivation of the tunneling rates calculation.

II. MODEL AND METHODS
We consider a superconducting island coupled to normal metallic electrodes, given by the Hamiltonian where the leads are described by with the electron energy, ξ νkσ , and the annihilation operator, c νkσ , referring to an electron in lead ν with momentum k and spin σ ∈ {↑, ↓}.We assume that each lead remains in internal equilibrium, described by a Fermi-Dirac distribution n F with chemical potential µ ν and temperature T .The superconducting island is described by where j labels a state with energy ε jσ , and the Bogoliubov-de Gennes (BdG) operators read Here, d and d † are annihilation and creation operators of electrons in the island, u j and v j are the BdG coefficients for the subgap states, φ is the superconducting phase operator, such that e −iφ annihilates a Cooper pair in the island, and ρ σ = ± for σ =↑, ↓.
The electrostatic repulsion term is given by where n g is the dimensionless gate-induced charge offset and N the (excess) electron number operator, accounting for the fermionic occupation of the states and the number of Cooper pairs on the island.Figure 1(a) sketches the device described by this model, in the case of only two leads and a single nanowire.Panels (b) and (c) represent the superconducting quasiparticle spectrum and the many-body energy levels, respectively.Finally, the tunneling between leads and superconducting island is described by the Hamiltonian where t is the tunneling amplitude.The corresponding tunneling rates are defined as Γ ν,j,σ = 2πρ F ν |t ν,j,σ | 2 , and we consider them k-independent (wideband limit), where ρ F ν is the density of states at the Fermi level of the lead ν.
In this work, we employ a zero-bandwidth approximation for the superconducting island, meaning that we consider only two relevant energy levels ε 0,↑ = ε and The first represents a low-energy spinpolarized subgap state whose spatial structure and degeneracy will be the main focus of our work.The state at energy ∆ is a spin-degenerate level and represents an effective description of the proximitized superconducting continuum which extends through the whole hybrid device.This approximation models a short nanowire where the quasiparticle spectrum above the proximity induced superconducting gap displays a quantization determined by the finite system size.The corresponding lowest energy excitation can mediate elastic cotunelling processes through the device.A possible alternative is to consider a number of degenerate states at energy ∆, each coupled to a single lead.Each of these states models a metallic quasiparticle continuum with a strong relaxation towards the bottom of the continuum band.With this choice, the states at energy ∆ can not mediate elastic cotunneling transport.Throughout the rest of the paper, we focus on the first scenario, where each nanowire hosts a single extended superconducting continuum, to avoid an excessive increase of the computational run time.However, this choice does not qualitatively influence the transport features at low bias, where the conductance is dominated by the subgap states.See Fig. 6 in App.B for a comparison of the two cases.

A. Formalism
In this work, we focus on the regime where Γ ν,j,σ T .In this weak coupling limit, a perturbation theory in the tunneling amplitude provides accurate results.Since we are not interested in the thermal excitations of the superconducting device, we additionally consider temperatures much smaller than the remaining energy scales.To compute the tunneling rates, we use the T -matrix formalism, describing the transition probability between two states where W if weights the rate through thermal distributions of the electrons in the leads, δ is the Dirac delta functional ensuring energy conservation between the initial and final states (E i = E f ) and which can be truncated at the desired order.In Eq. (7) and in the rest of the paper, we set = 1.The linear term in Eq. ( 8) describes the sequential tunneling; the higher order terms define, among other processes, the cotunneling contributions, which become progressively more important when the tunnel amplitudes t increase.In Sec.
A of the appendix, we provide expressions for the sequential and the (second-order) cotunneling rates, used in this work to evaluate the transport through the superconducting islands.
The quantum state of the island is described by |a =|N, N C , n , where N C is the number of Cooper pairs in the island and n is a vector representing the occupation of the subgap and continuum excited states.The time derivative of the occupation probability of a given state is given by the master equation Here, we assume that coherences between different states in the superconducting island are negligible; indeed, we consider either models without near degeneracies where coherences cannot develop or with no coupling between degenerate states, due to different parity or spatial support.In the stationary limit, the system does not evolve in time and Ṗ stat b = 0.These conditions, together with the normalization b P stat b = 1, form a linear system of equations for the occupation probabilities of the island states.Using the resulting stationary distribution, the current flowing from lead ν to the device is determined by the island transition rate to other states as where we take s = +1 (s = −1) for electrons tunneling in (out) of the island and we set the electron charge to unity.Here, δn(ν) is the net charge transferred from lead ν in a cotunneling process, and the sum runs over all the possible rates connecting the island state |b with any state |a (see Appendix A 2 for more detail on the cotunneling rates Γ cot ).In Eqs. ( 9) and (10) we are considering processes that change the island occupation by ±1 electron, while neglecting local and crossed Andreev reflection processes [71].This approximation is justified when the charging energy is the dominant energy scale, in particular E c ∆, meaning that there is a large energy penalty for changing the charge of the device by 2e.Throughout the rest of the paper, we will indeed consider this strongly Coulomb-blockaded regime.

III. SINGLE NANOWIRE GEOMETRY
In this section, we focus on resolving the spatial structure of subgap states in a single nanowire device coupled to two or three normal metallic leads.Regarding the spatial profile of the subgap state, we consider three situations: (i) a single subgap state delocalized at both ends of the nanowire but vanishing in its bulk, which we refer to as a "Majorana-like scenario", (ii) two degenerate subgap states each localized close to one end of the device, (iii) an extended subgap state that couples to all leads independently from their position.In all these cases, we consider that the subgap states are strongly spin-split, behaving as effectively spinless.This situation can be achieved by a strong magnetic field or by proximity with ferromagnetic materials, and corresponds to standard models to achieve the topological regime [72].Our goal is to distinguish between these possible scenarios based on the structure of the cotunneling steps and the nonlocal conductance.

A. Two-terminal device
We first focus on a two-terminal device, with leads coupled to the nanowire ends, which can host either a single subgap state or two degenerate ones located at the ends, see Figs. 2 (a) and (b).For simplicity, we choose a symmetric voltage drop at the contacts between the device and L and R leads by setting their chemical to µ L = +V /2 and µ R = −V /2, where V is the voltage bias.
If the subgap levels have zero energy ε = 0, which is the case for ideal Majorana zero energy modes, the Coulomb blockade structure is qualitatively similar in the two cases (data not shown).Instead, when ε = 0, the most evident difference is the strong reduction of the conductance peaks at V = 0 in the case of local states.This can be seen by comparing the conductance at V ∼ 0 and halfinteger n g in panels (c) and (d) of Fig. 2. Indeed, when the subgap states do not couple to both leads, transport is suppressed at small biases because the subgap states cannot support resonant tunneling of electrons.This results in a much weaker conductance for |V | < 2ε, ∆ − ε with respect to the situation where the subgap state connects both ends of the nanowire.However, zero-bias tunneling is still possible through elastic cotunneling, with an amplitude ∼ Γ L Γ R [71].
We also observe a clear difference in the structure of the cotunneling steps, highlighted by the colored tics in Figs. 2 (c) and (d).In the case of an extended state, two inelastic cotunneling steps appear at V = ∆ ± ε for the even and the odd valleys, corresponding to transitions to the lowest excited states (red tics) [60].In the even valley, it corresponds to the splitting of a Cooper pair, whose electrons end up in the subgap and the continuum of states.In the odd valley, the lowest excitation corresponds to promoting the electron in the subgap state to the continuum.There are two additional steps due to the breaking of a Cooper pair into two electrons that enter the quasiparticle continuum (blue tics).
In the case with two states with local support, each coupled to a single lead, the number of cotunneling steps increases and, in particular, the threshold for inelastic cotunneling in the even-parity valleys decreases substantially.Indeed, the lowest excited state with even parity has an energy given by the sum of the energies of the two subgap states [see panel (j) in Fig. 2].In the case considered in the figure with the two states having the same energy, this threshold is 2ε.We note that this process does not require the creation of a quasiparticle excitation above the superconducting gap, reducing the conductance step onset with respect to the nondegenerate situation.While the local conductance (through L and R) shows little dependence on the presence of a third lead and how it couples with the device (data not shown), the nonlocal conductance displays interesting features.First, notice that G M is an odd function of the voltage bias V , thanks to the symmetry R ↔ L we impose in our model: an inversion of the bias only changes the direction of the current flow between the L and R lead, while the cur-rent I M remains insensitive to the sign of V since the lead is always grounded.Hence, the differential conductance (G M = dI M dV ) is an odd function of the voltage bias.When the lead does not couple to the subgap state, the conductance at V = 0 is suppressed, as clearly seen by comparing left and right panels in Fig. 3.Even though this effect is somehow trivial because low bias transport only involves the lowest energy level in a gapped system, its consequence is of great importance: the middle lead probes the density of states inside the wire.Hence, the M lead can resolve the spatial profile of subgap states and detect whether they have a nonvanishing projection on a specific portion of the device.
The non-local conductance also shows a peculiar sign dependence on n g , illustrated by the sharp G M jumps in Fig. 3.These sign changes appear whenever there is a crossing between ground states or excited states with different parities, corresponding to the crossing of the parabolas in Fig. 1 (c) (see Appendix B for more details).Their dependence on V , instead, arises from the necessity of having a chemical potential on the L and R leads large enough to populate the excited states coupled also with the M lead.This behavior can be better understood by considering the energy dependence of the superconducting island states on n g and the sequential tunneling rates (see App.A 1 for their derivation).For instance, the rates describing the tunneling of an electron from the M lead to an empty quasiparticle state , where i and f label the initial and final states with energy E i and E f .Γ M |u j | 2 is the effective tunneling rate in the wideband limit multiplied by the local projector on the particlelike component of γ j .Since the temperature is small, the Fermi factor n F is a step-like function, activating the tunneling process mainly when E f − E i < 0. The opposite is true when we examine the process for extracting one electron from the occupied quasiparticle state γ j , meaning that it is activated when E f −E i > 0. Hence, the sign of the current in the M lead changes abruptly when n g tunes the energy levels of the device across the resonance E i = E f , as adding/removing electrons from M becomes more favorable.We show more details in Fig. 5.
In conclusion, a three terminal Coulomb-blockaded device has a highly tunable differential conductance alongside the standard transport suppression inside the Coulomb diamonds.As shown in panels (c)-(f) of Fig. 3, this is a robust feature that requires neither the presence of a zero energy subgap state nor a third lead coupled only to the superconducting continuum, even though the latter is useful to have an extended bias window where the M lead is effectively decoupled from the device.

IV. DOUBLE NANOWIRE
In this section, we extend the analysis of the nonlocal conductance to double nanowire setups; two parallel semiconductor nanowires are covered by the same ○ an extended state coupled to all the leads, and 4 ○ a local state that couples to R, while L and R couple to a different extended state.The parameters are the same as in Fig. 2.
floating superconducting island.The common superconductor allows the exchange of Cooper pairs between the wires, making them share a common charging energy, but it does not allow the tunneling of excitations above the gap from one wire to the other.In this context, we therefore introduce two separate states at energy ∆ to model the quasiparticle continua in the two wires.In contrast, we assume that subgap states can delocalize between the two wires.This setup for the device is inspired by recent experimental achievements [68,[73][74][75] and the scope of exploring the topological Kondo effect when the system couples to more than two leads [64][65][66].We consider that the nanowire ends couple to three different leads, as sketched in panels (a), (d), (g), and (j) of Fig. 4. As before, we consider that the L/R leads are symmetrically biased, while M is grounded.A different biasing condition is shown in Fig. 7 in the appendix.
Regarding the spatial structure of subagap states, we identify four possible situations: 1 ○ two subgap states that extend along a single nanowire, 2 ○ a subgap state localized at each lead-device interface, 3 ○ a common subgap state that couples with all three terminals, 4 ○ two subgap states localized respectively at the left and right ends of the double nanowire device.These scenarios are sketched in panels (a), (d), (g), and (j) of Fig. 4, while the panels below each sketch show the related local G L and nonlocal G M differential conductances.In all cases, we consider degenerate subgap states with energy ε = 0.3∆.
Subgap states with different energies would be easily detected by the different sizes of the Coulomb diamonds.The differences between the four cases are summarized in Table I.
Let us start from case 1 ○, where each nanowire hosts a single subgap state with support on both ends.In this case, there is no probe at the middle of any of the wires.Therefore, the transport features cannot distinguish between trivial extended states or nonlocal Majorana-like subgap states.The associated local conductance, G L = dI L dV , and the non local one, G M = dI M dV , are reported in Figs. 4 (b) and (c), respectively.Sequential tunneling processes contribute to the current between L and R at low bias, while they are suppressed for the M lead.A single sequential tunneling line is observed in G L at low bias voltages; in the odd valleys, the stationary distribution P stat has indeed a significant contribution from the configuration where a quasiparticle is frozen in the subgap state in the lower wire.This is due to the imbalance between the Γ M rates for incoming and outgoing particles at small temperatures.When the lower wire subgap state is occupied, the superconducting island cannot be excited by removing this quasiparticle via sequential tunneling through leads L or R, thus suppressing the leadingorder transport mechanism in G L .This is evident when comparing the odd valleys in Fig. 2(c) and Fig. 4(b): in proximity of the charge degeneracy points the sequential tunneling lines moving towards the odd valleys disappear and cotunneling becomes the dominating process at low bias.Also the weak conductance through the M lead at low bias voltage is due to cotunneling processes, where the device exchanges one electron with either L or R, and M , keeping the total charge on the superconducting island constant.
In case 2 ○, when a trivial subgap state localizes at the interface with each lead, the conductance is suppressed for |V | < ε.This is due to the local support of the subgap states, which cannot directly mediate sequential transport between different leads.Therefore, transport is dominated by elastic cotunneling at low bias voltages.The nonlocal conductance G M , Fig. 4(f), displays feint resonances in the even valleys.The reason behind is that tunneling processes between the superconducting island and the M lead involve (virtual) changes of the number of Cooper pairs in the island.The sign of the current depends on the charge difference between the ground and the lowest excited state.
Case 3 ○ corresponds to a single subgap state coupled with all three leads.Therefore, the same Coulomb structure is present in the local and the nonlocal conductance, although the latter is an odd function of the bias V and exhibits sign changes when there are level crossings in the many-body spectrum of the superconducting device.The sequential tunneling of electrons dominates transport, leading to strong conductance features at the charge degeneracy points for small bias voltages.Qualitatively, the signal of the M lead is equivalent to that shown in Fig. 3(f) for a single nanowire with an extended bound state.
The last possibility we analyse, case 4 ○, corresponds to two degenerate subgap states, localized one in the left part and one in the right part of the device.This situation is of particular interest, as it might correspond to a double nanowire geometry where Majorana modes localized at each device ends strongly hybridize due to their small spatial separation.The local conductance G L displays a parity dependent negative differential conductance (NDC) region at small bias, see Fig. 4(k).This is caused by a quasiparticle being trapped in the subgap state coupled to the R lead, blocking the current between the L and M leads.The appearance of the NDC region for either positive or negative bias is due to the electron-or hole-like excitation trapped in the right subgap state.In this regime, sequential transport through the right lead is suppressed, leading to a current reduction (not shown).For |V | > 2ε, the sequential transport channel to the right lead is open again, therefore making It also provides a relatively fast relaxation mechanism for the trapped quasiparticles in the right lead and the NDC region disappears.
The four cases show conductance steps inside the Coulomb valleys at finite voltage values.These steps appear when V matches the system excitation energy, opening the inelastic cotunneling channel: exchange of two electrons between the leads and the island, leaving it in an excited state (although keeping the island charge invariant).Depending on whether a bound state couples simultaneously to the L and R leads, cases 1 ○ and 3 ○, or not, cases 2 ○ and 4 ○, the number of cotunneling steps vary from 2 to 3 (colored tics in Fig. 4).A similar behavior has been described for the single wire situation, see discussion around Fig. 2.
Finally, we show in Appendix B another biasing situation, where L and M are symmetrically biased and R is grounded.We find that local and nonlocal transport can distinguish between the considered four situations also in that case.

V. CONCLUSIONS
In this paper, we analysed how the transport features of a multiterminal superconducting device with strong charging energy depend on the number and spatial structure of subgap bound states.In particular, we investigated the role of the spatial extent of the subgap states, which might not couple to all the leads attached to the device; we showed that the nonlocal differential conductance in multiterminal devices allows for a qualitative characterization of their spatial profiles for biases below the superconducting gap.We considered a zero bandwidth model to describe the low-energy features of a superconducting floating island consisting of a proximitized single semiconducting nanowire or a pair of parallel nanowires.These systems are indeed known to host subgap states that determine the transport properties of the device.We focused on a situation where the charging energy E C is the dominant energy scale, and we adopted a second-order master equations approach which allows us to characterize both the sequential tunneling of single electrons and the inelastic cotunneling events.The sequential tunneling signal, dominating for small leads-device tunneling strength, gives information on the energy and spatial structure of the subgap states mediating transport.When two leads are not coupled by the lowest-energy states, the zero bias conductance is strongly suppressed, leaving only a faint cotunneling feature.In this way, transport can determine whether a state has support on the two ends of the wire.In the same way, additional leads can be added to gain spatial resolution inside the wire.Therefore, the absence of sequential tunneling conductance at low bias is a way to discriminate between trivial extended states and possible topological states with only support at the ends of the wires.Moreover, the number and voltage of inelastic cotunneling steps allow us to determine the number and the energies of the subgap states.The cotunneling signal is therefore another indicator that can be used to characterize the spatial structure of the subgap states.The described features hold for an arbitrary number of discrete subgap states with energy |ε| > T, Γ.For ε < T, Γ, instead, different subgap state configurations may result in the same qualitative transport features, thus hindering their spatial characterization.
Finally, electron transport in multiterminal Coulombblockaded devices results in an interesting pattern of peaks with positive and negative differential nonlocal conductance, depending on the induced charge n g and on the voltage bias V .This allows to switch the direction of the current flowing in the grounded lead, or suppress it, without changing the potential difference between the source and the drain but only by tuning the induced charge on the whole superconducting island.
We presented results based on a perturbative analysis of the transport properties, valid when the temperature is larger than the coupling between the leads and the device.Complementary methods are needed to describe the low-temperature and strong-coupling regimes, where electron correlations effects are important, in the non-equilibrium situation [76][77][78].The master equation approach we presented, however, is less computationally intensive and provides a clear picture of the transport mechanisms as long as non-perturbative effects can be neglected. .2πρ F |t ν | 2 , with the lead density of states at the Fermi level ρ F .These rates induce transitions between charge states differing by one electron, denoted through the vectors n and n accounting for the fermionic occupation of the island states.

Cotunneling rates
Inside the Coulomb-blockaded region, inelastic cotunneling creates a series of steps in the differential conductance, where the voltage bias corresponds to energy differences between different states of the superconducting island.These are processes where one electron is transferred between two leads, leaving the island in an excited state.We also consider the elastic cotunneling that, instead, describes processes in which the island energy is conserved, and gives rise to the conductance background inside the diamond.However, we disregard local and crossed Andreev processes, where the island charge changes by 2e, as they are suppressed by the strong charging energy in the system.The corresponding rates are given by Eq. ( 7), where Here, H T (ν) describes the tunneling between the lead ν and the island, ν denotes a lead different from ν, and m 1,2 are virtual intermediate states.To derive this expression, we have imposed energy conservation, which leads to a function dependent on the energy of the tunneling electron from/to one of the leads, ω 1 .The cotunneling rate can be written as This expression for the cotunneling rates is divergent when the denominator in Eq. (A2) vanishes.To avoid the divergent behaviour, we regularize the divergences as explained in Ref. [79].The resulting integral can be formally solved analytically [80], which leads to a complicated expression involving special functions.We find that for T /E C 10 −3 , it is computationally more efficient to expand the Fermi distribution function into a sum of complex Matsubara-Ozaki frequencies [81] n where β α and r α are the approximated Matsubara frequencies and residues, respectively.Finally, Eq. (A3) can be evaluated using the residue theorem, yielding to a rather compact expression This sum can be truncated at α ≈ 100 Matsubara-Ozaki frequencies for the parameters used in the calculations.

Appendix B: Supporting results
In this appendix, we present additional results to support the chosen modelling of the considered devices and describe some details of the transport features presented in the main text.In particular we first focus on the processes describing the transport of the three-terminal single-nanowire devices, which determine the transport across the lead M ; then we address the comparison between different modelling of the quasiparticle states above the superconducting gap; finally, we consider a different choice of the voltage drop across the doublenanowire devices.To understand the three-terminal transport features in a Coulomb-blockaded superconducting device, it is useful to compare the conductance dependence on the voltage bias and on the induced charge with the low-energy many-body spectrum.In Fig. 5, we report the eigenvalue structure and the nonlocal conductance G M = dI M dV , close to the charge degeneracy point, for the device sketched in Fig. 3(a): a single Coulomb-blockaded nanowire hosting a subgap state that does not couple with the grounded lead M .The lowest energies of states with odd and even parity are represented by the continuous blue and orange lines, respectively, in Fig. 5(a)-(b).These states are connected via sequential tunneling involving a particle transfer to or from the subgap state, which does not contribute to the current flowing through the M lead.Hence G M is suppressed at low bias, as can be seen in Figs.5(c) and (d).The current in the M lead is activated only when V is large enough to populate excited states (via the L or the R leads) that can then relax to a lower energy state through a tunneling event between the superconducting continuum and the lead M .These processes correspond to transitions between a dashed line (i.e., states with excited quasiparticles in the superconducting continuum) and a continuous line with different parity (colors) in Figs. 5 (a) and (b).These events are highlighted by the oblique lines in panels (c) and (d), representing energy thresholds for populating states that can contribute to the current through M .The sign of The lines in the lower panels correspond to the excitation threshold denoted by the arrows in the upper ones.In panels (a) and (b), sequential tunneling via the M lead connects only states with different parity and charge occupation in the states at ∆.These processes always involve a relaxation from a higher energy state to a lower energy one, thus explaining the changes in sign of the differential conductance when two parabolas cross.The remaining parameters are the same as in Fig. 2 the current, and hence of the conductance, changes at the crossing between energy levels with different parity.This phenomenology is independent from the energy ε of the subgap states, whether it is zero, as in panels (a) and (c), or finite, as in panels (b) and (d).The latter situation only has a slightly richer structure of the Coulomb diamonds due to the energy difference between the first excited states in the even and odd sectors.
Another interesting effect shown in Fig. 3 is the sign mismatch between the sequential tunneling and the cotunneling contribution close to the lower edge of the Coulomb diamonds.This can be seen at V = ∆, where cotunneling (blurred signal) gives rise to negative conductance while sequential tunneling (sharp lines) contributes with positive conductance when n g approaches 5.5 from below.Both can be understood by considering the processes mediated by tunneling through the M lead for specific values of n g and V .When V > 0 and n g → 5.5, the lowest energy excitation that can relax through the cen-tral lead is the odd parity state with a quasiparticle in the superconducting continuum that couples with an incoming electron to create a Cooper pair.This process is thus associated with a positive (ingoing) current from lead M .Instead, the lowest energy inelastic cotunneling step in the odd valley corresponds to the excitation of a quasiparticle from the subgap state to the continuum, mediated by the virtual occupation of the even parity state with no quasiparticle present.In this second-order process, the only tunneling event through the M lead is the destruction of a Cooper pair into the high-energy quasiparticle and an outgoing electron, which carries a negative particle current and, thus, negative conductance.Similar arguments explain the sign change between sequential tunneling and cotunneling in other regions of the Coulomb diamonds.

One-vs two-state approximation to model the Bogoliubov quasiparticle continuum
Concerning the modelling of the quasiparticle continuum above the gap, in Fig. 6 we compare the twoterminal conductance for a nanowire with a single quasiparticle state at energy ∆ coupled with both leads and with two degenerate quasiparticle states at the same energy with finite supports on each end, as sketched in the upper panels.The nanowire hosts one or two subgap states at energy ε.As long as the voltage bias is smaller than ∆ − ε, transport is dominated by the properties of the subgap state and the spatial structure of the superconducting continuum does not affect the qualitative features in conductance of the device.When the bias is larger and the states at energy ∆ become accessible via sequential tunneling events, the "broken" continuum can trap a quasiparticle in one of the two edges of the device, suppressing transport and inducing a negative differential conductance region in the Coulomb diamonds.Notice, however, how this NDC is qualitatively different from that appearing in all data presented in the main text, as it is a sequential tunneling feature appearing in the local conductance of two-terminal devices.

Example of the multiterminal transport with different voltage drops
Finally, in Fig. 7 we present the conductance results for the double nanowire geometry with a different bias choice with respect to the data shown in Fig. 4 in the main text: here we bias symmetrically leads L and M , while lead R is left grounded.The main difference with respect to Fig. 4 is that now the biased leads are not connected by a single superconducting continuum at energy ∆.In particular, notice the different natures of the NDC regions appearing in panel (b) and (h): in the former, it is due to a quasiparticle trapped in the subgap state connected to lead M and, indeed, it is not reflected in G M .In the latter, it corresponds to a quasiparticle trapped in one of the two superconducting continua, in the upper or in the lower nanowire.This second case is similar to the data presented in Fig. 6(d), where the NDC is due to the broken superconducting continuum on a single nanowire, although the conductance in Fig. 7(h) is not symmetric in V because of the presence of the third terminal.○ an extended state coupled to all the leads, and 4 ○ a local state that couples to R, while L and R couple to another extended state.The parameters are the same as in Fig. 2. The conductance we show for lead M is GM = dI M dV R , hence it is always negative, with the exception of some small regions.

Figure 1 .
Figure 1.(a) Sketch of a Coulomb-blockaded hybrid superconducting nanowire device, hosting a subgap state (red).In this example, the nanowire couples to two leads (indicated with L and R).(b) Single-particle energies of the island, containing a subgap state at energy ε (black line), quasiparticle excitations above the gap ∆, and the lead states displaying Fermi-Dirac distributions (yellow) at different chemical potentials.(c) Lowest eigenergies in the presence of interactions as a function of the offset charge, ng = αVg with α being the lever arm.Blue (orange) lines correspond to states with odd (even) number of electrons in the island.We show an example where the nanowire hosts a spinless subgap state, where the solid and dashed lines correspond to the ground state and the excited state in each charge sector.

Figure 3 .
Figure 3. Nonlocal differential conductance GM = dI M dV through a grounded lead coupled to the middle of the wire.(a) and (b) respectively depict situations where the middle lead does not and does couple to the subgap state (red blobs).Panels (c) and (d) show the nonlocal conductance for a case where the subgap state energy is ε = 0, while panels (e) and (f) show cases for ε = 0.3∆.Panels (c) and (e) refer to the subgap state structure shown in (a) while panels (d) and (f) correspond to (b).The remaining parameters are the same as in Fig. 2 and all the tunnel couplings are considered to be equal.

Figure 4 .
Figure 4. Double wire system coupled via a trivial superconductor allowing the exchange of Cooper pairs, so the system has a common charging energy. 1 ○ One subgap state connects to L and R, while another couples to M (a).Panels (b) and (c) are the local conductance through the lead L and the nonlocal conductance through the grounded M lead, respectively.2 ○ Situation for three local states coupled to the different leads, 3○ an extended state coupled to all the leads, and 4 ○ a local state that couples to R, while L and R couple to a different extended state.The parameters are the same as in Fig.2.

4 ○
Sequential tunneling present in all terminalsSequential tunneling absent in GR 2 cotunneling steps 3 cotunneling steps Negative differential conductance in GL TableI.Summary of the main transport features for the four situations shown in Fig.4.

1 .
Features of the current across the lead M in the three-terminal nanowire

Figure 5 .
Figure 5. (a) Many-body eigenergies of the island for a spinless state at zero energy (∆ = 0.2EC).(b) Many-body energies for a subgap state at finite energy, ε = 0.3∆.(c) Differential conductance through a middle lead that does not couple to the subgap state [see the sketch in Fig. 3 (a)].Here we consider only the sequential tunneling contribution to the current.Panel (d) shows the corresponding results for ε = 0.3∆.The lines in the lower panels correspond to the excitation threshold denoted by the arrows in the upper ones.In panels (a) and (b), sequential tunneling via the M lead connects only states with different parity and charge occupation in the states at ∆.These processes always involve a relaxation from a higher energy state to a lower energy one, thus explaining the changes in sign of the differential conductance when two parabolas cross.The remaining parameters are the same as in Fig.2

Figure 6 .
Figure 6.Two-terminal conductance for different subgap and above gap states situations.Panels (a) and (c) show schematically the situation where a bound state extends from left to right (red), while the state at higher energy can be either extended (a) or local (c).Panels (b) and (d) show the corresponding local conductance.Panels (e) and (g) show the situation where the island hosts local subgap states, whose corresponding conductance is shown in (f) and (h).Notice the NDC region appearing in panel (d): they correspond to a particle blocked in one of the two state at energy ∆ and unable to tunnel out from the device on the other side, thus suppressing transport.Below the energy threshold required to populate the superconducting continuum (|V | = 2∆), transport is dominated by the properties of the subgap state alone: the low-bias peaks of panels (b) and (d) are almost identical and so are those of panels (f) and (h).

Figure 7 .
Figure 7. Upper left panels: sketch of the double wire system. 1 ○ one bound state connects L and R, while M connects to another one (a).Panels (b) and (c) are the conductance through the grounded M lead. 2 ○ three local states coupled to the different leads, 3○ an extended state coupled to all the leads, and 4 ○ a local state that couples to R, while L and R couple to another extended state.The parameters are the same as in Fig.2.The conductance we show for lead M is GM = dI M dV R , hence it is always negative, with the exception of some small regions.