Depinning and dynamics of vortices confined in mesoscopic flow channels

We study numerically and analytically the behaviour of vortex matter in artificial flow channels confined by pinned vortices in the channel edges (CEs). The critical current density Js for channel flow is governed by the interaction with the static vortices in the CEs. Motivated by early experiments which showed oscillations of Js on changing (in)commensurability between the channel width w and the natural vortex row spacing b0, we study structural changes associated with (in)commensurability and their effect on Js and the dynamics. The behaviour depends crucially on the presence of disorder in the arrays in the CEs. For ordered CEs, maxima in Js occur at commensurability w = nb0 (n is an integer), while for w ≠ nb0 defects along the CEs cause a vanishing Js. For weak disorder, the sharp peaks in Js are reduced in height and broadened via nucleation and pinning of defects. The corresponding structures in the channels (for zero or weak disorder) are quasi-1D n row configurations, which can be adequately described by a (disordered) sine-Gordon model. For larger disorder, matching between the longitudinal vortex spacings inside and outside the channel becomes irrelevant and, for w ≃ nb0, the shear current Js levels at ∼30% of the value Js0 for the ideal commensurate lattice. Around ‘half filling’ (w/b0 ≃ n ± 1/2), the disorder leads to new phenomena, namely stabilization and pinning of misaligned dislocations and coexistence of n and n ± 1 rows in the channel. At sufficient disorder, these quasi-2D structures cause a maximum in Js around mismatch, while Js smoothly decreases towards matching due to annealing of the misaligned regions. Near threshold, motion inside the channel is always plastic. We study the evolution of static and dynamic structures on changing w/b0, the relation between the Js modulations and transverse fluctuations in the channels and find dynamic ordering of the arrays at a velocity with a matching dependence similar to Js. We finally compare our numerical findings at strong disorder with recent mode-locking experiments, and find good qualitative agreement.

quasi-2D structures cause a maximum in J s around mismatch, while J s smoothly decreases towards matching due to annealing of the misaligned regions. Near threshold, motion inside the channel is always plastic. We study the evolution of static and dynamic structures on changing w/b 0 , the relation between the J s modulations and transverse fluctuations in the channels and find dynamic ordering of the arrays at a velocity with a matching dependence similar to J s . We finally compare our numerical findings at strong disorder with recent modelocking experiments, and find good qualitative agreement.

Introduction
The depinning and dynamics of the vortex lattice (VL) in type II superconductors is exemplary for the behaviour of driven, periodic media in the presence of a pinning potential [1]. Other examples range from sliding surfaces exhibiting static and dynamic friction and absorbed monolayers [2] to charge density waves (CDWs) [3], Wigner crystals [4] and (magnetic) bubble arrays [5]. Vortex matter offers the advantage that the periodicity a 0 of the hexagonal lattice can be tuned by changing the magnetic induction B. In addition, the effect of various types of pinning potentials can be studied. This pinning potential, arising from inhomogeneities in the host material, can be completely random, as in most natural materials, or can be arranged in periodic arrays using nanofabrication techniques [6,7]. In a variety of cases, correlated inhomogeneities occur naturally in a material, such as twin boundaries and the layered structure of the high-T c superconductors [8].
Depinning of the VL in a random potential generally involves regions of plastic deformations [9]- [13], i.e. coexistence of (temporarily) pinned domains with moving domains. For very weak pinning, the typical domain size can exceed the correlation length R c of the VL (see [11]) and the weak collective pinning theory [14] can be successfully used to estimate the critical current density J c [15,16]. However, as either the ratio of the VL shear modulus c 66 and the elementary pinning strength or the number of pins per correlated volume decreases, plastically deformed regions start to have a noticeable effect on J c . Recent imaging experiments [17] have shown directly that the rise in J c in weak pinning materials near the upper critical field B c2 , known as the peak effect [15,16], originates from such, rather sudden, enhancement of the defect density. This strong reduction of the VL correlation length is also accompanied by a qualitative change in the nature of depinning: for strong pinning, depinning proceeds through a dense network of quasistatic flow channels (filaments) such that the typical width of both static and moving 'domains' has approached the lattice spacing [9,10,18,19]. Depinning transitions via a sequence of static, channel-like structures have also been observed experimentally via transport experiments [20].
In superconductors with periodic pinning arrays (PPAs), matching effects between the lattice and the PPA become important. As shown first by Daldini and Martinoli [21,22], when the vortex spacing coincides with the periodicity of the potential, pronounced maxima can occur in J c , while at mismatch, defects (discommensurations) appear which gives rise to a reduced J c . In the last decade, many more studies of VLs in PPAs have appeared, both experimentally and numerically. Pronounced commensurability effects were found in films with 2D periodic pinning [6,7,23] for flux densities equal to (integers of ) the density of dots. In these systems, vortex chains at interstitial positions of the periodic arrays (e.g. at the second matching field of a square pinning array) can exhibit quasi-1D motion under the influence of the interaction with neighbouring, pinned vortices [24], as has also been observed in numerical simulations [25]. In addition, these simulations have revealed that, depending on the vortex interactions and the symmetry or strength of the PPA, a rich variety of other states and dynamic transitions can occur, often leading to peculiar transport characteristics.
Besides the above examples, the phenomenon of vortex channelling can also arise from the presence of grain boundaries in the sample. Historically, the 'shear' depinning of vortices in grain boundaries in low T c materials received considerable attention [26,27] because it could explain the quadratic decrease of J c near B c2 in practically relevant polycrystalline superconductors. More recently, the interesting issue of channelling of mixed Abrikosov-Josephson vortices in grain boundaries in high-T c materials was addressed in detail by Gurevich [28]. We also mention the recent interest in channel structures in superconductors as possible 'fluxon' rectifiers [29]. A system in which channel motion and its dependence on the structural properties of vortex matter can be studied systematically is that of narrow, weak-pinning flow channels in a superconducting film [30], see figure 1. The samples are fabricated by etching straight channels of width w etched 100 nm through the top layer of an a-NbGe/NbN double layer. With a magnetic field applied perpendicular to the film, vortices penetrate both the strong-pinning NbN in the channel edges (CEs) and the remaining NbGe weak-pinning channels. The strongly pinned CE vortices provide confinement to the vortices inside the channel, as well as the pinning (shear) potential which opposes the Lorentz force from a transport current J applied perpendicular to the channel. By changing the applied field H, one can tune the commensurability between the VL constants and the channel width, allowing a detailed study of the shear response and threshold for plastic flow as a function of the mismatch and the actual microstructure in the channel.
Phenomenologically, plastic flow in the channel occurs when the force density F = JB (with B/ 0 the vortex density) exceeds 2τ max /w, where τ max = Ac 66 is the flow stress at the edge (the factor 2 is due to both CEs) and w is the effective width between the first pinned vortex rows, defined in figure 1. Thus, the critical force density is given by The parameter A describes microscopic details of the system: it depends on lattice orientation, (an)harmonicity of the shear potential, details of the vortex structure in the CEs and the microstructure of the array inside the channel. Critical current measurements as a function of the applied field reflected this change in microstructure through oscillations of F s , shown in figure 2 for a channel with w etched ≈ 230 nm. Note that in such a narrow channel, the pinning strength due to intrinsic disorder in the a-NbGe is at the most 10% of F s (except for B 50 mT) and does not affect the oscillations. Since the natural row spacing of the VL is b 0 = √ 3a 0 /2, with a 2 0 = 2 0 / √ 3B, and in our geometry B µ 0 H, one can check that the periodicity of the oscillations corresponds to transitions from w = nb 0 to w = (n ± 1)b 0 with n integer, i.e. the principal lattice vector a 0 is oriented along the channel (figure 1). The envelope curve represents equation (1) with Brandt's expression for the VL shear modulus [31]: (b = B/B c2 is the reduced field and λ is the penetration depth) and a value A = 0.05. This value for A is close to the value u 2 /a 0 = 0.047 for the relative displacements at the crossover from elastic to plastic deformations as obtained from measurements on the peak effect [16,17]. This led to a qualitative interpretation of the reduction of F s at minima as being due to defects in the channel, which develop at incommensurability. However, recent developments [32]- [35] have shown that (strong) structural disorder may be present in the CE arrays, in which case the interpretation can drastically differ.
In this paper, we present numerical and analytical studies of the threshold force and dynamics of vortices in the channel system for various degrees of edge disorder. In an earlier paper [36], we studied the commensurability effects in the idealized case with periodic arrays in the CEs. In this situation, F s at matching (w = nb 0 ) is equal to the ideal lattice strength 2A 0 c 66 /w (the value A 0 = 1/(π √ 3) follows from Frenkels considerations [37]), while at mismatch dislocations develop, leading to A 0. The resulting series of delta-like peaks in F s versus matching parameter differed considerably from the experimental results, which could not be explained by thermal fluctuations or intrinsic disorder inside the channel. Therefore, we investigated the effect of positional disorder in the CE arrays on F s near commensurability (w ≈ nb 0 ) [34]. In this regime, the behaviour is dominated by the longitudinal displacements of vortices in the chains, i.e. quasi-1D, and F s is controlled by defects with Burgers vector along the channel. At weak disorder, we found a clear reduction of F s at commensurability caused by nucleation of defects at threshold, while the existing defects at incommensurability become pinned by disorder, leading to an increase of F s in the mismatching case.
The present paper first describes in detail these quasi-1D phenomena near commensurability and/or for weak disorder. Using a generalized sine-Gordon (s-G) model, we quantitatively describe how the structure and transport properties depend on the vortex interaction range and on weak disorder in the CEs. Besides the connection to our system, these results also provide a background for understanding quasi-1D vortex states and matching effects in artificial PPAs, including the effects of disorder which these PPAs may contain due to fabrication uncertainties.
The 1D model shows that, above a certain disorder strength, spontaneously (disorder) induced defects along the CEs dominate over incommensurability induced defects. The commensurability peak in F s is then completely smeared out with a value of F s at matching (w = nb 0 ) saturating at ∼30% of the ideal lattice strength. In the more general case of wider channels, the transverse degrees of freedom, especially away from matching (w/b 0 n ± 1/2), lead to new phenomena: under the influence of disorder, the channel array may split up into regions with n and n ± 1 rows, involving dislocations with Burgers vector strongly misaligned with the CEs. At sufficient disorder strength, such dislocations lead to a more effective pinning of the array than the 'aligned' dislocations around matching. F s then exhibits a smooth oscillation as a function of w/b 0 , similar to figure 2, with yield strength maxima occurring around mismatch. This behaviour resembles the classical peak effect, i.e. at mismatch, the enhanced ability of the arrays to sample configuration space allows better adjustment to the random CEs. In the last part of the paper, we show detailed simulations of both static and dynamical aspects of this behaviour, including a study of reordering phenomena at large drive. We find an ordering velocity of the arrays with a channel width dependence similar to that of the threshold force. Using a modified version of the dynamic ordering theory in [38], it is shown that such a behaviour can be explained by a reduction of the energy for formation of misaligned defect pairs away from matching. The numerical results at strong disorder are also in good qualitative agreement with recent mode-locking (ML) experiments on the channel system [32,33]. The outline of the paper is as follows. In section 2 the channel geometry and the simulation procedure are discussed. The first part of the paper deals with channels having hexagonal, ordered arrays in the CEs. In section 3 we present the s-G description and numerical results for a single 1D vortex chain in an ordered channel. In section 4 we show how the 1D behaviour extends to wider channels with multiple rows and ordered CEs. The second part of the paper deals with disordered channels. Section 5 describes the effects of weak CE disorder on the behaviour of a 1D chain, both analytically and using numerical simulations. The effects of weak disorder in wider channels are discussed in section 6. Section 7 describes the static and dynamic properties of wider channels in the presence of strong disorder, including an analysis of the reordering phenomena in this situation. A comparison with the dynamic ordering theory, the confrontation with experiments and a summary of the results are presented in sections 8 and 9.

