Coulomb drag and depinning in bilinear Josephson junction arrays

TL;DR

This paper investigates Coulomb drag and depinning in bilinear Josephson junction arrays, combining experimental observations, simulations, and analytical models to understand their interplay and predict new phenomena at higher voltages.

Contribution

It introduces a quasicharge model for these phenomena in bilinear Josephson junction arrays and derives analytical expressions, extending understanding beyond previous experimental results.

Findings

Numerical I-V curves match experimental data.

Analytical depinning-based expressions for current ratios.

Predictions of novel high-voltage phenomena.

Abstract

Coulomb drag and depinning are electronic transport phenomena that occur in low-dimensional nanostructures. Recently, both phenomena have been reported in bilinear Josephson junction arrays. These devices provide a unique opportunity to study the interplay of Coulomb drag and depinning in a system where all relevant parameters can be controlled experimentally. We explain the Coulomb drag and depinning characteristics in the I-V curve of the bilinear Josephson junction array by adopting a quasicharge model which has previously proven useful in describing threshold phenomena in linear Josephson junction arrays. Simulations are performed for a range of coupling strengths, where numerically obtained I-V curves match well with what has been previously observed experimentally. Analytic expressions for the ratio between the active and passive currents are derived from depinning arguments. Novel phenomena are predicted at voltages higher than those for which experimental results have been reported to date.

Equations

$$\mathcal{L}=\sum _{i,\zeta }\frac{C_J}{2}({\dot{\varphi }}_i^{\zeta }{)}^2+\frac{C_G}{2}({\dot{\psi }}_i^{\zeta }{)}^2+\frac{C_C}{2}({\dot{\chi }}_i{)}^2+E_J\cos ({\varphi }_i^{\zeta })$$

$$q_i^{\varphi ,\zeta }=C_J{\dot{\varphi }}_i^{\zeta };q_i^{\psi ,\zeta }=C_G{\dot{\psi }}_i^{\zeta };q_i^{\chi }=C_C{\dot{\chi }}_i.$$

$$\begin{array}{r}\mathcal{H}=\sum _{i,j}4e^2(n_i^{\zeta }+f_i^{\zeta }){\left[C^{-1}\right]}_{ij}^{\zeta {\zeta }^{\prime }}(n_j^{{\zeta }^{\prime }}+f_J^{{\zeta }^{\prime }})\\ -E_J\cos {\varphi }_i^{\zeta }-E_J\cos {\varphi }_i^{{\zeta }^{\prime }}\end{array}$$

$$\begin{array}{r}{m}_i^{\zeta }=-\sum _{j=1}^{i-1}{n}_j^{\zeta };F_i^{\zeta }=\sum _{j=1}^{i-1}{f}_j^{\zeta };Q_i^{\zeta }=\sum _{j=1}^{i-1}{q}_j^{\psi ,\zeta }+q_1^{\zeta };\\ X_i=\sum _{j=1}^{i-1}{q}_j^{\chi };n_i^{\zeta }=m_i^{\zeta }-m_{i+1}^{\zeta };q_i^{\psi ,\zeta }=Q_{i+1}^{\zeta }-Q_i^{\zeta }.\end{array}$$

$$q_i^{\varphi ,\zeta }$$

$$=Q_i^{\zeta }-2e(m_i^{\zeta }+F_i^{\zeta })+(-1{)}^{\zeta }{X}_i.$$

$$\begin{array}{r}\mathcal{L}=\sum _{i,\zeta }\frac{{\left(Q_i^{\zeta }+(-1{)}^{\zeta }{X}_i-2e(m_i^{\zeta }+F_i^{\zeta })\right)}^2}{2C_J}+\frac{{\left(Q_i^{\zeta }-Q_{i+1}^{\zeta }\right)}^2}{2C_G}\\ +\frac{{\left(X_i-X_{i+1}\right)}^2}{2C_C}+E_J\cos \left({\varphi }_i^{\zeta }\right).\end{array}$$

