Depinning of disordered bosonic chains

We consider one-dimensional bosonic chains with a repulsive boson-boson interaction that decays exponentially on large length-scales. This model describes transport of Cooper-pairs in a Josepshon junction array, or transport of magnetic flux quanta in quantum-phase-slip ladders, i.e. arrays of superconducting wires in a ladder-configuration that allow for the coherent tunnelling of flux quanta. In the low-frequency, long wave-length regime these chains can be mapped to an effective model of a one-dimensional elastic field in a disordered potential. The onset of transport in these systems, when biased by external voltage, is described by the standard depinning theory of elastic media in disordered pinning potentials. We numerically study the regimes that are of relevance for quantum-phase-slip ladders. These are (i) very short chains and (ii) the regime of weak disorder. For chains shorter than the typical pinning length, i.e., the Larkin length, the chains reach a saturation regime where the depinning voltage does not depend on the decay length of the repulsive interaction. In the regime of weak disorder we find an emergent correlation length-scale that depends on the disorder strength. For arrays shorter than this length the onset of transport is similar to the clean arrays, i.e., is due to the penetration of solitons into the array. We discuss the depinning scenarios for longer arrays in this regime.


I. INTRODUCTION
Depinning theory describes the onset of propagation in many different physical systems. Examples include electrical transport in charge density waves 1,2 , the critical current of magnetic flux lattices 3,4 in type II superconductors, the propagation of magnetic domain boundaries 5 and crack formation in strained materials 6 . It was recently shown that the onset of electrical transport in one-dimensional arrays of Josephson junctions is also determined by the depinning of the chargeconfiguration along the array 7 .
In this paper we consider a more general model, a discrete chain occupied by Bosons with a repulsive interactions that decays exponentially on long length-scales. In such a model the interaction between neighbouring islands can be expressed by introducing a continuous variable, quasi-charge/flux, whose value is determined by the distribution of bosons along the chain. Assuming that the continuous variable changes adiabatically, an effective model can be derived with the help of the Born-Oppenheimer approximation 7,8 . In the case that disorder is present in the system, depinning theory can be applied to find the critical driving force that leads to a steady boson transport through the system. The Josephson junction arrays studied in Ref. 7 represent a particular realization of this model. Alternatively, our results describe the dual system of quantum phase slip ladders. In the latter case (QPS ladders) the bosons are magnetic flux quanta that tunnel through quantum phase slip elements 9-11 that separate the loops in a ladder.
In voltage biased Josephson junction arrays, the depinning theory describes the transition from an insulating regime at low voltages to a transport regime at higher voltages. The critical voltage of the transition is referred to as the switching voltage. The insulating regime of the arrays is governed by an effective sine-Gordon-like quasicharge model.
In the study of the onset of transport in Josephson junction arrays, Ref. 7, the connection to standard depinning theory was established under the assumption of strongly disordered background charges (also referred to as the maximal disorder model). Under this assumption the disorder-term in the effective model is spatially uncorrelated, allowing one to apply the standard depinning theory. Additionally the mapping to the standard depinning theory assumes the array length N to be much larger than all other length-scales of the problem.
In this paper we study the depinning behavior of more general chain models in a parameter regime that is not found in Josephson junction arrays. Specifically we consider chains which do not meet one of the two aforementioned conditions, e.g. the chains are either short or only weakly disordered. In the case of short chains we find a saturation regime. In the case of weak disorder, we find spatially correlated long-range disorder in the effective sine-Gordon-like model. In the considered chain model the spatial correlation in the disorder decays approximately exponentially for large distances. The depinning process in systems with spatial correlations that decay with a power law has been studied in Ref. 12 with the help of functional renormalization group theory.

