Coulomb drag and depinning in bilinear Josephson junction arrays

When two systems with long-ranged Coulomb interactions are placed in close proximity, applying a voltage bias to only one of these systems can produce a current in both, even when there is no direct transfer of charge carriers. This effect, known as Coulomb drag [1], is of fundamental interest in condensed matter physics and nano-electronics, and has been observed in a wide variety of systems, including graphene bilayers [2, 3], 2D electron gases [4], edges states of fractional quantum Hall systems [5], quantum point contacts [6], nanowires [7, 8] and superconducting wires and films [9].

Another important transport phenomenon of interest in many-body systems is depinning [10], where transport in a driven system is inhibited by disorder 'pinning' the system in metastable states until the collective pinning force is overcome at some threshold driving strength. Depinning theory is crucial to the understanding of transport in systems such as charge density waves [11], magnetic domain boundaries [12] and flux-line lattices in type-II superconductors [13]. It has recently been shown that the onset of conduction in linear Josephson junction arrays at low voltages is also determined by the effects of depinning [14]. A unique system in which the interplay between depinning and Coulomb drag can be observed is the bilinear Josephson junction array.

Josephson junction arrays are ideal for studying low-dimensional electronic transport in a controllable system, as almost all parameters of the system can be tailored at will. Linear Josephson junction arrays have already been used to study depinning in a system where relevant parameters such as the amplitude of the pinning potential can be controlled [14]. To study both depinning and Coulomb drag, and their interplay, we focus on bilinear Josephson junction arrays. These systems consist of two linear chains of Josephson junctions coupled capacitively, so that the two arrays interact electrostatically but no direct transfer of charge carriers occurs between the two arrays. Experiments on these bilinear arrays have shown Coulomb drag and current mirror behaviour [15, 16], but the quantitative details of these effects have so far not been understood theoretically, nor modelled computationally.

Linear arrays are known to exhibit threshold behaviour when the charging energy of a single grain ${E}_{{\rm{C}}}={(2e)}^{2}/2{C}_{{\rm{J}}}$ is comparable to or greater than the Josephson tunnelling energy E_J [17]. In this regime the array acts as an insulator at small voltages, with a sudden transition to conducting behaviour at a switching voltage V_sw. Most work to date on these systems, both experimentally and theoretically, has described the threshold behaviour in terms of the minimum energy required to inject a single charge soliton into the array [17–20]. While this simple picture is elegant and attractive, it fails to account for effects arising from disorder in the system, which is always present in any realistic experiment. More recently, threshold behaviour in linear arrays was described in terms of the pinning of charge due to disorder in the array [21]. Such an approach has been able to achieve an excellent fit between the theoretical model and experimental data [14]. In this work, we extend this description to bilinear Josephson junction arrays.

Two different coupling geometries for bilinear arrays have been fabricated: straight coupling [16], where each site is coupled to only one site on the opposite array, and slanted coupling [15], where each site is coupled to two sites on the opposite array. Differences in the phase diagrams between these systems have been studied theoretically [22]. For concreteness, we will focus primarily of the case of straight coupling. However our methods can be readily applied to slanted coupling, and we will include results for both arrays where relevant.

In the bilinear Josephson junction array with straight coupling, depicted in figure 1, each site is coupled to two other sites on the same array via Josephson junctions, and to one site on the opposite array by a capacitance C_C. Each site also has a ground capacitance C_G. There is assumed to be no tunnelling of charge carriers between the two arrays.

Zoom In Zoom Out Reset image size

In a typical experiment, one of these arrays is held at a fixed voltage (assumed to be 0). We call this array the passive array. The voltage on the active array, however, is varied, so that an I–V curve can be determined. Both arrays remain insulating at low voltages. At the switching voltage V_sw both arrays become conducting, even though no voltage is directly applied to the passive array. The current through the active array is much larger, but the drag current in the passive array is clearly present and linearly proportional to the active bias.