$$\begin{array}{rl}\mathcal{H}=& \sum _i\{\frac{1}{2C_G}\nabla {\vec{Q}}_i^{\text{T}}\mathcal{M}\nabla {\vec{Q}}_i\\ & +\frac{1}{2C_J}{\left[\mathcal{M}{\vec{Q}}_i+2e\left({\vec{m}}_i+{\vec{F}}_i\right)\right]}^2-E_J\cos \left({\vec{\varphi }}_i\right)\}\end{array}$$

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

$$H_Q(m_i,{\varphi }_i)=\sum _i\frac{1}{2C_J}{\left[\mathcal{M}{\vec{Q}}_i+2e\left({\vec{m}}_i+{\vec{F}}_i\right)\right]}^2-\cos ({\vec{\varphi }}_i)$$

$$H=\sum _i\frac{1}{2C_G}\nabla {\vec{Q}}_i^{\text{T}}\mathcal{M}\nabla {\vec{Q}}_i+E_Q(\mathcal{M}{\vec{Q}}_i+2e\vec{F}).$$

$${{\rm Y}}^{\uparrow }=(1+\alpha )Q^{\uparrow }-\alpha Q^{\downarrow },{{\rm Y}}^{\downarrow }=(1+\alpha )Q^{\downarrow }-\alpha Q^{\uparrow }$$

$$H=\sum _i\frac{1}{2C_G}\nabla {\overrightarrow{{\rm Y}}}_i^{\text{T}}{\mathcal{M}}^{-1}\nabla {\overrightarrow{{\rm Y}}}_i+E_Q({\overrightarrow{{\rm Y}}}_i+2e{\vec{F}}_i).$$

$$\left(\begin{array}{c}{V}_{\text{rot}}^{\uparrow }\\ V_{\text{rot}}^{\downarrow }\end{array}\right)=\mathcal{M}\left(\begin{array}{c}{V}_{\text{lab}}^{\uparrow }\\ 0\end{array}\right)=\left(\begin{array}{c}(1+\alpha )V_{\text{lab}}^{\uparrow }\\ -\alpha V_{\text{lab}}^{\uparrow }\end{array}\right)$$

$$\left(\begin{array}{c}{I}_{\text{lab}}^{\uparrow }\\ I_{\text{lab}}^{\downarrow }\end{array}\right)={\mathcal{M}}^{-1}\left(\begin{array}{c}{I}_{\text{rot}}^{\uparrow }\\ 0\end{array}\right)$$

$$\frac{I_{\text{lab}}^{\uparrow }}{I_{\text{lab}}^{\downarrow }}=\frac{1+\alpha }{\alpha }.$$

$$\frac{I_{\text{rot}}^{\uparrow }}{I_{\text{rot}}^{\downarrow }}=-\frac{1+\alpha }{\alpha }.$$

$$I_{\text{lab}}^{\downarrow }=\frac{\alpha }{1+2\alpha }{I}_{\text{rot}}^{\uparrow }+\frac{1+\alpha }{1+2\alpha }{I}_{\text{rot}}^{\downarrow }$$

$$\frac{1}{C_G}\mathcal{M}\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,$$

$$\frac{1}{C_G}{\mathcal{M}}^{-1}\nabla {\overrightarrow{{\rm Y}}}_i+V_Q({\overrightarrow{{\rm Y}}}_i+{\vec{F}}_i)+\rho {\mathcal{M}}^{-1}{\dot{\overrightarrow{{\rm Y}}}}_i+L{\mathcal{M}}^{-1}{\ddot{\overrightarrow{{\rm Y}}}}_i=0$$

$$\mathcal{M}=\left(\begin{array}{cccccc}1+\alpha & -\alpha & 0& \dots & & \\ -\alpha & 1+2\alpha & -\alpha & & & \\ 0& \ddots & \ddots & \ddots & & \\ ⋮& & -\alpha & 1+2\alpha & -\alpha & \\ & & & -\alpha & 1+\alpha & \end{array}\right).$$

