The elastic depinning transition of vortex lattices in two dimensions

Large scale numerical simulations are used to study the elastic dynamics of two-dimensional vortex lattices driven on a disordered medium in the case of weak disorder. We investigate the so-called elastic depinning transition by decreasing the driving force from the elastic dynamical regime to the state pinned by the quenched disorder. Similarly to the plastic depinning transition, we find results compatible with a second order phase transition, although both depinning transitions are very different from many viewpoints. We evaluate three critical exponents of the elastic depinning transition. $\beta = 0.29 \pm 0.03$ is found for the velocity exponent at zero temperature, and from the velocity-temperature curves we extract the critical exponent $\delta^{-1} = 0.28 \pm 0.05$. Furthermore, in contrast with charge-density waves, a finite-size scaling analysis suggests the existence of a unique diverging length at the depinning threshold with an exponent $\nu= 1.04 \pm 0.04$, which controls the critical force distribution, the finite-size crossover force distribution and the intrinsic correlation length. Finally, a scaling relation is found between velocity and temperature with the $\beta$ and $\delta$ critical exponents both independent with regard to pinning strength and disorder realizations.


Introduction
The physics of disordered elastic systems is relevant for numerous systems, e.g. interfaces like magnetic [1] or ferroelectric [2] domain walls, contact lines [3], crack propagation [4], or periodic structures like charge-density waves [5] (CDWs), vortex lattices in type II superconductors [7], colloids [10] or Wigner crystals [11]. It is essential to understand the response of these systems to an external driving force, such as current-induced Lorentz force for vortices, electric field for charge-density waves, or magnetic field for domain walls. The competition between elasticity and disorder in all these systems leads to a great variety of phases. In particular the system remains pinned up to a critical driving force F c and starts sliding above F c . The transition from a pinned state to a sliding one is known as the depinning transition. It has been suggested by Fisher [12], in the context of CDWs, that the depinning transition could be regarded as a standard equilibrium critical phenomenon where the reduced force f = (F − F c )/F c and the velocity v would respectively act as the control parameter and the order parameter. This idea has proven useful in many other systems. This analogy suggests the existence of critical exponents, in particular β which describes the power-law response of the velocity at zero temperature v ∼ f β . However the analogy with standard critical phenomena is limited since the depinning transition is by nature a non-equilibrium transition. For example, unlike equilibrium phase transitions no divergent steady-state correlation length scale exists for an elastic line when approaching the depinning from below [13]. Connections between the depinning transition of elastic interfaces in random media and non-equilibrium phase transitions into absorbing states are reported in the literature [15]. The pinned phase where the dynamics is frozen may be considered as an absorbing state in which the system may enters but cannot leave. Two universality classes have been proposed for the interface growth. The linear growth of the interface is described in the quenched Edwards-Wilkinson equation within the random field Ising class [19,20,21], while Kardar-Parisi-Zhang nonlinearities in this equation make the depinning transition closely related to directed percolation [22]. Non-equilibrium phase transition into absorbing states has also been shown in periodic systems with quenched disorder like colloids: when periodically sheared a transition was observed from a diffusive liquid to a pinned state [24], while in the case of driven colloids the transition separates plastic flow from a pinned state [25]. The critical exponents found in both studies are consistent with aborbing state transitions in the universality classes of directed percolation or conserved directed percolation.
When the strength of the underlying disorder is strong, the theoretical description of the phenomenon is difficult and the question about the nature of the depinning transition remains an open problem. In periodic systems with a 2D displacement field of dimension (N = 2) (e.g., superconductor vortices, colloids and Wigner crystals), simulations [26,34,37] have shown that this situation leads to dislocations and plasticity: at the depinning threshold tearing of the lattice appears where regions of pinned particles coexist with particles flowing around them. When the pinning of the underlying sub-strate is weak the statics and dynamics are dominated by the system elasticity which favors order. After the work of Fisher, further analytical and numerical works have refined the analogy with equilibrium critical phenomena to define universality classes and critical exponents (see e.g. [38] and references therein). However in numerical simulations of periodic systems in d = 2 various exponents have been found. β ≈ 2/3 is found for CDWs [39], Wigner crystals [40], and colloids [37], while other works found β ≈ 0.5 [41] and β = 0.92 ± 0.01 [42] in colloids, and β = 0.35 [43] in stripe systems. For superconductor vortices in d = 3 the value β = 0.65 is found [44], while in d = 2 a surprising exponent β > 1 has been reported [45].
In this paper, we show a detailed study of the critical behavior of superconductor vortices above the elastic depinning transition, including a finite-size scaling analysis and the determination of three critical exponents. The elastic depinning transition investigated in this paper is very different from the plastic depinning transition studied in a previous paper [46]. In the plastic case tearing of the lattice at the depinning threshold is observed with moving vortices flowing along changing interconnected channels around pinned regions of vortices. On the contrary, at the elastic depinning threshold, all vortices depin together and flow along well-defined rough static channels in a topologically ordered structure. Furthermore, the velocity is periodic in time at the elastic threshold, while it is shown to be chaotic at the plastic threshold [34] with large broad-band noise at low frequency. The curvature of the velocity-force curve is also very different since the elastic depinning is caracterized by a velocity exponent β < 1 whereas β > 1 for the plastic depinning.
We perform large scale molecular dynamics simulations of 2D superconductor vortex lattices driven over weak random disorder. Our model belongs to the category of 2D periodic systems with N = 2 and short-range interactions. It could model 3D superconductors (either conventional or layered) in an effective 2D regime, i.e., when the vortex line tension is high enough for the lines to remain straight. The numerical model is developed in section 2. The depinning transition at zero temperature is detailed in section 3 for different lattice sizes and disorder realizations. Above the depinning threshold a power law response v F >Fc ∼ f β is found with β = 0.29 ± 0.03, and the transition seems to be governed by the correlation length divergence with an exponent ν = 1.04 ± 0.04. The same exponent governs the finite-size effects. In section 4 we study the depinning transition at finite temperature. A power law response for the velocity v F =Fc ∼ T 1/δ is found and the critical exponent δ −1 = 0.28 ± 0.05 is determined. A scaling relation for the driving force and the temperature is established. Varying system size, disorder realization and pinning strength, the scaling function and both β and δ critical exponents show some degree of universality. The critical exponent values of our study are discussed in section 5.