Model and numerical procedure
We consider straight vortices at T = 0 in the geometry as illustrated in figure 3 for the case of 1 row in the channel. The approximation T = 0 is well justified over a considerable range of experiments (see section 8). The CEs are formed by two semi-infinite static arrays. The distance between the first vortex rows on both sides of the channel is w + b 0 , with w the effective channel width. The vortices are assumed to be fixed by columnar pins in the CEs. The principal axis of the pinned arrays is along the channel direction x. A relative shift x is allowed between the arrays. In figure 3(a) the simplest configuration is shown, where CE vortices form a perfect triangular lattice. For x = 0, their coordinates are for m = 0 and frac(m/2) denotes the remainder of m/2. Disorder is incorporated in the model by adding random shifts d to the coordinates of the ordered arrays: as shown in figure 3(b). The amplitudes of the random shifts are characterized by disorder parameters x and y as follows: transverse relative displacements d y /a 0 are chosen independently from a box distribution [− y , y ]. The longitudinal shifts d x n,m are chosen such that the strain (d n+1 − d n )/a 0 along the channel is uniformly distributed in the interval [− x , x ]. The latter provides a simple way of implementing loss of long-range order along the CEs. For x and y , we study the following specific cases: x , y = 0 in sections 3 and 4, x ≡ , y = 0 in section 5 and x = y ≡ in sections 6 and 7.
To study the commensurability effects, the effective width of the channel is varied from a value w/b 0 ∼ 1-10. We assume that the vortex density inside and outside the channel are the same. The number of vortices in the channel is then given by N ch = (L/a 0 )(w/b 0 ), with L the channel length. In a commensurate situation, one has w = pb 0 and both the row spacing and (average) longitudinal vortex spacing in the channel match with the vortex configuration in the CEs. When w = nb 0 these spacings become different, leading to generation of topological defects. While this model differs from the experimental case where the applied field drives the incommensurability, the method offers a simple way of introducing geometrical frustration and study various (mis)matching configurations.
With a uniform transport current J applied perpendicular to the channel, the equation of motion for vortex i in the channel reads (in units of N m −1 ): V(r) is the vortex-vortex interaction potential, j labels other vortices inside the channel, the damping parameter γ is given by γ = B 0 /ρ f with ρ f the flux flow resistivity and f = J 0 is the drive along the channel. For films which are not too thin compared to the penetration depth λ and magnetic fields small compared to the upper critical field B c2 , the interaction V(r) is given by the London potential: where U 0 = 2 0 /2πµ 0 λ 2 and 0 is the flux quantum. In the simulations, we integrate equation (5) numerically for all vortices in the channel. We use a Runge-Kutta method with variable time steps such that the maximum vortex displacement in one iteration was a 0 /50. Distances were measured in units of a 0 (r = r/a 0 ), forces in units of U 0 /a 0 and time in units of γa 2 0 /U 0 . Following [38], the London potential was approximated by with a cut-off radius r c corresponding to r c 3λ/a 0 . We performed most simulations for r c = 3.33. Periodic boundary conditions in the channel direction were employed. For each w/b 0 , we relaxed the system to the ground state for f = 0. We found that this is best achieved by starting from a uniformly stretched or compressed n or n ± 1 configuration. For an initial configuration with N ch vortices distributed randomly in the channel, relaxation resulted in (slightly) metastable structures, even when employing a finite temperature, simulated annealing method. Some peculiarities associated with such structures are mentioned in sections 4 and 6.
After the f = 0 relaxation, the average velocity versus force (v-f ) curve was recorded by varying the force stepwise from large f max → 0 (occasionally f = 0 → f max → 0 was used to check for hysteresis). At each force we measured v(f ) = ẋ i i,t after the temporal variations in v became <0.5% (ignoring transients by discarding the data within the first 3a 0 ). In addition, at each force we measured several other quantities, e.g. the temporal evolution of r i and the time-dependent velocity v(t) = ẋ i i .

Single chain in an ordered channel
The first relevant issue for plastic flow and commensurability effects in the channel is to understand the influence of periodically organized vortices in the CEs (see figure 3(a)). The characteristic differences between commensurate and incommensurate behaviour can be well understood by focusing on a 1D model in which only a single vortex chain is present in the channel. Therefore, the CEs are assumed to be symmetric with respect to y = 0 (i.e. x = 0 in figure 3) and only the longitudinal degrees of freedom of the chain are retained. At commensurability, w = b 0 , the longitudinal vortex spacing a = a 0 . For w = b 0 , the average spacing a = 0 /(Bw) = a 0 b 0 /w does not match with the period a 0 in the edges and interstitials or vacancies developed in the channel. Their density c d is given by

Continuum s-G description
We first consider the interaction of a vortex in the channel with the periodic arrays in the CEs. As shown in appendix A, when B 0.2B c2 and λ a 0 the edge potential arising from this interaction is where k 0 = 2π/a 0 . For w = b 0 and y = 0, the associated sinusoidal force caused by the edge has an amplitude which we denote by µ: Next we consider a static chain of vortices inside the channel. The chain is most easily described in terms of a continuous displacement field u(x), representing the deviations of vortices in the chain with respect to the commensurate positions, i.e. u(ia 0 ) = u i = x i − ia 0 . The edge force is then given as f p = −µ sin(k 0 u). To describe the interaction between vortices within the chain, we assume that their relative displacements are small, ∂ x u 1. Then one can use linear elasticity theory. Taking into account that the interaction range λ > a 0 , the elastic force at x = ia 0 is Using the Fourier transform of V , the force due to a displacement u q (x) = Re(u q e iqx ) with wave vector q is Recasting this into a sum over reciprocal vectors lk 0 ± q and retaining only the l = 0 term, one obtains the following dispersive elastic modulus of the chain: For deformations of scale >2πλ, the elastic force is f el = κ 0 ∂ 2 x u with a long wavelength stiffness κ 0 = U 0 π(λ/a 0 ).
The equation of motion for u for a uniformly driven chain is obtained by adding the driving force f to the edge force and the intra-chain interactions resulting in γ∂ t u = f + f p + f el . Assuming for the moment that the long wavelength description is valid, the evolution of u is given by the following s-G equation: A useful visual representation of equation (13) is an elastic string of stiffness κ 0 with transverse coordinate u(x) in a tilted washboard potential (µ/k 0 ) cos(k 0 u) − fu. The s-G equation (13) has been thoroughly studied in different contexts (e.g. [39]- [41]). In the static case (f = 0), it has the trivial solution u = 0, corresponding to a commensurate chain, or kinked, incommensurate, solutions in which u(x) periodically jumps by ±a 0 , each jump representing a point defect in the channel. In the context of long Josephson junctions (LJJs, [39]), a kink corresponds to a Josephson vortex where the phase difference across the junction changes   (14) for λ/a 0 = 1, 2 and 9 (most extended line). Symbols: numerically obtained displacement field for an isolated interstitial for the corresponding r c . ----, defect shape (18) in the nonlocal limit.
by 2π. An isolated defect is represented by the familiar 'soliton' solution of the s-G model: where x c denotes the centre of the defect and the ± sign denotes a vacancy or interstitial (kink or anti-kink). The length l d represents the core size of the defect: with g the dimensionless ratio between the chain stiffness and maximum curvature of the pinning potential: as follows from equations (9) and (12). For λ/a 0 1, l d thus considerably exceeds the lattice spacing. The continuum approach is validated since ∂ x u d 2a 0 / l d 1. In figure 4 we have illustrated the characteristic defect shape (14), along with numerical data from a later section.
The long wavelength limit is only valid when l d considerably exceeds λ. Since l d grows only as √ λ/a 0 , the dispersion in the elastic interactions becomes important beyond a certain value of λ/a 0 . This value is estimated by setting λq d = 2πλ/ l d = 1 in (12), resulting in λ/a 0 9, in which case l d 54a 0 .
For a larger interaction range, one employs the following approach, first derived by Gurevich [28] for mixed Abrikosov-Josephson vortices in grain boundaries. Expression (10) for the elastic force can be written as an integral f el = (ds/a 0 )∂ 2 s V(s)u(x + s). For defects on a scale <λ, only the short distance behaviour of V has to be retained: Integrating the expression for f el by parts and adding the edge force and the drive, the equation of motion becomes: A static solution of equation (17) for a single defect is [28] k 0 u = π + arctan(±2πx/ l nl d ), with l nl d = 6π 2 a 0 , which is valid when λ > l nl d . The value for l nl d is nearly the same as the s-G core size l d for λ/a 0 = 9. This means that upon approaching the non-local regime, the increase of core size saturates at ∼60a 0 , while only the tails of the defect are affected according to equation (18), see figure 4. A more accurate calculation of the onset of the non-local field regime using Brandt's field-dependent vortex interaction (appendix A) shows that non-locality is only relevant for a channel in a superconductor with λ/ξ 50.
So far, we discussed isolated defects. For finite defect density, the repulsive interaction between defects of the same 'sign' causes a periodic superstructure in the chain. When c d grows to ∼1/ l d , the defects start to overlap significantly. For the (local) s-G model, explicit solutions for u have been obtained in terms of the Jacobi elliptic functions, for which we refer to [39,41]. Recently, also in the non-local limit where l nl d > λ, the 'soliton'chain has been described analytically [28], which we will not repeat here.