Previous theoretical work on these systems has focused on their static features such as their free energy [23] and phase transitions [24], or has focused on dynamics in the quasiparticle-dominated limit [25]. Here we are more concerned with the dynamical behaviour when each grain is deep in the superconducting regime, and we explicitly derive equations of motion for the charges on each island, which we use to numerically calculate I–V curves.

For clarity, we will first derive the Hamiltonian and equations of motion for the bilinear array with straight coupling. From this it is trivial to alter our model to the case of slanted coupling.

Across each capacitor and junction in the array there will be a difference in the phase of the superconducting condensate. We use ϕ, ψ and χ to denote the phase difference across the Josephson junctions (C_J), ground capacitors (C_G) and coupling capacitors (C_C) respectively (see figure 1). From the circuit diagram in figure 1(a), standard theoretical techniques [26, 27] can be applied to derive a Lagrangian for the system,

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \sum _{i,\zeta }\displaystyle \frac{{C}_{{\rm{J}}}}{2}{({\dot{\phi }}_{i}^{\zeta })}^{2}+\displaystyle \frac{{C}_{{\rm{G}}}}{2}{({\dot{\psi }}_{i}^{\zeta })}^{2}+\displaystyle \frac{{C}_{{\rm{C}}}}{2}{({\dot{\chi }}_{i})}^{2}+{E}_{{\rm{J}}}\cos ({\phi }_{i}^{\zeta }),\end{eqnarray} \tag{ 1 }$$

where the index i runs over different sites in an array and the index ζ labels the two layers, $\zeta =\uparrow ,\,\downarrow$.

Each phase in the Lagrangian has a conjugate charge variable

$$\begin{eqnarray}&&{q}_{i}^{\phi ,\zeta }={C}_{{\rm{J}}}{\dot{\phi }}_{i}^{\zeta };\qquad {q}_{i}^{\psi ,\zeta }={C}_{{\rm{G}}}{\dot{\psi }}_{i}^{\zeta };\qquad {q}_{i}^{\chi }={C}_{{\rm{C}}}{\dot{\chi }}_{i}.\end{eqnarray} \tag{ 2 }$$

These charges are related to the number of Cooper pairs on a superconducting island by charge neutrality and Kirchoff's laws, which allow us to write the Hamiltonian in terms of the number of Cooper pairs

$$\begin{eqnarray}&&{H}=\displaystyle \sum _{{i},{j},{\zeta }{\zeta }}4{{e}}^{2}({{n}}_{{i}}^{{\zeta }}+{{f}}_{{i}}^{{\zeta }}){[{{C}}^{-1}]}_{{ij}}^{{\zeta }{{\zeta }}^{{}^{{}^{{\rm{{\prime} }}}}}}({{n}}_{{j}}^{{{\zeta }}^{{}^{{}^{{\rm{{\prime} }}}}}}+{{f}}_{{\rm{J}}}^{{{\zeta }}^{{}^{{}^{{\rm{{\prime} }}}}}})-{{E}}_{{\rm{J}}}\cos {{\phi }}_{{i}}^{{\zeta }}{{i}}^{{\zeta }},\end{eqnarray} \tag{ 3 }$$

where ${C}_{{ij}}^{\zeta \zeta ^{\prime} }$ is the capacitance matrix which expresses Coulomb interactions between sites (and therefore depends on the particular coupling geometry) and ${[{C}^{-1}]}_{{ij}}^{\zeta \zeta ^{\prime} }$ is the (i, j) element of the inverse of ${C}^{\zeta \zeta ^{\prime} }$. f_i is the charge frustration on the ith junction due to, e.g., trapped charges or defects, and is modelled by a random number ${f}_{i}\in (-1,1)$, which we take to be evenly distributed across the entire interval (although other disorder distributions have been considered, such as weak or Gaussian disorder [21, 28]).