Numerical model
We study N v Abrikosov vortices interacting with N p random pins in the (x, y) plane. Periodic boundary conditions are used in both x and y directions. We consider the London limit λ L ≫ ξ s , where λ L is the penetration length and ξ s is the superconducting coherence length, i.e. we treat vortices as point particles. As in [47], the LAMMPS classical molecular dynamics code [48] is used to integrate the newtonian equations of motion with a velocity Verlet algorithm [49]. We use this parallel code in order to simulate large system sizes, which is crucial to study the elastic depinning since the Larkin domains are larger than in the case of plastic dynamics. However the LAMMPS code includes an inertial term proportionnal to the particle mass. The existence of a vortex mass has been discussed in the literature [50,55]. In general vortices are considered as massless objects although the vortex mass may be significant in the clean limit. In recent experiments [55], it is shown that Nb vortices close to T c ≈ 9K have an inertial mass which is in the intermediate range between the dirty and clean limit, showing that the inertial vortex mass should be a meaningful concept.
We describe the time evolution of a vortex i at a position r i with the equation where r ij is the distance between the vortices i and j located at r i and r j , r ip is the distance between the vortex i and the pinning site located at r p , and ∇ i is the 2D gradient operator acting on r i . F L = Fx is the Lorentz driving force in the x direction induced by an applied current, and η is the viscosity coefficient. The Langevin force F th i describes the coupling with a heat bath and can be expressed by a thermal gaussian white noise with zero mean and variance where λ, µ = x, y are the cartesian components of F th i and k B is the Boltzmann constant. The vortex-vortex pairwise repulsive interaction is given by a modified Bessel function and the attractive pinning potential is given by where R p the radius of the pins, and α v and α p are tunable parameters. We fix the strength of the vortex-vortex interaction by setting α v = 2.83 10 −3 λ L ǫ 0 where ǫ 0 is a characteristic energy per unit length. The number of pins is set to N p = N v and their radius is R p = 0.22λ L . The positions of the pins are chosen from a uniform random distribution. Each random sampling will be referred to as one "sample" or one "disorder realization". The average vortex distance is a 0 = λ L , and the vortex-vortex interaction is dealt with using a neighbor list method with a cutoff radius r c = 6.5 λ L . Since we are mainly interested in the permanent regimes we choose η/m = 0.1 (where m is the vortex mass) in the LAMMPS code such that the inertial term is small compared to the viscous term. Care has been taken to check that such ratio ensures that the long time behavior of the second order Newton's vortex dynamics is identical to the overdamped dynamics limit usually computed for the superconductor vortices. In particular we find identical vortex trajectories and identical values of the critical exponents using overdamped dynamics on smaller systems.
We use a unit system in which η = 1, λ L = 1, ǫ 0 = 1 and k B = 1. Depending of the relative strengths of the vortex-pin and vortex-vortex interaction, the dynamics can be either dominated by elasticity or disorder. Several values of the relative disorder strength α p /α v have been investigated. In sections 3.1, 3.2 and 4.1 we choose the relative disorder strength α p /α v ≈ 5.10 −3 , weak enough to enforce elastic behavior (plasticity is found above α p /α v ≈ 0.01). In section 4.2 we investigate several pinning strengths in the elastic regime, from α p /α v ≈ 2.10 −3 to α p /α v ≈ 7.10 −3 . Molecular dynamics simulations are performed for several system sizes, from N v = 270 up to N v = 20000 vortices, but keeping fixed vortex and pin densities. Moreover, various rectangular shaped basic cells of size (L x , L y ) have been investigated, from the almost square (L x , L y ) = (5, 6 √ 3/2)nλ L with n = 3, 8, 20 up to very elongated strips (L x , L y ) = (400, 20 √ 3/2)λ L . The strip geometry elongated in the driving force direction allowed the study of the critical depinning properties for large system sizes. Care has been taken to check that the basic cell anisotropy induced by the strip geometry does not alter the critical properties, and in particular the value of the velocity critical exponent β.

