Correlated transport through junction arrays in the small Josephson energy limit: incoherent Cooper-pairs and hot electrons

Linear arrays of Josephson junctions display surprisingly complex behaviour which belies their simple circuit diagrams. The earliest work on these circuits considered the strong, long-range interaction between charges within the array [1–5]. This Coulomb interaction between charges decays exponentially over a characteristic length that is set by the values of the circuit capacitances. During conduction, the repulsive interaction is counteracted by the voltage applied across the circuit, which 'pushes' the charges closer together. The equilibrium charge configuration is then reached when these two energies are balanced, resulting in a periodic charge distribution across the array. The movement of such a periodic charge state through the array in turn causes a periodic modulation of the current leaving the array. Through this process the statistical correlations in the charge distribution are converted into correlations in the current, i.e. correlated transport. The charges within a Josephson junction array (JJA) therefore display statistical correlations in both space and time, at the few electron level. The resulting correlated transport makes them a candidate circuit for applications in quantum metrology and current standards [6]. In contrast to the state of the art in charge-pump metrology circuits [7, 8], the use of JJA in metrology is currently limited by the precision of both the theoretical models and the experimental devices.

Although qualitative agreement between semi-classical theory and experiment has been seen in several devices [9–13], there are still many open questions involving hysteretic and thermal effects [14, 15], the role of quasiparticles and charge noise. The ultimate goal of a full quantum theory of these devices poses several challenges. Obtaining quantitative agreement at the level now possible in simpler Josephson devices [16–21] is difficult due to the interplay of charging energy, Josephson energy, charge noise and thermal effects. In addition, most experiments have focused on current-voltage (I-V) characteristics [6, 9] or microwave driven experiments [6, 22] which limits the amount of information available for comparison. An important exception are recent experiments [23–25] where the correlated transport of charges within a JJA was measured directly using a radio-frequency single-electron transistor (rf-SET). This provided detailed temporal information about the charge transport as well as the usual I-V characteristics. Of particular note is the observation of correlated transport for a range of magnetic fields, suggesting a smooth transition from normal- to super-current [24].

Using a relatively simple model, we demonstrate that correlated transport can be carried by incoherent Cooper-pair hopping in the limit of small ${{E}_{J}}$. We study the effect of equilibrium quasiparticles and incoherent Cooper-pair transport where the charge hopping is modelled using P(E) theory with the array's own intrinsic impedance (in contrast to [26] where an arbitrary, large impedance was assumed). Considering the effect of both temperature and magnetic field on conduction, we also show that the experimentally observed response is consistent with some level of self-heating of the array. This model opens the way forward for quantitative comparison between theory and future experiments and the development of a full microscopic model of transport in junction array circuits.

In the following section we derive our noise model due to the intrinsic array impedance and discuss details of the simulation method. In section 3 we consider the conduction through the array and show how correlated transport arises as a function of voltage. We then study the regime where Cooper-pair and quasiparticle transport coexists and results in a novel conduction regime (sections 4 and 5). Finally, in section 6 we consider the effect of magnetic field and the reduction of the superconducting gap on transport through the array.

Throughout this work, we consider a simple JJA circuit with a potential bias V applied across the array, see figure 1. Each Josephson junction has a Josephson energy ${{E}_{J}}$, an effective capacitance ${{C}_{J}}$ and a capacitance to ground ${{C}_{G}}$. The energy of the system is given by

$$\begin{eqnarray}&&H\left( \vec{q}\; \right)=\frac{1}{2}{{\overrightarrow{q\ }}^{T}}{{{\mathbf{C}}}^{-1}}\vec{q}+{{C}_{J}}V\vec{\delta }_{1}^{T}{{{\mathbf{C}}}^{-1}}\vec{q}\end{eqnarray} \tag{ 1 }$$

for a particular bias V and charge configuration $\vec{q}=e\vec{n}$. The source term ${{\vec{\delta }}_{1}}$ is equal to one for site 1 and zero everywhere else.

Zoom In Zoom Out Reset image size

In contrast to previous theoretical work on JJ arrays [14, 15], we take the conjugate phase to be across the ground capacitor ${{C}_{G}}$, which makes the charge variable ${{q}_{i}}$ simply the excess charge on the ith island, as is usually done when considering single-electron transistors. Although this choice has no additional observable consequences, it simplifies the implementation of the kinetic Monte–Carlo (KMC) algorithm for simulating the dynamics.

