Active Quasiparticle Suppression in a Non-Equilibrium Superconductor

Quasiparticle (qp) poisoning is a major issue that impairs the operation of various superconducting devices. Even though these devices are often operated at temperatures well below the critical point where the number density of excitations is expected to be exponentially suppressed, their bare operation and stray microwave radiation excite the non-equilibrium qp’s. Here we use voltage-biased superconducting junctions to demonstrate and quantify qp extraction in the turnstile operation of a superconductor–insulator–normal metal–insulator–superconductor single-electron transistor. In this operation regime, excitations are injected into the superconducting leads at a rate proportional to the driving frequency. We reach a reduction of density by an order of magnitude even for the highest injection rate of 2.4 × 108 qp’s per second when extraction is turned on.

In superconducting circuits it is important to minimize the number of non-equilibrium quasiparticles as they deteriorate the operation of various devices, such as the coherence of quantum bits based on Josephson junctions [1][2][3][4][5] or Majorana nanowires [6,7], cooling power of superconducting microcoolers [8,9], sensitivity of kinetic inductance detectors [10][11][12] and performance of superconducting resonators in other applications [13,14].In principle bringing the system to temperatures T much below the superconducting transition should reduce the number of excitations, as at k B T ≪ ∆ (here k B is the Boltzmann constant and ∆ the superconducting energy gap) their equilibrium number density n qp is suppressed exponentially.It has been demonstrated, however, that when the devices are operated, quasiparticle (qp) excitations are created, caused typically by the drive signals or by stray microwave photons from hotter stages of the refrigerator [1,[15][16][17][18], or by ionizing radiation [19].To overcome this "qp poisoning", several methods have been studied.These include introduction of normal metal traps [20][21][22], geometry optimization [15,23], vortex traps by magnetic field [13,[23][24][25], gap engineering by variation of the film thickness [26][27][28][29] or phonon traps [30,31].Recently, blockage of qps by a voltage filter-tuned superconducting gap has been demonstrated [32], and employing a lower gap superconductor as a qp trap has been analyzed in detail [29].Although the crucial role of low qp density was identified already in early studies of superconducting qubits [33][34][35] and related single-charge circuits [26,36], further understanding of the generation mechanisms and reducing n qp remains the topic of an ever increasing intense research activity [17,[37][38][39][40] as the effort to increase the coherence times of superconducting qubits continues.
For radio frequency (rf) driven superconductorinsulator-normal metal-insulator-superconductor (SINIS) turnstiles [41], operated as a source of quantized electric current [42], both drive-induced and background qps are the most severe limitation to reaching metrological current quantization accuracy [15,42].In this work we show that such a hybrid single-electron transistor (SET, see Fig. 1(a) for a typical sample) functions as a sensitive, practical, and quantitative detector of the qp density of its superconducting electrodes, that could be integrated with a variety of other mesoscopic superconducting devices to probe their n qp .Under rf drive, hybrid SETs present a turnstile for single electrons, a simple-to-operate candidate realization for a solid-state standard of electric current [41,43].For a wide range of parameters, the output current I = ef is determined only by the drive frequency f and electron charge e (see Fig. 1(d) for typical measurements in this regime).Here the non-equilibrium state results from qp injection (Fig. 1(b)) when the drive signal at frequency f is applied to the gate electrode of the transistor.Earlier work [15] demonstrated that the influence of drive-induced and environmental qps to a SINIS turnstile can be reduced by improving the geometry of the S electrodes and by shielding from residual microwave radiation, respectively.
Here we combine the turnstile with an independent, in-situ control -both extraction and injection -of qps.The qp poisoning can be reduced by suitably voltage biasing superconductor-insulator-superconductor (S 1 IS 2 ) junctions with different superconducting gaps (∆ 1 and ∆ 2 , respectively) where excitations are extracted from S 1 to S 2 as long as ∆ 1 < ∆ 2 [44][45][46][47].In an important experiment, this effect has been used to cool one of the leads of a single-Cooper pair transistor [48] in a similar manner as refrigeration by normal metal-insulator-superconductor junctions [49][50][51] has been used to cool down the normal lead.However, when using a single Cooper-pair transistor it was not possible to control the mechanism or rate of qp creation, making it more difficult to quantify the population reduction due to the biased S 1 IS 2 junction.On the contrary, hybrid single-electron transistors allow for quantitative control.To that end, here we demonstrate that direct S 1 IS 2 cooling of the turnstile leads offers promise to fully extract the drive-induced qps under typical pumping conditions, in particular when bulky or thick electrodes cannot be utilized.
The evacuation of qps is manifested by the stabilization of current to ef in the SET turnstile operation [15].The qp extraction can be tuned by varying the voltage biasing of the S 1 IS 2 junction.However, non-equilibrium qps are also subjected to recombination and diffusion processes along the electrode as depicted in panel (b) of Fig. 1.Due to the exponential dependence of n qp ≈ D (E F ) √ 2π∆k B T e −∆/kBT on ∆ the number of excitations is lower in the higher gap film.Here D (E F ) denotes the normal-state density of states at the Fermi energy and T is the electronic temperature of the superconductor.Therefore, providing sufficient energy to qps in the lower gap superconductor, biasing the junction such that the singularities in the superconducting densities of states align, will promote a transfer to the higher-gap superconductor (see the inset of Fig. 1(b)) [52].Although two superconductors with different energy gaps are required, the qp extraction in a voltage-biased tunnel junction is in sharp contrast to those gap engineering methods [26][27][28][29] where qps are passively trapped in lower-gap regions away from the relevant operational zones of the device.To implement the qp evacuation and the probing of n qp experimentally, we fabricate and measure a series of aluminum-based samples, cooled down in a dilution refrigerator reaching electronic base temperature T b ≈ 50 mK.
Here we show detailed results for two devices for confirmation of results, both with copper as the turnstile normal-metal island, separated by aluminum-oxide barriers from two aluminum leads.The samples were patterned by electron-beam lithography and metal deposition was done by multi-angle shadow electron-beam evaporation, see Supporting Information for further details on the fabrication process.One of the leads (left in Fig. 1(a)) is made wide to trap qps passively while the other (right) is long and narrow, this way promoting an excess qp population: due to the geometry, qp diffusion away from the turnstile junction and their subsequent relaxation is intentionally poor in this narrow lead [15].The samples were fabricated using standard electron-beam lithography and shadow mask techniques.The narrow lead is the S 1 part of a S 1 IS 2 Superconducting QUantum Interference Device (SQUID) whereas the fork-shaped electrode forms S 2 with a higher superconducting gap.The different superconducting gaps are achieved by depositing a thicker aluminum layer as the turnstile lead (d 1 ≈ 70 nm) and a thinner one for the rest of the SQUID, (d 2 ≈ 8 nm) [53].A magnetic field perpendicular to the plane of the sample is applied to produce a flux Φ = Φ 0 /2 (here Φ 0 = h/(2e) is the magnetic flux quantum) in the SQUID loop and therefore suppress the supercurrent in the S 1 IS 2 junctions so that subgap V SIS can be applied [48] (see also the Supporting Information).
We characterize the samples by sweeping the gate voltage V g between the "closed" (n g = 0) and "open" (n g = 0.5) states (with n g = C g V g /e, where C g is the gate capacitance), and stepping the SET bias voltage V b .The S 1 IS 2 junctions are first kept unbiased, named as the "cooler-off" case.The main parameters of the samples have been extracted from these measurements by fitting the measured current I to current-voltage curves calculated with a master equation approach taking into account sequential tunnelling and Andreev reflections [54] (see for example Supporting Information and Ref. 55).One of these fits can be seen in Fig. 1(c).For sample A (B) the total resistance of the tunnel barrier is R T = 159.9kΩ (63.0 kΩ), the charging energy of the island E c = 0.95∆ 1 (0.50∆ 1 ), the thickness of the N island d = 30 nm (40 nm) and the Dynes parameter [18] η = 1.0 × 10 −4 (7.5 × 10 −5 ).For both samples, the superconducting energy gap ∆ 1 of the leads is 180 µeV, and the island's lateral dimensions l × w are 1 µm × 100 nm.See Supporting Information for further details on the employed model where it is explicit that the island temperature is calculated according to power balance and that the area of a single conduction channel (A ch ) can be extracted from these data.Yet it is more precise to extract A ch from pumping measurements since the effects of the second-order tunnelling are more pronounced in that operation regime.Since qp diffusion is altered by the different lead shapes [15] the DC fits were made by leaving temperature of the long lead as a free parameter and assuming the base temperature of 50 mK for the wide lead, estimating an electronic temperature ∼ 150 mK for the long lead.The values of A ch , 9.5 nm 2 (10 nm 2 ) for sample A (B), are lower by about a factor of two compared to previously reported values [56]; yet the fits are relatively insensitive to the exact value of this parameter.We estimate these values from data with the cooler bias V SIS near to certain optimal point, where the Andreev effects are more noticeable, as will be seen later.By analysing the differential conductance of the SQUID it was possible to estimate ∆ 2 to be ∼ 235 µeV and ∼ 240 µeV ± 5µeV for sample A and B, respectively.
The strong cooling by the superconducting junctions is evident in the turnstile operation of the SET shown in Fig. 2. In this regime a sinusoidal signal with an offset equivalent to n g = 0.5 is applied to the gate electrode and the amplitude of this signal is swept while the SET current is measured.The presented data shows these measurements zoomed to the first current plateau (see Fig. 1(d) for a typical measurement in a wider range) for both samples.The plots display the behaviour with the SQUID at zero bias and at a finite bias close to where optimum cooling is achieved.Besides displacing the plateau level from the expected value of I = ef , an elevated qp density in the leads close to the junctions makes the level depend on the SET bias voltage and hence the values of the current at the plateaus will spread [15].In the two cases of Fig. 2 the "cooler-on" (V SIS = 0) curves show a smaller deviation from the expected current than the cooler-off ones.Additionally, the bias dependence of the plateau level is weaker with the cooler active.All this means that the density of excitations has been suppressed in the narrow lead by biasing the S 1 IS 2 junctions.For these particular devices, the remaining deviation is explained by the relatively low E c < ∆ which promotes Andreev tunnelling, as well as by leakage current [55] and remnant qp population as will be seen later.
To quantify the reduction of the qp density, the pumping data at the plateau shown in Fig. 2 are compared with a model (black solid curves) based on a master equation approach similar to that in the DC case.In the model the temperature of the island is varied at each amplitude value according to the proper power balance.Furthermore, since the used signal frequencies are high (f > τ −1 eff , τ eff being the effective qp relaxation time), it is suitable to model the system as if the temperature of the superconductors was constant throughout the operation cycle [16].Additionally, electron-electron relaxation is assumed to be fast enough to avoid branch imbalance.With these assumptions the comparisons with the data are done with the temperature of the narrow lead as the only free parameter (device parameters are fixed by the DC measurements).As noted before, n qp is related to the effective electron temperature T in the superconducting lead obtained by means of the thermal balance model.These simulations allow to deduce that the biasing of the S 1 IS 2 junctions effectively cools the narrow long lead.Furthermore, the use of the pumping operation of a hybrid SET as a sensitive detector of qps is justified by these calculations.Fig. 3(a) shows a specific example of sample B biased at V b = 100 µV and gate modulated with an amplitude of around A g = 3.5 mV.There it is evident that the plateau level grows with the excitation density in the lead.The behaviour is nearly linear at high population but flattens for low values.The behaviour of the qp density when the cooler bias changes is analysed by performing measurements where V SIS is swept while the gate drive amplitude A g is fixed to approximately middle of the first current plateau.Examples of such measurements are shown in Fig. 3(b) for device A, and in panels (c) and (d) for the B device.The plateau level and the spreading of them for different bias values of the SET, decrease until we reach eV SIS ≈ ±|∆ 2 − ∆ 1 | [57] (indicated by a dashed line in panels b-d).At this point the temperature of the lead is at a minimum and it changes only weakly until a threshold at which it starts to grow rapidly.Therefore we can assert that when the cooler bias voltage reaches the optimum point the excitation population in the lead is at a minimum and at higher biases it starts to grow up first at a low rate and then rapidly when there are peaks in sub-gap current of the S 1 IS 2 SQUID (see Supporting Information for further details).Similarly, biasing at V SIS < 0 affects the system in the same way but this time by extracting hole-like excitations (consider the inset in Fig. 1(b) with the lower singularities aligned).We can also understand this by looking at the modelled cooling power of the S 1 IS 2 junction [58] (see Supporting Information for details) when it begins to grow with V SIS towards a sharp peak, where the lead temperature is minimum, and then it decreases dramatically, corresponding to the zone with finite but less efficient cooling.Finally, at eV SIS ≈ ∆ 2 +∆ 1 qps are injected into the lead with the smaller gap.No sharp dip is observed in the measured I vs. V SIS curves at eV SIS = |∆ 2 − ∆ 1 |, similar to earlier works on S 1 IS 2 cooling [47,48].The singularity matching peak in the cooling power is likely washed out due to low-frequency noise in the S 1 IS 2 cooler bias voltage, finite sub-gap density of states, and local inhomogeneities in the superconducting gap.
The influence of the S 1 IS 2 cooler is qualitatively the same for both samples and independent of the operation frequency as well as of the SET bias, although the plateau levels are different, see the supplementary material for additional measurements to corroborate this fact.Note that this level approaches the ideal value but never reaches it even in the calculations for samples with low E c < ∆ (see Fig. 3(a)).In these devices with low E c we observe a clear crossover from excess qp-dominated to Andreev reflection-limited current quantization as the cooler is turned on.
We estimate the lead temperature for a wide range of operation frequencies.Figs.4(a) and (b) show that the lead temperature T extracted from the fitting procedure, and in turn n qp (Figs.4(c) and (d)) grow with increasing frequency also in the cooler-on case.However, in sample B the curve for the cooler-on case flattens at f > 40 MHz since the SET behaviour loses sensitivity to qp density in this operation regime (see Fig. 3(a)).For the near optimal cooler bias (V SIS ≈ 50 µV and 60 µV for A and B, respectively) the dependence of n qp versus the operation frequency is resemblant.This comes from the fact that qp transport is dominated by diffusion.There is dramatic drop in the qp density due to the biasing of the S 1 IS 2 junctions.By linearly extrapolating the densities to the f = 0 limit it is possible to see that there is a reduction between the "off" and "on" case by an order of magnitude also in this situation.We conclude that the cooler suppresses also the ever present excess qp population n qp,0 , generated by non-equilibrium environment [59][60][61], which in turn should diminish the subgap current in the DC regime.On the other hand it is seen that the biased S 1 IS 2 junctions cannot totally remove this excess qp population in the leads.We expect that even more efficient evacuation can be achieved by placing the S 1 IS 2 junction closer to the injection junction -or even under it -since the highest concentration of qps is near this point [21,51], and further diffusion, recombination and phonon pair breaking, among other phenomena can thus be avoided.
We measured further reference samples with identical aluminum lead geometry but without S 1 IS 2 SQUIDs and obtained the qp densities in them (for results see Supporting Information).The qp densities are similar to those shown in Fig. 4 for the cooler-off case, and thus there is no significant influence of the unbiased cooler junctions or the presence of their biasing circuit on the relaxation of the non-equilibrium excitations.Hence the presented cooler-off densities should correspond approxi-  mately to densities in the case if the S 1 IS 2 junctions were not present.Thus the diffusion model in [15], which is based on heat conduction, without considering normal metal traps, should hold.Within this model the qp density is given by where ℓ, ̟, ϑ are the dimensions of the lead (25 µm × 100 nm × 70 nm), P inj is the injected power which can be approximated by ∆ 1 f in the driving regime of the experiment [15].In Supporting Information we show that this is a good approximation.Notably, even at f = 100 MHz, the injected power for an aluminium-based device with ∆ 1 = 200µeV is only around 3.2 fW, extractable by a sub-micron S 1 IS 2 junction.In addition, ρ n is the aluminium normal state resistivity.Based on measurements of a number of separate test structures, we estimate ρ n ≈ 31Ωnm at 77 K (see the supplement for details).A few of these 4-probe structures were cooled down to 4.2 K where we estimate a further 10% decrease of ρ n below its 77 K value.The measured qp numbers are larger than those predicted by this model with no free parameters.Fitting Eq. ( 1) to the data using ρ n as a free parameter yields ρ n ∼ 90 Ωnm.These fits are shown as black lines in panels (c) and (d) of Fig. 4. Further experiments are needed to understand the discrepancy, observed also in [62] for devices with higher charging energy.At f = 0, V SIS = 0, we estimate n qp,0 ≈ 250 µm −3 for both samples, obtained from fits to dc IV charac-teristics of the turnstile.These background qp densities for the cooler-off case are roughly two orders of magnitude higher than those observed for highly-shielded SI-NIS turnstiles [15], or superconducting qubits and resonators [11,37,63].The excess in n qp,0 compared to Ref. 15 originates from the lack of a microwave-tight indium seal in the sample holders employed in these measurements and from the lack of any normal metal traps and restricted geometry -narrow, thin, long -of the S electrode with the S 1 IS 2 junctions.To bring down n qp,0 , a sample holder with a higher level of shielding can be utilized, and the S 2 electrode can be fabricated with separate lithography and deposition steps that do not limit the junction geometry as with the multi-angle shadow evaporation technique where both NIS and S 1 IS 2 junctions are created with a single mask.
The applicability of the qp extraction and detection techniques considered in this work extends beyond SINIS turnstiles and improving the accuracy of their pumped current.First, one can envision straightforward integration of S 1 IS 2 coolers for instance with absorbers of kinetic inductance detectors, or into superconducting resonators analogously to NIS coolers [64].When the bottom electrodes of the junctions are fabricated in a separate lithography and deposition step, their area and geometry can be adjusted at will.Furthermore, instead of tuning the gap by the film thickness, different superconducting materials can be utilized.Secondly, the SINIS turnstile, demonstrated here to function as a sensitive and direct qp probe, can be combined with various superconducting devices to measure the background qp density.This high-impedance, non-invasive probe can test the level of microwave shielding also in setups with sensitive nonsuperconducting devices.Finally, we have shown the driven hybrid turnstile to act as a highly controlled qp injector that could be used to investigate the qp sensitivity and qp trap efficiency of superconducting qubits and other devices.
In summary, we have been able to demonstrate active extraction of non-equilibrium excitations from a superconductor.A more than an order of magnitude reduction from values as high as 5.8 × 10 3 µm −3 to 260 µm −3 at the highest studied injection rate of 2.4 × 10 8 s −1 was achieved.Furthermore, in the limit of no qp injection (f → 0) we find a similar reduction although a finite population of environmentally created excitations remains.Our work shows that the qp density in the superconducting electrodes can be controlled by active qp traps, here demonstrated for the first time with in-situ control of qp injection.
Supporting Information.Calculations of the cooling power of S 1 IS 2 junctions are contained here, as well as a model for the power transfer in a driven system.Additionally, details on fabrication and measurement procedures are given.Further details on the theoretical background and numerical calculations of SINIS SETs as well as results for additional and reference measurements are provided.Finally, details on the aluminum resistivity are given.A simplified model of a resistively shunted superconducting junction can be developed for calculating the cooling power of one of the superconductors [1].The system consists of a Josephson junction composed of a tunnel barrier between two leads with superconducting gaps ∆ i and temperatures T i with i = 1, 2. The junction is connected in series to a resistance R S and fed by a bias V SIS .The system dynamics is given by where e is the elemental charge and I c is the critical current of the junction.For a single junction, I c can be approximated by [2] Here, f 1 denotes the distribution function of the superconductor 1 (often taken to be the Fermi-Dirac function), whereas R is the normal-state tunnel resistance.
In the experiments, instead of an individual S 1 IS 2 junction we use a DC superconducting quantum interference device (SQUID) with two junctions in parallel to be able to minimize I c by a magnetic flux Φ. Neglecting the inductance of the loop, the critical current is given by , where I c,L/R are the critical currents of the individual junctions and Φ 0 = h/(2e) is the flux quantum.We apply an external magnetic flux Φ = Φ 0 /2 through the SQUID loop so that I c = |I c,L − I c,R | is at a minimum, due to typical junction size asymmetry at a value around 0.5 nA.
Suppressing I c allows us to voltage bias the SQUID at subgap voltages, essential for the S 1 IS 2 cooling.In the simplest approximation, we neglect the cooling power of the SQUID junction that is further away from the SINIS turnstile, and model the SQUID as a single junction with the suppressed I c .
Combining Eqs.(1) one obtains V SIS (mV) In a Josephson junction φ = 0 when The power transferred from the lead i = 1 is written as P (1) (t) = P J (t), where P qp (t) is the power transmitted by quasiparticles (qp's) given by P (1)  qp In (5), n i denotes the density of qp states in lead i, modeled by the broadened BCS DoS Here, η is the Dynes parameters which models subgap leaks [3,4].Additionally, f i denote the distribution functions.
On the other hand the Josephson contribution [1] P J (t) can be written as P sin (t) sin ϕ, where When averaging the power in time the sine contribution vanishes and the power extracted from the lead i = 1 is When the junction is biased above critical current (V SIS > I c R S ) the cosine component vanishes in average.If ∆ 2 > ∆ 1 the power is positive and the average is peaked at eV SIS = ±|∆ 1 −∆ 2 |, as shown in Fig. S1 for various Dynes parameters.Notice how the peak gets smeared at large η.Additionally, observe that for large Dynes parameters the power becomes negative at biases closer to ∆ 1 , indicating injection of qp's to S 1 .For low leakage cases this heating does not start until close to eV SIS > ∆ 1 + ∆ 2 .Figure S2(a) shows the large-scale behavior of the current at the plateau against V SIS for sample A at f = 40 MHz and V b = 200 µV.As already shown in the main manuscript the lowest current is found around eV SIS = |∆ 2 − ∆ 1 | in agreement with the model.Here, the response is smeared likely by subgap leakage, local inhomogeneities of the gap, and noise in the bias voltage.Further discrepancies with the model such as the strong heating around ∆ 1 correspond to peaks in subgap current of the SQUID as can be seen in Fig. S2(b) and are not caught by this simple model.Furthermore, a much weaker heating peak at eV SIS ≈ ∆ 2 is present.Finally, observe that when the current-voltage curve starts to become linear at eV SIS > ∆ 1 + ∆ 2 there is a strong heating since qp's are now injected from S 2 to S 1 in accordance with the basic model.To simulate the subgap peaks, expressions for the cooling power should be developed and solved in a framework that includes multiple Andreev reflections, analogous to calculations of the electric current e.g. in Ref. 5. In a more detailed model also the temperature gradient along the S 1 electrode and the electromagnetic environment of the S 1 IS 2 junction should be considered.

POWER TRANSFER IN A DRIVEN SYSTEM
The process of transferring electrons from the island to the leads can be regarded as a jump stochastic process whose probability of jumping is governed by the equation The reverse process is assumed to vanish and Γ is the rate, which in the slow driving regime depends on the driving as Γ = Γ 0 e − β(1−ε) .For the present case, Γ 0 = ∆ e 2 R , with R the tunnel resistance and β = ∆ kBT (k B is the constant, T is the temperature and ∆ the superconducting gap).Furthermore, ε(t) = 2E c /∆ (n g (t) − 1/2) is the energy of the transferred electron normalized to the energy gap, where E c is the charging energy of the island and n g is the charge induced by the driving signal.
The solution of Eq. ( 9) is Now, notice that the power transfer is given by P (t) = Ė = ∆εΓ(t)p(t).The average transmitted power during one drive period of length τ = 1/f is then given by In the conditions of interest β ≫ 1 and Γ ≫ f , with f the frequency of the driving.The integral is thus ∼ 1 for a general periodic driving and P ≈ ∆f .A specific case of the dynamics of the instantaneous power transfer is shown in Fig. S3(a).Notice how it sharply peaks at a specific time making its integral nearly one and the energy transmitted during one cycle ∆.Fig. S3(b) depicts how the power evolves with the driven electron energy for a quarter of cycle.

Fabrication
The samples were fabricated on 4-inch silicon substrates covered by 300 nm thermal silicon oxide.The process is based on electron beam lithography (EBL, Vistec EBPG5000+ operating at 100 kV) and multiangle shadow deposition of metals in an electron-beam evaporator.First, a set of bottom gates (visible e.g. in Fig. 1(a) of the main text) and large ground plane electrodes were formed by depositing a 2 nm titanium adhesion layer, 30 nm gold, and another 2 nm Ti protection layer through a mask defined in a single layer positive resist (Allresist AR-P 6200).After conventional liftoff, the gate and groundplane electrodes are covered by a 50 nm insulating layer of Al 2 O 3 grown by atomic layer deposition (ALD).Next, another round of EBL and metal evaporation (2 nm Ti followed by 30 nm AuPd) is used to form bonding pads and coarse electrodes for connecting to the tunnel junctions, to be patterned in the third and final lithography step.For defining the S 1 IS 2 and NIS junctions, a Ge-based hard mask is used [6].It is composed of a 400 nm sacrificial layer of P(MMA-MAA) copolymer, covered by 22 nm Ge also deposited by e-gun evaporation, and finally a thin (approximately 50 nm) layer of PMMA on top.Before development of the Ge mask, the wafer was cleaved into smaller, typically 1 cm x 1 cm chips.The pattern defined into the PMMA resist is transferred to the Ge layer by CF 4 reactive ion etching (RIE).Subsequently, an undercut profile is formed in the copolymer layer by oxygen plasma etching in the same RIE system.
The tunnel junction deposition starts by evaporating a 13 nm layer of Al at a substrate tilt angle of approximately 55 • , resulting in an approximately 8 nm film that forms the fork-shaped S 2 electrode of the SQUID.Right after deposition, this film is oxidized in situ in the evaporator (static oxidation with typically 2 mbar of pure O 2 in the chamber for two minutes).A subsequent deposition of 70 nm aluminum at normal incidence angle (zero tilt) completes the S 1 IS 2 SQUID and forms the S 1 electrodes of the turnstile.Finally, the SINIS turnstile and the whole device are completed by a second oxidation (typically 1 mbar O 2 for one minute) and depositing 50 nm of copper under a tilt angle of 39.5 • (in the opposite direction compared to the S 2 layer), resulting in an approximately 40 nm thick Cu island.For measurements at cryogenic temperatures, a chip with an array of 3-by-3 devices is cleaved to fit a custommade chip carrier and electrically connected to it by Al wire bonds.Due to the limited number of measurement lines, typically 2-3 devices from the chip are selected to be characterized at low temperatures.

Measurements
Experiments were carried out in a custom-made plastic dilution refrigerator with base temperature of about 50 mK.Conventional cryogenic signal lines (resistive twisted pairs between room temperature and the 1 K pot flange, followed by at least 1 m Thermocoax cable as a microwave filter to the base temperature), connect the bonded chip to a roomtemperature breakout box.The rf line for applying the drive signal to the turnstile gate consists of a stainless steel coaxial cable down to 4.2 K, a 20 dB attenuator in the liquid helium bath, followed by a feedthrough into the inner vacuum can of the cryostat.Inside the cryostat, the rf signal is carried by a continuous superconducting NbTi coaxial cable from the 1 K stage down to the sample carrier.Magnetic field for the S 1 IS 2 SQUID was applied by current biasing a superconducting coil placed on the exterior of the inner vacuum can.Current and voltage biases were realized by programmable voltage sources and function generators.Current amplification was achieved by room-temperature low-noise transimpedance amplifiers (FEMTO Messtechnik GmbH, models DDPCA-300 and LCA-2-10T).The curves of the pumped single-electron current were typically iterated at least 10 times and averaged accordingly, neglecting those repetitions during which an offset charge jump had occurred.Amplifier offset currents were subtracted by comparing the pumping curves with their counterparts measured under source-drain bias of opposite polarity.

NUMERICAL CALCULATIONS
The current-voltage characteristics of a SET can be modelled from a stochastic master equation on the charging events n, where n designates the charge state of the island.The probability p (n) evolves as [7][8][9][10] where γ nn ′ is the total transition rate from the n state of the island to the n ′ state.Furthermore, in the steady state dp (n, t) /dt = 0 and since the charging events n are discrete one can express Eq. ( 12) as the matrix equation where Generalizing for two leads, as is the case of an SET, the transition rates are given by γ nn ′ = Γ l n→n ′ (δE) + Γ r n→n ′ (δE).Here, Γ is the transition rate between individual charging states n and n ′ , l designates events between the left lead and the island and r similarly stands for the right lead.Finally, δE is the energy change of the process given by δE ±,l/r 1e for single-electron (1e) and two-electron (2e) tunnelling, respectively.Here V l = κ l V b and V r = −κ r V b , κ l/r is the ratio between the junction capacitance and the total capacitance, V b is the bias voltage applied between the leads of the transistor.In Eqs. ( 14) and ( 15) n is the initial island excess charge, n g is the charge number induced by the gate voltage, E c is the charging energy and + (−) designates tunnelling to (from) the island.The explicit expression for the rates depends on the specific system and on the transition processes taken into account.Since here we consider transitions in NIS junctions and only single-electron and two-electron Andreev processes are taken into account, the rates are given by for 1e tunnelling, and for Andreev tunnelling, and δE is the corresponding energy change from Eqs. ( 14) or (15), respectively.In Eqs. ( 16) and ( 17), ∆ is the superconducting gap of the leads, R l/r is the tunnel resistance of the junction involved in the event either left or right and N is the number of conduction channels which can be written as A/A ch with A being the junction area and A ch is the area of an individual channel.The term δ takes into account the energy of the intermediate (single-electron tunnelling) state which has a finite lifetime [11].Additionally, f N is the Fermi-Dirac distribution of the normal-metal island, f S is that for the superconducting lead involved in the tunnelling event and n s is the superconducting density of states given by Eq. ( 6).Furthermore, Once the probability vector p is obtained, the current through the SET can be calculated as In order to get realistic results the temperature change of the island has to be taken into account.To do this one calculates the power transferred to the normal island in a single electron event as where N refers to normal metal and δE is again the related energy change cost from Eq. ( 14).The total power transferred to the island is calculated as Q = q • p with q i = QN i→i+1 + QN i→i−1 , where QN i→i±1 = QN,r i→i±1 + QN,l i→i±1 .The heat flow to the phonons in the N island is governed by Qe−ph = VΣ T 5 N − T 5 b where V is the volume of the island, Σ is the electron-phonon coupling constant (≈ 2.5 × 10 9 WK −5 m −3 for copper, close to previously measured values [12,13]), T N is the electron temperature of the island and T b is the temperature of the phonon bath which is considered to be the same as the cryostat mixing chamber temperature.An additional power transfer due to Andreev reflection is considered in the form of Joule heat, i.e.QA = I A V b where I A is the current due only to Andreev events and V b is the applied bias voltage.Finally, the heat balance is Qe−ph = Q + QA .This condition is used to solve for T N .
A similar method can be applied when a periodic gate drive is applied to the system [10].If the period of this driving is τ then it is valid to assume that the steady state probability satisfies p (t) = p (t + τ ).In order to solve for the probability the cycle is discretized in m intervals of length ∆t = τ /m, next the matrix A (k∆t) = A k is calculated for each interval as well as the b k vector.If we now build the propagator then the initial probability is given by since for a periodic driving p(0) = p(τ ).
A more useful form of the propagator is given by where In Eq. ( 23) 0 is an appropriate vector of zeros.This way Eq.( 12) can be reformulated in terms of the augmented rate matrix Ãk and a new probability vector p (t) = [p (t) q ] T , where q is the average charge transferred during one cycle.The final propagator can be decomposed as Then, Eq. ( 21) is reformulated in terms of this new propagator and the average charge is obtained as q = U b • p(0).The average current pumped is I = q /τ .To complement this calculation the temperature of the island is calculated self-consistently, the same power balance is applied and the average heat transferred is now calculated with a similar propagator.However, now b k has to be replaced by q k and therefore Q = Q /τ .

Sample with different ∆1 and ∆2
A third sample similar to A and B (see Fig. S4(a)) was measured.However, in this sample the turnstile leads are thinner (d 1 ≈ 30 nm) and therefore their superconducting gap is higher than in the main samples A and B. Additionally, the superconductor S 2 is thicker than in the main samples (d 2 ≈ 17 nm).The parameters of the SET, determined by comparison with simulations (see Fig. S4(b)), are ∆ 1 = 205 µeV, E c = 0.53∆ 1 , R l/r = 97 kΩ, A ch = 5.5 nm 2 , and η = 2.75 × 10 −5 .The NIS junctions of this structure turned out to be asymmetric with a capacitance ratio of C r /C l ≈ 0.2 where l is for the left and r for the right junction.
As expected, the narrow lead cools down when a bias voltage V SIS is applied to the SQUID (see Fig. S5) until a certain optimal point is reached.This bias is designated by the dashed black line in Fig. S5.As can be seen, the optimal cooling occurs at V SIS ≈ 25 µeV, a lower bias than for samples A and B due to the smaller difference |∆ 2 −∆ 1 |.
Notice that the behaviour shown in Fig. S5 is qualitatively the same as that observed in samples A and B.
In order to estimate the qp reduction in the narrow lead, the device was first operated in the cooler-off regime (V SIS = 0 µV).Using the same procedure as for samples A and B, the narrow lead temperature and excitation densities were determined in the turnstile operation at different gate signal frequencies by comparing these measurements with simulations.The results can be seen in Fig. S6.It can be seen that the densities are much lower than in the main measurements, which can be understood considering that qp's are more costly to generate due to the higher superconducting gap ∆ 1 .
In addition, the system was operated near the optimal cooling point (V SIS = 28 µV).It was observed that the qp density of the narrow lead varies only slightly throughout the different employed frequencies.The temperatures (and excitation densities) remain approximately at 180 mK (5.98 µm −3 ) for 10, 20 and 40 MHz, and there is a slight increase to 190 mK (12.32 µm −3 ) at f = 80 MHz.This could be due to an increased recombination rate that is only overcome at 80 MHz.Again, the qp densities have been reduced by an order of magnitude with respect to the cooler-off case.

Reference sample without S1IS2 cooler
In order to check if the unbiased S 1 IS 2 junctions affect the excitation density a reference sample without the SQUID was measured.As can be seen in Fig. S7(a), the geometry of the reference sample is otherwise the same as for samples A and B described in the main text.After comparing the DC measurements with simulations, it was determined that ∆ 1 = 190 µeV, E c = 0.63∆ 1 , R l/r = 209.5 kΩ, η = 8 × 10 −5 , and A ch = 20 nm 2 .
Following the same steps as with samples A and B, current measurements in the turnstile operation of the reference SET were done and compared to simulations as shown in Fig. S8.From these comparisons the temperature of the narrow lead is again determined at different driving frequencies (see Fig. S9(a)) and finally the qp density is extracted It is possible to see in Fig. S9(b) that the values of n qp are similar to those of samples A and B in the same driving frequency range and, also the zero frequency limit is in agreement with the main measurements.In general, the n qp values are slightly lower than for samples A and B since in the reference device the superconducting gap ∆ 1 is slightly higher.In conclusion, the bare presence of the S 1 IS 2 SQUID does not affect significantly the excitation population in the narrow superconducting lead.The resistivity of aluminum was measured on separate samples using the four probe method, see Fig. S10 for the experimental set up and a typical sample.The test wires were fabricated using a process identical to the full samples.The measurements were done with four different wire geometries, all with a thickness of 70 nm and lateral dimensions 25 µm × 100 nm, 25 µm × 200 nm, 50 µm × 100 nm, and 50 µm × 200 nm, which are resemblant in shape to the long leads of the actual devices.Two batches of test samples were fabricated with differing aluminum deposition rates, 2 Å/s and 15 Å/s, respectively.The measurements were done at room temperature (∼ 298 K) and by immersing the sample into liquid nitrogen (∼ 77 K).The results of these procedures are shown in Table S1 as the mean value (ρ) and the standard deviation (σ) of a distribution of measurements in units of Ωnm.
Since the aluminum deposition rate for the measured SINIS SETs was ∼ 2 Å/s, we estimate that the resistivity of the long leads in the normal state at 77 K is around 31 Ωnm.Furthermore, immersing a few test wires in liquid helium (∼ 4.2 K) allows us to expect a roughly 10% decrease in resistivity (ρ 4.2 K ≈ 28 Ωnm) compared to the values measured at 77 K.  S1.Measured resistivity data from the four probe method in aluminum nanowires.The mean (ρ) and the standard deviation (σ) for a series of measurements done at different conditions (see the text) are shown, together with the number of samples employed.The resistivity of bulk aluminum at room temperature is 27.09Ωnm [14].

FIG. 1 .
FIG. 1. Qp extraction probed in a SINIS SET.(a) Electronic scheme for turnstile operation of the hybrid SET with an integrated S1IS2 cooler, together with a false color scanning electron microscope image of sample A. Light red marks the normal metal and blue the superconducting parts, with light (dark) blue designating the lower (higher) gap superconductor.(b) Schematic sketch of the device.Quasiparticles are injected to the narrow lead where relaxation is inefficient.Qp excitations are extracted via biased S1IS2 junctions, and recombination and scattering processes are still present along the lead.The inset shows that when the S1IS2 junction is biased at eVSIS = |∆1 − ∆2| the singularities in the qp densities of states match and excitations are transferred from the lead with ∆1 to the reservoir with ∆2.(c) IV curves of sample B in the positive bias regime, for a large number of different gate voltages.The red lines are simulations for the envelope curves, and the blue dots are the experimental data.(d) Measurements for sample A in the turnstile operation at f = 10 MHz, V b = 120, 160 and 200 µV as blue, red and yellow symbols, respectively.Here the current is measured against the amplitude of the gate signal Ag and the red dashed line illustrates the ideal value of I = ef for the first plateau.

FIG. 2 .
FIG. 2. Improvement of the I = ef pumping plateaus by the S1IS2 cooler.Data for the cooler-off (open diamonds) and cooler-on (filled triangles) cases compared with a model based on a master equation approach (solid black lines).Blue, red and yellow curves correspond to V b = 120, 160, 200 µV, respectively.(a) Sample A operating at 80 MHz; in the cooleron case the SQUID is biased at 60 µV.(b) Sample B operating at 100 MHz; here in the cooler-on cases the SQUID is biased to 50 µV.

FIG. 3 .
FIG. 3. Sensitivity of current plateau to qps.(a) Calculated curves showing how the current at the plateau varies with the qp density.The parameters of sample B were used in these calculations along with a bias V b = 100 µV and a gate amplitude Ag = 3.5 mV.(b) Measured current at the plateau as a function of the bias applied to the S1IS2 cooler for device A at f = 40 MHz.From bottom to top the curves correspond to V b = 120, 160, and 200 µV, respectively.(c) As in (b) for B at 20 MHz.(d) As in (c) for 80 MHz.Dashed black lines indicate the optimal biasing point of the cooler.

FIG. 4 .
FIG. 4. Qp density variation as a function of drive frequency.(a) Estimated narrow lead temperature during the turnstile operation as a function of the gate signal frequency for sample A and (b) for sample B. (c) Estimated qp density at the junction as a function of the pumping frequency for sample A and (d) B, corresponding to the data in panel (a) and (b), respectively.The black solid lines are fittings to Eq. (1).

6 FIG
FIG. S1.Normalized cooling power of a S1IS2 tunnel junction for different Dynes parameters.The junction parameters used in the calculation are given in the figure.In addition, to demonstrate practical magnitude of P(1) we set R = 5 kΩ, which is typical for our devices, to obtain the scale on the right hand side.The inset shows calculations in a wider voltage range showing the onset of strong heating.The bias points corresponding to ∆1, ∆2 and ∆1 + ∆2 are indicated by arrows.
FIG. S2.(a) Pumped current of sample A at a constant drive amplitude on the I = ef plateau.It is measured at f = 40 MHz, V b = 200 µV and plotted as a function of the bias applied to the S1IS2 cooler over a wide range.(b) Sample A SQUID current-voltage curve at Φ ≈ Φ0/2.Dashed vertical lines indicate the bias voltages eVSIS = |∆2 − ∆1|, ∆1, ∆2, and ∆2 + ∆1.
FIG. S4.Sample with different superconducting gaps (a) SEM micrograph of the measured sample with smaller difference |∆2 − ∆1|.Blue designates superconductor and red normal metal.Different tones of blue mean different gap superconductors.(b) IV curve in the DC regime.Blue dots are measured data and red lines are calculated curves from which device parameters are estimated.
FIG. S6.(a) Narrow lead temperature as a function of the gate signal frequency.(b) Estimated qp density as function of the gate signal frequency.
FIG. S8.Current pumping regime in the reference sample.Open diamonds are measured data and solid black lines are simulated.Blue, red and yellow correspond to bias voltages V b = 100, 140 and 180 µV, respectively.The gate signal has a frequency of (a) 10 MHz and (b) 40 MHz.

3 )
FIG. S9.Qp density variation as function of drive frequency.(a) Narrow lead temperature as a function of the gate signal frequency for the reference sample.(b) Estimated qp density as function of the gate signal frequency, corresponding to temperatures in panel (a).