The role of rare events in the pinning problem

Type II superconductors exhibit a fascinating phenomenology that is determined by the dynamical properties of the vortex matter hosted by the material. A crucial element in this phenomenology is vortex pinning by material defects, e.g., immobilizing vortices at small drives and thereby guaranteeing dissipation-free current flow. Pinning models for vortices and other topological defects, such as domain walls in magnets or dislocations in crystals, come in two standard variants: i) weak collective pinning, where individual weak defects are unable to pin, while the random accumulation of many force centers within a collective pinning volume combines into an effective pin, and ii) strong pinning, where strong defects produce large vortex displacements and bistabilities that lead to pinning on the level of individual defects. The transition between strong and weak pinning is quantified by the Labusch criterion $\kappa \approx f_p/\bar{C}\xi = 1$, where $f_p$ and $\bar{C}$ are the force of one defect and the effective elasticity of the vortex lattice, respectively ($\xi$ is the coherence length). Here, we show that a third generic type of pinning becomes dominant when the pinning force $f_p$ enters the weak regime, the pinning by rare events. We find that within an intermediate regime $1/2<\kappa<1$, compact pairs of weak defects define strong pinning clusters that extend the mechanism of strong pinning into the weak regime. We present a detailed analysis of this cluster-pinning mechanism and show that its pinning-force density parametrically dominates over the weak pinning result. The present work is a first attempt to include correlations between defects into the discussion of strong pinning.