We have chosen to include only charge disorder in our model, neglicting the possibility of, for example, disorder in E_J and C_J. We expect disorder in these quantities to be much weaker in experimental systems, and to have the same qualitative effect as charge disorder (that is, to modify the potential energy profile).

The inverse capacitance matrix in equation (3) includes long-range electrostatic interactions between islands. In discussing threshold phenomena in linear chains, it has proven useful to move away from this highly coupled model to a so-called quasicharge model, which includes only nearest-neighbour interactions.

To this end, we introduce cumulative quasicharge variables

$$\begin{eqnarray}{m}_{i}^{\zeta } & = & -\displaystyle \sum _{j=1}^{i-1}{n}_{j}^{\zeta };\quad {F}_{i}^{\zeta }=\displaystyle \sum _{j=1}^{i-1}{f}_{j}^{\zeta };\quad {Q}_{i}^{\zeta }=\displaystyle \sum _{j=1}^{i-1}{q}_{j}^{\psi ,\zeta }+{q}_{1}^{\phi ,\zeta };\\ {X}_{i} & = & \displaystyle \sum _{j=1}^{i-1}{q}_{j}^{\chi };\quad {n}_{i}^{\zeta }={m}_{i}^{\zeta }-{m}_{i+1}^{\zeta };\quad {q}_{i}^{\psi ,\zeta }={Q}_{i+1}^{\zeta }-{Q}_{i}^{\zeta },\end{eqnarray} \tag{ 4 }$$

where the inclusion of ${q}_{1}^{\phi ,\zeta }$ in the definition of ${Q}_{i}^{\zeta }$ is an offset term fixing the number of charges on the first site.

The charge variables Q are defined as the polarization across the junctions rather than the number of Cooper pairs on each grain. As was shown in [20], the charge solitons in a Josephson junction array can not be thought of as a single Cooper pair dressed by a polarization cloud, but rather the soliton is the polarization cloud itself. In our model, it is this polarization that is our relevant charge variable, and not the number of Cooper pairs.

We follow the derivation given in [14, 21] and note that these variables are constrained by requirements of charge neutrality, such that

$$\begin{eqnarray}&&{q}_{i}^{\phi ,\zeta }=\displaystyle \sum _{j=1}^{i-1}{q}_{j}^{\psi ,\zeta }+{q}_{1}^{\phi ,\zeta }-2e({n}_{i}^{\zeta }+{f}_{i}^{\zeta })+{(-1)}^{\zeta }{X}_{i},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&=\,{Q}_{i}^{\zeta }-2e({m}_{i}^{\zeta }+{F}_{i}^{\zeta })+{(-1)}^{\zeta }{X}_{i}.\end{eqnarray} \tag{ 6 }$$

Using these relations we can rewrite the Lagrangian in equation (1) in terms of quasichagre variables

$$\begin{eqnarray}{ \mathcal L } & = & \displaystyle \sum _{i,\zeta }\displaystyle \frac{{({Q}_{i}^{\zeta }+{(-1)}^{\zeta }{X}_{i}-2e({m}_{i}^{\zeta }+{F}_{i}^{\zeta }))}^{2}}{2{C}_{{\rm{J}}}}+\displaystyle \frac{{({Q}_{i}^{\zeta }-{Q}_{i+1}^{\zeta })}^{2}}{2{C}_{{\rm{G}}}}\\ & & +\displaystyle \frac{{({X}_{i}-{X}_{i+1})}^{2}}{2{C}_{{\rm{C}}}}+{E}_{{\rm{J}}}\cos ({\phi }_{i}^{\zeta }).\end{eqnarray} \tag{ 7 }$$

