Charge filling factors in clean and disordered arrays of tunnel junctions

We simulate one-dimensional arrays of tunnel junctions using the kinetic Monte Carlo method to study charge filling behaviour in the large charging energy limit. By applying a small fixed voltage bias and varying the offset voltage, we investigate this behaviour in clean and disordered arrays (both weak and strong disorder effects). The offset voltage dependent modulation of the current is highly sensitive to background charge disorder and exhibits substantial variation depending on the strength of the disorder. We show that while small fractional charge filling factors are likely to be washed out in experimental devices due to strong background charge disorder, larger factors may be observable.

Scientific RepoRts | 5:17572 | DOI: 10.1038/srep17572 long devices maximal disorder is expected to completely suppress the U-dependence of the transport properties by self-averaging. The devices fabricated so far however are too short for full self-averaging as seen in standard condensed matter systems. It is therefore important to study the U-dependence in short arrays even in the case of maximal disorder. Improvements in fabrication technologies may reduce disorder in these devices. Therefore, it is important to understand the effects of varying levels of disorder and to which degree disorder would have to be reduced to see real improvements in the behaviour of these devices. We investigate charge filling factors by simulating a junction array model in the charging energy limit (negligible Josephson energy) with zero disorder. We then extend this model to include various strengths of background charge disorder and consider the effect disorder has on the charge filling structure. We specifically focus on the signatures of filling factors and background disorder in the current-voltage characteristics of the device.