I. INTRODUCTION
Broken-symmetry phases, as they appear in superconducting-, magnetic-, or density wave systems, exhibit physical properties on top of those originating from the underlying material. Typically, these ordered phases develop topological excitations (or defects) that govern the material properties, e.g., vortices in superconductors [1] or domain-walls in magnets [2 and 3]. Remarkably, it is the interaction between the material's and the topological defects that determines the static and dynamical properties of the latter, with pinning immobilizing vortices in superconductors guaranteeing the material's dissipationfree current transport [4 and 5] and fixing domain-walls in the magnet determining its coercive field [6]. On the fundamental side, pinning of topological defects constitutes a rich branch of disordered statistical physics with challenging phase-space and ergodicity properties, including the phenomenon of glassiness [7 and 8].
Traditionally, pinning in such systems was thought of as due to large ensembles of weak defects; the ensuing collective pinning theory [5,[9][10][11] has become a common framework for the description of pinning of superconducting vortices [7,[12][13][14], magnetic domain-walls [15][16][17], charge density waves (CDWs) [18 and 19] and other types of elastic media [20]. At the same time, an alternative viewpoint describing pinning due to a low density of strong centres was proposed early on, see Refs. [4] and [5]; recently, this strong pinning scenario has attracted increasing attention, particularly in studies of charge density waves [19, 21, and 22] and of magnetic flux-line lattices [23][24][25][26][27]. Although some effort has been made to qualitatively understand the crossover between the two regimes [28 and 29], a quantitative model describing this transition has not been developed so far. In this paper, we describe a new regime at the crossover between the two theories. We show that in a considerable part of the weak region, pinning is dominated by defect clusters co- Pinning mechanisms for flux lattices in type-II superconductors in the regime of low defect density npa0ξ 2 1. Shown is the critical force density Fc as a function of pinning strength κ ∝ fp/Cξ. For small κ → 0, pinning arises due to the collective action of a large number of defects within the Larkin volume Vc ∼ R 2 c Lc, resulting in the collective pinningforce density F coll ∼ (ξ/λ) 2 κ 3 fpnp(npa0ξ 2 ). For intermediate values 1/2 < κ < 1, pairs of defects in close proximity form strong pinning clusters that produce the cluster pinning-force density F clust ∼ (ξ/a0) 2 (κ − 1 2 ) 4 fpnp(npa0ξ 2 ); the latter dominates over the collective pinning result for a Labusch parameter κ > 1/2 + O[(a0/λ) 1/2 ]. For κ 1, pinning is strong, with the pinning-force density due to individual defects rising as (ξ/a0) 2 (κ − 1) 2 fpnp; the latter dominates over the cluster pinning when increasing κ beyond unity by the amount (npa0ξ 2 ) 1/2 . operating on short distances and forming strong pinning centers that are described with the tools of strong pinning theory. The dominance of these strongly-pinning small pairs over the weak collective ensembles can be traced back to the dispersive nature of the vortex elasticity. Pinning by rare events then interpolates between the strong pinning of individual defects and the random sum of weak pinning forces due to the many defects within the Larkin domains of collective pinning theory, as illustrated in Fig. 1.
The central problem arising in studies of pinned systems is the determination of the maximal driving (or critical) force density F c below which the system remains immobilized. This critical force is determined by the competition between the pinning centers characterized by their density n p and individual forces f p and the elastic properties of the manifold. In the present study, we focus on the vortex lattice [1] formed by flux lines or vortices, each carrying a superconducting flux quantum Φ 0 and characterized by a line energy ε 0 = (Φ 0 /4πλ) 2 (λ denotes the London penetration depth). The effective elasticityC ∼ ε 0 /a 0 , with a 0 the distance between vortices, captures the full elastic properties of the vortex lattice that combines the line tension and the interaction between vortices. The competition between pinning and elastic forces then can be quantified by the dimensionless Labusch parameter [4] κ ∼ f p /ξC, where ξ denotes the coherence length (or vortex diameter) in the superconductor. When κ increases beyond unity, individual pins change from weak to strong. The three scenarios, strong-, weak-collective-, and cluster-pinning as illustrated in Fig.  1 then provide different mechanisms and scaling laws for the critical force density F c .
The strong pinning paradigm rests on the assumption of a low defect density n p , such that κn p a 0 ξ 2 < 1 [28], and strong defects, i.e., κ > 1. In this setting, material defects act independently, resulting in a critical force density F c ∝ n p that is linear in the density n p of pinning centers. The task of calculating F c then simplifies considerably and even allows for a quantitative treatment: as defects act independently, the calculation of their contribution to the critical force density F c boils down to an effective single-particle problem where a strong defect interacts with an elastic manifold. The competition between potential and elastic forces does, however, add quite some complexity to the problem, with strong pinning inducing plastic deformations and bi-stable (pinned and free) states of the elastic manifold [5, 19, 21, and 30]. The non-symmetric occupation of these bi-stable solutions then generates a finite pinning force, with the critical force density derived from the maximally asymmetric occupation of metastable states given by F c ∼ (S trap /a 2 0 ) n p f p ∼ (κξ 2 /a 2 0 ) n p f p ; here S trap /a 2 0 defines the fraction of vortices falling into the defect trapping area S trap ∼ κξ 2 with longitudinal and transverse dimensions ∼ κξ and ∼ ξ [30 and 31].
Weak collective pinning instead, relies on the joint action of many defects, as individual weak pins with κ < 1 cannot hold the manifold. In the weak-collective pinning scenario, distant defects act with random forces on the manifold and their (random) addition within the Larkin volume V c ∼ λ 3 (λ/a 0 )/(κ 2 n p a 0 ξ 2 ) 3 (that contains a large number of pins) produces a critical force density F c ∼ [(ξ 2 /a 2 0 )n p f 2 p V c ] 1/2 /V c ∼ (ξ 2 /λ 2 )κ 3 (n p a 0 ξ 2 ) n p f p , where the factor ξ 2 /a 2 0 accounts for the fraction of defects within V c that overlap with the vortex cores. In fact, the collective force randomly accumulated in the Larkin volume V c produces an effectively strong pin [28] that satisfies the Labusch criterion κ(V c ) = 1.
In the present paper, we study the crossover between the strong-and weak-collective-pinning mechanisms near κ ∼ 1; this study leads us to the mechanism of pinning by rare events. Pairs of defects that reinforce one another appear with relative probability n 2 p and thus potentially compete with the force generated in the weakpinning scenario. In identifying suitable pairs, we find that closeby defects within the action volume ξ 2 a 0 of one defect define the relevant clusters; the density of such clusters then is given by (n p a 0 ξ 2 ) n p . Defects in one cluster act cooperatively rather then competitively. For defects with a pinning strength 1/2 < κ < 1, such neighboring pairs jointly produce a strong defect with 2κ > 1. Applying the strong-pinning formalism to these strong cluster-defects then produces a critical force density F c ∼ (ξ 2 /a 2 0 )(n p a 0 ξ 2 ) n p f p that is larger than the weak-collective force density by a factor (λ/a 0 ) 2 . This factor is a consequence of the dispersive nature of the tilt elasticity c 44 (k): while (non-dispersive) collective pinning involves the large Larkin scale R c > λ, clusterpinning appears on short distances below a 0 and hence involves the line rather than the bulk elasticity. Hence, we find a new transition region in the pinning strength κ where rare events, neighboring defects forming a strongpinning cluster, determine the critical force density F c .
The relevance of rare events has been pointed out before in the context of charge density wave pinning [32], where an analysis in D > 4 dimensions demonstrated the irrelevance of weak collective pinning. Instead, a finite but exponentially small (in the disorder strength) pinning-force density was found that originates from rare regions with anomalously coherent pinning. In our case, we deal with a D = 3 dimensional vortex lattice, where both types of pinning, weak collective and rare events contribute simultaneously, with the rare events identified as small defect pairs. The paper is organized as follows: in Sect. II, we discuss the formalism used in the description of vortex pinning for the generic case of an isotropic material and briefly present the main steps in the derivation of the pinning-force density F c in the strong-and weak-pinning scenarios and for the newly-introduced framework of pinning by close pairs of defects. In Sect. III, we first introduce the general two-defect problem for pairs of any size. In the overview section III C, we identify the strongly-pinning pairs and discuss their contribution to the pinning-force density as a function of the spatial sep-aration between the defects constituting the pair. We show that pairs of distant defects provide a smaller contribution, justifying the assumption of dominant pinning by rare clusters of close defect pairs. We proceed with a detailed analytical derivation of our results, involving an in-depth discussion of the effective anisotropic pinning potential of defect pairs (Sect. III D), of the effective Labusch parameter of defect pairs in Sect. III E, and the average pinning force of defect pairs in Sect. III F including a comparison to numerical results. Finally, in Sect. IV, we summarize our results and place them into context, including also some further directions of research.

II. VORTEX LATTICE PINNING
The pinning of a vortex lattice is an example of the (D+n)-dimensional random manifold problem; the latter describes a D-dimensional elastic manifold parametrized by ρ ∈ R D that is distorted with an n-dimensional displacement field u(ρ) ∈ R n due to a pinning potential ε pin (ρ, u). Assuming small distortions, the generic Hamiltonian describes this type of systems. Minimizing Eq. (1) yields the equation for the displacement field in the form, with the Green's function G(ρ); in reciprocal space, G(k) = 1/ck 2 . In the following, we first discuss the relevant properties of the real-space Green's function G(ρ) for our vortex problem and then turn to the peculiarities of the disorder potential ε pin (ρ, u) for the weak-and strong pinning situations.

A. Green's function
The vortex pinning problem considered here belongs to the class D = 3, n = 2 and the complex structure of the vortex lattice brings a number of modifications to the simple pinning model in Eq. (1). The Green's function for the vortex lattice (aligned along the z-axis) is in fact non-diagonal and features anisotropic and dispersive elastic moduli; focusing the discussion to isotropic superconductors and writing k = (K, k z ) with the transverse (K) and longitudinal (k z ) components of the reciprocal vector, it assumes the form [7] G αβ (k) = P αβ (K) with indices α, β ∈ 1, 2 and the projection operators P αβ (K) = K α K β /K 2 and P ⊥ αβ (K) = δ αβ − K α K β /K 2 .
The compression and tilt moduli c 11 (k) ≈ c 44 (k) ≈ (B 2 /4π)(1 + λ 2 k 2 ) −1 exhibit strong dispersion due to the long-range interaction between vortices; c 66 = BΦ 0 /(8πλ) 2 is the non-dispersive shear modulus (B ẑ is the magnetic field induced in the bulk of the superconductor). The corresponding real-space Green's function is obtained via standard Fourier transformation, with the integration over K restricted to the Brillouin zone of the vortex lattice, Of key importance will be the on-site Green's function G αβ (ρ = 0) = G(0) δ αβ . The integration in Eq. (4) then is dominated by transverse momenta near the Brillouin zone boundary K ∼ K BZ , and estimating the relevant longitudinal momentum by comparing the shear and tilt elastic energies c 66 K 2 ∼ c 44 (K BZ )k 2 z , we obtain the scaling result ]. The precise integration in Eq. (4) gives the result [23, 25, and 33] with a numerical factor ζ that depends on the chosen approximation for the elastic moduli.
To evaluate the spatial variations of the Green's function, we consider a simplified model of the vortex lattice elasticity: we drop the first term in Eq. (3) involving the large compression modulus c 11 (k) > c 66 and replace the projection operator in the remaining term by δ αβ , such that G αβ (ρ) = G(ρ) δ αβ . Our diagonal reponse function G[ρ = (R, z)] is characterized by a sharp and structured peak around the origin and a smooth decay ∝ 1/ρ at large distancesρ > λ, whereρ = (R, c 66 /c 44 (0) z) is the properly scaled distance due to the anisotropic elasticity of the vortex lattice. Going beyond the diagonal approximation does not change our strong pair-pinning results obtained below. Note that the function G(ρ) provides us with the displacement field u(ρ) = G(ρ)F due to a δ-force F δ(ρ) at the origin.
We first evaluate the Green's function in the nondispersive regime (large distances ρ), with the dominant contributions to the integration in Eq. (4) originating from small momenta λ 2 k 2 1, such that c 44 (k) ≈ c 44 (0). The anisotropy of the Green's function in Eq. (3) generates different decays along the directions longitudinal and transverse to the induced magnetic field, that is for ρ = (0, z) and ρ = (R, 0). To simplify the calculation, we remove this anisotropy by introducing the rescaled momentum vector q = (K, c 44 (0)/c 66 k z ) with c 44 (0)/c 66 = 16πλ 2 /a 2 0 , which leads to with c 44 (0) c 66 = (B 2 /16π √ π)(a 0 /λ) and the rescaled distanceρ = (R, c 66 /c 44 (0) z). Integrating over the Different domains analyzed in the evaluation of the scaled Green's function g(R, z) = G(R, z)/G(0). In the non-dispersive region (yellow) outside the ellipseρ 2 = R 2 + (a 2 0 /16πλ 2 ) z 2 ≈ λ 2 the Green's function decays ∝ 1/ρ, see Eq. (7). The value of the Green's function on the ellipse boundary is g(ρ ∼ λ) ∼ a 2 0 /λ 2 . Inside the ellipse (green), we find several regions characterized by different scaling results. For a0 z λ 2 /a0, the Green's function reads g(R, z) ∼ (a0/z)e − √ πR 2 /a 0 z + a 2 0 /λ 2 ; starting out at g(0, z) ∼ a0/z, it decays exponentially fast along R on the scale R ∼ √ a0z (region I, green) before saturating (ignoring slow logarithmic variations) at g ∼ a 2 0 /λ 2 [region II, light green]. For small longitudinal coordinates z a0, the Green's function evaluated on the z-axis is g(0, z) ≈ 1 − z/a h , with a h ∼ a0[ln(a0/ξ)] 1/2 the healing length, and its decay along the transverse coordinate R is governed by the same scale ∼ a h (region III, dark green). For z a0, R a0, the Green's function again saturates at the value g(R, z) ∼ a 2 0 /λ 2 . Increasing z at fixed R within the interval a0 < R < λ, the ratio g first increases and goes over a maximum when z reaches the value R 2 /a0 (red dashes); this feature produces a distinct ridge in the peak region of g. momenta q yields G(ρ) = 1/[4π c 44 (0)c 66ρ ] and the reverse transformationρ → ρ = (R, z) provides us with the result Eq. (7) describes the situation where the dispersion in the tilt modulus can be neglected, which is the case at large distances R 2 + (a 2 0 /16πλ 2 ) z 2 λ 2 , see the yellow region in Fig. 2; on the inner boundary (an ellipsoid with extensions R ∼ λ and z ∼ λ 2 /a 0 ), the Green's function assumes a constant value G ∼ (a 2 0 /λ 2 )G(0) and decays ∝ 1/ρ further out, see Eq. (7). Indeed, in order to drop the dispersion in c 44 , we require the q-integral in (6) to be cut by a large distanceρ (rather than the Brillouin zone), q 1/ρ, at values where qλ < 1 (rendering dispersion irrelevant), implying thatρ > λ. are not to scale. The dark-green peak in the center saturates to unity over a region ∼ a 2 0 ; at large distancesρ = R 2 + (a 2 0 /16πλ 2 )z 2 > λ, a smooth decay ∝ 1/ρ is observed (yellow). The peak at small distances (green) exhibits a dumbbell shape and gives way to a smoothly decaying background of elliptical shape at large distances (yellow); the two contours with g = 0.003 and g = 0.001 illustrate this change of shape from a dumbbell-to an elliptical form. Fig. 11 in Appendix A shows the detailed contour plot near the center of the structured peak, including the position of the ridge.
The evaluation of the Green's function at locations inside the ellipsoid requires proper integration both over small (k λ −1 ) and large (k λ −1 ) momenta; the full calculation is presented in Appendix A. For longitudinal distances z a 0 and arbitrary R, we find the interpolation formula with γ ≈ 0.577 the Euler-Mascheroni constant. This result provides us with various scaling regimes for the Green's function, as illustrated in Fig. 2. First, fixing R = 0 and going away from the origin along the longitudinal direction, the rescaled Green's function decays as G(0, z) ∼ (a 0 /z)G(0); the result in Eq. (8) matches the non-dispersive expression Eq. (7) at the crossover z ∼ λ 2 /a 0 .
Increasing R for z < λ 2 /a 0 , the Green's function is dominated by the first term in Eq. (8) that describes a Gaussian with height G(R = 0, z) ∼ (a 0 /z) G(0) and of width R ∼ √ a 0 z; with decreasing z this Gaussian peak becomes higher and narrows down to produce a dumbbell shape peak, see region I in the schematic Fig.  2 and the neck in the contour g = 0.003 in Fig. 3. In-creasing R beyond ∼ √ a 0 z, we enter region II in Fig. 2 where the second term in Eq. (8) dominates, interpolating smoothly between the peak and the non-dispersive result Eq. (7). This smooth interpolation through region II is of order G(R, z) ∼ (a 0 /λ) 2 G(0), with logarithmic corrections that become large at small values of z where the narrow dumbbell peak at the origin decays more rapidly. Note that beyond the point z ∼ λ 2 /a 0 where the decay length R ∼ √ a 0 z meets the ellipsoidal shell, the Green's function for R = 0 already assumes a value G(0, z ∼ λ 2 /a 0 ) ∼ (a 2 0 /λ 2 ) G(0) and no substantial variations with R are seen within the region I.
Increasing instead the longitudinal distance z at fixed R < λ, the Green's function first remains flat (region II), then steeply increases ∝ e − √ πR 2 /a0z upon entering the peak region I at z ∼ ±R 2 /a 0 , and then decreases smoothly ∝ 1/z, thus defining a maximum at z ∼ ±R 2 /a 0 . The resulting ridges located at the edges of the Gaussian peak are another manifestation of the dumbbell structure of the peak in G(R, z), see Figs. 3 and 11. The discussion has to be further refined in the regime of small z a 0 . As z → 0, the first term in Eq. (8) diverges for R = 0 and vanishes for R > 0, formally approaching the 2D delta-function ∝ δ 2 (R/a 0 ). In reality, accounting for the q-cutoff at the Brillouin zone boundary in Eq. (4) provides us with the finite result for the on-site Green's function G(R = 0, z = 0). An expansion in the longitudinal direction for small z a 0 then gives [34] with the healing length a h ∼ a 0 [ln(a 0 /ξ)] 1/2 . The decay length in the transverse direction at small z is affected by the single-vortex elasticity that becomes relevant near the Brillouin zone boundary [7 and 35]. Replacing the tilt modulus by c 44 (k) → c 44 (k) + (ε 0 /a 2 0 ) ln(a 0 /ξ) then entails a saturation of the decay scale R ∼ √ a 0 z in Eq. (8) at R ∼ a h ∼ a 0 (we ignore a factor ln(a 0 /ξ) in the scaling estimates) for z a 0 (region III). For R a h , we again cross over to the region II where the Green's function assumes the constant value ∼ a 2 0 /λ 2 , up to slow logarithmic corrections. The above analysis has been carried out for a simplified diagonal expression G αβ = G δ αβ . In a further step, one may replace the identity matrix δ αβ by the full transverse projector P ⊥ αβ (K), see Eq.
(3). Focusing on the non-dispersive regime, the q-integral in Eq. (6) picks up an additional angular dependence that depends on the geometry of the problem. For the component G xx evaluated in the xz-plane, we obtain the asymptotic dependence wherez = (a 0 /4 √ πλ)z is the scaled longitudinal length.
The result (10) then exhibits a modified anisotropy at large distances: the simple scaling G ∝ 1/ √ R 2 +z 2 in the expression (7) is replaced with G xx ∝ 1/2z wheñ z R and G xx ∝ 1/R at large R z. Finally, while G xy = 0, we find that G yy (R, 0, z) = (z/ √ R 2 +z 2 ) G xx . Note that G xx + G yy = G, as expected.
Having analyzed the elastic component in the pinning problem, we now turn to the discussion of the pinning potential ε pin (ρ, u) in Eq. (1) for the cases of strong pinning, weak collective pinning, and the pinning by rare clusters. Note that the smallest transverse scale R in the context of elasticity is the separation a 0 between vortices, while separations between defects as discussed below are considered small when R reaches the effective size ξ of defects. Hence, small lengths R take a different meaning when talking about the vortex lattice (elasticity) or the pinning landscape.

B. Strong pinning
We consider a lattice of flux lines or vortices aligned with the z-axis and described by the unperturbed vortex core positions R µ ∈ R 2 . The pinning force acts on the vortex cores and the pinning energy can be expressed in the form [with ρ = (R, z)] with ε µ pin [z, u µ (z)] the random pinning potential acting on the µ-th vortex line, Here, U pin (R, z) denotes the disorder potential generated by the material defects; assuming pinning due to pointlike defects located at r i = (R i , z i ), each with identical pinning energy e p , the disorder potential takes the form The factor p(R) in Eq. (12) describes the vortex form factor, e.g., for a δT c -type pinning mechanism [7], it reads p(R) = 1 − |ψ(R)| 2 , with ψ(R) the superconducting order parameter of the single-vortex solution to the Ginzburg-Landau equations. The simple Ansatz [36 and 37] |ψ(R)| = R/(R 2 +2ξ 2 ) 1/2 provides us with Lorentzian shape for the form factor, p(R) = 1/(1 + R 2 /2ξ 2 ). Combining Eqs. (12) and (13), we express the random pinning potential as with e p (R) = −e p p(R) the pinning potential due to a single defect; note that e p (R) is maximally negative for R = 0, i.e., pinning is maximal when the defect position R i coincides with the perturbed vortex position R µ + u µ (z). Substituting this result to Eqs. (2) and (11), we arrive at the equation for the displacement of the ν-th vortex in the form with the pinning force acting in the direction transverse to the field. The last relation above applies for the Lorentzian-shaped potential.
In Eq. (15), we sum over all interactions between defects and vortices. In practice, we assume that no more than a single vortex can be pinned by an impurity and neglect interactions of vortices with defects far away from the vortex core, The sum over the vortex index µ is then restricted to a single index µ(i) denoting the vortex closest to the impurity i. The relation (15) then allows to evaluate the displacement u µ (i) of the vortex µ pinned to the defect i at the position z i ; this is nothing but the vortex tip displacement of the µ(i)-th vortex, The set of equations (17) represents a system of N coupled non-linear equations for the displacements u i , with N the total number of defects.
Within the strong pinning paradigm, we assume that defects act independently, allowing for a further simplification of Eq. (17) where the displacement u i is ascribed exclusively to the action of the defect i; the summation over j in Eq. (17) then reduces to the term j = i, i.e., we neglect the force exerted by distant defects j = i on vortices µ(j) that contributes to the displacement u i via the non-local Green's function G(R µ(i) − R µ(j) , z i − z j ). It is exactly this simplification that will be dropped later on when considering strong pinning by pairs. The system of equations (17) then reduces to N independent equations with the effective vortex-lattice elasticity defined byC = 1/G(0), see Eq. (5).
The resulting pinning-force density is obtained by summing the forces from all pinning sites. Note that the solution u i in Eq. (18) depends only on the distance of the vortex from the pinning defect x i = R µ(i) −R i . The average pinning force density is thus where n p denotes the density of impurities and the average is taken with respect to the possible position vectors x; assuming a uniform distribution of relative distances x, the average then corresponds to a simple integration over x.
It turns out that the pinning force can be expressed as the gradient of the total pinning energy, , with e pin (x) involving pinning and elastic terms, If the solution u(x) to the on-site equation (18) is unique, implying a continuous evolution with x, the average pinning force density vanishes, as follows from a simple integration of f p [x + u(x)] over x, The single-defect Ansatz is thus meaningful only in the strong pinning regime where the solution for the on-site displacement is non-unique. In this case, different values (branches) for the total pinning energy e pin (x) describe pinned and unpinned vortex states, see Fig. 4. Proper averaging accounts for the occupation of these branches which is unsymmetric, resulting in a non-vanishing average pinning-force density. We perform this analysis for a radially symmetric potential with a force f p (r) =rf p (r).
The condition for the appearance of multiple solutions is provided by the Labusch criterion [4, 19, and 28] (we note that max f p (r) = f p (r m ) with the inflection point r m obtained from f p (r m ) = 0). Furthermore, we assume that the vortices are driven in the positive xdirection and we parametrize their trajectories x = (x, b) through the longitudinal vortex position x and an impact parameter b in the transverse direction (the distinction between the 'longitudinal' field direction along z and the 'longitudinal' direction of motion along x should be clear from the context). The resulting pinning force averaged over positions x then points in the negative x-direction and is evaluated in two steps: first, we perform an averaging over x at vanishing impact parameter b = 0 and then we average over contributions from vortex trajectories with finite impact parameters b = 0. For the case of a vanishing impact parameter b = 0, Eq. (18) can be reduced to one dimension,Cu = f p (x+u) (we have dropped the index i). This is equivalent to minimizing the total pinning energy e pin (x, u) = e p (x + u) + 1 2C u 2 with respect to u. Provided x falls into the bistability region, |x| ∈ [x − , x + ], there exist multiple solutions u f (x), u p (x) for the vortex tip displacement (denoting free and pinned vortex states) [23, 30, and 33]. Substituting these solutions to the total pinning energy provides multiple branches e f,p with the jumps in energy ∆e 1 pin = [e f pin − e p pin ] x=−x− and ∆e 2 pin = [e p pin − e f pin ] x=x+ occurring at the pinning (−x − ) and depinning (x + ) points, respectively.
The result (22) remains unchanged even for a nonvanishing impact parameter b = 0 [23 and 27], provided the vortex passes the defect within the pinning distance y p along the y-direction; for the radially symmetric case, it turns out that y p = x − and hence pinning occurs for impacts with |b| < x − . For |b| > x − , the pinning forces are small and multiple branches no longer exist, implying a vanishing average over x [23 and 27]. Finally, averaging the result (22) over y contributes a factor 2x − /a 0 and thus with −e x denoting the unit vector pointing in the negative x-direction.
Eq. (23) assumes different scaling forms for the regime of very strong pinning κ 1 and for moderately strong pinning κ − 1 1. In the first case, the jump sizes are related to the pinning potential depth via ∆e 1 pin ∼ e p and ∆e 2 pin ∼ κe p [23 and 27], together providing the estimate for the magnitude of the position-averaged pinning force f p x ∼ (κξ 2 /a 2 0 ) f p and a pinning force density This result is interpreted as a pinning force f p ∼ e p /ξ [see Eq. (16)] due to a single defect exerted within the trapping area [28 and 31] S trap = 2y p (x + + x − ) ∼ κξ 2 ; S trap /a 2 0 denotes the fraction of area occupied by trapped vortices.
For moderately strong pinning with κ close to unity (that is particularly relevant for the pinning by rare events), the energy jumps are evaluated by expanding the pinning force around the inflection point at r m , f p (r m ) = 0, where f p (r m ) = κC is maximally positive [27-29, and 33], In this situation, both jumps are identical and given by the expression [27 and 28] (note that f p (r m ) < 0) Using the scaling The pinning force density follows trivially, and vanishes at the Labusch point κ = 1, in accordance with the strong pinning criterion (21).

C. Weak collective pinning
When pinning is weak, κ < 1, individual defects fail to produce multi-valued solutions for the vortex displacement and the mean pinning force in Eq. (20) vanishes. Pinning then arises through the random action of defects within the collective pinning volume V c defined as the region where the spatial fluctuations of the vortex displacement u 2 (ρ) = [u(ρ) − u(0)] 2 remains bounded by the pinning scale, u 2 (ρ) ≤ ξ 2 . The displacement correlation function can be systematically evaluated from Eq. (2) using the disorder-averaged correlator of the pinning energy density Eq. (11) [7 and 11], with the correlation function k(u − u ) = d 2 R p(R − u)p(R − u ) related to the vortex form factor p(R).
A qualitative estimate for the displacement correlator is provided by summing up distortions originating from all defects within a finite volume. In the vicinity of a reference defect characterized by the pinning force f p ∼ e p /ξ, the distortion scale u 0 is given by the onsite Green's function, u 0 ∼ G(0) −1 f p . Expressing the on-site displacement through the effective vortex lattice stiffnessC = G(0) −1 and estimating the Labusch parameter as κ ∼ f p /Cξ provides us with u 0 ∼ κξ. Assuming small defect densities and hence large typical interdefect separations, the extension of the collective pinning volume falls into the non-dispersive regime of the Green's function, see Eq. (7). Defects located a distancẽ c , these displacements add up with a random sign, as the forces from different defects are randomly directed; furthermore, only the fraction ξ 2 /a 2 0 of defects that reside inside the vortex cores are directly attacking the vortices, resulting in a total squared dis- For small defect densities, as specified by the condition κ 2 n p ξ 2 a 0 1, the pinning radius R c λ indeed falls into the non-dispersive regime (note that κ 1). Finally, summing up the random force-contributions due to the active defects within the bundle volume

D. Pinning by rare events
The collective pinning scenario described above sums up small competing contributions to the vortex lattice distortions arising from typical fluctuations in the defect distribution, involving defects that lie far away from each other within the collective pinning volume. However, it does not account for the presence of rare clusters, where two (or more) weak defects act cooperatively, giving rise to an effectively strong pinning center; the latter then is supposed to produce a distortion exceeding the scale ξ of the pinning potential. In looking for promising candidate pairs, we consider Eqs. (7) and (8) that describe the decay of the Green's function; these imply that the vortex displacement is substantially suppressed beyond a distance ∼ a 0 away from the defect. Hence, two weak defects with 1 2 < κ < 1 can be combined into a stronglypinning object characterized by κ > 1 and producing a displacement u of order ξ only if they are at most a longitudinal distance z ∼ a 0 apart and pinning the same vortex core, i.e., they are separated by at most R ∼ ξ in the transverse dimension. This consideration then provides us with the density (n p a 0 ξ 2 ) n p of strongly-pinning pairs. With only the fraction ξ 2 /a 2 0 of those clusters being located within the vortex core area and each cluster exerting a pinning force ∼ f p , we arrive at the following estimate for the pinning force density due to defect pairs (with 1/2 < κ < 1 still close to unity), Assuming a magnetic field sufficiently above H c1 such that a 0 < λ, the pinning force due to such clusters dominates over the collective pinning contribution in Eq. (31) by a factor of (λ/a 0 ) 2 . This factor in fact arises due to the dispersion of the tilt elastic modulus; in order to trace its origin, we need to understand the explicit dependence of the quantities contributing to both pinning mechanisms on the elastic properties of the vortex lattice, a dispersive tilt modulus c 44 (k) and non-dispersive shear modulus c 66 . The collective pinning radius in Eq. (30) can be obtained by comparing the elastic energy and the pinning force density is estimated as We replace a factor f 3 p in (34) with κ 3 (Cξ) 3 , and On the other hand, our strongly-pinning pairs are small and the associated elastic scales R p and L p in the longitudinal and transverse directions are related by c 66 (u/R p ) 2 ∼ c 44 (K BZ )(u/L p ) 2 with the short scale elasticity c 44 (K BZ ); hence L p ∼ a 0 c 44 (K BZ )/c 66 ∼ a 0 , where we have chosen the smallest transverse scale R p ∼ a 0 of the lattice. The density of pairs then is given by n p (n p ξ 2 a 0 ) and the resulting pair pinning-force density is (assuming again κ ∼ O(1), cf. Eq. (32)) Finally, comparing the weak-collective-and clusterpinning force densities in Eqs. (36) and (37) provides us with that demonstrates that the cluster pinning dominates over the weak-collective pinning contributions due to the dispersion in the tilt modulus with its reduction ∝ (a 0 /λ 2 ) at the Brillouin zone boundary. The concept of pair-pinning described above can be extended to larger clusters, pushing the domain of pinning by rare events further down to smaller values of κ. Within the interval κ ∈ [1/n, 1/(n − 1)], n ≥ 2 and integer, n neighboring defects are required to form a strongpinning cluster with nκ > 1; the density of such clusters is given by (n p a 0 ξ 2 ) n−1 n p and the resulting pinning force density becomes F clust ∼ (ξ/a 0 ) 2 (n p a 0 ξ 2 ) n−1 n p f p . However, for pinning strengths κ ≤ 1/n with n ≈ 2 + 2[ln(λ/a 0 )]/[ln(1/n p a 0 ξ 2 )], the collective pinning dominates; given a low density of defects such that n p a 0 ξ 2 (a 0 /λ) 2 , this crossover lies close to n = 2, κ = 1 2 . The idea of pinning due to rare events has been previously touched upon in the context of charge density wave pinning in high dimensions, see Ref. [32]. In this case, the disorder-induced distortions accumulated over a finite-sized domain are not sufficient to induce pinning. This can be easily seen by considering the elastic Green's function G(ρ) ∝ ρ 2−D in D dimensions, yielding a total displacement accumulated within a pinning domain of size R that scales as u 2 (R) ∼ R 4−D , see Eq. (2). While for D < 4, the accumulated displacement will eventually exceed the threshold required for the existence of bistabilities at large domain sizes R, this is not the case for dimensions D ≥ 4. As noted by Fisher [32], this does not render the weak disorder irrelevant, since, although with exponentially small probability, one will always find rare domains with anomalously coherent pinning. The manifold is then pinned by such rare fluctuations rather than by the collective action of the disorder landscape. In our D = 3 vortex lattice, the situation is somewhat different: for D = 3, weak pinning is still active and competes with the pinning by rare events, which take the specific form of close-by defect pairs making up for a stronglypinning object. The latter mechanism dominates for defect strengths 1 2 < κ < 1 and a small density of defects. The dominance of pinning by such rare events is, however, not an inherent property of the pinning mechanism, but rather appears as a result of the specific, i.e., dispersive, elastic response of the vortex lattice.

III. TWO-DEFECT PROBLEM
We have seen in Sect. II B that the strong pinning paradigm assuming independent action of defects is meaningful only provided the Labusch parameter (21) satisfies κ > 1; in this case, the single-defect Ansatz gives rise to multi-valued solutions for the vortex displacement, what results in a finite averaged pinning force. Here, we consider the range of pinning strength 1 2 < κ < 1 and go a step beyond the single-defect ansatz by considering pairs of defects. Given N defects, of all possible N (N − 1)/2 pairings there will be a finite set of pairs that reach the strong pinning criterion, thus generating multi-valued solutions of vortex states; given that defects are dilute, these strong pinning pairs will be dilute as well and hence act independently.

A. Geometry
To find the relevant pairs, we consider two defects labelled by i = 1, 2 at positions (R i , z i ) and the associated vortices R µ(i) separated from the pins by x i = R µ(i) −R i ; the displacement fields u i at z i are solutions of the coupled equations (see Eq. (17)), The coupling g ∈ (0, 1] renormalizes the force at the site of the first impurity due to the action of the second impurity (and vice-versa) and reads While the impurity positions R i and displacements u 1,2 in (39) are continuous variables (with the small scale set by ξ), the vortex positions R µ (1,2) in (40) are restricted to the vortex lattice involving the scale a 0 . A coupling g of order unity implies that both impurities act with their maximal pinning force on the same vortex ; typical separations of such defects are below ξ in the transverse and below a 0 in the longitudinal direction. Hence, large couplings g are associated with close defect pairs lying within a volume a 0 ξ 2 . Small couplings g 1 refer to the situation where the impurities are separated far away from one another, of order several lattice constants a 0 , typically; in this situation, the defects act on different vortices and their mutual effect on the vortex pair is small.
As in Sect. II B, we consider a driving force applied in the positive x-direction and assume that the vortex lattice structure is preserved; under application of the drive, the vortices are displaced from their initial positions R 0 µ(i) by a constant shift of magnitude X along x, i.e, R µ(i) = R 0 µ(i) + X e x . It is convenient to reformulate the problem in terms of the mean vortex position x 2 (x1 + x2) relative to an effective defect characterized by an anisotropic (or non-radial) pinning potential e eff (g, ∆; r), with r = x + u and u = 1 2 (u1 + u2) the displacement of the effective vortex. For b = 0 the effective vortex passes through the center of e eff (r).
(relative to the defects) and the mismatch vector ∆, see Fig. 5. Note that the vortex positions x 1 and x 2 (relative to the defects) as well as the mismatch vector ∆ are restricted to the unit cell of the vortex lattice.
Pushing the vortex lattice a distance X along the xdirection, the mean vortex position x is parametrized as x = [x(X), b] with a fixed impact parameter b, while the vector ∆ remains constant. Since x 1 = x + ∆/2 and x 2 = x − ∆/2 (see Fig. 5), the vector ∆ is interpreted as the mismatch in the pinning by the two defects. If ∆ = 0, the defects are perfectly synchronized: for any X, the position of both vortices relative to the defects is the same, x 1 = x 2 , the pinning forces acting on both vortices are identical, and pinning by the defect pair is maximal. For a finite ∆ = 0, the two vortices are subject to different pinning forces, which reduces the total pinning strength. As shown in Fig. 5(b), the geometry can be reduced to one where an effective vortex at the position x impacts on an effective defect with a non-radial pinning potential e eff (g, ∆; r), with r = x+u and u = 1 2 (u 1 +u 2 ) the displacement of the effective vortex.

B. Averaging
Given the geometric layout of the strong pinning problem with two vortices and two defects, we have to find the associated pinning force density F pin by proper averaging. This averaging involves i) the averaging over trajectories x = (x, b) of vortex pairs with fixed mismatch ∆ and fixed coupling g, ii) the averaging over all possible mismatch vectors ∆, and iii) the averaging over couplings g in the pair pinning Eq. (39) that involves the relative distances between vortices R µ(1) −R µ(2) and the distance in elevation z 1 − z 2 of defects, see Eq. (40). The final result will provide us with a formula, Eq. (46), that expresses the pinning-force density F pin due to defect pairs in terms of the individual pair forces f pair [g(ρ), ∆] for defects separated by ρ and with a mismatch ∆ between them and the vortices. While this expression can be evaluated precisely using numerical techniques, here, we will discuss analytic results that are necessarily of approximate nature.
In the first step i), we fix ∆ and average over the vector x = [x(X), b] in a similar fashion as in Sect. II B. For each b, we average the aggregated pinning force exerted by the two defects while pushing X from large negative to large positive values and then take the average over the impact parameter b. This procedure provides us with the pinning force f pair of a defect pair at fixed mismatch vector ∆ and coupling g, where u 1 (x), u 2 (x) are the solutions to the coupled system, Eq. (39). As before, this average involves the jumps in energy between bistable solutions for the pair-pinning problem defined by Eq. (39).
In the next step ii), we average over mismatch vectors ∆ (the normalization a 2 0 follows from ∆ being restricted to the unit cell of the vortex lattice), where f pair = f pair · (−e x ) denotes the (negative) x-component of f pair ; note that for a vortex lattice pushed along the positive x-direction, the x-averaged pinning force points in the negative x-direction. The ycomponent of f pair vanishes after the ∆-averaging since it is compensated by the configuration with ∆ → −∆ and b → −b. Third, we determine the pinning force density F pin exerted by all defect pairs within a volume V by summing over pairs that are pinned at different separation (2) refers to the separation between the vortices and z 1 − z 2 is the distance between the defects along the z-axis. This final sum (or average) accounts for the dependence of the pair-force f pair (g, ∆) on the coupling g(ρ). Approximating the sum by an integral gives (note the factor of 1 2 to avoid double-counting of the defects) Carrying out one volume integral, we arrive at the final expression It remains to solve the coupled equations (39) and determine the resulting pinning force f pair (g, ∆) of defect pairs that enters the final expression (46) for the pinningforce density F pin . We first provide a qualitative overview of the results, before presenting the detailed derivations.

