Thermally induced subgap features in the cotunneling spectroscopy of a carbon nanotube

We report on nonlinear cotunneling spectroscopy of a carbon nanotube quantum dot coupled to Nb superconducting contacts. Our measurements show rich subgap features in the stability diagram which become more pronounced as the temperature is increased. Applying a transport theory based on the Liouville-von Neumann equation for the density matrix, we show that the transport properties can be attributed to processes involving sequential as well as elastic and inelastic cotunneling of quasiparticles thermally excited across the gap. In particular, we predict thermal replicas of the elastic and inelastic cotunneling peaks, in agreement with our experimental results.


Introduction
Due to proximity effects, a hybrid device made of a superconductor coupled to a mesoscopic normal conductor makes it possible to study a wide range of quantum phenomena. In particular, in the Coulomb blockade (CB) regime these include supercurrent transport carried by Cooper pairs [1][2][3][4][5][6], coherent electron transport in terms of multiple Andreev reflections [7][8][9][10], and quasiparticle transport [5,[11][12][13][14][15][16][17][18]. Andreev reflections lead to subgap structures with steps at bias voltage Δ ne 2 (  ∈ + n ) in the current-voltage characteristics [8-11, 13, 19, 20], which are smeared out by increasing the temperature [10,13]. In contrast, increasing the temperature favors quasiparticle transport, as it increases the probability of thermal activation of quasiparticles across the gap. The emergence of a zero-bias peak inside the Coulomb diamond by increasing the temperature [10,13] was explained in terms of resonant tunneling [11] of thermal quasiparticles. Recently the additional possibility of observing transport features due to sequential tunneling of thermally excited quasiparticles was theoretically proposed in [17] and experimentally confirmed in [18]. Such processes lead to thermal resonance lines within the Coulomb blockade region, parallel to the Coulomb diamond edges. Cotunneling processes due to quasiparticles, however, have so far been reported only for bias voltages above the superconducting energy gap [14,15]. In this work we present measurements in complete agreement with theoretical predictions regarding thermally excited quasiparticle transport in the cotunneling regime.
Cotunneling is a transport process in which the quantum dot (QD) is either excited (inelastic cotunneling) or kept in the same state as the initial state (elastic cotunneling) by means of events involving tunneling to an intermediate virtual state. Thus, for the inelastic case a bias threshold corresponding to the excitation energy is required to enable charge transfer [21]. In contrast with sequential tunneling processes, cotunneling at the lowest order is expected to be independent of the gate voltage.
We report on elastic and inelastic cotunneling spectroscopy on individual carbon nanotube (CNT) devices coupled to Nb superconducting leads. In the low-temperature limit, transport theory predicts for a CNT quantum dot superconductivity enhanced transport features at bias voltages Δ ± e 2 and Δ δ ± + e (2 ) m due to elastic and inelastic cotunneling of quasiparticles, respectively [14]. Here {δ m } is the set of excitation energies of the CNT from an N-particle ground state. With increasing temperature, we predict and observe the appearance of elastic and inelastic cotunneling features in the subgap region (i.e., for bias voltage amplitudes smaller than Δ e 2 ) due to thermally excited quasiparticles. In particular, the emergence of a zero-bias peak, corresponding to the thermal replica of the elastic cotunneling resonance, is expected. Our theoretical predictions are in good quantitative agreement with our experimental findings.
Individual single wall carbon nanotubes were grown on a highly p-doped Si/SiO 2 substrate by chemical vapor deposition [22]. The substrate acting as a global back gate is used to tune the electron occupation of the CNT. The source and drain electrodes were patterned on an individual single wall carbon nanotube by standard electron beam lithography and lift-off techniques. Here we report on measurements of two distinct samples. For sample A, figure 1, electrodes made of 3 nm Pd and 45 nm sputtered Nb with a spacing between electrodes on the order of 300 nm were used (see figure 1(a)); for sample B, figure A1 in the appendix, a metalization of 3 nm Pd and 60 nm sputtered Nb with a contact spacing on the order of 430 nm was applied (see figure A1(a)). To perform four-point measurements and as a resistive on-chip element, each superconducting electrode was connected to two leads made of AuPd to damp oscillations at the plasma frequency of the Josephson junction [23,24]. Low-temperature electrical transport measurements were performed inside a 3He/4He dilution refrigerator with a base temperature of 25 mK.
In both samples we observe regular CB diamonds over a large gate voltage range. Signatures of fourfold periodicity are observed in the measured gate range only for sample A. Figures 1(b) and A1(b) show the high-resolution measurements for the selected gate range for contacts in the superconducting state at temperature T = 25 mK and 30 mK, respectively. In both samples, lines of high conductivity are observed well inside the Coulomb diamonds; all these lines are horizontal, independent of gate voltage. To clearly identify them we restrict the gray scale for the differential conductance below the maximum conductance. Figure 1(c) shows a zoom corresponding to the region inside the diamond denoted ③ in figure 1(b) 3 . Horizontal lines are clearly visible and indicated by arrows in the conductance curve. T 300 mK and = T 1700 mK, respectively. The dashed white box corresponds to the range of gate voltages that is averaged to obtain the differential conductance curve shown on the left side of each figure.
One set of lines occurs at bias voltage ∼ ± V 0.52 SD mV (gray arrows). We ascribe it to elastic cotunneling processes at Δ = ± V e 2 SD . We extract Δ ∼ 0.26 meV for our superconducting film, compared with the expected value of Δ = 1.5 meV for bulk Nb. The mismatch of about a factor of 5 has already been reported in similar Nb-based devices [14,[25][26][27]. The reason for the gap reduction is still an open question. Possible explanations are the formation of niobium oxide, the thin composite of Nb and Pd, or the contamination of the lower Nb interface. For the deposited Nb/Pd strip a critical temperature of about 8 K was measured, where the resonant features remain present up to temperatures of about 4-5 K. Thus the transition temperature of the thin film is comparable to bulk Nb, that is, in contrast with the observed small value of Δ and the BCS relation Δ = k T 1.76 B . The inelastic part of the cotunneling spectra reveals excitations of the CNT quantum dot. Our data show a broad inelastic feature at a distance δ = 1.3 meV from the elastic line (black arrows). From additional stability diagrams for sample A, recorded at higher temperature and finite magnetic field to suppress superconductivity, we extract a charging energy Similarly, from the elastic and inelastic line, we can also extract Δ ∼ 0.23 meV and δ ∼ 0.11 meV for sample B. From additional stability diagrams in a regime in which superconductivity is largely suppressed, we identify a smaller charging energy ≃ E 3.2 C meV. The two samples have roughly the same superconducting gap Δ but differ in the charging energy E C , leading to different transport regimes. In both samples charging effects and the small coupling strength Γ Δ <  suppress Andreev processes, such that current is carried by quasiparticles. In sample Athe large charging energy further suppresses multiple quasiparticle processes. Thus the transport is dominated by sequential and cotunneling events. For sample B a simple description in terms of resonant tunneling of quasiparticles [11] may be conceived.
As the temperature is increased, new horizontal lines are observed. In sample A the novel lines arise for temperatures above ≈ T 600 mK at zero-bias and at bias voltage δ = ± V e SD . Figure 1(d) shows the same gate region as in figure 1(c) but now for the temperature = T 1.7 K. The additional lines, marked by stars, become more and more pronounced with increasing temperature. Andreev reflections do not give an explanation for the thermal behavior of such transition lines [4,[6][7][8]27]. Also, the Kondo effect cannot be the reason for the resonant peak at zero bias, as it has an opposite thermal behavior [28][29][30][31][32][33].
The feature of a zero-bias conductance peak is also supported by sample B, as shown in figure A1(c) in appendix A. The bias trace is taken in the middle of the Coulomb blockade valley at gate voltage ≈ − V 11.71 gate V. Upon increase of the temperature, one observes a rising conductance peak at zero bias and pairs of symmetrically displaced elastic and inelastic cotunneling peaks at finite bias. The feature at bias voltage Δ = ± V e SD and the thermal zerobias peak resemble data already reported in [10,13]. In analogous fashion, we expect these to be reproducible within the simple resonant model of [11]. The more complex behavior of sample A, where several cotunneling and sequential lines are observed within the CB diamond, clearly goes beyond the capability of the simple resonant picture that excludes Coulomb interaction. As shown hereafter, a full transport theory that includes all tunneling processes up to the second order in the strength Γ  of coupling to the leads can capture the experimental behavior in great detail.

Transport theory for S-CNT-S junctions
To understand the experimental observations, we consider a minimal model for a CNT quantum dot connected to two BCS-type superconducting leads. For the back-gated CNT we consider a single longitudinal mode incorporating orbital, m, and spin, σ, degrees of freedom. Coulomb interaction effects are considered within a constant interaction model, with U being the charging energy. The quadruplet CNT Hamiltonian thus reads (with = ± m 1, σ = ±1), where δ accounts for the breaking of the fourfold degeneracy of a longitudinal mode with energy ϵ d due to spin-orbit interaction and valley mixing [34].
The BCS superconducting leads are described by a conventional pair-interaction Hamiltonian on a mean-field level with respect to an offset energy E l 0 : This can be obtained by means of a particle-conserving Bogoliubov-Valatin transformation [35,36]: ⃗ v lk ( * ) [37,38]. Furthermore, the quasiparticles have an excitation energy 2 measured with respect to the electrochemical potential μ l . Finally, the BCS gap is defined by, where V | | characterizes the interaction potential between a pair of electrons. The connection with the superconducting leads is achieved by a single-particle tunneling where, for the sake of simplicity, the tunnel coefficient T l of lead l is considered to be spin, wave vector, and valley independent. The tunnel coupling strength can then be defined as , which is assumed to be energy independent.
We describe the time evolution of the system with the generalized master equation [39]: for the dynamics of the reduced density operator ρ red . This (still exact) equation allows a systematic perturbation expansion of the kernel superoperator τ K t ( , ) in powers of the coupling strength Γ  [40,41]. In the steady state limit and charge-conserved regime the master equation can be simplified further by applying the Laplace transform and its properties: i i i i can be represented in a diagrammatic language in which simple rules exist to directly obtain the corresponding analytical expression. In [42] these rules are derived and discussed in detail for the case of hybrid S-QD-S nanostructures.
An expression for the steady state current in terms of a perturbative expansion can be obtained in the same way. In particular, the net current of lead l is described by In the charge-conserved regime the reduced density matrix ρ χ χ′ i i is block diagonal (see appendix D). Thus the kernel element χ χ χ χ ′ ′ K f f i i up to the second order also represents the physical rate for processes transferring 0, 1, or 2 charges, depending on the charge difference between the states χ 〉 | i and χ 〉 | f . The problem of non-equilibrium hybrid superconducting-quantum dot junctions with an applied bias voltage is intrinsically time dependent. This can lead to time-dependent harmonic contributions to the stationary current associated with Andreev tunneling [9]. However, in the charge-conserved regime considered in this work, these harmonics are absent, and hence ρ → t˙( ) 0 red at long times. This is because the expectation values vanish since they break the conservation of total charge. Let us emphasize that, according to equation (4), we still have a finite superconducting gap and superconducting features (see appendix C for a detailed discussion).
Thermally assisted quasiparticle transport has up to now only been discussed in the context of sequential [17,18] and resonant [10,13] tunneling. Also responsible for the energy distribution of the fermionic quasiparticles, in addition to the BCS density of states (DOS), is the Fermi function. For high enough temperatures the Fermi function is thermally smeared, in the sense that quasiparticles can also occupy the high-energy branch of the DOS and thus can contribute to an additional transport channel. In the sequential tunneling regime, this gives rise to thermal replicas of the sequential tunneling transitions displaced by Δ ± e 4 in bias voltage (solid orange lines in figure 2(a)). When cotunneling processes are also taken into account, the number of expected thermal lines is greatly increased, as sketched in figure 2(a). In the figure we restrict ourselves to the exemplary Coulomb diamond denoted ③. Gate-dependent lines, induced by sequential processes, can be clearly distinguished from gate-independent cotunneling induced lines. Solid and dashed blue lines are transitions which are due to 'standard' sequential tunneling and cotunneling processes, respectively, i.e., contributions that are also present at low temperatures. Solid and dashed orange lines, in contrast, are due to thermally excited quasiparticles. Hence, they are present only at sufficiently high temperatures.
As already mentioned, standard elastic cotunneling lines are expected at bias Δ = ± V e 2 SD , and the inelastic cotunneling features occur at a bias

SD
, reflecting the excitation energy δ. Figure 2(b) visualizes the elastic cotunneling events in the many-body spectrum where the three-particle ground state is used as the reference energy. Choosing the center of diamond ③, corresponding to a certain gate voltage, the two-particle and four-particle ground states have the same energy. Thus, transitions from the three-particle ground state to the two-particle ground state and vice versa have the same probability as those from the threeparticle ground state to the four-particle ground state and vice versa, leading to elastic cotunneling. As shown following, thermal excitation of the lead quasiparticles yields thermal replicas at a bias Δ e 2 , smaller than for standard cotunneling features. We thus predict, in particular, the emergence of a cotunneling line at zero bias, being the thermal replica of the standard elastic lines at Δ ± e 2 . An exemplary contribution to elastic cotunneling in the diagrammatic language is shown in figure 3(a). Using the diagrammatic rules [41,42], the analytic expression is given by the kernel element including , and the energy difference δ = − ν χ E E E between the energy ν E of the virtual dot state ν〉 | and χ E of the dot state χ〉 | . Notice that in the example of figure 3(a), the state ν〉 | has one unit of charge more than state χ〉 | .The charges entering and leaving the dot carry the energies ω and ω′, respectively. An analysis of the double integral shows that, at low temperatures, it gives a pronounced contribution only in the case corresponds to the resonant case in which the highest occupied quasiparticle states in the source are aligned with the lowest empty quasiparticle states in the drain, such that elastic cotunneling onto and out of the CNT is possible. However, at higher temperatures thermally excited quasiparticles enable cotunneling transport also at zero bias. This mechanism is visualized in figure 3(b), where the numbers 1, 2, 3, 4 correspond to the tunneling events occurring at times τ ≡ ⩽ ⩽ ⩽ ≡ t t t t t 1 2 3 4 shown in figure 3(a). As seen in figure 3(b), if the thermally occupied quasiparticle states of the source are in resonance with the unoccupied quasiparticle states of the drain, elastic cotunneling through the dot can also occur at zero bias. The tunneling rate Γ ≡ for such a process is given by the expression in equation (8), adding the hermitian conjugated.
Mathematically, the condition for the onset of elastic cotunneling can be obtained from the analysis of the integrand ω ω′ I ( , ) of equation (8). This integrand is schematically depicted in  integrand ω ω′ I ( , ) in equation (8) is non-vanishing at only low temperatures in the blue region of the ω ω − ′ plane, depicted in figure 3(c). Upon increasing the temperature, the product is also non-vanishing along the orange stripes and on the orange spot.
In figure 3(c) figure 3(c). This corresponds to the gateindependent resonance at zero bias. In this simple resonance picture, we obtain the elastic forward cotunneling rate (see appendix E) in the middle of a Coulomb diamond by a first approximation of the integrand in equation (8): where we directly pointed out the bias dependence of the rate and introduced a degeneracy factor χ N depending on the state χ〉 | . Also including the backward process, the linear conductance is then approximated by in the limit Δ ≪ k T B . In particular, the former asymptotic characteristics indicate a transport property based on thermal excitation.
Analogously, subgap thermal replicas of the standard inelastic cotunneling lines are expected. We present a detailed analysis of the inelastic processes in appendix F and quote here the approximate result for the inelastic cotunneling rate: S D SD SD EC 2 similar to what was found in [14].

Comparison of theoretical and experimental predictions
In the following we use the BCS gap Δ, the excitation energy δ, and the charging energy E C extracted from the measured differential conductance plots to calculate the current through the CNT by means of the generalized master equation. Since the measured data revealed a relatively large critical temperature, we can assume a temperature-independent gap size in the considered temperature regime < T T 2 c . The calculations are performed by approximating . Figure 4(d) shows the corresponding experimental data for diamond ③. A short analysis of diamond ② is given in appendix B. In the bias and gate voltage range of figure 4(a) pronounced sequential tunneling lines and elastic and inelastic cotunneling features are seen. For better resolution we restrict the gray scale of the differential conductance below the maximum value. In figure 4(b) we focus on the Coulomb diamond denoted as ③. Beside the density plot we show the bias trace taken at the gate voltage marked by a white line, which supports the good quantitative agreement with the experimental data of figure 1(d). The standard cotunneling peaks (arrows) as well as their thermal replicas (stars) can be clearly recognized. The thermal behavior of the cotunneling features is illustrated in figures 4(c), (d), where the calculated and the measured differential conductance curves for different temperatures are presented. For the calculated curves we choose the same gate voltage as for the dashed white line in figure 4. For the experimental data we averaged a series of gate voltages marked by the box in figure 1(d). In both cases we emphasize that the standard cotunneling peaks are almost temperature independent, whereas the thermal replicas at zero bias and at δ = ± V e SD rise with increasing temperature.

Conclusions
In summary, we report on new cotunneling transport properties of a CNT contacted with two superconducting Nb leads based on thermally assisted quasiparticle tunneling. We observe the thermal replica of the elastic and inelastic cotunneling resonances with increasing temperature above 600 mK. These lead to an extra zero-bias peak and to an inelastic peak corresponding to the lowest excitation energy in the I V d d characteristics. To explain these non-equilibrium phenomena we derive a generalized master equation based on the RDM approach in the charge-conserved regime, applicable to any intradot interaction and finite superconducting gap. Modeling the CNT with a low-energy interacting spectrum, we find remarkable agreement with the experimental results concerning the thermal behavior of the additional cotunneling peaks. 4 We replace the Heaviside function Θ ω by a blurred step function. Despite γ being introduced empirically in this work, it can be shown that higher order processes involving quasiparticles lead to level broadening in the quantum dot and thus also to regularization of the divergence caused by the BCS density of states [11] similar to that provided by γ here.
We have in addition confirmed the prediction of a zero-bias peak due to thermally excited elastic cotunneling in another experimental setup. The description of sample B can be found in the main text. An atomic force micrograph of the studied quantum dot device is shown in figure A1(a). We observe regular Coulomb blockade diamonds over a large gate voltage range, also suggesting a defect-free CNT. In figure A1(b) we show high-resolution measurements for the selected gate range, including four Coulomb diamonds at temperature = T 30 mK. Inside the Coulomb diamonds we can identify gate-independent transition lines suggesting a symmetric coupling to the superconducting leads. To clarify the bias threshold of these horizontal lines, we take a bias trace of the region of interest at a fixed gate voltage ≈ − V 11.71 gate V pointed out by the dashed line. This enables us to observe the onset of a stable conductance peak for temperatures above ≈ T 600 mK which is more and more pronounced with increasing temperature. For that reason we assign the gate-independent conductance peak to a thermally assisted elastic cotunneling process. A detailed theoretical discussion follows in appendix D. However, a thermal replica of the inelastic cotunneling peak at bias voltage = ± V 0.11 SD meV cannot be clearly seen. This may be due to an overlap with the zero-bias peak.

Appendix B. Analysis of the Coulomb diamond ② of sample A
In figure B1 we show the bias trace of the measured differential conductance in the middle of the Coulomb diamond ② defined in the main text. It was obtained by the same averaging procedure as for the bias trace of diamond ③ explained in the main text. The curves for different temperatures include a richer peak structure than for the other diamond. We can identify the standard elastic cotunneling peaks at = ± V 0.55 ). The bias trace shows a zero-bias peak emerging at increasing temperature. To see the feature more clearly a conductance offset of about e h 0.03 2 was added systematically to each curve. The zero-bias peak is accompanied by the elastic and inelastic cotunneling peaks at negative and positive bias. the diamond ③ can be explained by a gate-dependent spin-orbit coupling in multielectron carbon nanotubes [44], as can also be seen in the overview figure 1(b). Besides, we can clearly recognize the rise of the thermal elastic cotunneling peak at zero bias with increasing temperature.
The additional features at bias voltage = ± V 1.2 SD mV cannot be explained by a singleshell model. Since the peak height is not temperature dependent, it must be a standard inelastic cotunneling feature. For our calculations we thus have to include a more complex excitation spectrum where the splitting to the next higher shell is smaller than δ. As was shown in [45], a two-particle ground state can lead to a rather complicated excitation spectrum where energetically close shells interact with each other, resulting in an effective shell splitting smaller than δ. Calculating such an effective Hamiltonian will remain a future task.
Also, the small peak at bias voltage ≈ − V 1.9 SD mV is almost temperature independent. By inspection of the stability diagram, we classify it as a cotunneling assisted sequential tunneling process (COSET). In such a COSET an excited state is populated by a preceding inelastic cotunneling process yielding to a gate-dependent sequential resonance peak inside the Coulomb blockade regime [41,46,47]. A more detailed discussion is left for future work.

Appendix C. BCS theory in the charge-conserved regime
In macroscopic superconductors with a large number of particles, the boson-like condensate is well described by a phase coherent state Φ〉 | with definite phase Φ. The presence of a relative phase between two weakly linked superconductors is at the origin of the Josephson effect [37,48,49]. In mesoscopic superconductors, charging effects due to Coulomb interaction break the degeneracy of states with a different number M of Cooper pairs. In such cases the phase Φ becomes uncertain and one has to project the state Φ〉 | onto a state 〉 M |2 with a fixed Cooper pair number M [50,51]. In this phase incoherent regime, the BCS Hamiltonian is thus properly diagonalized by means of the particle number-conserving Bogoliubov-Valatin transformation igure B1. Measurement of the differential conductance of sample A taken in the middle of diamond ② for different temperatures. More features are observed than for diamond ③, indicating a more complex excitation structure of the CNT spectrum. We can see two additional excitation transition resonances, indicated by arrows, not observed in diamond ③.
, as well as Cooper pair, Ŝ ( †) , operators. From the fermionic excitations described by the quasiparticle operators we demand γ γ δ δ = Moreover, the Cooper pair condensate and the quasiparticles are decoupled, i.e., By means of these commutator relations, we can further show that for the number operator N of the electrons it holds that i.e., the Cooper pair operator keeps the system in a state with a well-defined charge number: Together with equation (C.3), we conclude that the Cooper pair condensate is the vacuum state for the quasiparticles and that fermionic excitations can be described by In the phase incoherent regime, the equilibrium grand canonical density operator of the superconductor is given by  In the main text we claimed that the superconducting gap is not vanishing in the chargeconserved regime. The statement can be proved in the following way: the electron lead operators we need the Bogoliubov transform, equations (C.1) and (C.2). We then see also that Cooper pairs give contributions to the transport process through the acting operators Ŝ ( †) . We find: where in the last line we used the normalization condition . The second expectation value can be calculated in a similar way such that for the kernel component of the elastic cotunneling we obtain Figure E1. Diagrammatic representation of a relevant contribution to elastic (a) and inelastic (b) cotunneling. Necessary for the inelastic part is the energetic excitation of the final state χ′ in comparison with the initial state χ, with both states having the same amount of charge.
2 , ≡ ′ − ′ t t t 1 was applied in order to decouple the three time integrations. In addition, we expressed the energies ω, ω′ with respect to the electrochemical potential μ l .
To investigate the case when the double integral, and thus the kernel component, gives a relevant contribution to the transport dynamics, we analyze the integrand in detail. We are mainly interested in the region in the bias and gate voltage range in which the system is blocked into the ground state of the corresponding Coulomb blockade region in the sense of the sequential tunneling limit. If N electrons are trapped in the N-Coulomb diamond, the condition for a strong Coulomb blockade is μ . In our example the ground state energy for N charges is χ E , whereas ν E is the energy of the ( + N 1)-particle state ν〉 | . Hence, in the blockade regime is μ Δ − ≪ − is not vanishing. In figures E2(b) and E3(b) the roots of the denominators are represented by dashed lines. As explained in the main text, we are looking for the cases in which the roots meet the colored regions. In particular, the threshold for the onset of standard elastic cotunneling processes is obtained for those values of the bias voltage such that the ω ω = ′ root touches the corner of the blue region (see figure E2(b)). In that case the blue region includes the horizontal and the diagonal (ω ω ′ = ) zeros of the denominators. Thus the bias threshold Δ = ± V e 2 SD for the lowtemperature regime is obtained when the condition Δ μ ω ω Δ μ Note that for this bias voltage the diagonal zeros are located at the corner of the blue region, as seen in figure E2(b), where the product of density of state and the Fermi function has its largest value, resulting in a peak structure in the voltage characteristics. For higher temperatures additional scenarios have to be taken into account because the orange regions in figure E3(b) can no longer be neglected. Thus the condition for a strong Coulomb blockade has to be adapted to the low-bias regime, meaning μ Δ to prevent thermally excited sequential tunneling as shown in [18]. For our case this yields To see a rising of thermal elastic cotunneling, the diagonal ω ω = ′ root has to meet the orange region. Then one has to investigate the cases when the orange regions include the horizontal and the diagonal zeros of the denominators. In this situation we need only a minimal bias ⩾ V | | 0 SD . Thus, for large enough temperatures a remarkable contribution of the component of the kernel K EC for elastic cotunneling is always present in the bias-gate voltage range since the onset occurs at zero bias, as one can see by means of the condition Δ μ ω ω Δ μ − + = ′ = = − + D S . After analyzing the property of the integrand of the kernel element we can give a first approximation for the elastic cotunneling rate in the middle of a Coulomb diamond for the process shown in figure E1(a). As already explained we have to investigate only the integrand in the energy area, where ω ω ≈ ′. Furthermore we investigate the Coulomb diamond ③, in particular, as in the main text, where the center is placed at a gate voltage δ = + V E e ( ) .  as in equation (9). The result of a calculation of the linear cotunneling conductance in terms of the foregoing rate expression at zero bias is shown in figure E3(c).
Other diagrams contribute to elastic cotunneling. However, using the diagrammatic rules to evaluate their analytic expression, one realizes that they contain two different intermediate states, with one unit of charge more and one unit less than the state χ〉 | ; hence the two zeros of the corresponding denominator in the integrand are energetically far away from each other, resulting in a smaller contribution to the integral.
We also want to mention the dark orange dot in figure E3(b). In that region the kernel component K EC contributes at high temperatures only. This can be explained for the case where the diagonal zeros hit the area, resulting in a condition Δ μ ω ω Δ μ − + = ′ = = + D S for the bias threshold Δ = − V e 2 SD . The conductance peak in this bias region corresponds to an onset of a resonant charge current based on thermally excited quasiparticles in the drain producing unoccupied states in the low-energy branch of the BCS density of states even for high temperatures.