Clean array
We begin by considering an ideal model wherein the array is clean and homogeneous -there is no background charge disorder. We consider a one-dimensional linear array consisting of N = 50 islands each with a junction capacitance C J = 50 aF, an effective capacitance to ground C G = 2 aF and an effective electron temperature of T e = 30 mK. The array is driven by a symmetric voltage bias consisting of a small fixed bias V and varying offset voltage U, see Fig. 2. The current-voltage characteristics are non-monotonic functions of U due to charge filling factors within the array. We consider the normal conducting regime where the current is comprised entirely of electrons due to single electron tunnelling. The charge interaction length can be approximated by Λ = / C C J G , where Λ = 5 (or Λ = 15, where specified) for our simulations. The Hamiltonian consists of a charge interaction term and two source terms and is given by for a particular charge configuration  q. The two source terms δ  1 and δ N are equal to one for site 1 and N, respectively and zero everywhere else. The electrostatic interactions of the array are described by the N × N capacitance matrix where δE is the energy difference between initial and final charge configurations and f(E) is the Fermi function. In general, the influence of the circuit impedance on the transition rate can be included via P(E)-theory 18 . However, as C G /C J ≪ 1, we can consider the junction array to be a low impedance environment 19 . In this limit, P(E) ≈ δ(E), resulting in the expression given in Eq. 3 for the transition rates. These rates are used to select which transition occurs, based on a weighted probability. The total current is then computed based on the net movement, left and right, of charges over many Monte Carlo steps. To ensure the array has equilibrated, Monte Carlo steps of 10 5 -10 6 are used before collecting statistics with an additional 10 5 -10 6 steps and ensemble averaging over multiple data sets when considering disorder. For a detailed review of the KMC algorithm, we refer the reader to ref. 20 or ref. 21.
We can also calculate the analytical threshold voltage for a current to flow through the array (at T e = 0) 16,22 . In the limit of long interaction length (Λ ≫ 1) and assuming an initially empty array where e is the elementary charge. The threshold voltage and the current response are highly variable, depending on the characteristic value of the offset voltage U which sets the filling fractions within the array, as shown in Fig. 1. This behaviour is periodic in U, which is the focus of this paper, and is shown for half a period via the magenta line in Fig. 1.
The U-dependence can also be seen in Fig. 2 where we plot the current response as a function of U and V. The current is a periodic function of the offset voltage U for all simulated values of V, where the period of the current in U ∝ e/C G (the effective voltage on the ground capacitor). Therefore, the offset voltage required for an integer filled array (point of zero current) is . Due to particle-hole symmetry, the response is symmetric about the point U = e/2C G . The fine structure, seen at smaller voltages, is lost as V is increased and the response flattens out and becomes more rounded. It is this fine structure that is a signature of charge filling fractions.
In order to identify filling factors (p/q), we compute the time averaged charge-charge autocorrelation function 〈 n i (0)n i+m (0)〉 to measure the statistical correlations within the charge distribution. As the length of the array is not infinite, boundary effects are significant and the periodicity of the filling factors is strongest at the centre of the array.
From the current modulation in U, we can identify charge correlation patterns, i.e., peaks, troughs and dislocations in the patterns, as shown in Fig. 3. We can understand this structure by considering the filling factors of the array. Current blockade occurs at integer filling, i.e., the array is uniformly filled and it is difficult to move this pattern forward by injecting another charge or otherwise disrupting the pattern. In other words, integer filling factors correspond to states in which the charge configuration is rigid, therefore, the system requires significantly more energy to break this pattern and allow conduction to commence.
Fractional filling, on the other hand, allows current to flow more easily due to the presence of unoccupied sites between current carrying charges. For a perfectly periodic charge distribution that matches the boundary conditions of the array, we also see quenching. However, when the periodicity of the charges is disrupted, either via mismatch with the boundary or (near) degeneracy of two different charge configurations, current can flow more easily via interconversion of these charge states. As an example, periodic p/q filling fractions (e.g., 1/3 filling) are more difficult to move forward, resulting in current quenching. The combination of two different filling fractions (e.g., 1/3 and 1/2 filling) is easier to transfer through the array because the two factors can be interchanged. Therefore, as a function of U, we see current maxima and minima as aperiodic and periodic charge patterns, respectively. In Fig. 3, this effect can clearly be seen for the case of Λ = 15.
However, when the average spacing between charges approaches the interaction length Λ , we see a greater mixing of the filling fractions and the troughs in the current magnitude shift. As p/q → 1/Λ , we see peaks rather than troughs at rational values of p/q. This can be more clearly seen in Fig. 3 for the case of Λ = 5. When the mean charge separation is < 1/Λ , for example, 1/2, mixing of the charge configurations is weaker. Mixing at small Λ = 5 can also be seen in the autocorrelation functions in Fig. 3, where the fractions are no longer clear and become blurred, in contrast with the clarity of the same filling factors in the autocorrelation for large Λ = 15.
We observe splittings associated with defects in the dominant patterns, for instance, at Λ = 5, 1/2 filling in Fig. 3. The way in which these patterns match against the boundary conditions splits the equivalency of the different degenerate configurations. For example, at Λ = 5, 1/2 filling, we see that slightly above or below p/q = 1/2 gives a slightly lower current (the trough-peak-trough structure). This feature is not seen at Λ = 15, p/q = 1/2 because at large interaction lengths averaging over a greater number of sites washes out this effect. This is controlled by how easy it is to break the pattern, i.e., by adding an additional electron (or lack thereof) to make the pattern shift by half a period. As an example, the (metastable) state =  q 10101011010101 T is slightly harder to move along the array than the 'perfect' pattern of =  q 10101010101010 T . Therefore, depending on the matching of the boundary conditions and the presence of defects changes whether we see a peak or a trough.