C. Overview of results
Pinning is maximally strong if both defects are synchronized, i.e., ∆ = 0. In this case, u 1 = u 2 = u and Eq. (39) reduces to a single equation that is equivalent to the single-defect problem with renormalized pinning strength The condition κ eff (g, 0) > 1 for a strong-pinning pair then requires g(ρ) > g 0 (κ) with restricting the maximal separation between the defects, see Fig. 6. In order to arrive at an expression for f pair (g, ∆ = 0), it is convenient to express κ eff in terms of g and the critical value g 0 , With κ eff above but close to unity, we can make use of Eq. (27) and find that scales with (g − g 0 ) 2 ≤ 1. Defect pairs satisfying g(ρ) > g 0 (κ) can pin strongly at a finite mismatch vector ∆ as well. We then have to generalize the effective pair-pinning strength κ eff (g, ∆) to finite ∆ and the strong-pinning condition κ eff (g, ∆) = 1 will provide us with the ∆(g)-domain where pair-pinning is strong. The latter will allow us to determine the pairpinning force f pair (g, ∆). However, a quantitative evaluation of κ eff (g, ∆) for general g and ∆ is quite cumbersome, given the complex geometry of the problem.
Progress can be made by carrying out a perturbative analysis in ∆ ξ, which requires defect pairs to be at the verge of strong pinning, i.e., g(ρ) − g 0 (κ) 1. Such a calculation for the effective pinning strength κ eff (g, ∆) is carried out in Sec. III E; furthermore, it is shown there, that strong pinning with κ eff (g, ∆) > 1 is limited to small mismatches ∆ x < ∆ 0 x and ∆ y < ∆ 0 y with (see Eqs. (97) and (102)) see Fig. 7. An estimate for the ∆-averaged pair force can be obtained by combining the maximal pair force (51) at ∆ = 0 and the region in ∆ where pinning is strong, see Fig. 7, A more precise result for the average pair-pinning force f pair (g, ∆) ∆ is derived in Sec. III F, see Eq. (110). For g(ρ) − g 0 (κ) ∼ O(1), the action of a defect-pair can be analyzed only on a qualitative level. Such large coupling g implies that defects are closeby and the mismatch vector ∆ can go up to the vortex-core radius ξ. The averaging of the pair force f pair ∼ (ξ/a 0 ) 2 f p , see Eq. (27) for κ − 1 of order unity, over ∆ then gives Unlike in the previous case, a precise form for the dependencies of the effective pinning strength κ eff (g, ∆) and the pair-pinning force f pair (g, ∆) on the mismatch vector ∆ cannot be derived analytically. Quantitative insights can be made by a numerical treatment of the problem (see Sec. III F and Appendix B). Such proper averaging over the vector ∆ will only provide a numerical prefactor to the result in Eq. (54), that we do not consider here. Remarkably, Eq. (53) as originally derived from a perturbative analysis in small g(ρ) − g 0 (κ) 1 produces the correct scaling result Eq. (54) also for g(ρ) − g 0 (κ) ∼ O(1). On a qualitative level, we can thus extend the regime of applicability of the expression Eq. (53) to any coupling g within the interval [g 0 , 1].
Substituting f pair (g, ∆) ∆ from (53) into Eq. (46) yields the pinning-force density with the integration over distances ρ restricted through the condition g(ρ) > g 0 . We now proceed to discuss the resulting pinning-force density F pin due to the defect pairs. Of crucial importance in this discussion is the behavior of the coupling g(R, z) that is of order unity at distances R, z < a 0 and rapidly decays further away, see Sec. II A. We distinguish three cases, i) the limit κ → 1 2 , where g 0 → 1, implying that the coupling g > g 0 has to be close to unity. This condition demands that relevant distances R are far below a 0 , from what follows that both defects pin the same vortex. This regime, where defect pairs act on the same vortex, extends throughout all of the regime ii), where 1 2 < κ < 1/[1 + g(R = a 0 )]; in this regime, the coupling g is never large enough to produce strong pinning of different vortices that are always further apart than a 0 . Finally, in region iii), κ → 1 and g 0 → 0, hence, even a small coupling g is sufficient to establish a strong-pinning defect pair. In this case, distances R > a 0 become relevant and different vortices can get strongly pinned to separated defects.
Starting with i), we consider pinning strengths of individual defects close to the threshold κ = 1 2 . The expansion of Eq. (49) for κ − 1 2 1 gives the critical coupling g 0 (κ) ≈ 1 − 4(κ − 1 2 ) close to unity, and the condition g(ρ) > g 0 (κ) requires both defects to act on the same vortex. Using Eq. (9) further implies that g(0, z)−g 0 ≈ 4(κ− 1 2 )−z/a h > 0, hence the maximal longitudinal separation is limited by z 0 ≈ 4(κ − 1 2 )a h . Since for such range of coordinates we have g(ρ) the pinning-force density in Eq. (55) can be cast into the form where we have ignored the logarithmic factor in the expression for the healing length [see the definition of a h below Eq. (9)] and replaced a h ∼ a 0 . The result in Eq. (56) defines the onset of the pinning force density due to defect pairs for κ rising above κ = 1 2 ; it starts dominating over the collective pinning result, Eq. (31), as soon as the pinning strength surpasses the threshold κ = 1 2 by a small amount ∼ (a 0 /λ) 1/2 .
Increasing κ through region ii), i.e., staying below the threshold κ = 1/[1 + g(R = a 0 )], the critical coupling g 0 (κ) will remain finite, of order unity. Referring to Figs. 2 and 3, we note that in this situation, the separation ρ = (R, z) between defects must remain within the peak region III, since otherwise g rapidly decays to a value ∼ a 2 0 /λ 2 1 and the criterion g(ρ) > g 0 (κ) cannot be met. The integral in Eq. (55) is of order unity and the pinning-force density assumes the small-pair or cluster value Finally, when κ resides within the small interval 1 − a 2 0 /λ 2 κ 1, the critical coupling becomes small, g 0 a 2 0 /λ 2 , and the separation between the defects producing strong pair-pinning extends beyond region III into the regions I, II, and even the (yellow) non-dispersive region of Fig. 2. In this case, that includes region iii), we drop the small value g 0 (κ) against g(ρ) in Eq. (55) and rewrite the result for the pinning force density as The integration is dominated by small distances: the small-pair or cluster region (region III in Fig. 2) contributes the same estimate as in Eq. (57), A similar contribution arises from region I: we integrate over the transverse coordinate R < √ a 0 z and find the expression where the main contribution originates from the lower bound z ∼ a 0 (and hence also R ∼ a 0 ). The contributions of region II and the non-dispersive regime are smaller by a factor (a 0 /λ) 3 : In region II with transverse and longitudinal extension R ∼ λ and z ∼ λ 2 /a 0 , the Green's function assumes a constant value with g 7/2 ∼ (a 0 /λ) 7 and the integration gives a result In dealing with the non-dispersive region, we introduce the rescaled distanceρ = [R, (a 0 /4 √ πλ) z]; the nondispersive region is bounded below by |ρ| λ and the condition g(ρ) > g 0 (κ) translates to the upper boundary |ρ| a 2 0 /λg 0 . The integral over the non-dispersive region then takes the form that, given the large exponent 7/2, is determined by the lower boundρ ∼ λ. The total pinning-force density then sums up to In the limit κ → 1, we have g 0 (κ) → 0 and pairs of arbitrarily distant defects induce a finite pinning force. Indeed, the maximal distance between defects [providing the upper bound of the integral in Eq. (62)] diverges as ρ ∼ a 2 0 /λ(1−κ); this translates into maximal longitudinal and transverse separations z ≤ z 0 and R ≤ R 0 between defects with see also Fig. 6. However, the pair-pinning force in Eq. (62) is dominated by the lowerρ bound and contributions from distant defects are irrelevant, implying a finite integral even in the limit κ → 1. The origin of the power α = 7/2 can be traced back to the maximal size of the mismatch vector ∆ ensuring strong pinning, see Eqs. (52). The derivation of the latter requires a detailed quantitative understanding of the pinning mechanism due to defect pairs at the verge of strong pinning that is presented in sections III D, III E, and III F below.

