Tunable quasiparticle trapping in Meissner and vortex states of mesoscopic superconductors

Nowadays, superconductors serve in numerous applications, from high-field magnets to ultrasensitive detectors of radiation. Mesoscopic superconducting devices, referring to those with nanoscale dimensions, are in a special position as they are easily driven out of equilibrium under typical operating conditions. The out-of-equilibrium superconductors are characterized by non-equilibrium quasiparticles. These extra excitations can compromise the performance of mesoscopic devices by introducing, for example, leakage currents or decreased coherence time in quantum devices. By applying an external magnetic field, one can conveniently suppress or redistribute the population of excess quasiparticles. In this article, we present an experimental demonstration and a theoretical analysis of such effective control of quasiparticles, resulting in electron cooling both in the Meissner and vortex states of a mesoscopic superconductor. We introduce a theoretical model of quasiparticle dynamics, which is in quantitative agreement with the experimental data.


SUPPLEMENTARY NOTE 1. HYSTERESIS IN THE MEASUREMENTS UNDER FIELD
As pointed out in the main text, a remanent field of δH ≈ 2.5 mT and the asymmetry in the Meissner state are present at the sample location at zero applied field H = 0 due to the presence of superconducting parts in the sample-holder. The dc measurement presented Fig. 2(d) of the main text for Sample B has been reproduced (on the same sample) in a sample-holder that does not have the superconducting shield for which B = H. These measurements are shown in Supplementary Figure 1 for several bias current values I bias = 1 (blue triangles), 10 (red circles), and 100 pA (black squares), together with the theoretical model for the experimental data (solid lines of corresponding colors). The dc voltage V vs field B in this measurement is symmetric with respect to the zero applied field value (except for the vortex hysteresis intrinsic for the sample). The theoretical model presented in Supplementary Note 3 reproduces perfectly the experimental points.
Supplementary Figure 1. DC measurements without field distortion. Evolution of the voltage at a fixed bias currents I bias = 1 (blue triangles), 10 (red circles), and 100 pA (black squares) at the gate voltage C g V g /e = n g = 0.5 suppressing the Coulomb energy with the magnetic field for Sample B in a sample-holder that does not have the superconducting shield together with the theoretical model (solid lines of corresponding colors). The field B is swept from -30 mT to 30 mT, the direction is shown be the arrow.
In order to fit the theoretical model to our measurements versus field H performed in the sample-holder with the superconducting shield (i.e. with deformation of the field profile), a correction has to be done, as the field H applied through the coil differs from the effective "acting" field B seen by the sample. The field profile has been measured in the shielded sample-holder with a Hall sensor at 4.2 K and at 0.2 K. The first measurement has been done in the normal state when there is no magnetic shielding as a reference point (not shown). The second measurement of the effective field B versus the applied coil magnetic field H swept from -30 mT to 30 mT (at 0.2 K) is shown in Supplementary Figure 2 as black dashed lines. The arrows point out the direction of the sweep. For |H| > 20 mT, the sample-holder is fully normal and the effective field equals the applied one. For |H| < 20 mT, a nonlinear superconducting response from the sample holder is present leading to hysteresis. The red line in Supplementary Figure 2  (2)

SUPPLEMENTARY NOTE 2. HOMOGENEOUS APPROXIMATION IN A MESOSCOPIC SAMPLE
The samples studied experimentally are in the dirty regime, namely ≪ , where = ℏ /Δ the coherence length, D is the diffusion coefficient, l is the elastic mean free path, and ∆ 0 is the superconducting gap. The quasiparticle (QP) spectral characteristics in this case can be found from the Usadel equations (see, e.g., Ref. [1]) where g R = −g R* = cos θ(r) and f R = −f A* = −i sin θ(r) are the normal and anomalous Green's functions (superscripts 'R' and 'A' stand for 'retarded' and 'advanced').
Considering the experimental situation of a mesoscopic superconducting sample with the characteristic size R = 0.5 µm and the coherence length ξ ∼ 100−200 nm (depending on the diffusion coefficient), we have to verify if we can neglect the gradient terms in Sup. Eq. (1). For subgap energies the characteristic length scale of the function θ(r) can be estimated as follows: / 1 − /Δ . It is natural to assume the θ(r) function inhomogeneity to be small provided / 1 − /Δ > . This condition gives us the energy interval 1 − /Δ < / ~ 0.1, sufficient for the calculations of the electron-phonon heat flow and the thermal excitation leakage current δI for temperatures much lower than the superconducting gap. Indeed, the main contribution to and to δI is given by | /Δ − 1|~ /Δ ≃ 0.03 < / ≃ 0.1.
Certainly the above assumption is strictly valid only for the Meissner state: inside the vortex core the gap and the anomalous Green function turn to zero at the scale of the effective core radius rv which is of order of the coherence length ξ [2].
To avoid numerical solution of the Usadel equation we adopt in the main text the following approximate procedure. In the presence of vortices we assume that both the order parameter and θ(r) function vanish inside the vortex cores while outside the core regions we assume the θ(r) function to vary slowly and introduce, thus, its average θ over the region outside the vortex cores (omitting the spatial dependence in the notation). The deviations from the averaged order parameter ∆ beyond the cores also become small in this limit. Integrating now the above Usadel equation over the region outside the vortex cores we obtain Eq. (1) from the main text with the effective depairing parameter expressed through the superfluid velocity vs as Here the brackets 〈. . 〉 denote an average over the sample volume (over the central part of the sample A) with the excluded vortex core regions, φ is the superconducting order parameter phase and A is the vector potential determined by the external magnetic field B applied to the sample. The second Usadel equation in our approximation reduces to div v s = 0 and leads to the vanishing components of vs perpendicular to the sample boundary and to the boundaries of vortex cores. Here and further on we neglect the changes in the magnetic field B due to the screening currents flowing in the sample because of the smallness of the characteristic sample size R as compared to the effective screening length λ eff = λ 2 /d S . For our samples ≃ 230 nm [4] and d S = 20 nm, therefore ≃ 2.6 µm. Solution of the averaged Usadel equation gives us the expression for the hard gap E g in the density of states and for the order parameter ∆ as functions of Γ as 1,5-7 In the main text we focus on the case γ < 1 (Γ < ∆ 0 e −π/4 ), implying that the gap E g > 0 is non-zero.