Transport properties
With a uniform drive f , the transport properties strongly depend on the presence and density of defects in the channel. At commensurability (a = a 0 , c d = 0), a threshold force f s = µ is required, above which all vortices start moving uniformly. Their velocity is identical to that of an overdamped particle in a sinusoidal potential: v = f 2 − µ 2 [42]. The threshold µ coincides with the well-known relation between shear strength and shear modulus of an ideal lattice by Frenkel [37]: for a harmonic shear interaction, a value A = A 0 = a 0 /(2πb 0 ) = 1/π √ 3 is applied in equation (1). Identifying F s a 0 b 0 = f s = µ for w = b 0 , one finds: which coincides with the familiar expression for the shear modulus in the London limit: In appendix A we generalize the expression for the ordered channel potential to higher field and show that in that case also the potential is harmonic and that A = A 0 holds for a commensurate channel. At incommensurability, depinning of the chain is governed by the threshold force to move a defect. In the present continuum approach, such a threshold is absent. However, taking into account the discreteness of the chain, in which case equation (13) turns into a Frenkel-Kontorova (FK) model [36], a finite Peierls-Nabarro (PN) barrier exists to move a defect over one lattice spacing (see e.g. [41]). The magnitude of the PN barrier has been studied for a variety of cases, including FK models with anharmonic and/or long-range interactions [41,43]. For g < 1, f PN can amount to a considerable fraction of µ. Additionally, in this regime, anharmonicity may renormalize g [41,43] and cause pronounced differences between the properties of kinks (vacancies) and antikinks (interstitials). In our limit g 1, where ∂ x u 1 and harmonic elastic theory applies, these differences are small and the pinning force vanishes as f PN = 32π 2 gµ exp (−π 2 √ g). Hence, defects in an ordered channel give rise to an essentially vanishing plastic depinning current J s . 6 Considering the dynamics, for small drive f < µ, the motion of defects, each carrying a flux quantum, provides the flux transport through the channel. When defects are well separated, for c d < l −1 d , the mobility of the chain is drastically reduced compared to free flux flow and the average velocity v is proportional to the defect density: v = c d v d a 0 . Here v d is the velocity of an isolated defect at small drive. It can be calculated from the general requirement that the input power must equal the average dissipation rate: The last step arises from the space and time periodicity of u. Using ∂ t u = v d ∂ x u and the kink shape equation (14), one obtains the 'flux flow resistivity' at small defect density: where π 2 √ g/2γ = M d is the kink mobility in the s-G model [40]. For larger defect density, where defects start to overlap, this relation changes. The linear response for f µ may then be obtained from the solutions for u based on elliptic integrals [28,39,41].
For larger drive f µ, the 'tilt'-induced reduction of the (washboard) edge potential becomes important. This leads to an expansion of the cores of the sliding defects and causes a nonlinear upturn in the v-f curves. Exact solutions of equation (13) describing this behaviour do not exist. Therefore, we use a perturbative method similar to that in [22,44] which is able to describe the full v-f curve over a wide range of defect densities. It is convenient to define the displacements h(x, t) = u(x, t) − s(x, t), where s(x, t) = (q/k 0 )x + vt, with (q/k 0 ) = c d a 0 , is the continuous field describing the displacements of an undeformed incommensurate chain (i.e. straight misoriented string in the washboard potential) moving with velocity v. In terms of h, the equation of motion (13) and equation (20) can be written as The last term in equation (23) The reduced stiffness is and V ce is obtained from the full (λ-dependent) expression for V ce,0 in appendix A (equations (A.3) and (A.4)). The result is that the continuum approach is valid for λ/a 0 0.3. For the frequently occurring geometry of a vortex chain confined by vortices pinned in a square array, one can also estimate the continuum regime. In this case one replaces (w + b 0 )/2 in equation (8) by a 0 /2, yielding a periodic pinning force with µ 10 µ. Correspondingly, g λ/a 0 and lattice discreteness is unimportant for λ/a 0 2. The overlap of defects and the core expansion for f µ appears in the q and v dependence of h. Both effects cause a reduction of the relative displacements h. An approximate solution for h(v) is obtained by substituting equation (24) into equation (13), yielding the coefficients h m (details of the solution are deferred to appendix B).
The v-f relation equation (23) attains the form: The additional 'friction' force is represented in terms of the pinning frequency ω p = µk 0 /γ, the washboard frequency ω 0 = k 0 v and ω r = K 2 eff (c d )/γ which is the effective relaxation frequency for nonlinear deformations associated with a defect density c d = q/(2π), with K 2 eff (c d ) given in appendix B. At small v, the elastic relaxation time 1/ω r for the chain to relax is much smaller than the timescale 1/ω 0 between passage of maxima in the edge potential. This corresponds to the linear sliding response of the static structure of (overlapping) defects. For large v, 1/ω r 1/ω 0 meaning that the incommensurate chain is not given enough time to deform. This leads to expanded defects described by a sinusoidal variation of h with reduced amplitude (see appendix B).
The v-f curve then approaches free flux flow according to f − γv ∼ v −1 as for a single particle.
Recently, exact solutions describing the nonlinear dynamics and core expansion of mixed Abrikosov-Josephson vortices based on the nonlocal equation (17) have been derived in [28]. The resulting transport curves are very similar to those obtained from equation (25), see figure 5. We also note the similarity with the I-V curves obtained from a model for kinked Josepson strings [45] in high-T c superconductors with the field under an angle with respect to the insulating layers.

Numerical results
The simulations of symmetric channels ( x = 0) for w ∼ b 0 fully support the above findings. The interaction with the CEs for w = b 0 provides a maximum restoring force with a value 0.054, independent of the interaction cut-off r c used in the numerics. This value is in agreement with the dimensionless values for µ and c 66 in equations (9) and (19): a 0 µ/U 0 = 1/(6π) and ). The data points in figure 4 show the displacement field of a single defect (obtained by adding one vortex to a commensurate chain) for three values of the cut-off r c (i.e. various λ/a 0 ). We conclude that up to r c = 30 (λ/a 0 = 9), the s-G kink shape equation (14) forms a good description of a defect in the chain.
The data points in figure 5 show numerical results for the transport of a single chain in channels of various widths and r c = 3.33. The features discussed previously, i.e. the vanishing PN barrier and nonlinear transport, clearly appear in the data for incommensurate chains. We also plotted the results according to equation (25), with K 2 eff (c d ) evaluated using the results in appendix B for λ/a 0 = 1 and taking into account that µ slightly depends on w. The analytical treatment gives a very reasonable description of the data. Finally, we show in figure 6 the numerical results and analytical results of appendix B for the quasi-static and dynamic shape of the chain for w/b 0 = 0.97 (c d = 0.03/a 0 ). The numerical results closely mimic the analytic results, both for the kinked shape at small v and the core expansion with the associated reduction of h for large v.
To conclude this section we mention that, at incommensurability, due to the vanishing barrier for defect motion, the average velocity ẋ i (t) i has a vanishing ac-component. Only at commensurability, the washboard modulation is retained, the velocity at large drive being v(t) = v + (µ/γ) sin(ω 0 t).