The Lagrangian can be further simplified by introducing a dimensionless coupling parameter $\alpha ={C}_{{\rm{C}}}/{C}_{{\rm{G}}}$ and noting that, due to Kirchhoff's laws, ${X}_{i}=\alpha ({Q}_{i}^{\uparrow }-{Q}_{i}^{\downarrow })$. For notational clarity, we will also switch to a vector representation where ${\vec{Q}}_{i}={({Q}_{i}^{\uparrow },{Q}_{i}^{\downarrow })}^{{\rm{T}}}$ (and likewise for other quantities defined on both arrays), and introducing the discrete differential operator ∇ such that ${\rm{\nabla }}{\vec{Q}}_{i}={\vec{Q}}_{i+1}-{\vec{Q}}_{i}$. Finally, we introduce a coupling matrix ${ \mathcal M }$.

This gives us the Hamiltonian of the system,

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \sum _{i}\left\{\displaystyle \frac{1}{2{C}_{{\rm{G}}}}{\rm{\nabla }}{\vec{Q}}_{i}^{{\rm{T}}}{ \mathcal M }{\rm{\nabla }}{\vec{Q}}_{i}\ +\displaystyle \frac{1}{2{C}_{{\rm{J}}}}{[{ \mathcal M }{\vec{Q}}_{i}+2e({\vec{m}}_{i}+{\vec{F}}_{i})]}^{2}-{E}_{{\rm{J}}}\cos ({\vec{\phi }}_{i})\right\},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}{ \mathcal M }=\left(\begin{array}{cc}1+\alpha & -\alpha \\ -\alpha & 1+\alpha \end{array}\right),\end{eqnarray} \tag{ 9 }$$

where $\alpha ={C}_{{\rm{C}}}/{C}_{{\rm{G}}}$.

We now assume that the Q are slow-changing parameters compared with m and ϕ, which is always the case when there is a sufficiently large inductance (in [14] it was shown that the Bloch inductance of the array is sufficient so long as the system is driven adiabatically so as to avoid Landau–Zener transitions). In this limit, we may separate the time-scales of evolution of these parameters and apply a Born–Oppenheimer approximation. We take the portion of the Hamiltonian which depends on m_i and ${\phi }_{i}$

$$\begin{eqnarray}&&{H}_{Q}({m}_{i},{\phi }_{i})=\displaystyle \sum _{i}\displaystyle \frac{1}{2{C}_{{\rm{J}}}}{[{ \mathcal M }{\vec{Q}}_{i}+2e({\vec{m}}_{i}+{\vec{F}}_{i})]}^{2}-{E}_{{\rm{J}}}\cos ({\vec{\phi }}_{i}),\end{eqnarray} \tag{ 10 }$$

and take Q to be a constant, classical parameter. In this Hamiltonian, the coupling matrix ${ \mathcal M }$ only acts on the classical parameter Q, and not on the quantum operators m_i and ${\phi }_{i}$. Therefore, the Hamiltonian separates into the sum of single-site Hamiltonians, which can be easily diagonalized numerically. Diagonalizing H_Q for various values of Q, and taking only the lowest energy band, gives us an effective potential ${E}_{Q}({Q}^{\zeta },{F}^{\zeta })={E}_{Q}({Q}^{\uparrow },{F}^{\uparrow })+{E}_{Q}({Q}^{\downarrow },{F}^{\downarrow })$. The form of this potential is equivalent to the characteristic value of Mathieu's equation [29] $a({\nu }^{\zeta })$ with argument ${\nu }^{\zeta }=(1+\alpha ){Q}^{\zeta }-\alpha {Q}^{\zeta ^{\prime} }+2{{eF}}^{\zeta }$, where $\zeta ^{\prime}$ simply indicates the array opposite to ζ. This is a $2e$-periodic function which reduces to $a(\nu )\sim \cos (\nu )$ in the limit ${E}_{{\rm{J}}}\gg {E}_{{\rm{C}}}$. Applying this approximation, we get an effective semi-classical Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}\displaystyle \frac{1}{2{C}_{{\rm{G}}}}{\rm{\nabla }}{\vec{Q}}_{i}^{{\rm{T}}}{ \mathcal M }{\rm{\nabla }}{\vec{Q}}_{i}+{E}_{Q}({ \mathcal M }{\vec{Q}}_{i}+2e\vec{F}).\end{eqnarray} \tag{ 11 }$$