SUPPLEMENTARY NOTE 3. DC FITTING
Using the solution of the averaged Usadel equation, Eq. (1) from the main text, one can fit the IV characteristics shown in Supplementary Figure 1. Indeed, we consider a hybrid single electron transistor (SET), namely, a mesoscopic superconducting island tunnel coupled to the normal metal leads (NISIN). We apply a fixed bias current I bias through the normal leads and the constant gate voltage n g = C g V g /e = 0.5 to the gate electrode coupled to the island through the capacitor C g (see black and red lines in Fig. 2(c) of the main text) and measure the difference V of voltages V L,R = ±V/2 applied to the leads as a function of the magnetic field B seen by the sample.
In stationary state the current I bias flowing from one lead to another is equal in any cross section and it can be calculated in any of two junctions (for example, in the left one) as a sum over the island charge state k of the sequential tunneling rates Γ → (Γ → ) to (from) the island through the left junction. This sum is weighted with the probability p k of system being in this charge state, which is calculated using the standard rate equation for the balance of the probability fluxes 8-10 in the stationary case / = 0 with the tunneling rate Γ → ± = ∑ Γ → ± , and Γ → ± = Γ , ± given by Here , ± = ∓2 ( − ± 1/2) ∓ are the energies gained by the electrons tunneling to/from the island (being in the charge state k) through ith junction, R T /2 is the tunnel resistance of each junction. Here we focus on the magnetic field effects in the sample B (see Fig. 2 By fitting V(B) at I bias = 100 pA which is not affected by the vortex tail contributions one can extract the following values of fitting parameters α 1 = 0.38, α 2 = 0.438, and α 3 = 0.266 mentioned in the main text. Following [11] we attribute to all jumps in this plot with the change of the number of vortices in the sample and use the point of the first jump at B > 0 as the field of the first vortex entry B C = 14.4 mT. In this setup we don't see any transitions between vortex configurations with the constant vorticity like the transition to a giant vortex state (see, e.g., [12,13]). Using these parameters one can fit V(B) quite well at all bias current values with R/r v ∼ 1.7. The optimal value of the vortex core radius r v = 2.5 − 2.7ξ extracted from dc measurements in the Sample B is in perfect agreement with the previous theoretical works 2,3 . In subgap regime I bias = 1 and 10 pA the jump-like anomalies in V(B) become knee-like, but because of the abovementioned reasons we still associate each of them with the vortex entry or exit.

SUPPLEMENTARY NOTE 4. ELECTRONIC PUMPING
The electronic pumping of the Sample B at f = 5 MHz, when V bias ≈ 120 µV and n g ∼ 0.5, is shown in Supplementary Figure 3(a) with the field H swept from −10 mT to 2 mT. Contrary to what is observed in the Sample A, the increase of the magnetic fie ld increases the deviation from the current quantization I = ef , due to the effect of the screening current on the superconducting gap. This observation is in agreement with the theoretical model with the increasing number of QPs in S island with the field.
The electronic pumping of the Sample A at frequency f = 200 MHz, when V bias ≈ 250 µV and n g ∼ 0.5, is shown in Supplementary Figure 3

SUPPLEMENTARY NOTE 5. HEAT BALANCE EQUATION
In this section we describe the theoretical model of the relaxation of QPs in applied magnetic field by using the example of NISIN SET in the turnstile regime. By applying the constant bias voltages V L,R = ±V bias /2 to the normal leads and the periodic gate voltage n g (t) = C g V g (t)/e = n g 0 + A g sin(2πft) with a certain offset n g 0 , frequency f , and the amplitude A g to the gate electrode one can push electrons to tunnel through the system producing a time-dependent current I(t). This transport current I(t) flowing from one lead to another drives the NISIN turnstile out of the equilibrium by injecting nonequilibrium QPs into the S island. The power injected to the island increases with the frequency f and we model this increase in mean density of QPs in the superconductor by raising its electron temperature T relatively to the phonon bath temperature T 0 . Note that the quasiequilibrium Fermi distribution of electrons over energy f T (E) = [e E/k B T + 1] -1 is provided by the smallness of the inelastic electron-electron scattering time τ ee comparing to the operating time τ 0 = 1/f and the effective charging time e/I. 14 Due to the large electron-phonon relaxation length ≫ compared to the island size R we consider the heat balance equation 4  Within the optimal conditions of the proper turnstile shielding and optimized device geometry the averaged current I is close to its ideal value ef and the deviation δI = I − ef is mainly governed by nonequilibrium QP density in the S island near the junction The estimates for contributions in higher orders in small parameter ħ/e 2 R T are given in the following Note.

SUPPLEMENTARY NOTE 6. EXCESS CURRENT AS A FUNCTION OF ELECTRONIC TEMPERATURE
To calculate our main observable, the leakage current δI = I−ef in the NISIN turnstile we use the simplified version of the master equation given in Sup. Eq. (5)  . This expression contains the exponentially growing part with U which determines the dominant tunneling rate with maximal U for each time instant. We consider the offset n g 0 = 0.5 for simplicity and use the symmetry of the drive n g (τ 0 −t) = 1 − n g (t) focusing on the first half of the period with n g increasing from 0.5 − A g to 0.5 + A g . We assume that before the time instant t 1 the island is discharged k = 0 due to the domination rate Γ → among the others and the charging process is started at t = t 1 . The probability P 0 (t) to stay in the state k = 0 is decreasing with time t > t 1 as For typical frequencies the charging process occurs not far from n g = 0.5, therefore further we linearize the drive n g (t) ≈ 0.5 + 2πA g (ft− 1/4). As the island has been charged ( * ) = ≲ 1 (let's take ϵ = 1/2 for definiteness) the leakage current starts to flow. The number of excess electrons N l through the island can be written as the integral of the largest subleading rate Γ → ( ) governing the leakage current over the time interval t * < t < t 2 before this rate becomes the dominant one The leakage current can be calculated as follows ≃ 2 , where "2" accounts for the leakage during the second half of the period due to the symmetry k ↔ 1-k and L ↔ R.
for n g (t 1 ) = n g (t 2 ) = 0.5 at low enough electronic temperature ≲ 1 − ( /2 ) and ( ) = | | /2 − (1 − /T) + ln for n g (t 1 ) =1-n g (t 2 ) < 0.5 in the opposite case ≳ 1 − ( /2 ) . In this derivation we consider the operating frequency f to be small compared to the charging rate ~ (t * -t 1 ) -1 to avoid missing events. We neglect the relative corrections of order of | |/ (Γ → /Γ → and Γ → /Γ → for the case when the first term in Sup. Eq. (18) dominates for all rates). We don't take into account the factor ½ in N l during the time when Γ → ≫ (22) Γ → = Γ → ≃ Γ / , when the discharging occurs with the equal probability p L,R = 1/2 to left and to the right contact. We can do it, because during the integration of Sup. Eq. (20) , can go beyond the subgap range , < − suppressing the second term in Sup. Eq. (18) for Γ → exponentially ~ ( | , |)/ and keeping the rate Γ → to be the dominant one in the leakage current.
To avoid all these unimportant details we consider a certain A g -dependent numerical prefactor C ∼ 1 instead of the square brackets in Sup. Eq. (21) and come to Eq. (5) of the main text by using the assumption that the S gap near the junction (in the sample A) is close to ∆ 0 . Comparing Sup. Eq. (11) and Eq. (5) in the main text one can write down the following relation between the leakage current δI = I − ef and the QP density n qp near the junctions = ( ) / used in Fig. 4 of the main text to show the QP density scale.
Note that we neglect also the contributions of higher orders in the small parameter ħ/(e 2 R T ) like Andreev tunneling (see, e.g., [17,18]) due to rather large tunnel resistance of the sample contacts. Indeed, from the experimental side the attribute feature of Andreev tunneling is the additional peak in the beginning of each current plateau I = nef 17 which is not observed in all pumping measurements of this paper. From the theoretical side one can estimate the relative contribution δI AR /(ef ) of Andreev tunneling to the current as the ratio Γ AR /Γ[U ] of dc rates of sequential Γ[U ] = U/(e 2 R T ) and Andreev tunneling Γ ≃ ℏ /(4 ) in the above-gap regime. Here N = A/A ch is the number of channels in the tunnel junction, ≃ 6 • 10 nm is the area of the junction and A ch is the area of a single channel. Theoretical estimates given in [18] lead to A ch ∼ 2 nm 2 , while experimental observation 17 gives A ch ∼ 30 nm 2 . The upper bound estimate with A ch ∼ 30 nm 2 and R T = 577 kΩ for the sample A gives N ∼ 200 and / ( ) ~ ℏ/(4 ) ≃ 3. 10 which can be neglected comparing to the QP contribution.