D. Effective pinning potential e eff
In the following sections, we provide a systematic derivation of the results presented above. In a first step, we reduce the two-defect equation (39) to a single equation describing the interaction of a fictitious vortex with an effective pinning potential, see Fig. 5(b). We have seen that for vanishing mismatch ∆, the action of both defects is synchronized and the displacement of both vortices is identical, u 1 = u 2 . For a finite but small mismatch ∆, we reformulate the problem in terms of the mean position r of the displaced vortices relative to the defects and the internal relative position δr (we remind that x i = R µ(i) − R i is the unperturbed defect-vortex distance), Solving perturbatively for the internal coordinate δr will allow us to reformulate the two-defect problem in terms of a single equation for the mean 'fictitious' vortex tip position r. This reformulated problem then will involve an effective pinning potential e eff (g, ∆; r) exerting the pinning force f eff (r) = −∇ r e eff (g, ∆; r) on the 'fictitious' vortex.
In the above new coordinates, the two-defect problem of Eqs. (39) takes the form Expanding in δr to second order and subtracting one equation from the other provides us with an expression for δr, We write the gradient of the radial pinning force f p (r) = f p (r)r (withr = r/r the unit vector in radial direction) in terms of the projectors P ij =r irj and P ⊥ ij = δ ij −r irj , Using δ ij = P ij + P ⊥ ij , Eq. (68) is rewritten as Making use of the relation P α ij P β jk = δ ik δ αβ for the projectors with α, β ∈ { , ⊥}), Eq. (70) is easily inverted and provides a relation between the internal coordinate δr and the mean coordinate r of the vortex pair, Adding the two Eqs. (67) provides us with an equation for the mean vortex tip position, reminiscent of the single-defect case but with an effective pinning force Expanding to 2-nd order in δr (or ∆), the effective pinning force becomes f eff,k = (1+g)f p (r)r k + 1 2 (1+g)∂ i ∂ j f p,k (r)δr i δr j , (75) where the matrix of second derivatives of the pinning force can be expressed as ∂ i ∂ j f p,k (r) = γr irjrk + µ(r i δ jk +r j δ ki +r k δ ij ) with The sums ∂ i ∂ j f p,k (r)δr i δr j involve the expressions Combining Eqs. (75)-(79) gives the effective pinning force up to second order in the mismatch ∆, with the first three terms producing a radial force (along the vectorr), while the last term contributes a transverse force (along the vector ∆ − (r · ∆)r perpendicular tor). This effective pinning force can be written as the gradient of the effective pinning potential e eff (r), f eff (r) = −∇e eff (r), that has a much simpler form. Indeed, using ∇(r · ∆) = [∆ − (r · ∆)r]/r and fixing the integration constant by requiring e eff (r) → 0 as r → ∞, we find the effective pinning potential in the form Using polar coordinates ∆ = ∆(cos θ, sin θ) and r = r(cos ϕ, sin ϕ), the effective potential can be described in terms of the magnitude ∆ of the mismatch and the angle ϕ − θ enclosed by r and ∆, e eff (g, ∆; r) = (1+g)e p (r) −C ∆ 2 (1+g) Inserting the expressions for α ,⊥ , see Eq. (71), and expanding for large distances r ξ (i.e., small values of f p (r), f p (r)/r), we find that with the anisotropic terms appearing at finite mismatch ∆ proportional to f p (r), f p (r)/r ∝ (ξ/r) 4 vanishing faster than the isotropic term ∝ e p (r) ∝ (ξ/r) 2 .
In the limit g = 1 (with defects acting on the same vortex at the same height z), the effective potential becomes e eff (g = 1, ∆; r) = 2e p (r) + 1 4 (r · ∆) 2 f p (r) that corresponds, up to order O(∆ 2 ), to the simple superposition of two mutually shifted pinning potentials, e eff (g = 1, ∆; r) ≈ e p (r + ∆/2) + e p (r − ∆/2). In the limit g = 0, the two-defect problem (67) decouples and we can obtain independently each perturbed vortex tip position E. Effective Labusch parameter κ eff and strong-pinning range ∆ 0 We proceed with the calculation of the effective Labusch parameter κ eff (or pinning strength) for the anisotropic pinning potential of Eq. (82). For a single isotropic defect, the Labusch parameter κ is defined in Eq. (21) and involves the (maximal) potential curvature f p = −e p and the effective elasticityC. Going to the defect pair, the anisotropic potential (82) depends on the distance and arrangement of defects through the parameters g and ∆ = ∆(cos θ, sin θ). In the following, we consider a vortex with an asymptotic trajectory x = (x, 0) and determine the angular dependence (on θ) of the effective Labusch parameter κ eff (g, ∆) = κ eff (g, ∆, θ), i.e., we consider defect pairs with different angular arrangement relative to an asymptotically fixed vortex motion. Once the function κ eff (g, ∆) is known, the condition κ eff (g, ∆ 0 ) = 1 will provide us with the maximal misfit ∆ 0 (g) limiting the strong-pinning range. Trajectories with different angle of incidence and/or finite impact parameter b will be discussed later. Finally, the pinning strength in the vortex-defect system can be tuned by either changing the effective elasticityC or the energy scale e p of the defect potential-in the present discussion, we will tune κ eff via changingC.
We start our derivation of the pinning strength κ eff (g, ∆) by going back to its defining equation. For an isotropic pinning potential, this is given by (21), that in turn derives from the self-consistency equation [4, 5, and 28] Equation (85) allows us to connect incremental changes in the asymptotic and tip positions, δx = [1−f p (r)/C]δr, with jumps in δr occurring when 1 − f p (r)/C = 0. Combining this relation with the condition of its first appearance, f p (r) → max[f (r)], leads to (21). Finally, the maximum force derivative max[f (r)] is achieved at the inflection point r m defined via f p | rm = 0.
In the present anisotropic situation, Eq. (85) has to be generalized to its vectorial form (73) and incremental changes in asymptotic and tip positions of vortices are related via where is the Hessian matrix associated with the pinning energy landscape e eff (g, ∆; r). The Labusch criterion again marks the first appearance of an instability in the vortex tip position δr. We thus have to invert Eq. (86) and find the solution δr(δx)-a diverging result for δr then signals the presence of a jump in the vortex tip position. Approaching this divergence from the weak pinning domain, i.e., starting with a largeC and decreasing its value, the jump appears when the determinant in the matrix relation (86) vanishes, Evaluating the Hessian in cylindrical coordinates (r, ϕ), we obtain the matrix with β(r, θ) = −1 8 and functions γ(r, θ), δ(r, θ) that we do not need to calculate explicitly. Above, we have made use of the fact that the tip trajectory stays always close to the x-axis (up to corrections of order ∆ 2 ) and hence, we have set the angle ϕ in β(r, θ) to zero, β(r, θ) → β(r, θ). The condition of vanishing determinant (88) is equivalent to matching up the lower eigenvalue λ − (r, θ) < 0 of H withC, λ − (r, θ) +C = 0; furthermore, we need to find the location where this happens first, i.e., we have to determine the distance r eff m (θ) that generalizes r m to the anisotropic situation. Once this program is executed, the generalized Labusch parameter is given by that assumes unity at the weak-to-strong pinning transition and larger values on decreasingC further into the strong pinning region. Let us first consider the above generalized formulation of the Labusch parameter for the isotropic situation with ∆ = 0. Then H ij is already diagonal, with eigenvalues λ − (r) = −f p (r) < 0 and λ + (r) = −f p (r)/r > 0 close to the inflection point r m , where f p (r m ) = 0 and the maximum in −λ − = f (r) is realized. These results are fully in line with the previous discussion of the Labusch criterion (21).
The perturbative analysis of the anisotropic situation contributes corrections to order ∆ 2 that introduce an angular dependence of the results on θ. The eigenvalues λ ± (r, θ) of H, see (89), coincide, to order ∆ 2 , with its diagonal entries, since the off-diagonal terms only add a correction γ 2 ∆ 4 to the determinant appearing in their calculation. In particular, the lower eigenvalue assumes the form Following the definition in Eq. (90), we have to evaluate this expression at the generalized inflection point r eff m (θ). The latter remains close to r m , r eff m (θ) = r m + O(∆ 2 ), and using this Ansatz in (91), we find that the correction to r m is irrelevant since We thus arrive at the formal expression for the effective Labusch parameter Inserting the expression for β(r, θ) from above, we can rewrite this result into the convenient form When the vortex trajectory is oriented at a finite angle φ with respect to the x-axis, the angular dependence in (93) has to be replaced according to θ → θ − φ. Finally, the discussion for a finite impact parameter b can be easily reduced to the situation where the vortex approaches the defect center from an angle, see the discussion in Sec. III F below. Note that, while r eff m ≈ r m does not depend on angle to order ∆ 2 , the effective pinning strength κ eff (g, ∆, θ − φ) experienced by a vortex incident at an angle φ does. This implies that pinning is not equally strong when approaching the same defect from different direction. Rather, κ eff may be larger (or smaller) than unity when changing φ ∈ [0, 2π]. As a result, vortex trajectories crossing the same defect may undergo pinning and depinning jumps in some directions but not in other.
Next, we return to our vortex incident along x, substitute the anisotropic defect-pair potential Eq. (82) into the expression (93), and use κ = f p (r m )/C to find the explicit result Setting the pinning strength to its critical value, κ eff (g, ∆) = 1, we now can determine the misfit parameter ∆ 0 (g) below which pinning is strong. We first analyze the two special cases θ = 0 and θ = π/2, where the mismatch ∆ is parallel and perpendicular to the vortex trajectory, before generalizing the result to other angles θ.
For θ = 0, we obtain Since f p (r m ) < 0, the pinning strength decreases below its maximal value (1+g)κ as the mismatch ∆ is increased. By setting κ eff (g, ∆, 0) = 1 in Eq. (95), we can obtain the maximum longitudinal mismatch ∆ 0 x below which the pinning by the two defects is strong, Estimating f p ∼ f p /ξ 3 and f p /Cξ ∼ O(1) and dropping the factor 1 + g in the denominator, we can express the relevant g-dependence of ∆ 0 x in the parametric form For θ = π/2, we first rewrite the prefactor of the sin 2 θ term in Eq. (94). Introducing the notation β(r) = f p (r)/r, we find that Expressing the effective elasticityC through κ and β, C = κ −1 (βr) r=rm , the second factor in this expression simplifies to where in the last two steps we used (rβ) = f p (r m ) = 0 and β = −(2/r)β . Since f p (r m ) > 0, we have β > 0 at r = r m and hence the effective Labusch parameter κ eff (g, ∆, π/2) = (1+g)κ (100) again decreases with increasing ∆. Setting κ eff (π/2) = 1 then defines the maximal transverse mismatch for strong pinning, Using similar estimates as above, we obtain the result in parametric form It is interesting to compare the pinning strengths for different vortex-defect configurations where the vortex trajectories are either parallel or perpendicular to the mismatch vector ∆. Combining Eqs. (96) and (101) and using g 0 = 1/κ − 1, we obtain the following ratio for the maximal longitudinal and transverse mismatches ∆ 0 For a maximal coupling between defects, g = 1, such that the defect potentials directly add up [cf. Eq. (84)], the factor [· · · ] 3/2 in the above expression is unity. The ratio of the two scales then depends on the specific form of the pinning potential in the vicinity of the inflection point r m . For the Lorentzian pinning potential producing the pinning force in Eq. (16), the ratio of the two scales reads ∆ 0 y /∆ 0 x | g=1 = 3/2 and pinning is always stronger, i.e., κ eff (θ) is larger along the direction perpendicular to the mismatch vector. Furthermore, since β(r m ) < 0, the ratio ∆ 0 y /∆ 0 x grows as g is decreased away from unity towards its minimum value g 0 = 1/κ − 1. When g = g 0 , the term 1 − κ(1 − g) in the denominator of Eq. (104) takes the value 2(1 − κ), hence as κ → 1, pinning due to pairs of distant defects is always stronger in the direction perpendicular to the mismatch, regardless of the specific form of the pinning potential.
The maximal longitudinal and transverse mismatches ∆ 0 x and ∆ 0 y allow us to identify the region of applicability of the perturbative approach that we have used to derive the effective pinning potential. For a longitudinal mismatch with θ = 0, Eq. (72) tells us that δr ≤ ∆ 0 x /2α and using the bound α ≥ 1−κ(1−g) = κ(g+g 0 ), we find that δr ξ(g−g 0 ) 1/2 . The perturbative approach is valid provided that δr ξ which is the case for g − g 0 1. Similarly, for a transverse mismatch, we use α ⊥ ∼ O(1) and therefore δr ≤ ∆ 0 y /2α ⊥ ∼ ∆ 0 y . The perturbation δr then again remains small for g − g 0 1, except for the crossover to strong pinning where κ = 1, g 0 = 0, and ∆ 0 y ∼ ξ. In this case, the validity of the perturbative approach requires ∆ y ξ (see Eq. (72)), however, a comparison with the numerical results in Fig. 8 demonstrates that the perturbative approach still provides excellent agreement even for a mismatch ∆ y comparable with ξ.