A. The model
We consider a chain of islands, where N i are the discrete bosonic occupation numbers of the islands, N 0 is the average occupation number at equilibrium (chemical potential, positive charge background), f i are the random gate charges, b i , b † i are the bosonic annihilation and creation operators (N i = b † i b i ) and the bosonic tunnelling amplitude is given by t. In the limit We assume that the long-range behavior of the interaction matrix M i,j = M |i−j| is determined by an exponential decay on a length-scale Λ, One concrete example of this situation is a Josephson array, a chain of superconducting island coupled via Josephson junctions and with self-capacitances C 0 (capacitances to the ground) and capacitances C between the neighboring islands (junction capacitances). In this case E t = E J is the Josephson energy whereas the coupling matrix is given by This gives for the Fourier transform M (k) ≡ FT (M i−j ) and for M |i−j| the following expressions: where the junction charging energy is defined by E C ≡ (2e) 2 /2C. In particular the energy cost of a single charge in such an array is of the order M j,j ≈ ΛE C . Activated behavior with activation energy of order ΛE C was observed in Ref. 13. It has been realized long ago 14,15 that for Λ 1 the junction variables provide a more appropriate description than the island ones. We introduce, thus The resulting Hamiltonian has the form, with the modified coupling matrix Here θ i ≡ ϕ i−1 − ϕ i . One can easily check that m i and θ i are conjugate variables. A qualitative picture in the low energy, long wave length regime can be obtained from the Fourier transform of the coupling matrice, where k ∈ [−π, π]. The k → 0 behavior of M |i−j| is dominated by the exponential decay and for small k Λ −1 the Fourier transform of the interaction matrix M |i−j| is approximately constant, M (k) ≈ M (0), which leads to

B. Standard Villain transformation
The model (10) can be treated by the standard technique involving Villain approximation 15 . We omit all the details and only mention the fact that the spin wave part of the resulting model in the limit k Λ −1 , where U (k) ≈ M (0)k 2 is quadratic, corresponds to a Luttinger liquid 16 , where the Luttinger liquid velocity v and the Luttinger liquid parameter K scale with the original model param- The corrections to (15) due to vortices (phase slips) in this type of theories (see, e.g., Ref. 17) are of the form ∝ cos(2πp[m(x) + F (x)]) with p ∈ Z. The amplitude in front of this terms (fugacity of vortices) is predicted to be small in the limit Λ 1 14,15 so that without disorder the critical value of K is close to 2/π (it may be renormalized by disorder 16 ).
We assume that M (0) E t and therefore K 2/π. In this case the system is firmly in the charge density wave (CDW)-regime 16 and dominated by classical charge dynamics. The disorder F (x) enters the relevant phase slip terms ∝ cos(2π[m(x)+F (x)]) and can pin the charge density profile.

C. Alternative derivation
Here we generalize the derivation given in Ref. 8 for the case of a Josephson array described by the capacitance matrix (6) to the case of a more general matrix M i,j characterized by a screening length Λ. We start by rewriting the Hamiltonian (10) as where U 0 = U j,j and δU i,j = U i,j − U 0 δ i,j . Next we transform the third term with the help of the Hubbard-Stratonovich transformation, which introduces a new degree of freedom Q i , often referred to as the quasi-charge. This gives such that H = min Q [H{Q}]. Next where The Fourier transform reads For small wave vectors k Λ −1 we obtain U (k) U 0 and B(k) ≈ U (k) ≈ M (0)k 2 . Thus, assuming Q i changes slowly enough as a function of the coordinate i, i.e., changes on length scales longer than Λ, we can approximate For the Josephson arrays with the capacitance matrix (6) we obtain U 0 = 2E C and M (0) = 2Λ 2 E C = 2E C0 , where E C0 ≡ (2e 2 )/2C = Λ 2 E C . In this case the form of the third term of (22) is exact 8 .
The adiabatic dynamics of the model (22) without disorder was analyzed in Ref. 8. The inclusion of disorder is straightforward. The aim is to integrate out the degrees of freedom (m i , θ i ). For a given (adiabatic) trajectory Q i (t) the dynamics factorizes to independent dynamics of single junctions governed by the Hamiltonians The Born-Oppenheimer periodic potential is given by the ground state of the well known Hamiltonian (23), where Q i +F i serves here as the total quasicharge. The function E Q (Q) is periodic with period 1 as can be seen from (23).
Thus we obtain the effective potential energy of the whole array of the form This potential is supplemented by the kinetic energy.
The quadratic part of the Lagrangian T − U C gives again the Luttinger liquid with the parameters (16). Since we assume K 2/π, the periodic potential E Q (Q i + F i ) is relevant and pins the density profile. In what follows we investigate the charge pinning in this setup.