In this Hamiltonian, there exists only nearest-neighbour coupling between islands in the array. Coupling between arrays is mediated by the matrix ${ \mathcal M }$, and charge frustration in the array enters in the form of the disordered periodic potential E_Q.

In this model E_Q is responsible for pinning the system and preventing charge transport in the arrays. The presence of ${ \mathcal M }$ in the argument of the potential indicates that this pinning potential is highly coupled. To understand the depinning of each of the arrays separately, we introduce new charge variables which are decoupled in the pinning potential

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}^{\uparrow }=(1+\alpha ){Q}^{\uparrow }-\alpha {Q}^{\downarrow },\qquad {{\rm{\Upsilon }}}^{\downarrow }=(1+\alpha ){Q}^{\downarrow }-\alpha {Q}^{\uparrow },\end{eqnarray} \tag{ 12 }$$

which is the same as 'rotating' the vectors ${\vec{Q}}_{i}$ by the matrix ${ \mathcal M }$. In this rotated frame, the Hamiltonian is

$$\begin{eqnarray}&&H=\displaystyle \sum _{i}\displaystyle \frac{1}{2{C}_{{\rm{G}}}}{\rm{\nabla }}{\vec{{\rm{\Upsilon }}}}_{i}^{{\rm{T}}}{{ \mathcal M }}^{-1}{\rm{\nabla }}{\vec{{\rm{\Upsilon }}}}_{i}+{E}_{Q}({\vec{{\rm{\Upsilon }}}}_{i}+2e{\vec{F}}_{i}).\end{eqnarray} \tag{ 13 }$$

Coupling between the two arrays occurs only in the charging term, not in the pinning potential. The theory of depinning would thus lead us to expect that each array in this rotated frame has a separate pinning potential, and therefore a separate depinning transition. If a voltage were to be applied to only one array in this frame, only one array will depin, and there will be no drag current. However, when a voltage is applied only to the active array in the unrotated frame (as is the case experimentally), there is an effective voltage on each array in the rotated frame, given by

$$\begin{eqnarray}&&\left(\begin{array}{c}{V}_{{\rm{rot}}}^{\uparrow }\\ {V}_{{\rm{rot}}}^{\downarrow }\end{array}\right)={ \mathcal M }\left(\begin{array}{c}{V}_{{\rm{lab}}}^{\uparrow }\\ 0\end{array}\right)=\left(\begin{array}{c}(1+\alpha ){V}_{{\rm{lab}}}^{\uparrow }\\ -\alpha {V}_{{\rm{lab}}}^{\uparrow }\end{array}\right),\end{eqnarray} \tag{ 14 }$$

where the subscripts correspond to the voltage applied in the rotated frame and the physical laboratory frame. The voltage applied to the upper array in the rotated frame is larger than that applied to the lower array, so we expect this array to depin first. At this point, there will be a current flowing through the upper rotated array, but not the lower rotated array. This regime corresponds to the point in the lab frame where both arrays are conducting. The lab frame currents will be given by

$$\begin{eqnarray}&&\left(\begin{array}{c}{I}_{{\rm{lab}}}^{\uparrow }\\ {I}_{{\rm{lab}}}^{\downarrow }\end{array}\right)={{ \mathcal M }}^{-1}\left(\begin{array}{c}{I}_{{\rm{rot}}}^{\uparrow }\\ 0\end{array}\right),\end{eqnarray} \tag{ 15 }$$

which allows us to predict the ratio between the active and passive currents in the lab frame as

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{{\rm{lab}}}^{\uparrow }}{{I}_{{\rm{lab}}}^{\downarrow }}=\displaystyle \frac{1+\alpha }{\alpha }.\end{eqnarray} \tag{ 16 }$$