Ordered CEs and multiple chains
We now turn to the results for channels containing multiple vortex rows and ordered CEs. The simulations are performed with the full 2D degrees of freedom and r c = 3.33. We implemented an edge shift x(w) with a saw tooth shape (0 x a 0 /2). This ensures that, as we vary w, a perfect hexagonal structure is retained for w = pb 0 with p an integer. However, for w = pb 0 , the qualitative behaviour did not depend on x. Figure 7 shows v-f curves of commensurate channels, w/b 0 = n with integer n 2. In these cases the arrays are perfectly crystalline and have a shear strength f s = µb 0 /w, inversely proportional to the channel width and in accord with equation (1) . This is consistent with the fact that only the first mobile chains within a distance ∼b 0 from both CEs experience the periodic edge potential (see equation (A.3)) while the other chains provide an additional pulling force via the elastic interaction. This interaction brings an additional feature to the dynamics, namely shear waves (see also [46]). The shear displacements of rows n in the bulk of the channel can be described in continuum form, u n (t) → u(y, t), by the following equation of motion: At large v the CE interaction can be represented by oscillating boundary conditions. As shown in appendix C, this causes an oscillatory velocity component dh/dt with y-dependent amplitude and phase describing periodic lagging or advancing of chains with respect to each other: Here ω 0 is the washboard frequency k 0 v, f(y) = cos(y/ l ⊥,v ) cosh(y/ l ⊥,v ) and g(y represents the distance over which the amplitude and phase difference decay away from the CEs. Although, in principle, equation (27) is valid only for γv/µ 0.25, it provides useful qualitative insight into the dynamics at all velocities: at small velocity, l ⊥,v is large, meaning that for all rows the velocity modulation and phase become similar. Hence, for v → 0 the array may be described as a single vortex chain, which is the underlying origin of the fact that close to threshold the curves approach the 1D commensurate behaviour v = f 2 − f 2 s with reduced threshold f s = µ/n. At large velocity, l ⊥,v eventually becomes less than the row spacing. In that limit only the two chains closest to the CE experience a significant modulation. In appendix C we quantitatively analyse the friction force in this regime with the result: In the inset of figure 7, this behaviour is displayed for n = 9 by the dotted line. In the highvelocity regime, the result agrees well with the numerical data, at lower velocities equation (28) underestimates the true friction. Next we discuss the behaviour of incommensurate channels. The static vortex configuration for a channel of width w/b 0 = 3.92 is shown in the upper part of figure 8. A Delauney triangulation shows that the array consists of four rows with two pairs of 5-, 7-fold coordinated vortices at the CE constituting two misfit dislocations of opposite Burgers vector b and glide planes along x. Due to their mutual attraction, dislocations at the upper and lower CE are situated along a line with an angle of ∼60 • with x. The two edge dislocations thus form a 'transverse' stacking fault (TSF). In the lower part of figure 8, the structure for a channel with w/b 0 = 3.52 is shown. Here the density of stacks, given by c TSF = |(1/a 0 ) − (1/a)|, is enhanced. The dislocations at one side of the CE repel each other and are equally spaced, like the periodic superstructure for a single chain in figure 6. The slight misalignement between the 'upper' and 'lower' dislocations of a pair is due to the relative shift between the CEs: the exact orientation of the pairs is determined by the choice of x. For channel widths in the regime n < w/b 0 n + 1/2 with integer n, we find very similar structures but instead of TSFs consisting of vacancies, we now have TSFs consisting of interstitial vortices, again arranged in a periodic superstructure.  (29) for λ/a 0 = 1. Figure 9 shows the transport curves associated with these structures. As for the single chain, the presence of misfit defects causes an essentially vanishing threshold force. For a small drive, f < µb 0 /w, a low-mobility regime occurs associated with glide of the edge dislocation pairs along the CE. This allows for the elastic motion of a complete TSF, i.e. the vortices in the 'bulk' of the channel remain 6-fold coordinated.
It is interesting to study how the mobility due to the TSFs changes on increasing the number of rows. In figure 10 we plot the mobility per stack, M TSF = (dv/df ) f →0 /c TSF versus channel width. M TSF around each peak decreases with increasing c TSF . This is caused by overlap of the strain fields of the defects, in analogy with the behaviour for a single chain. The overall increase of the peak value of M TSF is related to a change in the size of an isolated TSF. An extension of the analysis in section 3 allows to describe this change quantitatively. For small n, the longitudinal deformations do not vary strongly over the channel width. This can be understood by considering shear and compression deformations related by the equation κ∂ 2 x u x + c 66 a 0 b 0 ∂ 2 y u x = 0. It follows that a longitudinal deformation on a scale l along the channel varies over a scale l ⊥ = l √ c 66 a 0 b 0 /κ perpendicular to the channel. In the case l ⊥ w, the transverse variation of u x (y) is small and can be neglected so that κ 0 in equation (13) can be replaced by an effective stiffness nκ 0 due to n rows. Similarly, the driving force is replaced by f → nf . This results in the same equation (13) with a rescaled edge force µ → µ/n. Accordingly, the longitudinal size of a defect (TSF) is given by l TSF = 2πa 0 √ ng and the mobility of an isolated TSF by (compare M d below equation (21)): As shown by the drawn line in figure 10, this form gives a reasonable description of the data up to n = 3. Working out the condition l ⊥ w given above for the validity of equation (29), one obtains w (l d /2) √ c 66 a 0 b 0 /κ 3b 0 , in agreement with the data. At larger n, 'bulk-mediated' elasticity [47] leads to decay of the longitudinal deformations towards the channel centre. We also note that, due to the increase of l TSF with n, the density c TSF for which defects are non-overlapping, decreases on increasing n.
In the v-f curves of figure 9, we observe at larger velocity, features very similar to the transport of the 1D chain: for f µ/n, the effective barrier is reduced, leading to core expansion of the TSFs. Accordingly, the curves approach free flux flow behaviour. As in the commensurate case, this approach is initially slower than f − γv ∼ 1/v due to additional oscillating shear deformations in the channel for f µ/n. Figure 11 summarizes the behaviour of the shear force f s , taken at a velocity criterion v ≈ 0.01µ/γ, versus the matching parameter. At integer w/b 0 = n, the threshold is f s = µ/n, but we note that it can be reduced due to a finite edge shift x. At mismatch, f s is essentially vanishing. Near 'half filling', w/b 0 n ± 1/2, where the arrays switch from n to n ± 1 chains, a small enhancement of f s is observed. In this regime, the static (f = 0) structure was obtained by annealing from a random initial configuration, attempting to determine the exact switching point. This results in metastable structures with coexisting n and n ± 1 row regions (or longitudinal stacking faults, LSFs) bordered by dislocations with misoriented Burgers vector, see [36]. The increase in f s is caused by the finite barrier for climb-like motion of these dislocations, by which an LSF can move as 'giant' defect through the channel. For sufficiently large drive, (part of ) the LSFs are annealed which may result in hysteresis for up/down cycled v-f curves. We will discuss these 'mixed' n and n ± 1 structures in more detail in sections 6 and 7 in the context of disordered CEs. We also note that the integer chain structures with TSFs away from half filling differ from the results in [36]. The structures there, obtained from a random initial configuration, contained point defects unequally distributed among rows, yielding 'gliding' dislocations within the channel. Such structures are also slightly metastable but the conclusion of vanishing f s for incommensurate, integer chain structures, drawn in [36], remains unaltered.

Single chain in a disordered channel
We will now consider the influence of disorder in the CE arrays on transport in the channels, focusing in this section on the characteristics of a single chain for w/b 0 ∼ 1 with only longitudinal degrees of freedom. The CE disorder is implemented with longitudinal random shifts as described in section 2. We note that both CEs remain 'in phase'; the effect of quenched phase slips or dislocations between the CEs will be treated in the discussion in section 8.

Disordered s-G equation
First we consider the form of the channel potential in the presence of a weak disorder. To that end we generalize equation (A.1) in appendix A and express the CE potential at r 0 = (x, y = 0) in terms of the vortex density ρ e in the CEs: with ρ e (k) the Fourier transform of ρ e . For weak disorder (∇ · d 1), this density can be expressed in terms of the displacement field d in the CE as follows [48]: ρ e (r e , d) ( where δρ e = i cos [K i (r e − d(r e ))] represents the microscopic modulation due to the lattice (K i spans the reciprocal lattice) while ∇ · d reflects density modulations. As described in appendix D, this decomposition of ρ e leads to two contributions to the potential: where in the second term δµ(x)/µ = π √ 3∂ x d. The term V l represents long-range potential fluctuations and is smooth on the scale ∼a 0 . Its correlator l (s) Assuming that ∂ x d has short-range correlations (on the scale ∼ a 0 /2) and a variance (∂ x d) 2 = 2 /3 as in the simulations, l can be written as The exponent α and the prefactor C α depend on the disorder correlations between rows in the edge: α = 2 when the strain ∂ x d(x) is identical for all rows and α = 1 when the strain is uncorrelated between rows. The term V p in equation (31) is the quasi-periodic potential arising from δρ e of the vortex rows nearest to the CEs. The amplitude fluctuations δµ/k 0 are characterized by (see appendix D): To obtain the energy of the vortex chain and the equation of motion, the vortex density inside the channel, ρ c , is decomposed similar to ρ e : where u is the displacement field of the chain. As shown in appendix D, in the limit λ > a 0 , the resulting interaction with the CEs can be written as represents the s-G functional for an ordered channel, and H a , H s are the disorder contributions due to amplitude fluctuations and random coupling to the strain: The contains contributions from local and non-local strains. The latter dominates for λ > a 0 (see appendix D). Hence s (s) = V s (x)V s (x + s) l (s). The model described by H = H SG + H a + H s is also used to describe LJJs or commensurate CDWs with weak disorder, however with different disorder correlations. In the former case, the term H a in equation (34) describes local variations in the junction critical currents [49,50]. For CDWs, a disorder contribution of the form H a arises from the so-called backward scattering impurities, while the term H s originates from 'forward' scattering impurities [51]. We also note that our model differs from the usual Fukuyama-Lee-Rice model for CDWs [52], in which commensurability is ignored either due to strong direct random coupling to u (δµ(x) µ) or due to large mismatch.
In principle, the equation of motion for the chain is given by γ∂ t u = −δH/δu. However, it has been shown in previous studies [47], [53]- [55] that in the moving state a convective term −γv∂ x u should be included. While irrelevant for the depinning process, such a term can be important for the dynamics, and for completeness we include it. 7 The resulting equation of motion is In writing (35) we have assumed, for simplicity, that the elastic deformations in the presence of disorder can be described by the long-wavelength stiffness κ 0 . Ignoring the last term, (35) describes the transverse displacements u(x) of an elastic string in a tilted 'washboard' potential with random amplitude µ(x)/k 0 and random phase φ(x) = − x −∞ dx V s (x )/κ 0 . The latter represents a u-independent random deformation of the chain.

Numerical results
The influence of disorder on the threshold force and the dynamics of the chain are directly visible in numerical simulations. The simulations were performed using r c = 3.33a 0 and channels of length L 1000a 0 . Figure 12 shows several v-f curves for channels with 0.93 < w/b 0 < 1.1 at a disorder strength = 0.025. We first focus on the result for the commensurate case w = b 0 . The disorder leads to a significant reduction of the threshold f s with respect to the pure value µ = 0.054. The reduction is enhanced on increasing , as shown for = 0.05 in the inset. The origin of the reduction is that disorder lowers the energy barrier for the formation of vacancy/interstitial (kink/antikink) pairs in the chain. Figure 13(a) shows the time evolution of the displacements u i = x i − ia 0 upon a sudden increase of f to a value f = 0.049 > f s at t 1 . For t < t 1 , u is 'flat' and the 2D crystal formed by the chain and the CEs is topologically ordered. At t = t 1 , the motion starts at an unstable site (at x/a 0 500) by nucleation of a vacancy/interstitial pair visible as steps of ±a 0 in u. We henceforth denote the force at which this local nucleation occurs by f n . The defects are driven apart by the applied force and when their spacing becomes ∼l d , a new pair nucleates at the same site. This process occurs periodically with rate R n , leading to the formation of a domain with defect density c d = R n / v d and a net velocity v = c d v d a 0 = R n a 0 with v d the average defect velocity. In the present case of weak disorder, v d is essentially the same as for = 0. For a further increase of the force to f = 0.053, an increase of the nucleation rate is observed. In [34] we showed that in larger systems, coarsening occurs in the initial stage of depinning due to a distribution of unstable sites. However, after sufficiently long times, the stationary state consists of one domain around the site with the largest nucleation rate (smallest local threshold f min n ) with vacancies travelling to the left and interstitials to the right. It is interesting to compare this to a study of CDWs with competing disorder and commensurability pinning [53]. Using a coarse-grained version of equation (35), it was found in [53] that in the pinning dominated, lowvelocity regime, the so-called interface width W(L) = (u(x) − u ) 2 x grows linearly with the system size L. The mechanism of defect nucleation, which we observe naturally, explains this phenomenon. In addition, we found that, at depinning, the average velocity v = R max n a 0 can be described by R max n ∝ (f − f min n ) β with a depinning exponent β = 0.46 ± 0.04, similar as previously reported for 1D periodic media [56].
The defective flow profile does not persist up to an arbitrary large force. In the commensurate v-f curve in figure 12 and its inset, a small kink, is observed for f µ. Associated with this kink we find a transition to a much more ordered state. We have illustrated the temporal evolution of vortex displacements in this state in figure 13(b) for f = 0.08. The 'staircase' structure has vanished and the relative vortex displacements are greatly reduced. In fact, in the above-mentioned study [53], a very similar transition in the CDW dynamics was found, and was shown to be of first order. We leave the precise dependence of this transition on disorder and vortex interactions in our channels for future studies.
We now turn to the incommensurate case.
The v-f curves with w = b 0 in figure 12 all exhibit a finite threshold instead of the vanishing threshold for the incommensurate channels without disorder (figure 5). With disorder, the defects that are present in the channel, couple to the disorder in the CE, which causes a pinning barrier f d . This barrier has a distribution along the channel {f d } and maximum value f max d . We now focus on the curves with small defect density c d 1/ l d for which the defects are individually pinned. In this regime, the threshold force f s satisfies f s f max d . The precise threshold behaviour depends on the distribution of barriers {f d }, similar to those for LJJ and CDW systems [49]. As an illustration, we show in figure 14(a) the evolution of displacements for a channel with w/b 0 = 0.99 for a force just above threshold f(t > 0) > f s . The static configuration at t = 0 (f < f s ) shows that the disorder breaks the periodicity of the 'soliton' chain. For t > 0, depinning starts with the defect at x 270a 0 and proceeds via a 'collision-release' process between the moving defect and its pinned neighbour. Thus, for f f s strong local variations in the defect mobility exist and the overall chain velocity depends n . For f µ, the structure of both defects disappears again. The resulting dynamic state resembles that of the large velocity profile shown in figure 6, but with additional 'roughness' due to the weak CE disorder.
The v-f curves in figure 12 at w/b 0 = 0.93 and w/b 0 = 1.1, for which the defect density in the chain is larger, exhibit a smaller threshold force. In this regime the interaction between defects starts to become important and f s is determined by collective pinning of the defects. This situation was studied analytically for the case of Josephson vortices in a disordered LJJ in [50]. We will not consider this situation explicitly but we note here that, as the disorder and the typical pinning force on the defects increases, the onset of the collective pinning regime shifts to larger defect density, where defect interactions are stronger [50].
In figure 15 we show the dependence of f s on channel width, both for the weak-disorder regime treated above and for larger disorder. The data at = 0.025 exhibit a sharp peak at w = b 0 , reflecting the gap between minimum nucleation threshold and maximum defect pinning force. Larger disorder however rapidly smears the peak, being eventually completely suppressed for 0.15. The origin of this behaviour is a spontaneous nucleation of defects in the static chain at larger disorder. This is conveniently illustrated via the changes in the distribution of the individual defect pinning force {f d } and that of the nucleation threshold {f n } for increasing disorder strength, shown in figure 16. The data were obtained by simulating hundreds of short channels (L = 100a 0 ) both with one vacancy, yielding {f d }, and without 'geometric' defects, yielding {f n }. While for = 0.025, the distributions are separated (formally, in infinite systems, such separation only exist for bounded disorder, see the next section), for larger disorder they start to overlap and they become nearly identical for = 0.075. This implies that nucleated defect pairs at w = b 0 can remain pinned, while at incommensurability defects may be nucleated before the 'geometrical' defects are released. In other words, regardless of the matching condition, the static configuration contains both kinks and antikinks.
While for bounded disorder, static defects first appear at a disorder strength defined as c , the complete collapse of the peak in f s is associated with the presence of a finite density of disorderinduced defects, of the order of the inverse kinkwidth ∼l −1 d . We define the disorder strength at which this occurs as * , here * 0.15. An example of the displacement fields for this disorder strength is shown in figure 17 for w = b 0 . The lower panel shows the displacements for f f s , relative to the displacements in the CE. Clearly, the static configuration has numerous defects. In general, the approach to the critical pinned state occurs by avalanches in which local nucleation and repinning, i.e. non-persisting nucleation events, drive the rearrangements. The upper panel displays the evolution of displacements above threshold, revealing a growth of 'mountains' due to persistent nucleation, superimposed on a disordered background.
The effect of large disorder on the shape of the v-f curves is shown in the inset to figure 15. All curves now exhibit essentially a linear behaviour, 8 except in a small regime 8 A detailed study of the v-f curves at strong disorder showed that the friction force f − γv exhibited, besides the usual decrease on increasing v above threshold, a shallow minimum and then a weak increase when further increasing v. These features did not depend on system size or simulation time step and varied little with commensurability. The observed behaviour is in marked contrast with analytical and numerical studies in standard CDW models, which find a correction f − γv ∝ v −1/2 at large velocity for a 1D system, see [57]. The difference could arise from the additional u independent random force or the remaining (relatively small) commensurability pinning which are present in our system, but absent in the previous studies. This issue deserves further investigation. just above f s . We note also that, in this disordered regime, a gradual transition to a smoother displacement field occurs at larger forces, similar to the dynamic transition found for CDWs. Going back to figure 15, we should also mention the overall asymmetry of f s with respect to w/b 0 = 1 and the slight decrease of f s on increasing w/b 0 at large disorder. These effects are unrelated to the competition between commensurability pinning and disorder discussed so far, but simply reflect the overall decrease of the edge potential for larger width.

Analysis of pinning forces and crossover to strong disorder
Using the results of section 5.1, we now analyse in more detail the dependence of the pinning force on disorder and the vortex interaction range. We focus on the average pinning strength of isolated defects, which we derive here in a semi-quantitative fashion (the formal calculation is deferred to appendix D). Our analysis applies to the case of weak disorder, i.e. we assume that the defect shape is unaffected by disorder [50]. Extrapolation to larger disorder provides a useful estimate for the crossover value * at which the commensurability peak vanishes. We conclude the section with a summary of previous results [34] for the threshold forces f min n and f max d in the special case of bounded disorder.
The disorder correction equation (34) to the energy of the vortex chain consists of a term H a , due to amplitude fluctuations in the periodic potential, and a term H s , due to random coupling to the strain. We first evaluate the typical pinning energy of a defect E 2 a due to the amplitude fluctuations. The local fluctuations are assumed to be uncorrelated on a length scale a 0 , and have a variance (δµ/k 0 ) 2 . Hence, for a defect in the chain, which extends over a range l d , the resulting random potential has a variance E 2 a (δµ/k 0 ) 2 (l d /a 0 ) µ 2 2 l d a 0 /4. The typical pinning force on a defect is then given by E 2 a / l d which reduces to The typical pinning energy E 2 s of a defect due to the term H s in equation (34) is estimated in a similar way: the mean-squared energy due to coupling of a single fluctuation in V s to the strain of a defect is ∼ s (0)(r d /a 0 ) 2 (2a 0 / l d ) 2 , where r d is the range of s , given below equation (34), and a 0 / l d represents the strain. On the scale of a defect, there are l d /r d such fluctuations. Thus the associated random potential for a defect has a variance E 2 s ∼ s (0)(r d /a 0 ) 2 (a 0 / l d ) 2 (l d /r d ). Taking s (0) l (0) and using equation (32), in which case r d = λ, yields E 2 s 2C α (U 0 λ) 2 (λ/a 0 ) α /(l d a 0 ) (the factor 2 comes from the refined calculation in appendix D) The typical pinning force E 2 s / l d due to random coupling to the strain is then given by The ratio between the two characteristic defect energies is E 2 2 , which shows that, particularly for increasing λ/a 0 , the dominant pinning is due to random coupling to the strain. Henceforth we use only this contribution. We next estimate the disorder strength where the commensurability peak vanishes. As mentioned before, this collapse occurs when the density of disorder-induced defects becomes ∼ l −1 d , in other words, when the typical energy gain of a defect due to disorder becomes of the same magnitude as its bare elastic energy (dx/a 0 )(κ 0 /2)(∂ x u) 2 µa 0 √ g. This leads to * ∝ C −1/2 α g −(2α+1)/4 . For the particular case of random strains that are identical for all rows (α = 2), * is given by * This can be compared to the numerical data in figure 15. Even though those results were obtained for λ/a 0 = 1 (g = 3π), formally outside the regime of validity of our analysis, the predicted value * 0.18 is in reasonable agreement with the data. In the particular case of bounded random strains in the CEs, the distribution of the nucleation force {f n } at commensurability is bounded from below by f min n and that of the defect pinning force {f d } is bounded from above by f max d (at weak disorder). For completeness we give here the previously derived results [34] for these extremal values: both occur due to disorder fluctuations on the same length scale as that on which the displacement field u(x) varies. For a defect, this naturally corresponds to l d . The associated maximum defect pinning force is f max d /µ ∝ g 3/2 (for uniform strains [34]). For nucleation, at f µ, the appropriate length scale is l san , the extent of a so-called small amplitude nucleus [40]. Due to the nonlinearity of the pinning force, l san itself depends on the force, i.e. l san (f ) > l d and it diverges for f → µ. As shown in detail in [34], this leads to a minimum nucleation threshold given by 1 − (f min n /µ) ∝ [g 3/2 ] 4/3 . From the condition f max d = f min n , one then obtains the disorder strength c at which pinned defects can first appear spontaneously in the system: c g −3/2 < * .

Wide channels with weak disorder
We now consider how channels of larger width, in which vortices have the 2D degrees of freedom, behave in the presence of weak-edge disorder. Close to commensurability, the effects we find are similar to that for a single chain. However, around 'half filling' the importance of the transverse degrees of freedom of the channel vortices become apparent.

Behaviour near commensurability
For the commensurate case, w/b 0 = n, weak CE disorder causes a reduction of the threshold with respect to the ideal value f 0 s = µ/n, see the data in figure 18. The reduction originates from defect formation at threshold, as illustrated for w = 3b 0 in figure 19. Three snapshots of the displacements of individual rows inside the channel are displayed in figure 19(a). The first snapshot is for f < f s 0.7f 0 s and yields the 'flat' profile. The subsequent snapshots for f > f s reveal simultaneous nucleation and motion of a pair of 'oppositely charged'TSFs, each terminated at the CEs by a pair of edge dislocations (see figure 19(b)). The macroscopic, stationary motion of the array is governed by periodic repetition of this process at the least stable nucleation site. In figure 19(c), we show the vortex trajectories during nucleation of the TSF shown in figure 19(a) and (b). Very similar images were obtained in decoration experiments at the initial stage of VL depinning in NbSe 2 [11], implying that even for weakly disordered VLs, defects may nucleate at depinning (see also [12]). We also note that in simulations of a rapidly moving 2D VL [58], the same nucleation mechanism as in figure 19, but relative to the co-moving frame, was identified as source of velocity differences between chains.
For incommensurate channels close to commensurability, the TSFs which are caused by the mismatch are pinned by the disorder. For the weak disorder strengths considered here, the

Behaviour around 'half filling'
The dependence of f s on the channel width for = 0.05 is shown in figure 20 for two disorder realizations. For larger channel widths, smaller channel lengths were used with L 1000/(w/b 0 ). The data around matching (w/b 0 n) reflect the nucleation and pinning of TSFs as discussed above: f s at commensurability is reduced compared to the pure value µ/n and the commensurability peak is considerably broadened, particularly for larger n, due to the pinning of TSFs. The apparent discontinuity in the peak may be an artefact of the finite channel length (see the discussion in section 5.2 on the distribution of nucleation sites). A new, and robust feature, however, is that around 'half filling' distinct maxima in f s appear. The origin of these maxima is illustrated by the static structure for w/b 0 = 3.48 in figure 21. The triangulation shows that, in addition to aligned dislocations with b x, also misaligned dislocations appear with Burgers vector at an angle of about ±60 • with the CEs. These misaligned defects are locally stabilized by the disorder in the CEs and thus also pinned by the disorder which leads to the increase in the threshold force. In addition, the projection of the driving force along the glide direction is always smaller for misaligned dislocations than for aligned dislocations. It is seen that the misaligned dislocations separate regions with n rows from regions with n ± 1 rows. Either of these regions may thus be considered as an LSF. In section 4 we mentioned that in the absence of disorder, LSFs are metastable, the structure with a single integer number of rows and a regular distribution of TSFs is energetically somewhat more favourable. In presence of weak disorder, however, LSFs are stable and have the lowest energy.
As for the dynamics at f > f s , the dislocation structure and flow pattern are generally different from the static pattern of LSFs. We illustrate this in figure 22 where the vortex trajectories at various (increasing) driving forces are shown. For small drive ( figure 22(a), f 0.25µ/(w/b 0 )), a region of plastic motion within the channel is seen at about the location   where the static pattern shows a four-row structure. Vortex transport through these fault zones occurs by repeated nucleation and annihilation of misaligned defects. At a fixed driving force, an LSF remains at the same position, although its boundaries fluctuate over a distance of at most 2 − 3a 0 . This contrasts the situation in the absence of disorder where LSFs can move along the channel via a 'climb'-like process (see [36]).
For different driving forces, the location and amount of either n or n±1-row regions or fault zones is different, as shown in figure 22(b)-(d). In this particular segment, the n = 4 region expands on increasing the drive but at other locations the reverse can occur. Moreover, after cycling the force, a different structure can occur at the same drive, i.e. no unique structure exists at a given force. This may also lead to small hysteresis in the v-f curves. Overall, we see that the transverse degrees of freedom in the channel, in combination with disorder, give rise to an important new mechanism of yield strength enhancement. In the following section, we analyse these structures in more detail in the context of strong disorder.

Wide channels with strong disorder
The behaviour of f s versus w/b 0 in figure 20 still exhibits considerable discrepancies with the experimental data in figure 2. Clearly, the CE disorder underlying these experimental data 9 differs from the type of weak CE disorder considered so far. Motivated by recent imaging experiments [35], which showed glassy vortex configurations in the NbN edge material, we now consider the case of strong CE disorder. As will be shown, in this regime the effect of transverse degrees of freedom and the presence of misaligned defects provide the main mechanism for the critical current oscillations.
The simulations we discuss here in detail were performed at a large disorder strength = 0.2. The only remaining order in the CE arrays is their preferred orientation with the principle lattice vector along the CEs. The system sizes were typically such that wL 1500a 0 b 0 . We also allowed for quenched defects between the two CEs (uncorrelated longitudinal strains in the upper and lower CE arrays), but at these large disorder strengths this is not essential.

Static structures, yield strength and depinning
In figure 23 we show triangulations of the static vortex configuration for channels of width w/b 0 = 3, 3.5 and 4. As seen in figures 23(a) and (c), even in the matching case, dislocations are present due to the large random stresses from the CE. While some of these may originate from quenched 'phase slips' between upper and lower CEs, we checked that, also without these 'phase slips', the matching structure at this disorder strength always exhibits twisted bonds (defined by two adjacent dislocations of opposite charge) at the CEs as well as oppositely 'charged' TSFs terminated by defect pairs at the CE. While most defects are thus situated along the CE with b x, occasionally they can be located inside the channel and have misoriented Burgers vector.
Turning to the mismatch case in (b), it is seen that a region of four rows (left) coexists with a three-row region (right). In between these regions, there must be dislocations with misaligned Burgers vectors. In addition, numerous other misaligned defects are visible, rendering local destruction of the chain alignment with the CEs (although less frequent, the latter can also occur in the matching state, see figure 23(c)).
To further characterize the disorder in the f = 0 structures, we analysed, for the regime 2.5 < w/b 0 < 3.5, the average length of domains without misaligned dislocations, L , as well as the total fraction, L tot /L, of regions with n = 3 rows. The results are shown in figure 24. As observed, both quantities are maximum at w/b 0 = 3 and decay considerably away from matching: for |w/b 0 − 3| > 0.3 the average length of 'correlated' n = 3 regions in the static structure is no more than 10a 0 and they make up less than ∼50% of the channel. The remaining fraction 1 − (L tot || /L) contains misaligned dislocations and small regions with 2 (w/b 0 2.5) or 4 (for w/b 0 3.5) aligned chains.
In figure 25 we show the behaviour of the threshold force for = 0.2. The modulation of f s with channel width is still present but it has changed considerably compared to the weak disorder case. The sharp maxima at integer w/b 0 have vanished, very similar to the case of the 1D chain at strong disorder, see figure 15. Instead, we now observe smooth oscillations, with maxima in f s for w/b 0 n + 0.65 and minima for w/b 0 n + 0.15. The maxima, although occurring slightly above 'half filling', are of similar nature as the local maxima at w/b 0 n ± 1/2 for weak disorder: they are related to numerous misaligned defects present in the structure around mismatch. They enhance the flow stress compared to that of the structures around matching with predominantly aligned defects.
The differences in threshold force are also reflected by the vortex trajectories at small velocity. In figure 26 we show these trajectories for channel widths w/b 0 = 3.1, 3.6 and 4.1, close to the extrema in f s . The first thing to notice is that the trajectories at mismatch ( figure  26(b), w/b 0 = 3.6) are densely interconnecting, on a scale ∼a 0 , i.e. the motion is fully plastic and creation and annihilation of misaligned defects occurs over nearly the full channel length. For the 'near matching' cases in (a) and (c), the motion occurs mainly in the form of integer chains. However, for the small velocity considered here, the dynamics still exhibits a considerable amount of plastic motion. This partly occurs due to vortices which remain stuck at the CEs (also  visible in (b)) and partly due to narrow interconnecting regions. In fact, comparing the average length of the regions without inter-row switching in (a) with the data for the static structure in figure 24, it is seen that near matching the structure at small v is more disordered than the corresponding static structure, reflecting nucleation of misaligned defects in regions which were free of such defects at f = 0.

Analysis of dynamical properties
We now consider in more detail the properties of the moving structures at f > f s . We first show in figure 27(a) two characteristic v-f curves associated with a minimum (w/b 0 = 3.1) and a maximum (w/b 0 = 3.6) in flow stress. It is observed that the enhanced threshold in the latter case also translates in a larger dynamic friction, f − γv, of the driven structure. In addition, the latter curve exhibits a small positive curvature in the small-velocity regime v 0.01. This nonlinearity is related to the strong plastic nature of the motion in this regime. A convenient way to characterize the amount of plasticity is to calculate the mean-squared displacement of vortices from their centre-of-mass (M) positions [60]: where α = x − x M , y − y M denote longitudinal and transverse displacements, respectively. As shown by Kolton et al [60] for a 'bulk' 2D VL, W α (t) can be characterized by W α = R α t ξ α . For example, when ξ y = 1 we have normal transverse diffusion (caused by 'random' switching of vortices between chains) with R y = D y the diffusion coefficient. However, in the channels W y will become bounded at long times (large x M ) due to the finite channel width. In the inset (i) to figure 27(a), we have illustrated this behaviour for w/b 0 = 3.6 and f = 0.021. For x M 15a 0 (the point x M = 0 was chosen in the steady state after transients had disappeared), W y increases linearly as in usual diffusion, but for larger displacements (times) W y levels off and eventually saturates. In addition, even when all vortices remain in their chain (no transverse diffusion), W y initially increases to a value W y,c due to finite chain 'roughness'. Such a behaviour is observed for x M 2a 0 in the inset (ii) to figure 27(b), where W y is shown for a more coherent flow situation at w/b 0 = 3.1 and f = 0.041. In practice, we found that W y,c is always reached for x M < 3a 0 , while for the long time (large distance) behaviour significant levelling of W y occurs only when x M 10a 0 . The appropriate regime used to characterize real diffusion is therefore given by 3a 0 < x M < 10a 0 (i.e. x M = 7a 0 ). Further, D y itself does not directly reflect the density of chain-switching points along the channel. Defining the 1D density of such switching points as ρ 1D sw N sw /L, the rate of switching increases linearly both with ρ 1D sw and with the average velocity: sw v. Hence, the diffusion constant is given by D y a 2 0 /t sw = ρ 1D sw va 2 0 . Being interested in ρ 1D sw , we therefore divide out the intrinsic velocity dependence of D y and calculate sw versus f are shown in figure 27(b) for the two cases w/b 0 = 3.1 and 3.6. Around the minimum in flow stress (•), the overall value of ρ 1D sw is considerably smaller than around the maximum (•), similar to what is seen in the trajectories in figure 26. For both cases, ρ 1D sw is clearly reduced on increasing the force. This reflects both a decrease in the number of fault zones in the moving structure as well as a suppression of switch events (called 'transverse freezing' in [60]) within regions already organized in n moving rows. Around matching, ρ 1D sw smoothly vanishes at f ∼ µ ≈ 0.05, 10 indicating complete dynamic ordering into an n = 3 row structure without transverse wandering. For w/b 0 = 3.6, an ordering transition is also observed but it occurs at much larger drive (f ∼ 3µ, results not shown) and the array orders into an n = 4 row configuration, with a reduced spacing b < b 0 between the chains and average vortex spacing a > a 0 within the chains. At the end of this section, we will show how in this large-drive regime the number of rows changes with w/b 0 .
The insets to figure 27 also show the longitudinal mean-squared displacements. For strongly plastic flow at mismatch (inset (i)), W x is large and ξ x is close to 2, as expected when some vortices remain stuck at the CEs. For the more coherent situation in (ii), where transverse switching has nearly ceased, W x is smaller and ξ x 1. We always find an exponent 1 < ξ x < 2, similar to the results for 2D VLs in [60]. Interestingly, even without transverse wandering (ρ 1D sw = 0), W x increases indefinitely (with ξ 1), indicating that the moving integer chain structure still exhibits slip events and (local) velocity differences between the chains.
We illustrate some more aspects of the structures at large drive in figure 28, where the vortex trajectories and triangulations of a single snapshot are displayed for w/b 0 = 3.1 and 3.55, both at f 2µ. For w/b 0 = 3.1, we do not see any transverse wandering in figure 28(a). Figure 28(b) shows that in this case the dynamic structure exhibits only dislocations with b x, mainly located at the CE but also occasionally between chains inside the channel. The latter dislocations are possibly dynamically nucleated. They are non-stationary in the co-moving frame [58,61] and lead to slip events and the growth of W x as was discussed above. Turning to the mismatch case (figures 28(c) and (d)), it is seen that the dynamic structure consists of n = 3 and n = 4 row regions coexisting in the channel. At the driving force considered here, ρ 1D sw has a finite but small value a 0 ρ 1D sw 0.007, which is solely due to switching of vortices in the fault zones separating the 3 and 4-row regions. Within these regions transverse wandering is absent. At yet larger forces the minority 3-row regions vanish and complete ordering into four rows occurs, as for the 10 The decay of ρ 1D sw (f) is approximately exponential, regardless of (mis)match, but generally shows a small discontinuous jump to zero at large velocity where a 0 ρ 1D sw < 0.001. However, the exact force or velocity at which this jump occurs and the value of a 0 ρ 1D sw at the jump showed considerable scatter for different disorder realizations or when varying the system size. For practical purposes we therefore study the ordering transition defined through two specific criteria for a 0 ρ 1D sw > 0.001. case w/b 0 = 3.6. For a given driving force, the n-row regions again remain quasi-static during motion. The triangulation in (d) exhibits the expected misaligned dislocations at the fault zone but in general also aligned dislocations are present between the chains within an n-row region (not shown in the figure). It is also interesting to compare the velocities in the two coexisting regions with n and n ± 1 rows. Denoting the vortex velocity in an n-row region by v n and the longitudinal vortex spacing there by a n , flux conservation implies that nv n /a n = n v n /a n . We checked that, in both regions, the average flux density 1/(a n b n ), with b n = w/n the row spacing, was equal to 1/(a 0 b 0 ) within ∼4%. Therefore, a n = (n /n)a n and consequently the average vortex velocities are equal, v n = v n . However, the local washboard frequency, ν n = v n /a n , is different in both regions. Indeed, the spectrum of the velocity fluctuations in channels with dynamic coexistence of n and n ± 1 rows showed two shallow fundamental peaks at frequencies ν n /ν n = n /n. We however note that, both around matching (where a single peak occurs at ν ∼ v/a 0 ) and around mismatch, the amplitude of the washboard peak(s) decays on increasing the channel length L. In addition, for large velocities, the mixed n/n ± 1 structures ultimately anneal into a single n or n ± 1 domain, causing the collapse of one of the peaks.
The simulations also allow to explicitly show the influence of the transverse degrees of freedom on the modulations of the dynamic friction force (and, ultimately, the critical current). Generalizing the expression for the friction force equation (20) in section 3.2 including transverse fluctuations leads to where f x fric and f y fric denote the contribution to the total friction due to longitudinal and transverse fluctuations, respectively. Figure 29(a) displays the numerical results for these quantities, normalized by the total friction force f − γv as obtained from the v-f curves. The sum of the two data is 1 as it should be, confirming the correct numerical evaluation of these quantities. The total friction is also shown for clarity (figure 29(b)) and is essentially the same as the data shown in figure 25. Clearly, the relative contribution of longitudinal fluctuations to f fric decreases on approaching a mismatch situation, while f y fric /f fric increases accordingly. These qualitative features remain present also at larger velocities, where permeation modes between chains are being suppressed.
Finally, we consider the continuous modulation of structural (dis)order in the moving arrays when varying the channel width. As was shown in figure 27(b), for a fixed velocity, the density of switching points ρ 1D sw is maximum around mismatch while for fixed channel width it decreases with velocity. We can then define an ordering velocity v c operationally as the point at which ρ 1D sw is reduced below a certain threshold. In figure 29(c) the behaviour of v c for channel widths w/b 0 > 2 is shown for two criteria. The lower curve (v c,1 ), with a 0 ρ 1D sw ≈ 0.01, corresponds to a situation with ∼70% of the channel length 'transversely frozen' into integer chain regions of length 10a 0 . The upper curve (v c,2 ), with a 0 ρ 1D sw ≈ 0.002, corresponds to nearly fully annealed arrays. As observed, v c,1 increases smoothly by about an order of magnitude between a matching and mismatching situation, while v c,2 increases by a factor 5 (except for w/b 0 > 5). Regardless of the criterion, the amplitude of the modulation of v c is considerably larger than the amplitude of the f s modulation, which we will further discuss shortly.
For v v c,2 the arrays all completely anneal into a single n-row structure without permeation modes. Figure 29(d) displays the number of rows n of these structures versus w/b 0 . The switching from n to n + 1 rows is seen to occur at half integer channel widths for w/b 0 > 3 but the transitions 1 → 2 and 2 → 3 take place below these points, at w/b 0 1.35 and w/b 0 2.4, respectively. Around the transition points, there are regions in which n and n ± 1 rows coexist at smaller velocity v v c,2 . The widths of these regions were found to be w/b 0 = ±0.05.

Discussion
The results for f s versus w/b 0 at strong disorder, in the previous section, show strong resemblance to the measurements of the critical current versus field B in the experimental channel system, (see figure 2). For a more detailed comparison, we combine both data in figure 30, represented in terms of the parameter A (see section 1): for the experiment, A exp is determined from A exp = F s w/(2c 66 ) and is plotted versus B 1/2 , for the simulations, A = (f s /µ)(w/b 0 )A 0 . 11 As observed, the shapes of the oscillations are in very reasonable agreement, although the overall value of A from the numerics is larger than the experimental value. 12 The important conclusion from the simulations √ 3Bρ 0 )), where ρ f 0.2Bρ 0 /B c2 was used, with ρ 0 the normal resistivity and B c2 the upper critical field. Taking λ = 1.1 µm, B 0.012(w/b 0 ) 2 T for the channel system under investigation, B c2 = 1.55 T and ρ 0 = 2 µ m, yields (v/a 0 )(γa 2 Occasional checks using the latter criterion however yielded results for f s and A/A 0 essentially the same as in figures 25 and 30. 12 There are several possible explanations for the discrepancy between A exp and the numerical result for A. The experimental system consists of channels with walls of finite height at the CEs. Screening currents along these walls here is thus that the maxima in the measured critical current do not correspond to traditional commensurability peaks, but are caused by enhanced plastic motion and transverse deviations in mismatching channels with strong edge disorder. Previously, we indeed obtained experimental evidence for this mechanism via Introduce at first occurrence ML experiments [32] in which an rf-drive is superimposed on the dc-drive (see also [63]). Interference between the former and collective modes of deformation in the moving array leads to plateaux in the current-voltage curves. The plateau (ML) voltage then directly yields the number of moving rows n and transitions from n → n ± 1 rows were observed to coincide with maxima in critical current. 13 More detailed ML experiments provided a wealth of additional information on the dynamics of the arrays [33]. In particular, a minimum velocity was required to observe the ML phenomenon, which we identified as ordering velocity. This ordering velocity exhibited a strong upturn away from matching, similar to the behaviour of v c in the simulations. As proposed in [33], v c can be estimated from a modified version of the dynamic ordering theory of Koshelev and Vinokur (KV) [38]. Firstly, instead of the 2D random potential in [38], for the channels v c is related to the short-range-correlated random stress from the CE, with rms amplitude ∼ε ce c 66 and ε ce the random strain. Secondly, when thermal fluctuations can be ignored compared to typical defect energies (see below), v c is inversely proportional to the energy k B T p for creation of small defect pairs (see also [65] , with k BTp A p c 66 a 2 0 /2π (per unit vortex length). The typical pairs referred to are those with misaligned Burgers vector, which are responsible for breakup of the chain structure. Hence, the increase in v c away from matching implies a decrease in the pair formation energy, which is accounted for by including the parameter A p 1 in k BTp , while setting A p ≡ 1 at matching. Experimentally, A p decreased to ∼0.1 at mismatch. The behaviour of v c and A p in the simulations is analysed using the dimensionless form of the above formula for v c : The data for v c,1 around matching (where A p = 1) yield as measure for the random strain ε ce 0.19 and, near mismatch, a reduction of the defect pair-creation energy by a factor A p 0.1. The latter is in very reasonable agreement with the experiments, and the reduction of k B T p near mismatch also qualitatively agrees with the large number of disorder-induced fault zones observed in the static structures near mismatch (figure 23). However, the result for ε ce is considerably larger than the value ε ce = 0.025 found experimentally. This is also manifest in the fact that in the experiments the pinning frequency, defined as f s /(γa), always exceeds the ordering frequency v c /a, while in the simulations f s γv c . Within the modified KV theory, the pinning frequency and the ordering frequency are directly related via [33]: v c /a = τ(f s /γa) 2 may cause the average distance between the first mobile row and the first pinned rows to be larger than that in the simulation (where it is b 0 at matching), leading to an overall reduced value of the shear interaction. In addition, the precise amount of disorder in the experimental system is unknown (imaging studies only exist on related geometries and outside the relevant field regime, see [62]. Further, differences in the amount of longitudinal and transverse positional disorder may exist. Finally, compared to the experiments, the simulations were performed for relatively small ratio λ/a 0 . 13 In the experiments we found that the switching point n → n ± 1 together with the maxima and minima in J s can be shifted to different magnetic fields by changing the field history. This implies a non-trivial relation between the effective width w and the physical channel width w etched , depending on the amount of screening currents along the channel walls, see [64]. Experimentally, τ was found to be independent of the matching condition. The decrease in the defect pair creation energy (∝A p ) then relates to the increase in the pinning frequency (yield strength) away from matching (∝A) as A p ∝ 1/A 2 . The numerical results for A p versus A behave similarly and can be fitted by A p ∝ 1/A ς with ς ∼ 2-3, but data collected over the full range of w/b 0 show too much scatter to make a more detailed comparison. Nevertheless, the KV model qualitatively accounts for the enhanced amplitude of the v c modulation compared to that of the yield strength.
As for the discrepancy between the experimental and numerical values for ε ce or v c , one should keep in mind that experimentally the ordering velocity is determined from the onset of an n-row ML plateau, i.e. it corresponds to the velocity at which coherent n-row regions first appear, while incoherent regions may still exist in other parts of the channel. In the simulations, this may thus correspond to v c determined using a larger criterion for the density of switching points ρ 1D sw . In addition, the superimposed rf-drive in the experiments may assist reordering of the structure, also leading to smaller values of v c . For future studies it is interesting to directly compare numerical simulations of channels with mixed rf-dc drive with the experiments and test which criterion best represents the reordering phenomenon.
Additionally, the experimental results indicated that in the large-drive regime (where the ML amplitude saturates at a constant value), the coherently moving fraction of vortices does not exceed ∼40% (at matching), while it was reduced on approaching mismatch. This feature appears at odds with the simulations, where eventually the arrays all order into a completely transversely frozen n-row domain at large drive, regardless of the matching condition. At present we do not have a good explanation for this discrepancy.
Finally, we comment on the T = 0 approximation made throughout this study. When thermal fluctuations are important, not only does one expect the dynamic ordering behaviour to be changed (see [38,66]), also the sharp threshold behaviour in the v-f curves will be smeared and activated flow may occur. The relevant energy scale for both phenomena is again the energy for formation of small defect pairs k B T p . In the low-magnetic field experiments in [33], this energy was ∼2 orders of magnitude larger than k B T and on comparison with these results the T = 0 approximation was justified. To compare with the experiments near the melting field in [66], it would be required to include thermal fluctuations in the simulations.

Summary
We have presented a detailed study of the properties of vortices confined to narrow flow channels with pinned vortices in the CEs. In the experimental system which motivated this work [30,32,33], the threshold force (yield strength) shows pronounced commensurability oscillations when the natural vortex row spacing is varied through integer fractions of the channel width. The analysis and simulations presented in this paper show that in a mesoscopic channel system, the dependence of threshold on commensurability as well as the dynamics of vortices in the channels drastically vary with the amount of disorder in the confining arrays. At zero or weak disorder, the system behaves similar to 1D LJJ systems and defects at the CEs reduce the yield strength. At large disorder, the behaviour involves transitions from quasi-1D to 2D structures, where an increase in the amount of plastic deformations enhances the yield strength similar to the situation in the classical peak effect in superconductors.
We first presented a generalized s-G description for a 1D vortex chain in an ideally ordered channel. In this case (or for channels with multiple chains near commensurability), the threshold force has sharp peaks at commensurate widths, whereas it is essentially vanishing at incommensurability due to easy glide of 'aligned' defects, i.e. defects with Burgers vector along the CE. The model was then extended to study the effects of weak disorder in the confining arrays. Simulations and analytical results showed that this disorder causes the sharp maxima in the threshold force at matching to be lowered and broadened due to nucleation (at matching) and pinning (at mismatch) of edge defects. Apart from these defects, the arrays respond elastically in this regime, both near threshold and at large drive. We studied numerically the relevant edgedefect dynamics and, using the s-G model, we analysed the crossover to strong defect pinning on increasing the disorder strength.
For large disorder in the CEs, matching between the longitudinal vortex spacings inside and outside the channel becomes irrelevant and the peaks in threshold force around matching completely vanish with a 'saturated'value for f s of about 30% of the ideal lattice strength. Around mismatch however, the arrays become susceptible to the formation of defects with Burgers vector misoriented with the channel direction. Such defects either locally break up the integer chain structure or exist at the boundaries of n and n ± 1-row regions, coexisting in the channel. At large disorder, they are strongly pinned and cause the threshold force to exceed that around matching. Approaching a matching condition, the density of misaligned defects is reduced and a smooth modulation of f s results, with minima near matching. The depinning transition always involves plastic deformations inside the channel, but the amount of plasticity drastically increases away from matching. Using the density of transverse switching points (obtained from the transverse diffusion in the moving structures) as dynamic 'order parameter', we study the evolution of the moving structures by changing the channel width and the drive. The arrays reorder (partially) into transversely frozen n-row regions at a velocity v c that shows a similar modulation with commensurability as the threshold force. Finally, we compared the modulations of f s and v c at strong disorder with the available experimental results and with the dynamic ordering theory in [38] and found good qualitative agreement.
for the reduction of superfluid density with field and an additional attractive interaction due to overlapping vortex cores. The Fourier transform of this interaction reads: where b = B/B c2 , λ = λ/(1 − b) 1/2 and ξ = ξ C (2−2b) 1/2 the effective coherence length with C∼1. It is convenient to split equation (A.1) for the total potential in terms of the contribution V m of row m (see figure 3). Integrating over k y and using Poisson summation yields where k 0 = 2π/a 0 and the prefactors are First of all, we neglect in equation (A.3) the l = 0 terms which represent uniform (x-independent) interaction. Secondly, when λ a 0 , as practically encountered in films, k 0,λ ≈ k 0 . Then the |l| > 1 terms can be neglected, resulting in a sinusoidal potential |l| = 1. Further, in the summation over m only the contributions from the m = ±1 terms are significant. Next we employ the relation ξ 2 /a 2 0 = b √ 3/4π and rewrite k 0,ξ = (1/a 0 ) 4π 2 + 8π(1 − b)/( 3C 2 b). The resulting expression for the total edge potential is . For a channel with w = b 0 , the maximum µ(b) of the sinusoidal pinning force −∂ x V ce,0 at y = 0 is then given by where we used e −π √ 3 1/24π 2 and It can be checked that the associated shear modulus c 66 = π √ 3µ(b)/2a 0 is very similar to the interpolation formula of Brandt equation (A.2). Additionally, the edge potential equation (A.5) is harmonic for all fields. Hence the ideal flow stress of a commensurate, ordered channel, is characterized by Frenkels value A 0 = 1/π √ 3 for all fields. 14 In the low-field limit b 0.2 (and λ/a 0 1/π), the above expressions reduce to equations (8) and (9) in section 3.1.
Using the field-dependent vortex interaction equation (A.2), one can derive the parameters in the s-G description of section 3.1, generalized for higher fields. The equation for the chain stiffness becomes The reduced stiffness is obtained from Taking into account these refinements in equations (10) and (17), the defect width in the non-local regime becomes One can obtain the typical crossover field b nl (or typical λ/a 0,nl ) at which non-local behaviour sets in for a chain in an ordered channel by equating equation Hence, the non-local regime is absent for a channel in a material with λ/ξ 50 (and thickness d λ). For λ/ξ 60, non-local behaviour occurs for b > b nl with b nl 0.2. This is to be compared with the estimate a 0,nl < λ/3π resulting from a simple London interaction, see section 3.1.
The amplitudes h m are obtained by inserting equation (24) into the equation of motion (13) with the wavelength-dependent elasticity κ(q) from equation (12). Since h can become of the order of a lattice spacing a 0 , one expands the 'sin' term in (13) up to second order in h. Furthermore, we assume that h m decays rapidly upon increasing m and we keep only the three lowest order contributions of h m with m 3. Collecting terms of equal wave number, one obtains the following set of approximate equations: where K m,q = m 2 q 2 κ(mq) and h * denotes the complex conjugate of h. At small v, the real components of h vanish and the amplitudes describing the shape of the quasi-static deformations are given by However, the most important effect of the coupling is that above a characteristic velocity v * (see below), |h 2 | and |h 3 | decrease rapidly as |h 2 | ∼ |h 1 |/v and |h 3 | ∼ |h 1 | 2 /3v, so that only the first Fourier mode survives. This mode is described by (v v * ): Since equation (B.8) can be obtained by neglecting the elastic force terms κ n,q , it describes the deviation from the average velocity of a single particle in a periodic potential. This has the correct small v behaviour (where higher modes play a role) and large v behaviour (where h 2 ≈ h 3 ≈ 0). First-order perturbation using only h 1 to first order in equation (B.2) yields the same functional form but with K 2 eff replaced by K 2 1,q , supporting the analytical interpolation made in obtaining equation (B.12). Finally, using equation (B.1) one arrives at equation (25) in section 3.1. The solution (C.3) describes a periodic velocity modulation ∂ t h of each chain, with a y-dependent amplitude |h | and phase shift. The latter reflects periodic lagging and advancing of the inner rows with respect to the outer ones. Both the decay of |h | away from the CE and the phase shift strongly increase with velocity through l ⊥,v . In figure C.1 we have illustrated this behaviour for small velocity γv/µ = 0.4 (a) and larger velocity γv/µ = 2 (b). Although the result in (a) is actually out of the regime of validity of the high-velocity expansion, the 'in-phase' behaviour at small v is a feature that qualitatively agrees with the simulations. For comparison we show in figure C.2 numerical results for the oscillating velocity component for w/b 0 = 9 and small drive, f 2f s µ/5 (γv/µ ≈ 0.16). In addition to the small phase shift, one observes that the modulation is highly non-sinusoidal, which is not captured by our approximate solution. Based on the above, one can also evaluate the dynamic friction force f − γv. Using equations (20) and (23), the expression for the elastic-continuum u(y, t) would be f − γv = (γ/vw) dy dt[∂ t h(y, t)] 2 . However, in the velocity regime where our solution applies, the length l ⊥,v 2b 0 . Since we have a discrete number of vortex rows, it is the two outer rows at y = ±w/2 which give the dominant contribution. Evaluating f − γv using equation (C.4) then leads to equation (28) in section 4.

Appendix D. Disordered channel potential and pinning of defects
In this appendix we derive the disorder corrections to the channel potential and the effect on pinning of the chain. For the vortex interaction, we use the London potential of equation (6). Starting from equation (30)  where k ± = k 0 ± k and the second line is valid for k 0.4k 0 . The local strains thus also provide variations in the potential height. A similar conclusion holds when transverse shifts in the CE are included. In the simulations, the mean-squared strain is (∇ · d) 2  The term A y,y (x) = (2π) −1 dk V(k, y)V(k, y ) cos(kx) can be approximated by model and the remaining terms reflect the corrections due to disorder. Finally, we denote the first correction, which is due to the amplitude fluctuations, as H a and the second, coupling to the strain, as H s . For weak disorder, we can now calculate the effect of disorder on a defect in the chain by assuming that the shape of a defect at x, u d (x − x) = 2a 0 arctan [exp(±2π(x − x)/ l d )]/π, is unchanged by disorder [50]. The pinning energy of a defect due to the term H a is: The effect of the coupling to the strain is given by E s (x) = (a 0 ) −1 dx V s (x )∂ x u d (x − x) which has the following correlator: where the term C l 2C α (λ/a 0 ) α in square brackets is due to the non-local contribution V l and the factor 4 arises from the term ∼ ∂ x d in V s . Hence, for large λ/a 0 the non-local term dominates in the coupling to the strain.