F. Average pinning force
In this section, we use the findings on the effective Labusch parameter to calculate the average pinning force due to two defects coupled by g. Restricting first to vortex trajectories x = [x(X), 0] passing through the center of the effective pinning potential, we evaluate the xaveraged pinning force for the longitudinal component (along the x-direction), see Eq. (22), where ∆e eff,1 pin and ∆e eff,2 pin denote the jumps in the effective pinning potential defined as in Eq. (19), i.e., e eff pin (x) = e eff [r(x)] + 1 2C [x − r(x)] 2 . Provided that κ eff (g, ∆, θ) > 1, the jumps ∆e eff,i pin , i ∈ {1, 2} in the pinning energy are given by Eq. (26) with the replacements f p (r m ) → f eff (r m ) and κ → κ eff ; otherwise, for κ eff (g, ∆, θ) < 1, the effective pinning force vanishes after the averaging. This gives the two-defect pinning force Using Eq. (103) for the effective Labusch parameter and noting that f eff (r m ) = (1+g)f p (r m )+O(∆ 2 ), we obtain The pinning force thus decays with ∆ from its maximal value ∼ f p (g − g 0 ) 2 at ∆ = 0 to zero as the mismatch increases to ∆ 0 x and ∆ 0 y along and perpendicular to the vortex trajectory, respectively, see Fig. 7. In Fig. 8, we compare this analytic formula to the numerical results (see Sec. B) in the regime g − g 0 1 for different angles θ. Note that, while the (perturbative) analytic results assume a small mismatch ∆ ξ, they remain applicable for angles close to θ = π/2, where ∆ becomes comparable to ξ.
The geometric complexity arising at finite impact b, see Fig. 9, produces interesting new features, e.g., asymmetric pinning and depinning jumps or even trajectories with only one of the pinning/depinning jumps real- ized. In the generic situation, the vortex tip associated with the trajectory x = [x(X), b] undergoes a jump every time the position r hits the distance r m + O(∆ 2 ) from the center of the effective pinning potential. At the instance of the jump, the vortex asymptotic position is x = r − f eff (r)/C [see Eq. (73)], i.e., at a distance x p = r m − (1 + g)f p (r m )/C + O(∆ 2 ) from the center of the pinning potential. Pinning and depinning then occur at the asymptotic positions x = (±x p cos φ, x p sin φ), with φ ≈ arcsin(b/x p ). Note that for the single-defect pinning in Fig. 4 corresponding to a large κ, the energy jumps associated with pinning and depinning appear at different asymptotic distances x − , x + from the defect center; for marginally strong effective pinning with κ eff (g, ∆) − 1 1, we can neglect the difference x + − x − ∼ ξ(κ eff − 1) 3/2 (see Refs. [27 and 33]) and set The angles enclosed between x and ∆ at the pinning and depinning events are θ p = θ + φ and θ dp = θ − φ, see Fig. 9. The jump size at the pinning transition does not depend on the direction of the vortex motion and thus can be evaluated from the vortex trajectory passing directly through the center of the effective pinning potential but at an angle θ + φ, i.e., ∆e 1 eff (g, ∆, θ, b) = ∆e 1 eff (g, ∆, θ + φ, 0). Similarly, ∆e 2 eff (g, ∆, θ, b) = ∆e 2 eff (g, ∆, θ − φ, 0) at depinning. Hence, the pinning and depinning jumps assume different values at finite impact b. The pinning force is then expressed as f pair (g, ∆, θ, b) = 1 2 [∆e 1 eff (g, ∆, θ + φ, 0) + ∆e 2 eff (g, ∆, θ − φ, 0)]/a 0 . Furthermore, since for marginally strong pinning and b = 0 trajectories, the pinning and depinning jumps in energy are equally-sized, we express the resulting pinning force in terms of the forces exerted on the b = 0 trajectories, With the pinning and depinning jumps no longer equal, we may encounter situations where one of the jumps is absent. This is the case for misfits ∆ with θ = 0 and b = 0; e.g., when κ eff (g, ∆, θ+φ) > 1 but κ eff (g, ∆, θ−φ) < 1, the vortex undergoes a pinning jump in energy when its asymptotic trajectory passes the circle of radius x p for the first time but does not undergo any depinning jump when the asymptotic trajectory crosses the circle a second time, see Fig. 10(b) (note that the corresponding trajectory with opposite impact parameter y = −b will undergo a depinning jump but will not jump upon pinning).
The pinning force averaging is done through integra-  A vortex passing at the transverse distance b = 2.5 ξ experiences an effective Labusch parameter κ eff (g, ∆, φ) = 1.011 (with φ = arcsin(b/xp)) and the vortex undergoes pinning and depinning jumps (green dots) when its asymptotic position crosses the circle of radius xp. The blue arcs denote the sections on the circle of radius xp where the effective pinning is strong; for g − g0 1 and κ approaching unity, pinning is stronger along the direction perpendicular to the mismatch vector ∆ (that is for φ = π/2) than along the direction parallel to the mismatch, see the discussion below Eq. (103). The tip position (right plot) jumps at the pinning and depinning events, as illustrated by the two pairs of green dots. (b) Mismatch vector enclosing an angle θ = π/4 with the vortex trajectory. At this angle, pinning is strong, κ eff (g, ∆, θ) = 1.005, and a vortex with vanishing impact parameter b = 0 undergoes both pinning and depinning jumps. However, for the vortex passing at b = 2.5 ξ, the corresponding effective Labusch parameters read κ eff (g, ∆, θ+φ) = 1.014 and κ eff (g, ∆, θ − φ) = 0.996 and the vortex undergoes a jump only upon pinning. The right plot shows the corresponding tip trajectories with its jumps. tion over the mismatch ∆ and the impact parameter b, The integration over ∆ is restricted to the stronglypinning region, i.e. κ eff (θ) > 1. Rewriting the integration in terms of ∆ = (∆ x , ∆ y ), we compute the factor and obtain the averaged pinning force Keeping systematically the corrections in the jump position r m + O(∆ 2 ) is equivalent to replacing b = x p sin φ + O(∆ 2 ) in the integration leading to the Eq. (109). Carrying out the integration of the additional O(∆ 2 ) term would contribute with a quartic correction ∝ ∆ 4 to Eq. (110) which we ignore here. With x p ∼ ξ and the results for the maximum mismatches ∆ 0 x , ∆ 0 y in Eqs. (97) and (102), we obtain a parametric estimate for the pinning force originating from defect pairs in the form, with the following interpretation: the defect pair induces the pinning force ∼ f p rescaled by the factor (g − g 0 ) 2 that accounts for the distance between defects, a factor ξ 2 /a 2 0 ∼ S trap /a 2 0 that represents the areal fraction where vortices are trapped, and the additional factor ξ 2 /a 2 0 together with the distance-dependent factor (g − g 0 )(g + g 0 ) 1/2 that derive from the constraint on the mismatch ∆. The result Eq. (111) is the basis for the evaluation of the pinning force due to defect pairs at different separations in Sec. III C.