D. Edge bias
From now on we employ the terminology of Josephson junction arrays and put 2e = 1. As we are primarily interested in the transport properties of the chains, we introduce a bias V at the edge: To simplify the treatment in terms of the depinning theory we transform the system from a boundary biased situation to a spatially homogeneous driving by introducing a parabolic shift in Q and F , and the corresponding potential part of the Hamiltonian with a driving force V /N , In this formulation the problem corresponds to the discrete version of the well known depinning problem in onedimension 2 . The elastic energy of the field Q i is determined by the elastic constant C 0 . The elastic field is pinned by the random pinning potential E Q (Q i +F i ) and driven by the homogeneous driving force V /N . In the pinned regime the applied force V /N is not strong enough to overcome the potential barrier imposed on the elastic object Q i by the random pinning potential.
In the case that no driving force is applied, V = 0, the form of the elastic object is determined by a competition between the elastic term (Q i − Q i+1 ) 2 and the pinning term On small length-scales, where the elastic energy term dominates, Q i is approximately constant. The field Q i changes on large length-scales where the pinning potential dominates. The crossover between the two regimes happens at the length scale L p , which was first determined by Larkin for a flux line lattice in type II superconductors 4 . The length L p goes by many names depending on the physical systems that are pinned. In type II superconductors it is called Larkin length, in ferromagnets with domain boundaries Imry-Ma length 5 and for charge density waves it is called Fukuyama-Lee length 1 . In this work we use the term Larkin length.
Once the driving force V /N exceeds a critical force V cr /N , the pinning potential is overcome and the elastic object starts to move through the disordered medium. An intuitive argument to find the value of the critical driving force can be found by comparing the driving force to the pinning force at the Larkin length 4 . The distribution of Q is rigid on length-scales up to the Larkin length. The elastic object can only start to move when the driving force exceeds the collective pinning force on a segment with length L = L p .

E. Strong disorder
We first consider the strongly disordered model for which the results of the standard depinning theory 1,2 are directly applicable. To make the connection to these results the difference between the original disorder (f i ) and the effective quasi-disorder before (F i ) and after (F i ) the parabolic shift in the quasicharge is important. In terms of the original disorder the strongly disordered model is defined by, where p(f i ) is the probability distribution of the disorder f i and Θ H is the Heaviside Θ-function. This model corresponds to the strongest possible disorder in the considered chain model. A frustration f i with an absolute value larger than 1/2 is compensated by placing an additional (anti)-boson on the i-th island of the chain. The disorder is bounded by ±1/2 and a box-distribution of disorder-charges inside these boundaries is the maximal disorder. While f i itself is not spatially correlated, in the effective quasi-charge model, the quasi-disorderF i is correlated between different islands i and j, At first this seems to deviate from the normal situation in depinning theory where the disorder in the system is not spatially correlated 2 . However, in the depinning theory, only correlations in the pinning potential are relevant to the behaviour of the system. The potential E Q is a function of the quasi-charge with a periodicity of 1.
Since the disorder f i is box distributed in an interval that corresponds to the periodicity of the potential, the offsetF i can be absorbed into another uncorrelated boxdistributed disorder term f b i , From the point of view of the potential E Q , the quasi-disorderF i is equivalent to a spatially uncorrelated disorder term f b i in the maximally disordered model. Another way to determine whether spatial correlations inF i affect the quasi-charge model, is to calculate the disorder-averaged correlation function of the pinning potential: where the correlation function R(Q) is given by, Since the correlator of the pinning potential is proportional to a Kronecker delta, the pinning potential is not spatially correlated.
We have now seen that for the maximal disorder model we arrive at an effective model that conforms with the standard assumptions of depinning theory. In this case the Larkin length and the critical driving force are well known (see for example Ref. 1).
The approximate value of the Larkin length in onedimensional systems is given by 1 , The relevant parameters of the chain are the energy E C , the tunnelling amplitude E J , the chain length N and Λ.
To express the Larkin length in terms of these parameters we have defined the functionR, where E max Q is the amplitude of the random pinning potential E Q (Q i +F i ). The correlation functionR is a function of the dimensionless ratio of the tunnelling matrix element and E C . The function needs to be determined numerically only once for all possible values of chain length and C 0 .
Similarly the depinning force can be expressed in terms ofR and is given by 2 , Further corrections to this intuitive approach can be obtained from renormalization-group-approaches 2,19,20 . We use the approximation Eq.40 in this work. In most systems, where depinning theory is applicable, the system size is much larger than the Larkin length and it is a good approximation to assume infinite system size. We now turn our attention specifically to short finite chains. From Eq. 41 we see that the critical driving force decreases with increasing Λ. At the same time the Larkin length increases, In finite chains the Larkin length becomes equal to the system size N when Λ reaches the value, Increasing Λ further while keeping E C constant only increases the elastic energy E C0 coupling neighbouring islands. The Larkin length, the length-scale on which Q is approximately constant, should increase. However in this limit the Larkin length is equal to the system size and therefore the field Q is approximately constant along the whole array. For Λ Λ N the critical driving force is independent of the interaction length as long as E C is kept constant. A lower boundary for V cr is approximately given by which is the critical driving V cr one finds for Λ = Λ N . In reality V cr saturates for smaller Λ, when N is of the same order of magnitude as L p (for comparison see the numerical simulations in Sec. III). This leaves the principal behaviour of Eq. 44 unchanged and contributes a prefactor of order one in the expression for the critical driving force.