As we increase the voltage in the lab frame, the voltage on the lower array in the rotated frame will increase until it too reaches its threshold. As can be seen in equation (14), the effective voltage felt by the lower array will be negative, so that in the rotated frame the two currents will flow in opposite directions. In the rotated frame, the ratio between currents is

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{{\rm{rot}}}^{\uparrow }}{{I}_{{\rm{rot}}}^{\downarrow }}=-\displaystyle \frac{1+\alpha }{\alpha }.\end{eqnarray} \tag{ 17 }$$

Switching back to the lab frame, the current on the passive array is

$$\begin{eqnarray}&&{I}_{{\rm{lab}}}^{\downarrow }=\displaystyle \frac{\alpha }{1+2\alpha }{I}_{{\rm{rot}}}^{\uparrow }+\displaystyle \frac{1+\alpha }{1+2\alpha }{I}_{{\rm{rot}}}^{\downarrow },\end{eqnarray} \tag{ 18 }$$

which, upon using equation (17) to express this in terms of only one of the rotated currents, gives us ${I}_{{\rm{lab}}}^{\downarrow }=0$.

Here we have seen a separation of the pinning of each array. When the first threshold is reached, the active array depins and we enter the Coulomb drag regime. As voltage is increased beyond this first threshold, we predict that there will exist a second threshold voltage, where the passive array will depin. After the second threshold, Coulomb drag behaviour will cease and the passive current will drop to zero (provided there is no voltage directly applied to the passive array).

From the Hamiltonian in equation (13) we can derive equations of motion for the variable ${\vec{Q}}_{i}$,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{C}_{{\rm{G}}}}{ \mathcal M }{\rm{\nabla }}{\vec{Q}}_{i}+{V}_{Q}({ \mathcal M }{\vec{Q}}_{i}+{\vec{F}}_{i})+\rho {\dot{\vec{Q}}}_{i}+L{\ddot{\vec{Q}}}_{i}=0,\end{eqnarray} \tag{ 19 }$$

which in the rotated frame becomes

$$\begin{eqnarray}&&\displaystyle \frac{1}{{C}_{{\rm{G}}}}{{ \mathcal M }}^{-1}{\rm{\nabla }}{\vec{{\rm{\Upsilon }}}}_{i}+{V}_{Q}({\vec{{\rm{\Upsilon }}}}_{i}+{\vec{F}}_{i})+\rho {{ \mathcal M }}^{-1}{\dot{\vec{{\rm{\Upsilon }}}}}_{i}+L{{ \mathcal M }}^{-1}{\ddot{\vec{{\rm{\Upsilon }}}}}_{i}=0,\end{eqnarray} \tag{ 20 }$$

where ${V}_{Q}(Q)=\partial {E}_{Q}(Q)/\partial Q$, ρ is the resistance per site of the array and L is the inductance. ρ and L are phenomenological and have been including to ensure numerical convergence (it should be noted that ρ is not necessarily related to the normal-state resistance R, but is rather the sub-gap resistance that is observed even in the Cooper-pair limit [30]). Since these terms couple to $\dot{Q}$ and $\ddot{Q}$ respectively, these terms do not influence the location of the threshold voltages. Because ρ and L are strictly phenomenological, their values are chosen for numerical convenience. This model is therefore not able to quantitatively predict the absolute magnitude of the currents in the system. Solving these equations of motion should, however, correctly return qualitative structure of the I–V curve, as well as quantities independent of L and ρ, such as the ratio between currents ${I}^{\uparrow }/{I}^{\downarrow }$.