The capacitance matrix of the system is given by

$$\begin{eqnarray}{\mathbf{C}}=\left( \begin{array}{ccccccccccccccc} {{C}_{G}}+2{{C}_{J}} & -{{C}_{J}} & 0 & \cdots \\ -{{C}_{J}} & {{C}_{G}}+2{{C}_{J}} & -{{C}_{J}} & \cdots \\ 0 & -{{C}_{J}} & {{C}_{G}}+2{{C}_{J}} & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array} \right)\;,\end{eqnarray} \tag{ 2 }$$

where the inversion of this matrix can be computed very quickly numerically, but can also be expressed analytically [27–31], a fact we employ when deriving the interaction energy and threshold voltage. For this circuit, we define an interaction length $\Lambda$ between charges in the array, where

$$\begin{eqnarray}&&\frac{1}{\Lambda }={{{\rm cosh} }^{-1}}\left( 1+\frac{{{C}_{G}}}{2{{C}_{J}}} \right)\end{eqnarray} \tag{ 3 }$$

and in the limit of small ground capacitance $\Lambda \approx \sqrt{{{C}_{J}}/{{C}_{G}}}$.

To understand the origin of correlated transport in a JJA, we first derive the effective interaction energy between charges within an array. Taking the large array limit $N\gg\Lambda >1$, we find

$$\begin{eqnarray}&&U\left( {{n}_{i}},{{n}_{j}} \right)=\frac{1}{2{\rm sinh} \left( {{\Lambda }^{-1}} \right)}\left[ \frac{n_{i}^{2}}{2{{C}_{J}}}+\frac{n_{j}^{2}}{2{{C}_{J}}}+\frac{{{n}_{i}}{{n}_{j}}}{{{C}_{J}}}{{{\mbox{e}}}^{-|i-j|/\Lambda }} \right]\end{eqnarray} \tag{ 4 }$$

which is the energy for two charges at positions i and j in an infinite length, zero-biased array. In the limit of ${{C}_{G}}\ll{{C}_{J}}$, we can write ${{C}_{J}}{\rm sinh} {{\Lambda }^{-1}}\approx {{C}_{J}}/\Lambda \approx \sqrt{{{C}_{J}}{{C}_{G}}}$. The first two terms in equation (4) are the respective charging energies and the third term is the interaction energy between the charges, which decays exponentially over the interaction length $\Lambda$.

We can also compute the threshold voltage for quasiparticles or Cooper-pairs to be injected from the voltage source,

$$\begin{eqnarray}&&V_{{\mbox{th}}}^{\left( \kappa \right)}=\frac{\kappa {\mbox{e}}}{2{{C}_{J}}\left[ {\rm exp} \left( 1/\Lambda \right)-1 \right]}\end{eqnarray} \tag{ 5 }$$

where $\kappa =1,2$ for quasiparticles or Cooper-pairs respectively. Both the interaction energy and threshold voltage expressions are derived from pure energetic considerations. The details of the array conductance and correlated transport depend critically on the form of the hopping rates for charges within the array, which in turn depend on the charge carrier and the influence of temperature, magnetic field and the specific noise model.

We are specifically interested in the small ${{E}_{J}}$ limit so we assume that coherent oscillations are completely suppressed and that the evolution of the system can be described using P(E) theory [32–34]. Typically when computing the influence of noise on smaller circuits, the dominant source is the impedance of the leads and measurement circuitry. In the case of JJ arrays however, this leads to the conclusion that noise dominates at the edges of the array but within the array its influence decays exponentially over the length scale $\Lambda$. Although in principle this could lead to coherent behaviour deep within the array, we assume (as others have done [26]) that within the array, locally generated noise dominates. We specifically consider the noise generated by the internal impedance of the array itself, however additional contributions from background charge fluctuations [35] and two-state defects [16, 19, 20, 36–38] within both the junctions and the substrate may also play a role.

To derive the internal impedance seen by a particular junction we assume that every other junction within the array has an impedance ${{Z}_{J}}$ given by a capacitor with capacitance ${{C}_{J}}$ and a resistor with resistance ${{R}_{T}}$ in parallel. Generally we assume that the resistance of each junction is large as compared to the resistance quantum, although the explicit value of ${{R}_{T}}$ proves to be unimportant. We neglect the Josephson inductance in the limit of small ${{E}_{J}}$.