Zero temperature
For several disorder realizations we start from a perfect triangular lattice at high velocity and the driving force is slowly decreased down to the sample-dependent critical depinning force F sample c below which the system is permanently pinned. We compute the mean critical depinning force F c where the overline indicates average over disorder realizations. We show in Fig.1a the evolution of F c with respect to the relative disorder strength α p /α v for N v = 8000 vortices in a basic cell of size (L x , L y ) = (400, 20 √ 3/2)λ L . The rapid increase in the mean critical depinning force indicates a crossover from elastic dynamics dominated by elasticity to plastic dynamics dominated by disorder: such a sharp increase was for example found in a previous work on a similar system [37]. In Fig.  1b the typical trajectories of the vortices at the elastic depinning threshold are displayed for α p /α v ≈ 5.10 −3 . The vortices flow in elastically coupled rough static channels and the structure is topologically ordered: all vortices depin together with the same mean velocity, which means that all vortices keep the same neighbors as they move. The dynamics is jerky and the velocity of the center of mass is periodic in time where the period corresponds to the time for each vortex to replace its preceding neighbor in the same channel.

Velocity-force response
For each system size (L x , L y ) we define the reduced velocity, is the time averaged longitudinal velocity of the center of mass of the vortices and F sample c is the critical force, both measured for a given disorder realization, and the overline is an average over all disorder realizations for a fixed value of f . This approach has been used for example in numerical simulations of elastic interfaces in a disordered medium [56].
We set the relative disorder strength α p /α v ≈ 5.10 −3 and we plot in Fig. 2 the reduced velocity v with respect to the reduced force f for two different system sizes: (L x , L y ) = (100, 50 √ 3/2)λ L averaged over N = 21 samples in Fig. 2a, and (L x , L y ) = (400, 20 √ 3/2)λ L averaged over N = 14 samples in Fig. 2b. Three regions appear in Fig. 2. Region I is the manifestation of the finite size effects in the system whose signature is the single-particle regime [46], where v ∼  hysteresis decreases when the system size increases, which therefore confirms that such hysteresis phenomenon is just a finite size effect. In region II a power law regime v ∼ f β with β < 1 is measured, which we identify with the critical regime of the continuous depinning transition. Finally, in region III the system is far above the critical depinning threshold and approaches the asymptotic linear behavior v ∼ f obtained in the absence of disorder. For each system size we compute the depinning exponent β from a power law fit in the region II of the disorder averaged v−f curves (see Fig. 2). Fig. 3a shows the evolution of β with respect to the transverse size L y for a fixed longitudinal size L x = 100λ L . We see that the value of β for L y ≥ 18 √ 3/2λ L becomes independant of the transverse size L y . In particular β does not depend on the basic cell anisotropy since we measure identical values in square basic cells (L x , L y ) = (100, 120 √ 3/2)λ L . In Fig. 3b we show the evolution of β with respect to the longitudinal size L x for various transverse sizes L y ≥ 18 √ 3/2λ L . It can be seen that β has reached a constant value for L x ≥ 100λ L . Taking the mean value of β computed for both L x ≥ 100λ L and L y ≥ 18 √ 3/2λ L , we obtain the result β = 0.29 ± 0.03. Note that taking the mean of individual values β i for each disorder realization, i.e. without averaging the velocity over disorder, gives a similar value β = 0.27 ± 0.04.