$$\begin{array}{rl}\mathcal{H}=& \frac{1}{2C_G}\nabla {\vec{Q}}^{\text{T}}\mathcal{M}\nabla \vec{Q}+E_Q(\mathcal{M}\vec{Q})\\ =& \frac{1}{2C_G}\nabla {\overrightarrow{{\rm Y}}}^{\text{T}}{\mathcal{M}}^{-1}\nabla \overrightarrow{{\rm Y}}+E_Q(\overrightarrow{{\rm Y}})\end{array},$$

Full text

I Introduction

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 Narozhny and Levchenko (2016), is of fundamental interest in condensed matter physics and nano-electronics, and has been observed in a wide variety of systems, including graphene bilayers Gorbachev et al. (2012); Li et al. (2016), 2-D electron gases Zheng and MacDonald (1993), edges states of fractional quantum Hall systems Orgad and Levit (1996), quantum point contacts Khrapai et al. (2007), nanowires Glazman et al. (2006); Klesse and Stern (2000) and superconducting wires and films Duan and Yip (1993).

Another important transport phenomenon of interest in many-body systems is depinning Brazovskii and Nattermann (2004), 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 wavesFukuyama and Lee (1978), magnetic domain boundaries Imry and Ma (1975) and flux-line lattices in type-II superconductors Larkin and Ovchinnikov (1979). 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 Vogt et al. (2015). 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 Vogt et al. (2015). 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 Shimada and Delsing (2000); Shimada et al. (2012), 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_{C}=(2e)^{2}/2C_{J}$ is comparable to or greater than the Josephson tunnelling energy $E_{J}$ Haviland and Delsing (1996). 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 arrayHaviland and Delsing (1996); Fedorov et al. (2011); Homfeld et al. (2011); Hermon et al. (1996). 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 arrayVogt et al. (2016). Such an approach has been able to achieve an excellent fit between the theoretical model and experimental dataVogt et al. (2015). In this work, we extend this description to bilinear Josephson junction arrays.

Two different coupling geometries for bilinear arrays have been fabricated: straight coupling Shimada et al. (2012), where each site is coupled to only one site on the opposite array, and slanted coupling Shimada and Delsing (2000), where each site is coupled to two sites on the opposite array. Differences in the phase diagrams between these systems have been studied theoretically Choi et al. (1998a). 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 Fig. 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.

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 Hu and O’Connell (1996) and phase transitions Choi et al. (1998b), or has focused on dynamics in the quasiparticle-dominated limit Walker and Cole (2013). 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.

II Theoretical Model

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 $\phi$, $\psi$ and $\chi$ to denote the phase difference across the Josephson junctions ($C_{J}$), ground capacitors ($C_{G}$) and coupling capacitors ($C_{C}$) respectively (see Fig. 1). From the circuit diagram in Fig. 1 (a), standard theoretical techniques Devoret (1997); Zagoskin (2011) can be applied to derive a Lagrangian for the system,

[TABLE]

where the index $i$ runs over different sites in an array, and the index $\zeta$ runs over the two arrays.

Each phase in the Lagrangian has a conjugate charge variable

[TABLE]

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

[TABLE]

where $f_{i}$ is the charge frustration on the $i$th 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 disorderVogt et al. (2016); Walker et al. (2015)). $C_{ij}^{\zeta\zeta^{\prime}}$ is the capacitance matrix which expresses Coulomb interactions between sites (and therefore depends on the particular coupling geometry). The inverse capacitance matrix 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

[TABLE]

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 Ref. Hermon et al. (1996), 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 Vogt et al. (2016, 2015) and note that these variables are constrained by requirements of charge neutrality, such that

[TABLE]

Using these relations we can rewrite the Lagrangian in Eq. 1 in terms of quasichagre variables

[TABLE]