We can now formulate the effective circuit model (figure 1(b)) for the junction impedance ${{Z}_{J}}$ and the impedance of the rest of the array $Z\left( \omega \right)$. Using the theory of continued fractions [39] we can derive the array impedance $Z\left( \omega \right)$ for an arbitrary length array, however the resulting expression is unwieldy. To obtain a more useful expression, we derive an effective impedance by considering the resistive and capacitative response of the array individually. Computing the resistance of half an array as ${{R}_{{\mbox{HA}}}}=N{{R}_{T}}/2$ and the effective half-array capacitance ${{C}_{{\mbox{HA}}}}=\left[ {{C}_{G}}+\sqrt{{{C}_{G}}\left( {{C}_{G}}+4{{C}_{J}} \right)} \right]/2$, we then approximate the total impedance of the half-array as that of the resistive and capacitative components separately (see figure 1(c)). The total impedance seen by a junction is then given by ${{Z}_{t}}\left( \omega \right)={{\left[ i\omega {{C}_{J}}+Z{{\left( \omega \right)}^{-1}}/2 \right]}^{-1}}$ where

$$\begin{eqnarray}&&Z{{\left( \omega \right)}^{-1}}=i\omega {{C}_{{\mbox{HA}}}}+\frac{2}{N{{R}_{T}}}\end{eqnarray} \tag{ 6 }$$

is the contribution from each half-array (assuming the resistive and capacitive responses are decoupled). In the large impedance regime we consider (${{R}_{T}}\gg{{R}_{K}}=h/{{{\mbox{e}}}^{2}}$, ${{C}_{G}}\ll{{C}_{J}}$) this is a good approximation when compared to the (exact) continued fraction solution.

Integrating $\operatorname{Re}\left[ {{{\mbox{Z}}}_{t}}\left( \omega \right) \right]$ over all frequencies, we obtain the amplitude for a delta-function approximation to the array impedance,

$$\begin{eqnarray}&&\operatorname{Re}\left[ {{{\mbox{Z}}}_{t}}\left( \omega \right) \right]\approx \frac{2\sqrt{2}\pi \delta \left( \omega \right)}{\sqrt{\left( {{{\mbox{C}}}_{G}}+4{{{\mbox{C}}}_{J}} \right)\left( {{{\mbox{C}}}_{G}}+2{{{\mbox{C}}}_{J}}+\sqrt{{{{\mbox{C}}}_{G}}\left( {{{\mbox{C}}}_{G}}+4{{{\mbox{C}}}_{J}} \right)} \right)}}\end{eqnarray} \tag{ 7 }$$

which in the limit of ${{C}_{G}}\ll{{C}_{J}}$ gives $\operatorname{Re}\left[ {{{\mbox{Z}}}_{t}}\left( \omega \right) \right]\approx \left( \pi /{{{\mbox{C}}}_{J}} \right)\delta \left( \omega \right)$ as expected [33]. Approximating for short times, we can then evaluate the integral over the time correlation function [33] and obtain an expression for the P(E) function,

$$\begin{eqnarray}&&{{P}_{\kappa }}\left( E \right)=\frac{1}{\sqrt{4\pi {{\kappa }^{2}}{{E}_{c}}{{k}_{B}}T}}{\rm exp} \left[ \frac{-{{\left( E-{{\kappa }^{2}}{{E}_{c}} \right)}^{2}}}{4{{\kappa }^{2}}{{E}_{c}}{{k}_{B}}T} \right]\end{eqnarray} \tag{ 8 }$$

where the effective charging energy ${{E}_{c}}={{{\mbox{e}}}^{2}}/2{{C}_{A}}$ is now that of the array impedance,

$$\begin{eqnarray}&&{{C}_{A}}=\frac{1}{2\sqrt{2}}\sqrt{\left( {{C}_{G}}+4{{C}_{J}} \right)\left( {{C}_{G}}+2{{C}_{J}}+\sqrt{{{C}_{G}}\left( {{C}_{G}}+4{{C}_{J}} \right)} \right)}.\end{eqnarray} \tag{ 9 }$$