Many of these results are readily extended to the case of slanted coupling. One must replace the 2-component quasicharge vector ${\vec{Q}}_{i}={({Q}_{i}^{\uparrow },{Q}_{i}^{\downarrow })}^{{\rm{T}}}$ with the $2N$-component vector ${({Q}_{1}^{\uparrow },{Q}_{1}^{\downarrow },\ldots {Q}_{k}^{\uparrow },{Q}_{k}^{\downarrow },\ldots {Q}_{N}^{\uparrow },{Q}_{N}^{\downarrow })}^{{\rm{T}}}$, which we shall denote by $\vec{Q}$ as it should always be clear from context whether we are referring to the 2-component vector or the $2N$-component vector. Then the 2 × 2 coupling matrix ${ \mathcal M }$ is replaces by a $2N\times 2N$ matrix

$$\begin{eqnarray}{ \mathcal M }=\left(\begin{array}{ccccc}1+\alpha & -\alpha & 0 & ... & \\ -\alpha & 1+2\alpha & -\alpha & & \\ 0 & \ddots & \ddots & \ddots & \\ \vdots & & -\alpha & 1+2\alpha & -\alpha \\ & & & -\alpha & 1+\alpha \end{array}\right).\end{eqnarray} \tag{ 21 }$$

This leads to a Hamiltonian

$$\begin{eqnarray}{ \mathcal H } & = & \displaystyle \frac{1}{2{C}_{{\rm{G}}}}{\rm{\nabla }}{\vec{Q}}^{{\rm{T}}}{ \mathcal M }{\rm{\nabla }}\vec{Q}+{E}_{Q}({ \mathcal M }\vec{Q})\\ & = & \displaystyle \frac{1}{2{C}_{{\rm{G}}}}{\rm{\nabla }}{\vec{{\rm{\Upsilon }}}}^{{\rm{T}}}{{ \mathcal M }}^{-1}{\rm{\nabla }}\vec{{\rm{\Upsilon }}}+{E}_{Q}(\vec{{\rm{\Upsilon }}}),\end{eqnarray} \tag{ 22 }$$

where, as with straight coupling, $\vec{{\rm{\Upsilon }}}={ \mathcal M }\vec{Q}$. This gives us an equation of motion identical to equation (20), but with the slanted definitions of $\vec{Q}$, $\vec{{\rm{\Upsilon }}}$ and ${ \mathcal M }$.

The form of the coupling matrix for the slanted array is the same as that of the capacitance matrix for a linear Josephson junction array (with the identifications ${C}_{{\rm{J}}}\to \alpha$ and ${C}_{{\rm{G}}}\to 1$), and so it can be inverted analytically using the same methods [23, 31].

Much of the analysis of the array with straight coupling depended on the fact that the Hamiltonian separates neatly into a sum of two-site terms This is not the case with slanted coupling, and we are not able to obtain the same analytical results.

The equation of motion equation (20) enables us to numerically determine I–V curves for the bilinear array with either slanted or straight coupling. In figure 2 we present such I–V curves for arrays with straight coupling calculated at a variety of different couplings α using parameters given in the caption to figure 2, which were chosen for numerical speed and convenience. The existence of two distinct pinning thresholds can clearly be seen, and the predicted current ratios hold. The qualitative form of the I–V curves at low voltages (i.e. before the second threshold) are in good agreement with the published experimental results of [16].

Zoom In Zoom Out Reset image size

I–V curves for the case with slanted coupling are very similar in form, but with lower threshold voltages. This is in contrast to the results reported in [15], where qualitatively different results were observed. This is not surprising, however, as the I–V curve presented in that work was obtained at finite magnetic field and with a finite voltage applied to the passive array. Furthermore, the array investigated in [15] consisted of SQUID-loops rather than simple Josephson junctions and are therefore susceptible to the effects of low-frequency flux noise and other complications not considered in the present work.

The exact location of the threshold voltages depends on the disorder realisation in the system. Different experimental systems will have different disorder realisations, and therefore there are a range of possible threshold voltages. In any numerical simulation one disorder realisation must be used, and therefore no single numerical I–V curve will exhibit exactly the same threshold voltage as a corresponding single experimental I–V curve unless both the theoretical and the experimental systems have exactly the same disorder realisation. General qualitative trends, however, such as the dependence of the thresholds on the coupling α, remain true across different disorder realisations.