Finite-size scaling
Finite-size systems used in numerical simulations lead to rounding and shifting effects in second-order phase transitions that can be analysed within the finite-size scaling theory [57]. In addition to providing a clearer picture of the thermodynamic limit of the transition from finite size observations, this approach allows to extract further information about diverging lengths in the system. From a general theorem for disordered systems [58] a finite-size scaling length can be defined from the statistical properties of a large number of finite-size samples. The finite-size scaling exponent ν FS characterizing the divergence of such a length at the threshold should satisfy the inequality ν FS ≥ 2/d.
We first examine the critical force distribution for each system size, and determine their width ∆F c (L x ) for the relative disorder strength α p /α v ≈ 5.10 −3 . In Fig. 4a we display ∆F c (L x ) versus L x which scales as [59], ∆F c (L x ) ∼ L −1/ν FS x with ν FS = 1.09 ± 0.07, in agreement with the inequality ν FS ≥ 2/d.
Similarly, for each system size we identify the appearence of finite size effects with the crossover force F cross between the critical regime and the single-particle regime (see regions I and II in Fig. 2). We plot in Fig. 4b the width ∆F cross (L x ) of the crossover force distribution, which behaves like ∆F cross (L x ) ∼ L −1/ν FS x with the same exponent ν FS = 1.09 ± 0.02.
Furthermore, we expect the correlation length ξ to behave in the critical regime as ξ ∼ (F − F ∞ c ) −ν where ν is the intrinsic correlation length exponent, and F ∞ c is the critical force for the infinite system. At the crossover force F cross the correlation length becomes of the order of the system size L x , so that the crossover force F cross varies with the longitudinal size L x as   From a power law fit shown in Fig. 5, we extract the value ν = 1.04 ± 0.04 which is consistent with ν FS computed previously. Therefore, a diverging length above the depinning threshold with an exponent ν = ν FS seems to appear. Our results suggest that this length coincides with the intrinsic correlation length ξ. We show in Fig. 6 the direct computation of the correlation length ξ extracted from the velocity correlation function fit with the functional form < v x (0)v x (x) >∼ C 0 x −κ e −x/ξ (see [60]). We find that both fit parameters C 0 and κ are almost constant when F − F ∞ c varies. For a given disorder realization and (L x , L y ) = (100, 18 Fig. 6 shows the evolution of ξ as a function of F − F ∞ c . As expected, when F → F + cross finite size effects actually appear when ξ ∼ L x with L x = 100λ L . Furthermore the correlation length ξ diverges as ξ ∼ (F − F ∞ c ) −ν for F > F cross in the critical regime, and the power law fit exponent gives ν = 1.2 ± 0.2 which nicely agrees with Fig. 5. Therefore, our results strongly suggest that, unlike CDWs [59,61,63], the intrinsic correlation length controls both the critical force distribution and the finite-size crossover force.
We can now write a scaling relation between the velocity, the driving force and the system size L x at T = 0 as, v ∼ L −β/ν x g(f L 1/ν x ) This scaling relation is relevant only close enough to the depinning transition and at large enough L x .
We take the disorder averaged response v(f ) for different system sizes L x = (40, 100, 200, 400)λ L and various transverse sizes from L y = 18 √ 3/2λ L to L y = 120 √ 3/2λ L , and for the relative disorder strength α p /α v ≈ 5.10 −3 . In Fig. 7 we plot vL β/ν x versus f L 1/ν x using the value of β = 0.29 ± 0.03 and ν = 1.04 ± 0.04, retaining only the points of the region I and II. A good collapse of the data onto a single curve is obtained. The behavior of the scaling function g(x) is g(x) ∼ x β for x ≫ 1, and in contrast to ordinary phase transitions g(x) → 0 for x → 0. This unusual property has already been observed in numerical simulations of the depinning transition of driven interfaces [64]. It is due to the single-particle regime which can partially hides usual finite-size effects.

Critical exponent δ
In this section we study the finite temperature effects on the depinning transition. For F < F c at T > 0 the vortex mean velocity does not vanish since thermal fluctuations provide sufficient energy to overcome local energy barriers: this results in a rounded depinning transition. Using the analogy with standard critical phenomena, a power law response v F =Fc ∼ T 1/δ can be defined at F = F c . We show in Fig. 8 the typical reduced velocity-temperature response that we obtain. Two different behaviors can be seen in the v(T ) response. Below the critical force, v goes to zero faster than a power law resulting in concave v(T ) curves, while above the critical force v approches a nonzero limit as T goes to zero resulting in convex v(T ) curves with an horizontal asymptote on the left. The effective critical force F * c is defined as the force at which the convexity changes. In agreement with a second-order phase transition, we can extrapolate at F * c a power law response v F =F * c ∼ T 1/δ (dashed line in Fig. 8) from which we measure the critical exponent δ −1 = 0.28 ± 0.05, where the error bar is deduced from the different lines one can draw to extrapolate the power-law behavior.
We would like to emphasize that this new determination of the disorder averaged critical force is consistent with the value of F c analysed at T = 0 for the same samples in section 3. Going back to the v(f ) curves at T = 0 and using the effective critical force F * c , we find β = 0.28 ± 0.05 which nicely agrees with the value β = 0.29 ± 0.03 extracted from the zero temperature simulations of section 3. This larger uncertainty comes from the smaller number of samples simulated at finite temperatures.
Finally, note that we observe an identical mean value δ −1 = 0.28 when taking the mean of individual values δ −1 i obtained from the velocity-temperature curves for each disorder realization.
Because of the power law responses of the velocity with T and f , we assume that the velocity v is a generalized homogenous function of driving force f and temperature T . We want to express the relation between v, f , and T in terms of dimensionless quantities and to make this relation independent of the prefactors in the two power law responses. We define v 0 and T 0 as and we introduce the dimensionless velocity v = v/v 0 and the dimensionless temperature T = T /T 0 . We can then define the following scaling ansatz, where S(x) is the scaling function with two branches S + and S − corresponding to f > 0 and f < 0, respectively. The power law dependencies v f >0,T =0 = f β and v f =0 = T 1/δ imply that S(x) has the following asymptotical behavior: The S + branch goes asymptotically to the horizontal axis for T → 0 defining a driving dominated regime, while the S ± branches have both an oblique asymptote with slope δ −1 for f → 0 defining a temperature dominated regime. The change of variables (v, T ) → ( v, T ) is equivalent to choosing the intersection of the asymptotes as the origin of coordinates.
To show the existence of this scaling behavior, we compute the velocity without averaging over the disorder, i.e. we compute v i for each of the N = 14 samples that we have studied, with i = 1, .., N. This choice allows to better see if all data scale on a single curve. Different system sizes L x = (100, 400, 1000)λ L i.e. N v = (5000, 8000, 20000), different pinning strengths 10 3 × α p /α v ≈ (2, 3, 5, 7) and different realizations of the random position of the pinning centers have been investigated. In Fig. 10 we plot v i |f | −β with respect to T |f | −βδ and a good collapse of all individual samples to a single curve with two branches is obtained. The values of F c used to compute f are sample specific. The two prefactors v 0 and T 0 are not universal and their value varies for each sample, but the scaling functions S ± remain unchanged. We verified that averaging over disorder as previously defined, we obtain the same scaling functions.
The collapse of our data implies that a true critical regime is observed with power law responses for the velocity with respect to the driving force and to the temperature. Moreover the critical depinning exponent β and the thermal exponent δ −1 do not depend on the strength and the positions of the pins. Theses results suggest some degree of universality of our model for the elastic depinning of vortices.

Discussion
In our study we find a continuous elastic depinning transition for weak pinning strength. Both the depinning exponent β and the thermal exponent δ, and the scaling function linking velocity v, temperature T and driving force f show some degree of universality with respect to disorder. A large variety of β values has been found in the literature as cited in section 1. Our value is close to the value β ≈ 1/3 measured for interfaces [60] in a space dimension d = 2. However, our result should rather be compared to other periodic systems. For example, our result is close to the value β = 0.35 recently measured in anisotropic periodic systems (stripe systems) [43] for which the displacement field dimension is N = 2 in a d = 2 dimensional space. Much studies appear in the literature for the case N = 1, in particular for CDWs where numerical simulations [39] in d = 2 give β = 0.65 ± 0.05 and ν = 0.5 ± 0.1, whereas our results give very different values β = 0.29 ± 0.03 and ν = 1.04 ± 0.04. Also, we find ν ≈ 1 for the intrinsic correlation legnth exponent while ν = 1/2 is expected for CDWs above threshold. Moreover, our results of section 3.2 suggest ν = ν F S in contradiction with CDWs. Actually, the poor agreement between our results and CDWs may originate from the fact that N = 1 for CDWs whereas N = 2 in our vortex simulation. We note that inserting ν = 1 in the scaling relation ν = 1/(2 − ζ) which originates from the statistical tilt symmetry [20] gives ζ = 1. When inserting both values in the hyperscaling relation β = ν(z − ζ) [19,20,63] (where z and ζ are respectively the dynamic and roughness critical exponents) together with z = 5/4 [65], we obtain β = 1/4 which is close to our result β = 0.29 ± 0.03 but may not be satisfactory since ζ = 0 is expected to hold in the hyperscaling relation.
Finally, going back to the connection between the depinning transition and absorbing-state phase transitions, we note that our result β = 0.29 ± 0.03 is close to the value found in d = 1 directed percolation or conserved directed percolation, and that ν = 1.04 ± 0.04 is close to the value of ν ⊥ in d = 1 directed percolation.

Conclusion
In this paper we have numerically studied the 2D superconductor vortex dynamics in random media. A crossover from plastic dynamics dominated by disorder to elastic dynamics dominated by elasticity is found. We investigated the depinning transition in the elastic regime. Above the depinning threshold all the vortices depin together and flow in elastically coupled rough static channels. Three dynamical regimes are found: the first one is controlled by the finite-size of the simulation box leading to the so-called single-particle regime, the second one is the critical regime from which a power law response v ∼ f β at T = 0 can be extracted with β = 0.29 ± 0.03, while the third one is the recovery of the ohmic regime. A finite temperature analysis gives the thermal rounding exponent δ −1 = 0.28 ± 0.05 defined by the power law response v ∼ T 1/δ at the threshold force. A finite-size scaling analysis at zero temperature suggests (in contrast with CDWs) the existence of a unique diverging length above the depinning threshold with an exponent ν = 1.04±0.04, which controls the critical force distribution, the finitesize crossover force distribution and the intrinsic correlation length. A scaling law for the temperature and the force dependence of the velocity is found, consistent with a second order phase transition. Both critical exponents and the scaling function are independent of the disorder (strength and position of the pins) in the range of parameters we studied, indicating some degree of universality. Nevertheless the comparison with other similar periodic systems may suggest that large universality can not be found for 2D elastic depinning. We expect that such elastic depinning transition occurs in relatively weak pinning regimes of superconductor vortices. Within transport measurements where the voltage response to an applied current is equivalent to the velocity response to an applied force, the critical exponents β and δ may be accessible from finite temperature measurements at different driving currents.