As this P(E) function describes the contribution due to noise generated within the array itself, this contributes to hopping rates for all the processes within the array. In principle, close to the edges of the array, this noise contribution will be modified by the external impedance but for simplicity we ignore this extra contribution, given that the length of the array considered is ten times the interaction length.

To study the transport of (correlated) current in the superconducting limit, we must consider the contribution of both Cooper-pairs and equilibrium quasiparticles. Although non-equilibrium quasiparticle distributions are possible and have been found to be an important contribution to loss processes in small JJ circuits [40–42], we ignore these contributions as the assumption of thermal equilibrium is adequate when considering the transport through JJ arrays at moderate effective electron temperatures. For similar reasons, we do not consider co-tunnelling processes in this analysis as the first order process dominates in all parameter regimes of interest.

The rate for an equilibrium quasiparticle to move between two charge states which differ by total energy $\delta E$ is given by

$$\begin{eqnarray}&&{{\Gamma }_{{\mbox{eqp}}}}\left( \delta E \right)=\frac{1}{{{{\mbox{e}}}^{2}}{{R}_{T}}}\mathop{\int }\limits_{-\infty }^{\infty }\int d\epsilon \;d\epsilon \prime \;\frac{\mathcal{N}\left( \epsilon \right)}{\mathcal{N}\left( 0 \right)}\frac{\mathcal{N}\left( \epsilon \prime +\delta E \right)}{\mathcal{N}\left( 0 \right)}f\left( \epsilon \right)\left[ 1-f\left( \epsilon \prime +\delta E \right) \right]{{P}_{1}}\left( \epsilon -\epsilon \prime \right),\end{eqnarray} \tag{ 10 }$$

where $\mathcal{N}\left( \epsilon \right)/\mathcal{N}\left( 0 \right)$ is the normalized BCS density of states, $f\left( \epsilon \right)$ is the Fermi function and the change in energy $\delta E$ also takes into account the chemical-potential shifts associated with charges entering or leaving the leads. This rate scales linearly with $1/{{R}_{T}}$ above a characteristic gap whose size is set by both the superconducting gap $\Delta$ and the P(E) function which is centred around ${{E}_{c}}$.

In the limit of incoherent Cooper-pair transitions, the Cooper-pair hopping rate is simply given by

$$\begin{eqnarray}&&{{\Gamma }_{{\mbox{cp}}}}\left( \delta E \right)=\frac{\pi }{2\hbar }E_{J}^{2}{{P}_{2}}\left( \delta E \right).\end{eqnarray} \tag{ 11 }$$

The asymmetry of this function provides a compelling argument as to why the noise should be dominated by the large internal impedance of the array itself. If one was to derive a hopping rate assuming vanishing impedance, this would result in a peak centred about zero and therefore approximately diffusive charge movement within the array. Whereas in the large impedance limit, the lossy transport derives directly from the ease with which the array can absorb energy from the individual charges. In addition, we find that the results obtained using the high-impedance environment are in better agreement with the experimentally observed current voltage characteristics [24, 43].

To simulate the time evolution of a JJA we use the KMC algorithm, whereby charge transitions are chosen stochastically based on the relative weight of the relevant transition rates [26, 44–47]. This results in an output consisting of successive charge configurations and the time between each transition, from which average quantities can be easily computed. For all results we use a combination of initialization and measurement phases in the Monte–Carlo simulation. We initialize the system in the 'empty' state, where there are zero charges on each island. The system is then time-evolved according to the KMC algorithm for sufficient time that the charge distribution within the JJ array has equilibrated (typically ${{10}^{5}}$–${{10}^{6}}$ Monte–Carlo steps). The quantities of interest are then computed over a further ${{10}^{6}}$ time steps such that we observe variances of order 1% for the average charge on a particular site within the 50 site array.