F. Weak disorder
In the weak disorder case the bare disorder f i is not evenly distributed in the interval [−1/2, 1/2]. One particular system we have in mind is a ladder of quantum phase slip (QPS) junctions 10,11 . In such a system, superconducting wires are arranged in a ladder configuration, such that a one-dimensional chain of superconducting loops is formed. The superconducting wires that are shared by neighbouring loops contain a very thin section that forms the QPS-junction. Magnetic flux quanta in the loops assume the role of the bosonic particles. The QPS-junctions between the superconducting loops provide the hopping matrix element and the coupling matrix M i,j is the inverse inductance matrix of the system. In a ladder configuration of superconducting wires the inductance matrix has the exactly the previously mentioned tridiagonal form Eq. 6, as long as the kinetic inductance dominates over the geometric inductance.
Due to the lack of large magnetic dipoles in the vicinity of such a system, a weak disorder limit is more likely to be realized than in Josephson junction arrays.
We consider two models of weak disorder: (i) the weak box disorder with the disorder strength γ < 1; (ii) Gaussian disorder with a standard deviation σ < 1/2. In the weak disorder models the spatial correlation in the quasi-disorderF i+1 can not be neglected in the argument of the effective potential E Q (Q). The maximal value of the disorder f i is smaller than the periodicity of the potential E Q and the long range correlation in the quasi-disorderF i can not be absorbed in the potential. The correlation function of the pinning potential therefore acquires a long range correlation component. We decompose the correlation function into short and longrange components, with the δ-correlated component R(Q) and the long range correlation function R 2 (Q, i, j). Due to the long range correlations the intuitive picture of the depinningtransition is not valid anymore. For a long range correlation function, that decays with a power law, the problem has been approached with the functional renormalization group method (FRG) in Refs. 12 and 21.
It has been shown 6,22 that these long-range correlations lead to the emergence of a new length-scale in the pinned system, the typical correlation length L corr . The roughness function w(x) of a pinned system shows a different behaviour, namely a variation in the roughness exponent ζ rough , depending on whether the system is probed at length-scales smaller or larger than the correlation length 6 . We derive typical correlation lengths for the two weak disorder models under the assumption that E Q can be approximated as a cosine-potential, To calculate the correlation function of the pinningpotential of two different chain sites j and k we set, without loss of generality, j < k. The correlation function in the weak box-disorder model is given by an integral over the disorder, wherep(F j ) is the probability distribution of the quasidisorder F j . Expressing the cosine in terms of exponentials one finds that the absolute value of the correlation function R 2 is bounded by an envelope function R E , The long-range correlation function decays exponentially with the distance k−j and the correlation of the pinningpotential decays on the length-scale, As expected the correlation length goes to zero in the limit of the maximal disorder and diverges in the clean limit without disorder, For a Gaussian distribution of the bare disorder f i , the correlation function is bounded by the exponential function R G , The correlation length is determined by the standard deviation σ of the bare disorder, We can again test the limits of infinitely broad and nondisordered distributions, In the broad limit the Gaussian disorder shows the same asymptotic behaviour as the box-disorder distribution when approaching the maximal disorder limit. The maximal disorder limit is consistent with a very broad bare Gaussian distribution. In the opposite limit the Gaussian distribution corresponds to a homogeneous shift in the definition of the quasi-charge and the correlation length diverges. The correlation length L corr marks the crossover between a disorder free and a strongly disordered array. On length-scales smaller than the correlation length the value of the disorder F i is approximately constant and constitutes a mere shift in the field Q. If the weakly disordered system is probed on these length-scales it behaves like a clean chain. On larger length-scales the value of the disorder changes significantly and the system behaves like a disordered chain. This transition is shown in the next section with the example of the dependence of the threshold voltage on the length of the chains.