This scheme can be easily generalized to the case where the voltage on both arrays is non-zero. In this case we obtain an I–V surface, with the current being a function of both ${V}^{\uparrow }$ and ${V}^{\downarrow }$. To study the depinning thresholds and Coulomb drag behaviour change as both voltages are varied, we calculate a conduction diagram, as depicted in figure 3. Slight asymmetries in the diagram arise due to asymmetry in the disorder realisation. The conduction diagram for the (physically unrealistic) completely clean bilinear array ($\vec{F}=0$) has no such asymmetry.

From the conduction diagram we can see that there exist regions where, in the absence of coupling, a current would flow in one direction, but due to overpowering Coulomb drag it instead flows in the other direction. These regions of overpowering Coulomb drag are marked in figure 3 by red arrows. The conduction diagram for slanted coupling, figure 4, displays even larger regions of overpowering Coulomb drag given the same parameters. This can be understood as being due to a larger effective coupling between the two arrays with slanted coupling, as in there are more coupling capacitors.

In numerical simulation, if one wishes to observe the second depinning threshold it is vital that the voltage is applied adiabatically and that the system retains memory of past charge configurations. This is to be expected as the microscopic details of depinning theory depend heavily on the existence of metastable states [10]. When the system is driven diabatically we see a current not because the insulating state is unstable against the formation of current, but because the system has been driven sufficiently fast to push it away from stable configurations. Simulations in which the voltage was increased slowly displayed second threshold behaviour (figure 2). In figures 3 and 4, however, we are interested only in the conduction behaviour near the first threshold. These simulation were performed with a diabatic switch-on of voltage, so the second threshold could not be seen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Experimentally, to observe the second depinning threshold one should ensure that the this threshold lies below the quasiparticle excitation gap, $V\approx 2N{\rm{\Delta }}/e$. The second threshold can be kept at a relatively low voltage if one ensures that the coupling to ground C_G is small, as this ensures that both the dimensionless coupling α and the effective interaction length Λ will be large.

One must also ensure that the coupling is not too strong. Our entire theoretical approach is based around an understanding that the system is in the insulating regime of the quantum phase diagram. It has been shown that strong coupling between the arrays leads to an insulator-superconductor quantum phase transition in bilinear JJ arrays, even when each array is individually in the insulating state [22]. We therefore only expect our theory to hold when $\sqrt{{E}_{{\rm{J}}}{C}_{{\rm{C}}}/{(2e)}^{2}}\lt 2\sqrt{2}/\pi \approx 0.9$ for straight coupling and $\sqrt{{E}_{{\rm{J}}}{C}_{{\rm{C}}}/{(2e)}^{2}}\lt 2/\pi \approx 0.6$ for slanted coupling (as is the case in all of our simulations).

Bilinear Josephson junction arrays provide an ideal situation in which to study Coulomb drag and depinning in a controllable system. The fact that each array experiences a separate pinning force leads to two depinning thresholds in the I–V curve, and the interplay between depinning and Coulomb drag leads to novel behaviour. Conclusions made with respect to Josephson junction arrays may be applicable to other systems in which both Coulomb drag and depinning effects are present.

I–V curves of bilinear arrays for both straight and slanted coupling display the same threshold behaviour, however both threshold voltages are lower in the case of slanted coupling. Furthermore, we see that from the conduction diagrams that the array with slanted coupling has much larger regions of overpowering Coulomb drag.

We thank T Duty and H Shimada for useful discussions. This work was supported in part by the Australian Research Council under the Discovery and Centre of Excellence funding schemes (project numbers: DP140100375 and CE170100039). Computational resources were provided by the NCI National Facility systems at the Australian National University through the National Computational Merit Allocation Scheme supported by the Australian Government.