To study the correlated transport, we consider the charge–charge correlations within the array as this proves to be more computationally efficient than computing the current–current correlations directly. A key observable for JJA is the average charge on the jth site, $\left\langle{{n}_{j}}\right\rangle$. Computing the charge–charge autocorrelation function $\left\langle{{n}_{j}}\left( \tau \right){{n}_{j}}\left( 0 \right)\right\rangle$ on a particular site then gives a measure of the extent to which the charge distribution within the array is statistically correlated with itself at some later (or earlier) time τ. A further quantitive of interest is the Fourier transform of this correlation function, $\mathcal{F}\left[ \left\langle{{n}_{j}}\left( \tau \right){{n}_{j}}\left( 0 \right)\right\rangle \right]$, which gives a direct measure of the spectral response that would be observed in experiments [23–25, 44]. Although the Fourier transform of the charge–charge autocorrelation $\left\langle{{n}_{j}}\left( \tau \right){{n}_{j}}\left( 0 \right)\right\rangle$ can be computed directly [44], linearly sampling the resulting charge vector at a high bandwidth $>10$ GHz and then taking the autocorrelation and fast Fourier transform directly proves to be more efficient and less sensitive to numerical noise.

We are interested in the conduction processes of arrays in the superconducting regime, with ${{E}_{J}}\ll{{E}_{c}}$ and non-negligible superconducting gap. This means that transport with both possible charge carriers, Cooper-pairs and quasiparticles, is strongly suppressed and therefore the experimental observation of not only conduction, but correlated conduction in this limit is somewhat surprising. In addition, we find theoretically that over a large range of temperatures and voltages, correlated conduction in such circuits is largely suppressed which implies a relatively small part of parameter space in which these correlations can be observed. To focus on an experimentally relevant case, we consider a parameter regime (in terms of array length, ${{C}_{J}}$ and ${{C}_{G}}$) which corresponds approximately to that studied in [24].

A key feature of the parameters chosen in our analysis is that the incoherent Cooper-pair rate has its maximum ($4{{E}_{c}}$) at the edge of the low temperature equilibrium quasiparticle gap ($2\Delta$). Therefore we choose the charging energy to be ${{E}_{c}}=170\;\mu {\mbox{eV}}$ and the superconducting gap $\Delta =340\;\mu {\mbox{eV}}$, see figure 2. In this regime, the position of the P(E)-function within the quasiparticle gap, means that incoherent Cooper-pair hopping cannot be considered a resonant process but in fact must compete directly with the equilibrium quasiparticle rates, especially for large energy differences (larger than $\Delta$).

Zoom In Zoom Out Reset image size

Figure 2 illustrates the hopping rates for quasiparticles and Cooper-pairs as a function of energy difference at an effective electron temperature of T = 50 mK and T = 350 mK, with ${{E}_{J}}={{E}_{c}}/20=8.5$ μ eV. The effect of the P(E) function is to broaden the Cooper-pair rates as well as increasing the gapped region in the quasiparticle rates and smoothing the edge of this gap. The net result is that over a large range of energies both Cooper-pairs and quasiparticles can hop with approximately equal rates, especially at the typical transport voltages applied to this system.

Throughout our analysis, we model an array of N = 50 junctions with values of the junction capacitance ${{C}_{J}}=412.6$ aF and ground capacitance ${{C}_{G}}=25.9$ aF such that $\Lambda =4$. As we investigate the effect of varying ${{E}_{J}}$ while keeping the ratio $\Delta /{{E}_{c}}=2$, we consider junction resistances varying between ${{R}_{T}}=5{{R}_{K}}=129\;{\mbox{k}}\Omega$ and ${{R}_{T}}=125{{R}_{K}}=3.23\;{\mbox{M}}\Omega$ according to the standard relationship [48],

$$\begin{eqnarray}&&{{E}_{J}}=\left( \frac{\hbar }{2{\mbox{e}}{{R}_{T}}} \right)\left( \frac{\pi \Delta }{2{\mbox{e}}} \right){\rm tanh} \left( \frac{\Delta }{2{{k}_{B}}T} \right).\end{eqnarray} \tag{ 12 }$$

The absolute magnitude of the I-V characteristics of the array depends on both ${{E}_{J}}$ and temperature. For our chosen parameters we find that conduction is sporadic with a large variance for electron temperatures less than $T\approx 150$ mK as the system can become trapped in meta-stable charge states which do not easily decay due to the large quasiparticle gap. These trapping states are inconsistent with the observed reproducible and smooth response seen in experiments [24] and their absence implies an effective electron temperature which is higher than the nominal base temperature of the experiments.