Effects of disorder
Experimentally, imperfections or a degree of disorder is always prevalent in the system (e.g., in the substrate and within the junctions themselves), therefore, the inclusion of disorder makes our analysis more applicable to experimental observations. Furthermore, an important question is which features of this analysis are preserved in an inhomogeneous array. In this section, we extend our model by including various background charge disorder and investigating the weak and maximal disorder limits. We simulate boxed disorder (a uniform distribution) by adding a random fractional charge to each site, q ∈ [− ηe, + ηe]. The maximum/minimum value of disorder on any given site is given by ± ηe, where η indicates the width of the disorder distribution. This type of potential disorder simulates static background charges and is often used to study mesoscopic systems 9,12,23 .
The effect of the background on the current-voltage response can be seen even at nominal disorder (e.g., η = 0.01). While the structure remains largely unchanged, the current magnitude steadily decreases, as shown in Fig. 4. As the magnitude of static offset charges increases, the response becomes noisier and eventually only the most prominent structure remains (e.g., integer and 1/2 integer filling). However, even at very small current magnitude (strong disorder) periodic structure in U remains. The current scales as η because disorder constricts the conduction channel, making it more difficult for charge to flow (as opposed to a clean array). At strong disorder, a larger V is required to overcome the random static background charges and force current through the array.
Maximal disorder should be reached at η = 0.5 as this value of disorder corresponds to the point where any value larger than 0.5 can be negated by one additional charge tunnelling onto the island. Indeed, at η ≈ 0.5, our model shows maximal disorder and all of the structure is lost and charge filling factors are no longer seen. In Fig. 5, we show that the current is approximately constant across U at η = 0.5, indicating that maximal disorder has been reached (for this value of V) and a lack of U-dependence. Then at η = 0.6-0.9, the current oscillates again, with opposite curvature (on average) to the η ≤ 0.4 currents, with the current once again becoming approximately constant at η = 1. In order to understand this behaviour, we consider three different values of disorder; less than, equal to and greater than maximal disorder, depicted in Fig. 6. When the disorder distribution is broader than maximal disorder (η = 0.5), adding an individual whole charge causes the distribution to wrap back around (or fold back) into the interval ± 0.5. This folded disorder distribution can be mapped back to the original box-distribution form with an additional 'base' contribution spanning the whole [− 0.5, 0.5] interval by a shift of 0.5e (corresponding to the opposite current curvature observed in Fig. 5).
In Fig. 7, we plot the current response as a function of η for three different values of Λ at both the maximal and minimal I-U points. The applied voltage is scaled such that V = 0.8V th for each value of Λ , where V th is calculated via Eq. 4. In the clean array limit, the maximum current depends on Λ because the current response in the transport regime is not exactly Ohmic, but shows a linear current-voltage dependence with a non-zero current offset. The maximum/minimum currents as a function of U for all studied values of Λ equilibrate at approximately η = 0.5, i.e., maximal disorder. This is because disorder suppresses the current variations as a function of U, seen in Fig. 5. The oscillatory behaviour of the current at 0.5 < η < 1 can also be seen, with the currents at η = 0.5 and η = 1 equal for each value of Λ due to maximal disorder folding. Note that throughout this discussion Λ < N.
To better understand the role of the background charge disorder, we introduce the concept of a disorder correlation length, the distance over which charges on adjacent sites are still correlated. Beyond this length, the background disorder has the effect of randomising the electrostatic environment that the charges see and therefore, islands separated by more than this length can be considered independent. Assuming that the cumulative effect of charge disorder can be thought of as a random walk style process 24 , we define the disorder correlation length as where η is the value of charge disorder. When η ≥ 0.5, L corr ≤ 1, which corresponds to a delta-correlated disorder since the minimal length-scale of the system is the inter-site-distance L = 1. Such a delta-correlated disorder corresponds to the maximal disorder model 13,24,25 . Conversely, η → 0 corresponds to L corr → ∞, i.e., the clean array limit. When V is constant, L corr = 20, 35, 50 have near identical I-U characteristics as a function of N (data not shown). In contrast to the N-dependence of the threshold voltage 26 , the system exhibits similar behaviour whether L corr is equal to, less or greater than N. At N = L corr , the onset of transport changes from a boundary to a bulk effect and consequently, the dependence of V th on N changes. Here, however, Points of maximal and minimal current at zero disorder are shown (dashed lines). The system reaches maximal disorder at η = 0.5. When 0.5 < η < 1, the current exhibits oscillatory behaviour, however, the system reaches maximal disorder again at η = 1. The mean current across U and the variance in the disorder realisations for 0.1 ≤ η ≤ 1 is shown via error bars which depict one standard deviation. For very broad distributions, η ≫ 0.5, the η-dependence will abate due to decreasing 'base' to 'top' ratio.
we are considering a tunnelling array that is already in the transport regime and the modulation of the current due to a change in filling factors is naturally a bulk effect.

Discussion
In conclusion, we have studied charge filling factors as they manifest in a (unrealistic) clean tunnel junction array. These effects are predicted to be difficult to observe experimentally as the intrinsic disorder in experimental devices likely approaches or exceeds the limit for their observation. However, while small fractional filling factors are expected to be unobservable in present day devices, our model predicts that observation and reproducibility of larger factors may be possible (excluding maximally disordered arrays). For example, current blockade near integer filling should be quite robust to static offset charges. Signatures of charge filling factors could also be observable in other one-dimensional systems. The clean array analysis is applicable to arrays of quantum phase slip elements 27 where the effects of localised flux disorder should be much weaker and should display longer disorder correlation lengths. The maximal disorder analysis is more applicable to Josephson junction arrays, where background charge offsets dominate the behaviour.