The Lagrangian can be further simplified by introducing a dimensionless coupling parameter $\alpha=C_{C}/C_{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})^{\textrm{T}}$ (and likewise for other quantities defined on both arrays), and introducing the discrete differential operator $\nabla$ such that $\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,

[TABLE]

[TABLE]

where $\alpha=C_{C}/C_{G}$.

We now assume that the $Q$ are slow-changing parameters compared with $m$ and $\phi$, which is always the case when there is a sufficiently large inductance (in Vogt et al. (2015) 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}$

[TABLE]

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 McLachlan (1947) $a(\nu^{\zeta})$ with argument $\nu^{\zeta}=(1+\alpha)Q^{\zeta}-\alpha Q^{\zeta^{\prime}}+2eF^{\zeta}$, where $\zeta^{\prime}$ simply indicates the array opposite to $\zeta$. This is a $2e$-periodic function which reduces to $a(\nu)\sim\cos(\nu)$ in the limit $E_{J}\gg E_{C}$. Applying this approximation, we get an effective semi-classical Hamiltonian

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 Eq. 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

[TABLE]

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

[TABLE]

which, upon using Eq. 17 to express this in terms of only one of the rotated currents, gives us $I_{\textrm{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 0 (provided there is no voltage directly applied to the passive array).

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

[TABLE]

which in the rotated frame becomes

[TABLE]

where $V_{Q}(Q)=\partial E_{Q}(Q)/\partial Q$, $\rho$ is the resistance per site of the array and $L$ is the inductance. $\rho$ and $l$ are phenomenological and have been including to ensure numerical convergence (it should be noted that $\rho$ 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 Ågren et al. (2001)). Since these terms couple to $\dot{Q}$ and $\ddot{Q}$ respectively, these terms do not influence the location of the threshold voltages. Because $\rho$ 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 $\rho$, 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})^{\textrm{T}}$ with the $2N$-component vector $(Q_{1}^{\uparrow},Q_{1}^{\downarrow},\dots Q_{k}^{\uparrow},Q_{k}^{\downarrow},\dots Q_{N}^{\uparrow},Q_{N}^{\downarrow})^{\textrm{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\times 2$ coupling matrix $\mathcal{M}$ is replaces by a $2N\times 2N$ matrix

[TABLE]

This leads to a Hamiltonian

[TABLE]

where, as with straight coupling, $\vec{\Upsilon}=\mathcal{M}\vec{Q}$. This gives us an equation of motion identical to Eq. 20, but with the slanted definitions of $\vec{Q}$, $\vec{\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_{J}\rightarrow\alpha$ and $C_{G}\rightarrow 1$), and so it can be inverted analytically using the same methods Hu and O’Connell (1996); Cole et al. (2014).

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.

III Results of numerical simulations

The equation of motion Eq. 20 enables us to numerically determine I-V curves for the bilinear array with either slanted or straight coupling. In Fig. 2 we present such I-V curves for arrays with straight coupling calculated at a variety of different couplings $\alpha$ 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 Ref. Shimada et al. (2012).

I-V curves for the case with slanted coupling are very similar in form, but with lower threshold voltages.

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 $\alpha$, 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 Fig. 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 Fig. 3 by red arrows. The conduction diagram for slanted coupling, Fig. 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 Brazovskii and Nattermann (2004). 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 (Fig. 2), but the second threshold was not seen in simulations where the voltage is applied in sharp steps (i.e. diabatically).

Experimentally, to observe the second depinning threshold one should ensure that the this threshold lies below the quasiparticle excitation gap, $V\approx 2N\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 $\alpha$ and the effective interaction length $\Lambda$ 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 stateChoi et al. (1998a). We therefore only expect our theory to hold when $\sqrt{E_{J}C_{C}/(2e)^{2}}<2\sqrt{2}/\pi\approx 0.9$ for straight coupling and $\sqrt{E_{J}C_{C}/(2e)^{2}}<2/\pi\approx 0.6$ for slanted coupling (as is the case in all of our simulations).

IV Conclusion

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.

V Acknowledgements

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.