A fundamental assumption of P(E) theory is that the energy associated with hopping of charges is dissipated within the circuit itself and in the quasiparticle degrees of freedom. It is this excess energy which can lead to higher effective electron temperatures, an effect which has been observed in qubit [40, 49, 50] and SET experiments [35]. The effective temperature seen by the charge carriers can depend in general on the bias conditions, junction properties, the superconducting gap and even the timescale over which the experiment is conducted [35, 51, 52]. For simplicity in the remainder of this analysis, we assume an electron temperature of T = 200 mK.

Figure 3 shows the I-V response at T = 200 mK as a function of Josephson energy, demonstrating a characteristic $E_{J}^{2}$ dependence (see inset to figure 3). This indicates that the dominant processes are Cooper-pair hopping events, equation (11). More interesting is the behaviour of the threshold for conduction, which goes through two distinct transitions. Initially conduction rises rapidly once the voltage is greater than the threshold for injection of quasiparticles is reached (although there is some thermal broadening of this threshold). At this point the simulations clearly predict that conduction is possible, however it is at extremely low currents which are well below the detection threshold of existing experiments. A second distinct transition to observable current levels is observed in the range 3–4 mV, which one could associate with the injection of Cooper-pairs. This however is not the case as the injection threshold for Cooper-pairs is only a factor of two larger than that for quasiparticles, which is well below the observed threshold. In fact, this second transition to a higher current state is associated with the interplay of quasiparticles and Cooper-pairs and the 'filling' of the array, as will become apparent in section 5.

Zoom In Zoom Out Reset image size

We now turn our attention to the correlated transport of charges through the array. To observe such correlations, the low voltage (close to threshold) regime is usually considered, such that the average charge density is low and the charges are separated by approximately $\Lambda$, see for example [25, 53]. This conventional correlated transport regime is illustrated in figures 4(a) and (b), where we plot the charge distribution within the array as a function of time for two different voltages within the 'low current' state. At very low voltage (just above threshold), the conduction is partially diffusive due to the relative high electron temperature (T = 200 mK). As the voltage increases (1 mV $\lesssim V\lesssim 4$ mV), the charge distributions become progressively more correlated and more closely mimics the 'conventional' correlated transport in normal arrays [53]. We see a characteristic periodic distribution of charges within the array a fixed point in time, due to the interplay between the applied voltage and the repulsive interaction between charges. As the charges move systematically along the array, this spatial correlation between the charges also manifests as a temporal correlation of the type observed in experiments [25].

Zoom In Zoom Out Reset image size

To more clearly see this temporal correlation between the charges, in figures 4(c) and (d) we show the charge occupancy autocorrelation function as a function of delay time τ. The damped oscillations illustrate that the charge occupancy on a particular site is correlated with occupancy at some characteristic time later. This characteristic time corresponds to the mean spacing between charges as they travel along the array. We see that the amplitude of the oscillations is stronger at large voltage biases as the charges are packed more tightly together, reducing the thermal 'jitter' observed just above threshold.

This temporal correlation also varies with position within the array, being strongest in the middle, weaker at the start (as the voltage source is uncorrelated) and non-existent at the end of the array (for this voltage bias configuration). This variation in correlation strength has been studied previous for different bias configurations [53] and arises due to two effects. First and foremost, for this bias configuration the average charge distribution is approximately linear across the array, approaching zero at the end which is connected to ground. As each charge spends a vanishing amount of time on the final site, the corresponding autocorrelation function does not show the temporal correlations. There is also a more subtle effect which arises at low voltage bias, due to the periodic distribution of charges being partially 'pinned' by the boundary conditions in the array. This results in a slight but periodic modulation of the average charge distribution. These effects have been studied in depth for the case of normal conduction in [53]. In this analysis, we are specifically interested in the interplay of normal and superconducting processes and therefore limit ourselves to the charge occupancy statistics measured at a fixed point in the middle of the array.

To investigate correlated transport at higher voltages, we plot the spectral response (Fourier transform of the charge–charge autocorrelation function) measured in the middle of the array, see figure 5. The peaks, corresponding to temporally correlated charge states, form almost immediately once the circuit transitions to the higher current state and are still observable up to at least V = 13 mV. The average charge density at these voltages is $\left\langle n\right\rangle \approx 1.5$ which is well above the regime usually considered. For comparison with earlier work [53], in this model the charge distribution has already reached a linear drop throughout the array well before the transition to high current.

Zoom In Zoom Out Reset image size