III. SIMULATIONS
We obtain the critical driving force V cr by numerically solving the equations of motion of the field Q i that can be obtained from the Hamiltonian (28), The function V Q is the pinning force given by, To guarantee numerical convergence we have introduced a mass M and a linear dissipative term with a dissipation constant α R . Similar numerically simulations of the switching voltage in arrays of normal tunnel contacts have been conducted in Ref. 23. The critical driving force V cr is determined by adiabatically applying the boundary force V and determining whether a stable solution for the field Q i can be found. Although V is increased slowly, the switch-on time of the driving force V in the numerical simulation is finite. The phenomenological dissipative term has to be included to compensate the small transport velocityQ i introduced by the switch-on of V . The introduction of a phenomenological term is also a standard tool in the derivation of the depinning force V cr in renormalization-group-treatments of pinned systems 2 .
The mass M and the dissipation parameter α R both affect the dynamical properties of the system, however they have no influence on the breakdown of the static state. In the example of a Josephson junction array, the mass M corresponds to an inductance and α Q corresponds to an Ohmic resistance. In a quantum phase slip ladder M corresponds to a capacitance. In all simulations we choose the tunnelling amplitude and the coupling energy to be equal, E J = E C , so that the potential As long as the chain is more than twice as long as Λ, Vcr is independent of the length N and proportional to Λ. The critical driving force has the value predicted by an analytic estimate by Haviland and Delsing 24 . In the region where the chain is shorter than Λ the system is in the zerodimensional limit. The critical driving force is proportional to N and does not depend on Λ.
E Q is close to a cosine potential. The length of the chain N and the parameter Λ are varied.

A. The clean chain
We first simulate the clean model that has been used as the default model in a number of experimental papers on Josephson junction arrays 24,27,28 . While this model does not take into account charge disorder 29 , it might be more relevant for quantum-phase-slip-arrays than Josephson junction arrays as the former lack the strong charge disorder that can be found in the latter.
In the clean case, a simple argument to determine the critical driving force V cr can be found in Ref. 24. In the continuum limit (Λ 1) for long chains (Λ N ) the effective model of the clean chain is equal to the sine-Gordon model with a modified potential. The solutions of the standard sine-Gordon equation of motion are the well known solitons 8,26 , with the soliton velocity v. The spatial derivative of a static soliton v = 0 has a maximal value of The boundary driving force V takes the form of a boundary condition on the spatial derivative at x = 0, This can be used to estimate the maximal boundary force V for which a static soliton can exist at the array ends, In the Josephson junction arrays this force corresponds to the switching voltage at which the array switches from insulating to transport behaviour. The critical driving force does not depend on the array length and is proportional to the interaction length Λ. Both features are confirmed by the numerical simulations in Fig. 1 and Fig. 2.
In the limit Λ > N the spatial dependent field Q i takes the same value on all islands of the chain, Q i → Q and the coupled equations of motion simplify to a single equation of motion, The one-dimensional clean chain model reduces to a zerodimensional model. The critical driving force increases linearly with array size and is independent of Λ (Fig. 1). When the interaction length Λ is comparable to the inter-site distance Λ < 2 we are no longer in the continuum limit and the analytic approximation Eq. 70 is not valid. The switching-voltage is proportional to the length N and the Λ-dependence can be fitted to an exponential behaviour, as seen in Fig. 2. Only one set of fitting parameters β, γ is used for all four simulated chain-lengths. The change of the switching voltage behaviour can be understood in the following way. The interaction length Λ is a measure for the ratio of the elastic coupling between neighbouring islands and the depth of the pinning-potential. For small interaction lengths Λ < 2, Q-excitations can be created at the driven end of the chain without leading to complete depinning. During the adiabatic increase of force V the whole chain is filled with non-propagating Q-excitations. The depinning transition of these Q-excitations is determined by the Peierls-Nabarro-potential 25,26 . This give rise to a Λ-dependence of the form of Eq.73. In the context of Josephson-junction-arrays this was discussed by Fedorov et al. for the depinning of a single 2e-charge-excitation 26 .