IV. SUMMARY AND CONCLUSION
We have extended the strong pinning paradigm into the weak pinning domain by accounting for correlations between defects. The most relevant correlations arise from defect pairs-they reduce the critical Labusch parameter (or pinning strength) κ ∼ −e p /C for strong pinning from its standard value κ c = 1 to κ c,pairs = 1/2. When decreasing the individual defect's pinning strength κ towards the critical value κ c , the strong pinningforce density vanishes as F pin ∝ n p (κ − 1) 2 and strong pair-pinning takes over. Upon a further decrease of κ, the pair-induced strong pinning-force density F pin scales with n 2 p and vanishes at κ = 1/2 according to F pin ∝ n 2 p (κ − 1/2) 4 . The contributions to F pin from higher-order correlations between n defects scale as (n p ) n and quickly become irrelevant, with the dominant contribution to the pinning force density being taken over by weak collective pinning with F pin ∝ n 2 p κ 3 as κ drops below 1/2.
The origin of the pair-induced strong pinning condition κ > 1/2 is easily understood-for two defects that overlap in position, their joint pinning strength doubles and they reach the strong pinning criterion 2κ > 1. The substantial enhancement in pinning strength remains in place for small defect pairs that are separated at most by ∼ ξ and ∼ a 0 in transverse and longitudinal (field-) directions and act on the same vortex. However, this is not the full story: with κ > 1/2 approaching unity from below, pairs separated by transverse distances beyond a 0 can constitute a strong-pinning pair as well. These pairs, rather than pinning the same vortex, will pin different vortices. The interaction between these two defect-vortex entities is transmitted by the elastic properties of the vortex lattice, specifically, the static nonlocal Green's function G αβ (R, z), and determines the effective pinning strength of the extended pair which is smaller than the one of a small pair.
The Green's function G αβ (R, z) describes the displacement field u(R, z) for a δ-force acting at the origin and hence the distortion at the site of a second defect that is positioned a distance (R, z) away from the first defect. In our analysis, we have simplified the expression for the Green's function and considered its diagonal, transverse part G(R, z); the result, shown in Fig. 3, exhibits a sharp asymmetric and structured peak in the shape of a dumbbell. This complex real-space structure has not been considered before and is expected to be present in the full expression for the response matrix G αβ (R, z) as well.
For our extended pairs, the effective Labusch parameter or pinning strength κ eff , rather than simply doubling κ, scales as κ eff ∼ (1 + g) κ, with g = g(R, z) = G(R, z)/G(0, 0) < 1; hence, the partner defect contributes to the strong pinning with a reduced weight. Extended pairs within a distance determined by the condition g((R, z) > 1/κ − 1 ≡ g 0 (κ) thus potentially contribute to strong pair-pinning; geometric considerations refine this analysis and produce an effective Labusch parameter κ eff (g, ∆) that depends on distance (through g) and on the misfit ∆ between the defect-pair and the vortex lattice, with a finite ∆ further reducing the effective pinning strength κ eff , see Eq. (103).
The effective Labusch parameter κ eff (g, ∆) exhibits a non-trivial angular dependence encoded in the direction of ∆. While for isotropic single-defect pinning, strongpinning jumps appear near the inflection points arranged in a circle of radius r m , for an anisotropic potential as in Eq. (83), strong-pinning jumps appear on arcs that grow with decreasing elasticityC or increasing pinning strength e p as illustrated in Figs. 10. The direction away from the defect center where these arcs make their first appearance depends, besides the direction of ∆, on the detailed shape of the pinning potential, see Eq. (104).
With contributions to strong pair-pinning arising both from small pairs pinning one vortex and extended pairs pinning two separated (and relatively misfitted) vortices, the question arises about their relative total weight. It turns out, that the extended-pair force decreases with the scaled distanceρ = [R 2 + (a 2 0 /16πλ 2 )z 2 ] 1/2 as ∝ρ −7/2 , that makes the small-pair contribution (originating from pairs in a small volume ξ 2 a 0 ) dominate the strong-pair pinning-force density F clust in Eq. (32).
As follows from the above discussion, the elastic properties of the vortex lattice take an important role in the calculation of the strong-pair pinning force. Furthermore, they also define the dominance of strong-pair pinning over weak-collective pinning that is reduced by the factor (a 0 /λ) 2 . This reduction is a consequence of the non-local interaction between vortices producing a dispersion in c 44 (k). While collective pinning in the nondispersive regime (with a Larkin length R c > λ) involves a 'stiff' lattice with c 44 (0) = B 2 /4π, the small defect pairs involve the soft lattice with c 44 (k) = B 2 /4πλ 2 k 2 ; with the relevant k ∼ K BZ ≈ √ 4π/a 0 , the lattice is softer by a factor ∼ a 2 0 /λ 2 and hence pair-pinning is stronger. The large factor λ 2 /a 2 0 also guarantees, that strong-pair pinning is larger than weak-collective pinning deep in the dispersive regime where the lattice becomes softer.
It is interesting to compare the situation described above with the one studied by Fisher [32]: Focusing on weak defects and large dimensions D > 4, it turns out that weak-collective pinning is ineffective due to the fast spatial relaxation of the manifold's distortions. Pinning then is exclusively due to rare configurations appearing in the random pinning landscape. In our analysis, we start from the opposite limit, strong defects that pin the elastic manifold (here, vortices) individually. Upon decreasing the pinning strength κ below unity, we loose the strong pinning of individual defects and would expect weak-collective pinning to take over in D = 3, where distortions decay slowly, proportional to the inverse distance. Instead, due to the non-local interaction between vortices producing a dispersive elastic response, we find that specific rare events, small defect pairs, take over and produce the leading pinning mechanism. In a nondispersive elastic medium, the stiffening at large scales is absent and the two types of pinning, rare and collective, come with equal (parametric) weight.
While weak collective pinning arises from typical fluctuations in the defect distribution, strong pair-or clusterpinning arises from rare fluctuations. In reality, both types of fluctuations coexist and hence simultaneously contribute to the pinning force density F pin . Similar to the addition of resistivities arising from different scattering mechanisms in the Matthiessen rule describing metallic transport, the pinning-force densities from different pinning mechanisms should be added up to the total pin-ning force F pin ≈ F coll + F clust when describing the transport in a superconductor. However, given the inductive response of a superconductor, this corresponds to an addition of (critical) currents rather than voltages. Microscopically, comparing the distance between small pairs, d pairs ∼ [n p (n p a 0 ξ 2 )] −1/3 , with the size of the collective pinning volume V c ∼ L c R 2 c ∼ (λ/a 0 )R 3 c , one notes that V c contains many pairs. Hence, when dragging a vortex system slowly over the pinning landscape, one should observe a complex stick-slip type motion where small slips of individual vortices depinning from defect pairs combine with large slips of collectively pinned vortex bundles. It would be interesting to observe the motion of such a pinned vortex system in a numerical simulation.
Another future topic of interest is the further investigation of the real-space structure of the Green's function G αβ (R, z), both theoretically as well as experimentally. In particular, it would be interesting to come up with a proposal for an experiment that is sensitive to the nontrivial dumbbell structure of the peak in the response function.
(even though only the cases with zero, one, or two jumps have been discussed in the analytic part of this paper).
As for the single-defect case, it turns out that the pinning force exerted on the vortices can be expressed as the gradient of the total pinning energy. Indeed, taking the gradient of e pair pin (x) ≡ e pair pin [x, ∆; u 1 (x), u 2 (x)] defined through Eq. (B1), we find where the partial derivative ∂e pair pin /∂x is taken at fixed u 1 , u 2 . The term in Eq. (B2) involving the u-derivatives vanishes since ∂e pair pin /∂u 1,2 = 0 at the minimum. Taking the partial x-derivative in Eq. (B1) provides us with −∇ x e pair pin = f p [x + ∆/2 + u 1 (x)] + f p [x − ∆/2 + u 2 (x)], which is precisely the pinning force exerted by the two defects on the distorted vortices. The x-averaged pinning force along the trajectory x = (x, b) is then written as with ∆e pair,j pin = lim ε→0 [e pair pin (x j − ε, b) − e pair pin (x j + ε, b)] quantifying the energy jump at the position x j . Integrating ∂e pair pin /∂y would in general give a non-vanishing contribution to the pinning force in the y-direction; it is however compensated by the configuration with ∆ → −∆ and b → −b after averaging.
The result of this numerical evaluation is compared with the analytic result Eq. (107) in Fig. 8 for the case of marginally-strong pair-pinning g − g 0 1 and shows a very good agreement with the analytic result even at large values of the mismatch ∆ of order ξ at angles θ close to π/2, in which case the theoretical estimates made in Sec. III E do not guarantee the validity of the perturbative approach. For a Lorentzian shape potential and parameters used in Fig. 8 with θ = π/2, we find that the scaling factor in Eq. (101) assumes a value [(g − g 0 )/(1 + g)(g + g 0 )] 1/2 ≈ 0.17 and the prefactor contributes a factor ≈ 4.46, such that ∆ 0 y ≈ 0.78 ξ.

Numerical minimization
The main challenge in the minimization of the twodefect pinning energy in Eq. (B1) is to properly track the local minimum representing the current state of the vortex (the occupied branch) and to ensure that the minimization algorithm does not overshoot to another minimum.
We define u = (u 1,x , u 1,y , u 2,x , u 2,y ) and minimize the function e pair pin (u, x), see Eq. (B1), with respect to u. For a fixed asymptotic position x, we use Newton's method to iterate u, as long as the matrix of second derivatives H αβ = ∂ 2 e pair pin /∂u α ∂u β evaluated at u i remains positive-definite. The parameter γ is chosen to bound the step size |u i+1 − u i | < δu max to a pre-defined maximum value δu max . The method converges if the initial guess u 0 lies close to a local minimum and the new minimum is used as an initial guess for the minimization after changing the parameter x.
The appearance of at least one negative eigenvalue of H(u) during the minimization signals the disappearance of the local minimum and triggers the jump to another minimum. Minimisation through this region is performed using the Nelder-Mead algorithm with the initial simplex size set to δu max . Once the positive-definite region in the neighborhood of the new minimum is reached, the minimisation procedure switches back to Newton's method.
Obtaining the size of the energy jumps to the desired accuracy requires precise determination of the jump points x j = (x j , b) where the currently occupied local minimum disappears. This is achieved by repeated interval-halving: assume that for x = (x 0 , b) the Newton minimization converged to a local minimum u 0 , that is used as a starting point for the minimisation at the next position x = (x 0 +δx, b). The appearance of a region with a negative eigenvalue of the Hessian during this minimisation indicates the presence of the jump point x j in the interval (x 0 , x 0 + δx). Another minimisation is thus performed for x = x 0 + δx/2 that reduces the interval either to (x 0 , x 0 + δx/2) (if the negative eigenvalue region appears during the minimisation) or (x 0 + δx/2, x 0 + δx). The further iteration of this procedure locates the jump point to the required precision.