At these higher voltages, the correlations derive from a more complicated mechanism involving both single-charge (quasiparticle) and double-charge (Cooper-pair) excitations, which is not seen in a normal conducting array (without superconductivity). This additional transport mechanism is responsible for the transition to a higher current state for $V\gtrsim 5$ mV (see figure 3). Here the voltage is large enough to push the charge excitations closer together so that the interplay between quasiparticle and Cooper-pair tunneling, with the corresponding island charges, is energetically favourable and charges can flow freely again through the array.

The exact mechanism driving this interplay of odd and even charge states can be seen in figures 6 and 7. Odd charge states (n = 1) form which are relatively stationary, on the time scale of the current flow, and can sit on neighbouring sites due to the high bias voltage. These quasiparticle states form a quasi-static background charge distribution throughout the array. Transport then proceeds via movement of incoherent Cooper-pairs where the moving charge states are either $n=+2$ when there is a background charge of n = 0, or $n=-1$ when the background state is $n=+1$. Figure 7 illustrates the conduction pathway for this transport mechanism where the charge distribution across 10 sites within the array is depicted for a series of example Monte–Carlo steps. Although we see a combination of charge states $n=-1,0,1$ and 2, current is always carried by movement of Cooper-pairs from left to right.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The stability of this transport mechanisms relies on the high average charge density and the stationary quasiparticle states (due to the gapped hopping rate, see figure 2). It is only when both conditions are met, as a result of both high voltage and vanishing magnetic field, that we see this transport mechanism, where incoherent Cooper-pair hopping accounts all of the current flow through the array, once the system has equilibrated.

This new transport mechanism has a lower effective resistance as it is dominated by the Cooper-pair hopping rate ${{\Gamma }_{{\mbox{CP}}}}\propto E_{J}^{2}$ rather than simply the normal resistance ${{\Gamma }_{{\mbox{CP}}}}\propto 1/{{R}_{T}}$, resulting in a much higher current state than is seen at lower voltages. The interplay between $n=+2$ states and $n=-1$ states means that periodic spacing between current carriers need only be maintained for a few sites before conversion between charge states, which results in a more even distribution of correlations throughout the array. For even higher voltages ($V\gg 10$ mV), this effect is again washed out due to higher occupancy charge states and the correlations disappear.

This new transport regime also proves to be more robust in the presence of disorder, particularly in comparison to the zero-bias conductance regime [54]. Performing the same simulations with background disorder modelled as an initial random fractional charge state [55] allows us to mimic how the system would respond to random static variations in the local electric potential on each site. We initialize the system with charge states $n\in \left( -1,1 \right)$, we then evolve the system for enough time to reach equilibrium before collecting statistics. As the transport mechanisms considered only move charge from one site to another, the fractional offset due to the initial disorder is preserved modulo $1e$.

This method of modelling disorder results in no significant change in the I-V characteristics or correlation response in the high current regime. This is because the high current regime relies on the formation of a random distribution of static 'background' charges, the exact distribution of which varies from run to run. Including disorder simply selects one particular distribution (or a subset thereof) preferentially to minimize the total energy of the system. The transport mechanism proceeds in exactly the same fashion as without disorder, via incoherent Cooper-pair hopping through the background charge distribution. It should be noted that disorder does have a significant effect on the initial (quasiparticle dominated) low voltage regime, although as stated earlier, this regime is not experimentally resolvable for these parameters.

A key motivation for this analysis is the experimentally observed correlated transport as a function of magnetic field. At this point we have demonstrated that correlated transport can occur in the zero magnetic field, small ${{E}_{J}}$ limit. We now consider the characteristics of the current flow as a function of an applied magnetic field. We assume the Josephson energy displays the usual magnetic field dependence [48],

$$\begin{eqnarray}&&{{E}_{J}}\left( B \right)=\left( \frac{\hbar }{2e{{R}_{T}}} \right)\frac{\pi \Delta \left( B \right)}{2{\mbox{e}}}{\rm tanh} \left( \frac{\Delta \left( B \right)}{2{{k}_{B}}T} \right),\end{eqnarray} \tag{ 13 }$$