B. The maximally disordered array
Here we present the critical driving force obtained from numerical simulations of the maximally disordered model. In Fig. 3 we compare the dependence of V cr on the parameter N with analytic estimate in Eq. 41. At large N , where the array is longer than the Larkin length N > L p , we find that the numerical simulations fit to the expected linear dependence on the system length.
For small system lengths the switching voltage does not increase linearly with N , as expected in the saturation regime where the Larkin length is comparable to the system size (Eq. 44).
The numerically determined dependence of V cr on Λ is shown in Fig. 4 and Fig. 5. For small Λ the inter- site distance is comparable to Λ and the continuum limit of the standard depinning-picture does not apply. For large Λ the Larkin-length is comparable to the chainlength N and we observe a saturation of V cr with Λ. The saturation sets in for, where α sat is of order of one. Comparing the analytic estimate Eq. 44 with the saturation points we expect α sat in the range 2.5 ≤ α sat ≤ 3.5.
In this intermediate regime we expect a power-law behaviour with an exponent of − 2 3 (Eq. 41). Fitting the numerical data to a power-law we obtain the exponents −0.49 ± 0.05 (N = 150) and −0.56 ± 0.03 (N = 195). However this is limited by the numerically accessible chain-lengths and we can not obtain a robust confirmation of the value of the exponent of Λ from the numerical simulations.

C. Weak disorder and emergent correlation length
To validate our analytic model of the introduction of a new length-scale by weak disorder, we also simulate the depinning-transition of the weakly disordered chain. We choose the disorder strengths γ = 0.25, L corr (γ = 0.25) ≈ 10 and γ = 0.125, L corr (γ = 0.125) ≈ 40. In Fig. 6 it is shown that the system undergoes a transition when the array-length is equal to the correlation length, N = L corr . Below N < L corr the chain is described by the clean chain model (γ = 0). Above the transition the critical driving force increases linearly with N . The N -dependence of V cr matches the maximally disordered model γ = 1. When the correlation length is significantly larger than the array size we can approximate all correlated disorder terms F i by a single value F i ≈ F . The perfectly correlated disorder term F can be absorbed into the definition of the quasi-charge and the system is equivalent to the clean array without disorder F i = 0.
When the length of the chain exceeds the correlation length one has to distinguish between two cases. The case when the correlation length is smaller than Λ and the case where it is larger. The first case requires a careful treatment to map the weakly disordered case to an effective strongly disordered model. Here we limit ourselves to the simpler second case. In this case one can understand the behaviour of the critical driving force with the following simple argument.
The typical length of a soliton Λ is smaller than the correlation length and static solitonic solutions of the field Q can exist in the chain. On the one hand this leads to the creation of a boundary soliton at the edge of the driven system that corresponds to the boundary soliton in the clean case. This gives rise to an offset critical driving force V offset cr . Since the chain is longer than L corr it can be subdivided into domains of length L corr . To switch into the conduction regime, the applied driving force needs to overcome the transport threshold in each domain, where the transport threshold is proportional to the critical driving force in the clean chain Eq. 70, This mechanism explains the linear increase in V cr seen in Fig. 6 for N > L corr .

IV. CONCLUSIONS
In this paper we have studied the depinning behaviour of discrete bosonic chain models that can be described by an effective Hamiltonian in the adiabatic limit that is similar to the sine-Gordon model. The most experimentally relevant realization of this model are linear arrays of Josephson junctions, however another possible realization is a ladder configuration of superconducting wires with quantum phase slip elements separating neighbouring superconducting loops.
We used analytical considerations and numerical simulations to determine the critical driving force required to overcome the pinning of bosons in the chain. In the parameter regime that corresponds to experimentally studied arrays we reproduce the recently observed behaviour 7 . Going to new parameter regimes, namely short chains and weakly disordered chains, we see a saturation regime in short chains where the Larkin length exceeds the system length and the critical driving force is independent of the decay length Λ of the repulsive interaction. In the weak disorder regime we observe the emergence of a new correlation length-scale L corr . Both effects show good agreement between the analytic results and the numerical simulations