therefore both the superconducting gap and ${{E}_{J}}$ will contribute a magnetic field dependence to the conduction properties. As the magnetic field increases, the closing of the superconducting gap increases the rate of quasiparticle hopping while the Cooper-pair rate is dropping due to the reduction in ${{E}_{J}}$. For our parameters, this results in non-monotonic variation in the current for given voltage and the disappearance of the second (high current) transition for $\Delta \left( B \right)/\Delta \left( 0 \right)\lesssim 0.5$. In figure 8 we plot the relative fraction of Cooper-pair hopping compared to all hopping events as a function of both voltage and magnetic field. To avoid numerical noise from the finite simulation time, we only include data for currents greater than 0.5 aA. Looking at the relative contribution of Cooper-pair to quasiparticle events, we see that for $\Delta \left( B \right)/\Delta \left( 0 \right)\gtrsim 0.5$ Cooper-pairs dominate in the region of interest (4 mV $\lesssim V\lesssim 15\;{\mbox{mV}}$). For lower voltages, quasiparticle transport plays a more important role as the charge density within the array is closer to zero. Quasiparticles also dominate at large magnetic fields (small superconducting gaps) as one would expect.

Zoom In Zoom Out Reset image size

We now consider the spectral response of the charge–charge correlations at fixed current (I = 1 pA) as a function of superconducting gap size. Figure 9 illustrates the spectral response as a function of gap size where correlated transport is observed in both the small and large gap limits, although the correlations disappear in the intermediate regime. The frequencies where we observe a peak in the spectral response in the low and high gap regimes differ by a factor of two, due to the dominant charge carriers being either quasiparticles ($Q=1e$) or Cooper-pairs ($Q=2e$) where Q defines the magnitude of the effective charge carrier. This is consistent with magnetic field dependence of the relative fraction of Cooper-pairs plotted in figure 8 as well as the experimentally observed relationship between frequency response and current in [24].

Zoom In Zoom Out Reset image size

Although the distribution of single occupancy excitations in the high current regime is approximately static on the time scale of the current flow, long time fluctuations are observed which correspond to random rearrangement of the quasiparticle background. These slow fluctuations can be seen as fluctuations of the current response in the high current transition regime in figure 3. This slow background charge rearrangement leads to a non-negligible zero frequency component at low magnetic field, see figure 9. As the magnetic field is increased, this contribution to the spectra broadens as the hopping rate for quasiparticles increases. This transition regime corresponds to where the high current transport regime is no longer correlated due to strong fluctuations of the quasi-static background but the overall charge density is still too high to observe conventional (normal) correlated conduction. The exact details of this transition depend strongly on the characteristics of the array (effective impedance, temperature, Josephson energy). At sufficiently high magnetic field, conduction is simply dominated by quasiparticles and we see a single frequency peak corresponding to correlated transport.

This model reproduces the observation of correlated transport in both the high and low magnetic field limits [24], however several important questions remain. A key experimental observation is the continuous transition from $Q=1e$ to $Q=2e$ which is not reproduced in our simulations. Although broadening due to the finite response of the rf-SET would explain the continuous transition, this would still require coexistence of both the $Q=1e$ to $Q=2e$ dominated phases which is not observed here. However, given that both background disorder and the bias position of the rf-SET itself can stabilize particular quasi-static charge configurations, a plausible explanation is that the system fluctuates between both transport states, the result of which is then time averaged during rf-SET detection. Secondly, we have not considered the effects of variations in the effective electron temperature and in the impedance seen by the charges as a function of time and/or position within the array.

We have developed a model for transport due to both incoherent Cooper-pairs and equilibrium quasiparticles which explicitly includes the ability of the array impedance itself to dissipate energy. Such a model fundamentally modifies the response of the charge carriers and leads to qualitatively different predictions for the array response which help to explain several experimentally observed effects. We observe correlated transport due to incoherent Cooper-pair hopping which can occur at higher temperatures and higher voltages than one would expect for normal (correlated) conduction. This new correlated transport phase is carried via movement of double excitations (Cooper-pairs) through a random distribution of stationary single excitations (quasiparticles). The breakdown of this phase as a function of applied magnetic field is consistent with the observation of correlated transport in both the low and high field regimes, even for small Josephson energy.

The authors wish to acknowledge useful discussions with A Shnirman, T Duty, P Delsing, K Walker and N Vogt. This work was supported by the Victorian Partnership for Advanced Computing (VPAC). We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of Karlsruhe Institute of Technology.