Nonsteady relaxation and critical exponents at the depinning transition

We study the nonsteady relaxation of a driven one-dimensional elastic interface at the depinning transition by extensive numerical simulations concurrently implemented on graphics processing units. We compute the time-dependent velocity and roughness as the interface relaxes from a ﬂat initial conﬁguration at the thermodynamic random-manifold critical force. Above a ﬁrst, nonuniversal microscopic time regime, we ﬁnd a nontrivial long crossover towards the nonsteady macroscopic critical regime. This “mesoscopic” time regime is robust under changes of the microscopic disorder, including its random-bond or random-ﬁeld character, and can be fairly described as power-law corrections to the asymptotic scaling forms, yielding the true critical exponents. In order to avoid ﬁtting effective exponents with a systematic bias we implement a practical criterion of consistency and perform large-scale ( L (cid:2) 2 25 ) simulations for the nonsteady dynamics of the continuum displacement quenched Edwards-Wilkinson equation, getting accurate and consistent depinning exponents for this class: β = 0 . 245 ± 0 . 006, z = 1 . 433 ± 0 . 007, ζ = 1 . 250 ± 0 . 005, and ν = 1 . 333 ± 0 . 007. Our study may explain numerical discrepancies (as large as 30% for the velocity exponent β ) found in the literature. It might also be relevant for the analysis of experimental protocols with driven interfaces keeping a long-term memory of the initial condition.


I. INTRODUCTION
The depinning transition of an elastic interface driven in a random medium is a paradigmatic example of universal out-of-equilibrium dynamical behavior in disordered systems [1,2].From a purely theoretically point of view, considerable progress has been made in the last years thanks to a fruitful interplay between the very powerful analytical [3][4][5] and numerical techniques specially developed for studying equilibrium [6,7], depinning [8][9][10][11][12], and creep [13,14] of strictly elastic manifolds.From the experimental viewpoint, on the other hand, the understanding of this particular problem is directly relevant for various experimental situations where the elastic approximation is well met, such as magnetic [15][16][17][18] or ferroelectric domain walls [19][20][21][22], contact lines in wetting [23,24], and fractures [25,26].It has also been useful for making a connection between laboratory friction experiments and the observed spatial and temporal earthquake clustering [27].In spite of this progress, classifying the above systems in depinning universality classes remains, in most cases, a challenge.In this respect, we show in this paper that a detailed understanding of the experimentally relevant non-stationary dynamics at the depinning transition can be useful.
The physics of an elastic interface in a random medium is controlled by the competition between quenched disorder (induced by the presence of impurities in the host material), which promotes the wandering of the elastic object, and the elastic forces, which tend to make the object flat.One of the most dramatic manifestations of this competition is the response of these systems to an external drive, particularly the depinning transition phenomenon.In the absence of an external drive, the ground state of the system is disordered but well charac-terized by a self-affine rough geometry with a diverging typical width w ∼ L ζeq , where L is the linear size of the elastic object and ζ eq is the equilibrium roughness exponent.When the external force is increased from zero, the ground state becomes unstable and the interface is locked in metastable states.To overcome the barriers separating them and reach a finite steady-state velocity v it is necessary to exceed a finite critical force, above which barriers disappear and no metastable states exist.For directed d-dimensional elastic interfaces with convex elastic energies in a (D = d + 1)-dimensional finite space with disorder, the critical point is unique, characterized by the critical force f = f c and its associated critical configuration [28,29].This critical configuration is also rough and self-affine such that w ∼ L ζ with ζ the depinning roughness exponent.When approaching the threshold from above, the steady-state average velocity vanishes as v ∼ (f − f c ) β and the correlation length characterizing the cooperative avalanche-like motion diverges as ξ ∼ (f −f c ) −ν for f > f c with a typical diverging interevent time ξ z , where β is the velocity exponent, ν is the depinning correlation length exponent, and z is the dynamical exponent [5,[30][31][32].At finite temperature and for f ≪ f c , the system presents an ultra-slow steadystate creep motion with universal features [4,33,34] directly correlated with geometrical crossovers [13,35].At very small temperatures the monotonic increase of the correlation length with decreasing f below f c shows that the naive analogy breaks, and that depinning must be regarded as a non-standard phase transition [13,14].The transition is then smeared-out, with the velocity vanishing as v ∼ T ψ exactly at f = f c , with ψ, the so-called thermal rounding exponent [36][37][38][39][40].
The non-stationary dynamics at depinning is a different and interesting manifestation of the competition between elasticity and disorder, but it has received con-siderably less attention than the steady-state dynamics.Near the threshold the time needed to reach the driven non-equilibrium steady state can be very long, since the memory of the initial condition persists for length scales larger than a growing correlation length ℓ(t) ∼ t 1/z .Being limited only by the divergent steady-state correlation length ξ and the system size L, the resulting non-steady critical regime is macroscopically large, t ξ z , L z .It is hence relevant for experimental protocols, where it is in general difficult to assure history independence, as the memory of the initial condition is erased only by this slow transient process.Analogously to non-driven systems relaxing to their critical equilibrium states [41], the transient dynamics of a driven disordered system displays interesting, though different, universal features [42].
In this paper we study the non-steady relaxation of an elastic string in a random medium by extensive numerical simulations, extending the study of Ref. [43] in several ways.We first show that the thermodynamic critical force, although strictly non-universal, can be defined unambiguously and has a unique value; this is the force that drives the system into the so-called short-time-dynamics (STD) scaling form yielding the steady-state critical exponents.Second, we study corrections to the dynamical scaling form which appear at intermediate ("mesoscopic") times, but well above the expected non-universal microscopic time regime set by the microscopic disorder correlation range.We show that these corrections are rather robust under changes of the microscopic disorder, including its random-bond or random-field character.Having done this, we implement a practical criterion for separating these corrections and exploit the parallelism of graphics processing unit (GPU) computing to perform, to the best of our knowledge, the largest-todate molecular dynamics simulations of an elastic line in a random-medium described by a continuous quenched Edwards-Wilkinson (QEW) equation.This allows us to achieve long times in the non-steady critical regime and get, with the STD method, un-biased and precise estimates for the critical exponents for this depinning universality class: β = 0.245 ± 0.006, z = 1.433 ± 0.007, ζ = 1.250 ± 0.005, and ν = 1.333 ± 0.007.

Organization of the paper
In Sec.II we define the model and quantities of interest.In Sec.III A we describe the general scaling forms for the non-steady critical relaxation that are used for the analysis.In Sec.III B we briefly discuss the numerical methods, leaving all computational details for the Supplemental Material [44].In Sec.IV we show all the results.We first analyze the geometry of the critical configuration and the thermodynamic limit of the critical force in Sec.IV A. Then we study the non-steady relaxation of different observables in Sec.IV B. We characterize in Sec.IV C the deviations from the expected scaling forms of Sec.III A, and in Sec.IV D we analyze in detail these corrections to scaling, which explain the existence of biased effective exponents in the problem.In Sec.V we summarize and discuss all the results and in Sec.VI we conclude and give perspectives.

II. MODELS, PROTOCOLS, AND OBSERVABLES
We study a one-dimensional elastic line described by a single-valued function u(x, t).Here, u(x, t) measures the line's instantaneous transverse displacement from the x axis at time t.In the continuum, the overdamped interface dynamics is described by the QEW equation of motion, where f is a uniformly applied external force and F p (u, x) is the random pinning force.Without lack of generality we can set the friction coefficient η = 1 and the elastic constant c = 1.The pinning force has zero average and correlator The overbar represents the average over the disorder realizations and ∆(u) is a short-ranged function, of range r f = 1.We consider two cases for ∆(x).In the socalled random bond (RB) case the elastic line moves in a random potential such that F p (u, x) = −∂ u U (u, x) with U (u, x) bounded, and thus u ∆(u) = 0.In the random field (RF) case, the random potential U (u, x) is unbounded and diffuses as a function of u, with diffusion constant u ∆(u) > 0 [4].In turn, the pinning potential or forces can be sampled from different distributions.
Here we consider Gaussian and uniform (constant) distributions.Analyzing all the above cases separately will allow us to detect possible departures from the expected universal behavior and corrections to scaling at relatively short time and length scales.It is a convenient and safe procedure to discretize the interface displacement in the x-direction, keeping u(x, t) as a continuum variable.Doing so, we define the center of mass velocity for an interface of size L as, which, given Eq. ( 1) and η = c = 1, is nothing else but the spatial average of the instantaneous total forces acting on the line.From the instantaneous center of mass displacement, we can define the instantaneous quadratic width of the interface as, The geometrical properties of the line as a function of length-scale can be conveniently described using the averaged structure factor where q = 2πn/L, with n = 1, . . ., L − 1.We are interested in the time dependence of the above disorder-averaged observables with the elastic line initially flat, i.e., u(x, t = 0) = 0. Other, more complex, initial conditions may also be considered, but they do not provide more information on the dynamics of relaxation as long as they are not correlated with the random potential.Choosing an initially flat line allows us to easily detect the existence of memory depending on the length-scale, since correlations in this system always develop a roughness with self-affine length regimes described by positive exponents.For a globally self affine line we have S(q) ∼ q −(1+2ζ) , thus yielding the roughness exponent ζ.Such behavior is expected to hold in the socalled random-manifold regime; a regime where q is small compared to both the Larkin wavevector q c ∼ ℓ −1 c , [with w(ℓ c ) ∼ r f ] and the discretization wave vector q d ∼ 1.For our parameters, we have q c ∼ q d of order O(1).

III. METHODS
We describe here the scaling forms expected for the non-stationary relaxation of an elastic line in a random medium at large times, and introduce the numerical methods used to simulate the dynamics.

A. Universal non-steady relaxation
It was originally shown [45] using renormalization group techniques for a type A model, and then widely extended empirically to other models, that the short-time behavior of the order parameter's n-th moment close to a critical point follows the universal scaling relation Here t is the time, h = (H − H c )/H c is the reduced driving field of the transition associated with the order parameter (such as magnetic field, temperature, etc.), L is the system size, b is a scaling factor, m 0 is the order parameter's initial value, µ is a universal exponent associated with the short-time behavior, while β, ν, and z are the usual critical and dynamical exponents.This relation was analytically shown to be valid in the limit of small initial order parameter m 0 ≪ 1, with t larger that some non-universal microscopic time and smaller than the equilibration time τ eq ∼ L z .Nevertheless, in agreement with numerical simulations [46][47][48] and analytical approximations in mean-field models [49], it was shown that the following homogeneity relation is valid in the short (but macroscopic) time regime when starting from an ordered condition (commonly m 0 = 1).
In our system, and making the usual analogy with standard phase transitions, we consider the velocity v as an order parameter and the adimensionalized force (f − f c )/f c as the reduced driving field.While it is not clear yet how to implement a protocol yielding an initial condition equivalent to m 0 ≪ 1 and testing the complete relation Eq. (7), the implementation of an ordered equivalent initial condition is very simple, and a relation like Eq. ( 8) has been numerically checked in this model [42,43].In fact, the completely ordered initial condition should correspond to the infinite-velocity configuration.Since at large velocities the effect of the quenched disorder mimics a small thermal noise which vanishes with increasing v, this initial condition should correspond to the completely flat condition (i.e., we make a "quench" from an infinite to a finite force).By analogy, we therefore expect the velocity short-time behavior where h = |f − f c |/f c and the function ṽ± has two branches depending on the sign of f − f c .By choosing b = t 1/z in (9), Here, the functions ṽ+ and ṽ− are such that for h > 0 and L ≫ h −ν , v + → const and v − → 0 in the large t limit.It is worth noting here that while h −ν can be associated with the geometrical correlation length ξ diverging in lim f → f + c in the steady state, this is not true in lim f → f − c .We can, however, find a divergent correlation length ℓ relax ∼ (f c − f ) −ν , not observed in the steady-state geometry, but associated with the deterministic part of the avalanches that are produced in the steady-state dynamics of the f < f c lowtemperature creep regime [13,14].Hence, the interpretation of Eq. ( 10) is slightly different from the one derived from Eq. ( 8) for the Ising model (and other similar standard phase transitions), where divergent equilibrium correlation lengths do exist on approaching the critical point from both sides.At depinning, the relaxational dynamics described by Eq. ( 10) is valid in the short-time regime, after which the velocity reaches a steady-state value if f > f c , while for f < f c the velocity is blocked in a metastable state with memory of the initial condition and v → 0 [43].When T > 0, and for f < f c , activated processes permit the system to continue relaxing towards a steady-state with finite velocity [39].
Exactly at the critical point (h = 0) and in the limit L → ∞ we expect a power-law behavior for the velocity, Equation (11) can also be heuristically derived from a very simple physical picture, by assuming that the nonsteady relaxation at the critical point is controlled by a single growing correlation length ℓ(t) ∼ t 1/z , and that in the macroscopic time limit, some critical steady-state relations are "instantaneously" obeyed, with ℓ(t) playing the role of ξ, when ξ ≪ L. Hence, the steady-state relation v ∼ ξ −β/ν translates into v ∼ ℓ(t) −β/ν ∼ t −β/νz , leading to Eq. (11).The physical picture behind this "mapping" is that length scales below ℓ(t) are expected to be steady-state quasi-equilibrated, while above it the memory of the initical condition is still kept.This suggests that the structure factor should be described, above some microscopic scale and for ℓ(t) ≪ L, by S q (t) ∼ q −(1+2ζ) S(qℓ(t)) ∼ q −(1+2ζ) S(qt 1/z ), (12) with S(y) ∼ const for y ≫ 1, and S(y) ∼ y 1+2ζ for y ≪ 1.This is confirmed by numerical simulations [42,43].Note that this simple scaling is also obeyed in the absence of disorder but in the presence of thermal noise [50].For the depinning transition S(y) is expected to be a much more complicated function however, as the disorder in principle couples all Fourier modes.The above scaling relation also implies for the global width, or roughness.Note that this simple result depends on the choice of a flat initial condition.
Otherwise we would have a memory contribution, which can be included as w(t) ∼ t ζ/z w(t 1/z L −1 ), with w(y) ∼ const for y ≤ 1, and a function depending on the initial condition for y ≪ 1.The flat initial condition is thus convenient since w(y) = const for all y.Since Eq. ( 11) is valid exactly at h = 0, by evolving the system at different driving forces from a flat initial condition it is then possible to determine the critical force and extract the steady-state critical exponents.For the correct application of this STD method, we need, however, a good criterion to decide the time range in which to fit the exponents accurately.This is particularly tricky if scaling corrections are important and also described by power laws, since they can lead to effective power-law decays, yielding biased incorrect exponents.As we show later on, for our case, we can exploit known depinning scaling relations and combinations of observables to find a good criterion.

B. Numerical simulation
The numerical simulation protocol is roughly the same as the one used in Ref. [43], but here, in addition, we simulate different kinds of disorder, and implement massively parallel algorithms.Below we briefly summarize the method.We refer the reader to the Supplemental Material [44] for further details; in particular, for the description of our GPU-based parallel implementations of the algorithm which allow us to computationally reach very large system sizes, running simulations with speedups above 300x.
In order to numerically solve Eq. (1) the system is discretized in the x direction into L segments of size δx = 1, i.e., x = 0, . . ., L − 1, while keeping u(x, t) as a continuous variable.Computing the forces at each time-step, the integration of Eq. ( 1) is done using the Euler method.
To model the continuous quenched random potential we can either read from a precomputed array, or generate dynamically, uncorrelated random numbers with a finite variance at the integer values of u and x and use interpolation to get F p (u, x).In this work we consider both cubic-spline and linear interpolation between the random potential values sampled at the integers.This changes only the shape of the microscopic force correlator, but not the universality class.By generating directly the force field at integer positions from zero-mean uncorrelated random numbers in both directions we can get a RF disorder, while by doing the same for the potential field and deriving the force we can get a RB disorder.In both cases Eq. ( 2), with a short-range force correlator, is automatically satisfied.
For our numerical simulations we have used periodic boundary conditions in the longitudinal direction, so that u(0, t) is elastically coupled with u(L − 1, t).In the case where we construct continuous splines for the disorder potential or we limit ourselves to read a precomputed force field, we enforce periodic boundary conditions also in the transverse direction, thus defining an L × M system.On the other hand, when a dynamically generated disorder is used, the system size in the u direction is virtually infinite, although in our implementation it can be forced to be periodic as well [44].
A critical configuration u s c (x) and a critical force f s c can be unambiguously defined for a periodic sample of size L × M and with a given disorder realization.They are defined from the pinned (zero-velocity) configuration with the largest driving force f in the long time limit dynamics.They are the real solutions of such that for f > f s c there are no further real solutions (pinned configurations).Middleton theorems [28] can be used to devise an efficient algorithm which allows the critical force f s c and the critical configuration u s c (x) to be obtained for each independent disorder realization iteratively, without having to solve the actual dynamics Eq. ( 1) nor directly invert the non-linear system of Eq. ( 14).In the following section we use such an algorithm [10,11] to study the finite-size effects on f s c and in particular to obtain the appropriate thermodynamic limit f c = lim (L,M)→∞ f s c controlling the universal nonsteady relaxation.

A. The critical configuration and the thermodynamic critical force
We start by analyzing the geometry of the critical configuration and the critical force in periodic samples of size L × M , using the exact algorithm of Refs.[10,11], for a cubic-spline RB potential.This allows us, on one hand, to get a very precise estimate for the exponent ζ directly from steady-state solutions for individual samples u s c (x) and use it to disentangle the combinations of exponents appearing in the non-steady universal part of the relaxation.On the other hand, it allows us to get the value f c for the thermodynamic critical force, and show, from an anisotropic finite-size analysis, that it is unambiguously defined.In the next section we show that the results obtained in this section are fully consistent with the non-steady dynamics: The long-time geometry tends to the self-affine one of the critical configuration and, in particular, f c is exactly the force for which the velocity relaxation asymptotically displays a power-law decay in time, as long as the growing correlation length remains smaller than the system size L We determine the critical configurations u s c (x) and the associate critical force f s c (L, M ), for several samples of size L × M .To determine ζ we calculate the timeindependent averaged structure factor S q of the critical configurations and fit the expected S q ∼ q −(1+2ζ)  behavior.The advantage of using S q over the global width w(L) ∼ L ζ is that it allows us to better detect the possible presence of undesirable crossover effects, and the characteristic associated length-scale.In particular, when determining ζ, one should be cautious about the crossover to the random-periodic regime, with exponent ζ RP = 1.5 occurring at the length-scale In our case we have thus chosen M L ζ for extracting ζ.In Fig. 1 we observe that S q ∼ q −(1+2ζ) holds to high accuracy in the long wavelength limit (discretization effects appear at very short wavelengths) for a system of size L = 2048 and M = 13777.From a power-law fitting giving 1 + 2ζ = 3.50 ± 0.01 we deduce ζ = 1.250 ± 0.005, with better precision but still in very good agreement with previous estimates obtained with the structure factor at finite velocities in the steady-state [52] and with the scaling of the global width w(L) ∼ L ζ [8,29].
In order to study the thermodynamic limit of the disorder averaged critical force, f s c (L, M ), we use M = kL ζ , with k a finite control parameter.Any finite value of k leads to a parametrized family of universal critical force distributions, ranging from the Gaussian to the Gumbel distribution [53].For our present purposes, it is important to know how the thermodynamic limit L → ∞ depends on the aspect-ratio parameter k.In Fig. 2 we show f s c against L, for different values of k.As we can see, the size-dependent average critical force f s c tends to a unique value in the limit L → ∞ when keeping k fixed.Furthermore, at finite L, f s c is smaller than f c for k 2, and larger than f c for k 2, where f c therefore represents the thermodynamic critical force.This behavior is consistent with the crossover from a random-periodic depinning (k → 0), with a thermodynamic critical force smaller than the random-manifold one, towards the extremestatistics random-manifold depinning with an infinite thermodynamical critical force (k → ∞), with fluctuations described by the Gumbel distribution [53,54].What is important to remark here [55] is that for the whole range of finite values of k that we have analyzed, the curves slowly converge to the same thermodynamic limit 001.This value is in very good agreement with the value obtained for the same microscopic model with other methods, such as quasi-statically pulling the string with a soft spring [56].
Therefore, we see that although f c is not universal, for fixed microscopic disorder correlations it attracts the infinite family of sample geometries that is compatible with the random-manifold setting M = kL ζ .Therefore, this setting alone removes the ambiguity of defining the critical force of a periodic sample in the thermodynamic limit.Since the short-time relaxation of the string is controlled by a single growing correlation length ℓ(t) ∼ t 1/z describing the slow development of correlations from microscopic scales, it can not be sensitive to the dimensions or aspect ratio of the computational box, provided that ℓ(t) ≪ L and ℓ(t) ζ ≪ M .It is then natural to expect that f c , which is sensitive only to the microscopic pinning correlator for the random-bond family M = kL ζ , is the force that must be tuned to target the universal non-steady relaxation regime.

B. Short-time-dynamics scaling
We now turn to the study of the relaxation dynamics.We start by reproducing previous results [43], with exactly the same model, but using larger systems, in order to make visible the effects that were neglected before.
In Fig. 3 we observe the averaged string velocity behavior against time.When the applied force f is greater than the thermodynamic critical force f c , v(t) saturates to a finite value at a time which increases as we approach f c from above.On the other hand, for forces smaller than f c , v(t) goes to 0 after a transient, which is longer as we approach f c from below.Exactly at f = f c we expect v(t) ∼ t −β/νz after a microscopic time regime.One can in principle use this as a criterion to determine f c .Of course, it is very difficult to hit exactly f c , but what we can do is to bound it from above and below.As we approach f c it takes longer and longer times (and requires better averages) to determine if a velocity curve for a given force is saturating or vanishing.Big system sizes help to reduce noise, since v(t) self-averages.A detailed inspection of the simulations results presented in Fig. 3(a) permits us to determine f STD c = 1.915 ± 0.002 using the STD method, in complete agreement with the extrapolation to the thermodynamic limit of the exact critical force obtained in Sec.IV A. This shows the consistency of the two methods, and from now on we refer to f STD c simply as f c .With the value of f c , from Eq. ( 9) one can test a joint scaling form for the force and time-dependent velocity.Considering L → ∞ and using b ∼ t 1/z in Eq. ( 9) we get   where ṽ± are universal functions, and the ± sign indicates whether the critical point at f c is approached from above or from below.This scaling relation is tested in Fig. 3(b).In order to appreciate the difference from previous studies of the same problem the following values for the exponents have been used: z = 1.5 [43], β = 0.33 [52], and ν = 1.33 [43].It can be observed that the curves collapse into two sets: on one hand those with f < f c , with ṽ− (x) going to zero, and on the other hand those with f > f c with a saturation to a finite value of ṽ+ (t).
In Fig. 4 we show the string averaged structure factor at different times during the evolution of the system in the presence of the applied force f = f c = 1.915.As can be seen in Fig. 4(a), at large length scales (small values of q) the system keeps memory of its initial state, which we have chosen to be flat.At small length scales though (but not small enough to explore the lattice effects), the geometry of the string becomes approximately self affine, being characterized by a behavior S q ∼ q −(1+2ζ) where ζ is the roughness exponent.In Fig. 4(b) we test the scaling hypothesis of Eq. ( 12), by using z = 1.5 [43] and ζ = 1.25 [29,52].
As we can observe in Figs.critical exponents β = 0.33, ν = 1.33, ζ = 1.25, and z = 1.5 produce rather good collapses.We will argue however that these exponents are effective and have an appreciable bias, which causes the small deviations appreciated in those collapses.This can be perceived only by simulating large systems and reaching several orders of magnitude in time.In Fig. 5 we show v(t) and w(t) vs t at the critical force f c .Above the microscopic time regime of order t mic ∼ 1 we observe in v(t) a crossover at a time t cross ∼ 50, between two approximate powerlaw decays with exponents ≈ 0.16 for short but mesoscopic times, and ≈ 0.13 for times larger than a few tens [see Fig. 5(a)].A similar crossover is seen in Fig. 5(b) where we observe a crossover between two approximate power-law decays, from ≈ 0.67 to ≈ 0.84.This behavior is unexpected from the universal non-steady relations v(t) ∼ t −β/νz and w(t) ∼ t ζ/z , which generally assume a microscopically very short transient regime before reaching the truly universal macroscopic time regime.In order to understand this behavior, in the next section we analyze this crossover in detail.

C. Robustness of the crossovers
The first question we can ask is how robust are the observed crossovers and in particular, how much they depend on specific details such as the intensity or the shape of the microscopic disorder correlator.
For a given disorder distribution, we have first checked, by simulating v(t) at the corresponding critical force, that the crossover time (t cross ∼ 50) does not change appreciably by changing the nature of the disorder.
In order to study the dependence of the crossover on the shape and nature of the microscopic disorder correlator, we have fixed ∆(0), defined in Eq. ( 2), and analyzed four different cases: • Case 1: Gaussian-distributed disorder potential with cubic spline interpolation between the integers (RB), • Case 2: Uniformly-distributed disorder potential with cubic spline interpolation between the integers (RB), • Case 3: Gaussian-distributed disorder potential with linear spline interpolation between the integers (RB), • Case 4: Gaussian-distributed random forces (RF) without interpolation between the integers.
Case 1 corresponds to the usual Gaussian-distributed disorder with cubic spline that we used in Sec.IV B. Case 2 consists in using a uniformly distributed disorder instead of the Gaussian-distributed numbers of the first case.In Case 3 we determine the continous potential using linear interpolation instead of a cubic one, so the pinning force as a function of the position has discontinuities at the integers.These first three cases correspond to random-bond disorder, as the quenched potential fluctuations Case 4 helps us in particular to answer the question of whether the observed crossovers may be related to the crossover from RF to RB expected for depinning [5].
In Fig. 6 we show the time dependent velocity for the four cases, each one at its corresponding thermodynamic critical force, determined with the STD method.Apart from a factor depending on f c , the qualitative behavior of v(t) appears unaltered for times longer than some microscopic threshold t mic ∼ 1.This time is controlled by the microscopic correlator range [43], which is the same in the four cases, r f = 1.If we fit pure power laws for all curves in the region t t cross ≈ 50 we find, within the error bars, consistent values for the exponent β/νz ∼ 0.129, as shown in Table I.By plotting vt 0.129 vs t, we clearly see that corrections to scaling at intermediate times (t mic < t < t cross ) are present.They have an approximate algebraic decay, and seem to have a weak dependence on the microscopic disorder.We thus conclude that the observed value for β/νz is universal as expected.Moreover, the corrections to scaling appear to be robust: they do not strongly depend on the intensity, shape, or nature (RB vs RF) of the short-ranged correlated microscopic pinning force.
The crossovers in Fig. 6 can be thus possibly ascribed to leading order power-law scaling corrections to the universal non-steady dynamics, similar to what is observed in the Ising, Potts, and XY models for instance [57][58][59].
As is customary, one can try here to fit a "corrected" formula in order to get an unbiased value for β/νz: However, with the present data, we find that this procedure is inaccurate as it involves several parameters introduced through an ad hoc formula for the correction.This takes us to our next step, the search for a different practical criterion to directly estimate t cross .We start by observing that in the truly universal non-steady macroscopic critical regime, the quantity should reach a time-independent constant value as a consequence of the scaling relation β = ν(z − ζ): since at f = f c in the macroscopic regime, v(t) ∼ t −β/νz and w(t) ∼ t ζ/z , then γ(t) ∼ const.This criterion was used to determine the depinning threshold in the long-range discrete elastic model [12], since any departure from f c would make γ vanish, or increase with time at very long times.For our present purposes, we note that if the measured γ(t) depends on time, it means that we have not reached the critical scaling regime yet.Therefore, a constant γ(t) is, at least, a necessary condition.In Fig. 7 we show the behavior of γ(t) with time for Case 1.After a microscopic regime we can see a decreasing behavior of γ(t) with time, which implies an "effective" unbalanced relation [ν(z − ζ) − β] eff < 0. The system is then in what we call a "mesoscopic" regime where the critical dynamics scaling [Eq.( 8)] is still not valid.Nevertheless, it is possible to fit a reasonable power law for v(t) in this regime, if larger times are not available.We hence see the importance of having a good criterion to decide where to fit the exponents: if the time is not large enough, the resulting exponents will be effective and have a systematic bias.Moreover, using γ(t) as a method to determine f c may also suffer from these mesoscopic time effects, as its initial decay may be incorrectly ascribed to having an applied force larger than f c .Af- ter that long crossover, we can observe that γ(t) develops a plateau starting at t ∼ 1000.This plateau is cutoff by finite-size effects, which introduce an L-dependent characteristic time.We also observe that scaling corrections in v(t) or w(t), shown in Figs. 5 and 6, are important in the mesoscopic time regime, clearly visible in the time-dependence of γ(t).We thus argue that within the plateau in γ(t), scaling corrections are already negligible, and the true critical exponents can be consistently extracted with the STD method.As we see in Fig. 7, this criterion is strongly limited by finite size effects, reducing the time range of the plateau.We observe that in order to obtain a sufficiently wide time window to fit the exponents accurately (i.e., a few decades in time), we really need to attain big system sizes, much larger than L = 8192.

D. Large-scale simulations
With the numerical method implemented so far, in periodic samples of size L×M , we are limited [60] to system sizes not much bigger than L = 8192.This is because we need sufficiently large M in order to avoid undesirable periodicity effects, i.e., a crossover to the random-periodic class [61].We hence move to our "memory-free" numerical implementation, where we dynamically generate the underlying disorder instead of reading it from a previously stored array, similarly to what was done in Ref. [36] but with concurrent computations.As we explain in the Supplemental Material [44], in this implementation it is simpler to work with linear splines for interpolating the potential, therefore corresponding to Cases 3 (RB) and 4 (RF) of Sec.IV C, for which we have already shown that scaling corrections are present.
In Fig. 8 we present the behavior of γ(t) with time for a string of size L = 2 25 = 33554432, for Case 3. Simulating such large systems give us the benefit of the self-averaging of the velocity.Indeed, in Fig. 8 the smooth curve for f c was averaged only over five samples, while the other two correspond to only one sample.Second, notice that, compared with Table I, we are now able to improve the determination of the thermodynamical critical force with the STD method, pushing the uncertainty to the fourth decimal place, f c = 1.5652 ± 0.0003 for Case 3. At this value of f c we obtain a well defined plateau in γ, lasting approximately two decades, starting at t ∼ 1000.Ignoring the microscopic time regime t < t mic ∼ 1, we can now perform a good fit of the data at f = f c with the expression Here we have introduced the fitting parameters G 0 , t g and α g , obtaining G 0 ≃ 1.143, t g ≃ 0.142, and α g ≃ 0.56.This allows us to quantify the development of the plateau, and to confirm the power-law shape of the scaling correction in γ(t).We are now prepared to consistently fit both the scaling corrections in w(t) and v(t) and the combination of critical exponents ζ/z and β/νz respectively.In Figs. 9 and 10 we show w(t) and the velocity v(t) as functions of time, for L = 2 25 = 33554432 and for an applied force f c = 1.5652, corresponding to a RB linearly interpolated random-potential generated with Gaussian numbers.In Fig. 10, we have also included the behavior of the string for two other forces, just above (f = 1.5655) and just below (f = 1.5649) the estimated depinning force f c .For w(t) the corresponding three curves are indistinguishable, so we show only the curve corresponding to f = f c .We can see very nice asymptotic power-law behaviors of both quantities, spanning several orders of magnitude.Nevertheless, as noticed for γ(t) in Fig. 8, a proper power-law fit can be accomplished only after deciding a starting point for the scaling regime.With  the criterion of a plateau in γ(t) such a starting point could be decided up to some arbitrariness, but it can be included in the final uncertainty of the fitted exponents.In addition, we propose corrected expressions for w and v that can be used to fit the data during all the non-steady evolution, excluding the microscopic regime 0 ≤ t < t mic ∼ 1.The corrections to scaling in these quantities then read The fits give W 0 ≃ 0.839, t w ≃ 0.485, and α w ≃ 0.712 for the width, and V 0 ≃ 0.721, t v ≃ 0.110, and α v ≃ 0.715 for the velocity.The insets of Figs. 9 and 10 confirm that the corrections can be fairly described by powerlaws, as observed earlier in Sec.IV C for smaller systems.
The fitted parameters are to some extent sensitive to the choice of t mic but this arbitrariness can be included in the error bars, and the proposed form of the corrected critical behavior is robust.
For completeness, we also present the corrections to scaling for the RF case.In Fig. 11 we can see the behavior of the string velocity with time when we apply the critical force for the RF case f c = 0.9713 ± 0.0002 and different forces close to f c .The complete picture of the RB case is repeated, but now the expression Eq. ( 20) gives back slightly different exponents for the corrections to scaling.The fit gives V 0 ≃ 0.53, t v ≃ 0.382 and α v ≃ 0.927 in this case.Nevertheless we cannot confirm whether or not the correction exponents are distinguishable between the RB and RF cases, since they are very sensitive to the choice of the arbitrary quantity t mic , the time from which we start the fitting.For example for the RF case α v changes between 0.69 and 0.99 when t mic is varied from 2 to 10.
However, for the RF case we find a β/νz which is indistinguishable from the one extracted for the RB case, as detailed in the following section.With the precision we have, we can not rule out, however, that these exponents are actually the same for the RB and RF cases, and thus maybe universal.

E. Critical exponents
We are now in position to accurately and consistently determine the critical exponents β and z from the nonstationary dynamical behavior of the system.In Fig. 12 we show our determination of β/νz and ζ/z [presented as (1 − ζ/z) since these are expected to be equal] extracted from v(t) and w(t), respectively.To check the consistency of the fitted values we have used three different fitting procedures: (a) First, as is customary, we use the corrected expressions (19) and (20) to extract (ζ/z) and (β/νz), respectively.The fit should be performed above the microscopic regime t > t mic .Since the microscopic  In black (full lines) we show the time tmic on and after which we fitted using the corrected expressions (19) and (20).In red (dotted lines) is shown the time tplat from which we can consider γ(t) to be constant, and we fitted w(t) and v(t) with pure power laws v ∼ t −β/νz and w ∼ t ζ/z , without corrections.In green (dot-dashed lines) we show the exponents as a function of t0, with t0 such that we fitted the data with pure power laws in the window [t0, 4t0].
time t mic is not precisely defined we let it vary in a reasonable range t ∈ [2,20] and compare the values fitted for (ζ/z) and (β/νz) according to ( (b) Second, we use the "plateau in γ" criterion to fit our exponents.That is, we first observe the behavior of γ(t), and propose a time t plat , above which we can consider the plateau well developed.From t > t plat we fit separately w(t) and v(t) with the pure power-law expressions w = w 0 t ζ/z and v = v 0 t −β/νz , respectively.Again, since t plat is not precisely defined, we allow it to move in a wide window t plat ∈ [500, 40000] and observe how the fitted exponents (ζ/z) and (β/νz) behave.
(c) Third, we perform pure power-law fittings of our data, as in the previous case, but now we perform fits within time windows t ∈ [t 0 , 4t 0 ] where t 0 is arbitrarily chosen in a wide range of values.This is to show how the resulting exponent changes depending on the selected time window and to mimic the situation where long-time runs of large systems are not available.
Figure 12 shows that the three fitting procedures lead to the same result β/νz ≈ 1 − ζ/z ≈ 0.128, for sufficiently large t mic , t plat and t 0 .Interestingly, the values fitted for 1 − ζ/z and β/νz in the third procedure are very sensitive to t 0 , and show a clear tendency to produce a positive bias on both at small times.This explains the previous observations of larger effective values for (β/νz) [43].On increasing t 0 , while keeping the window selection as [t 0 , 4t 0 ] (other choices different from the factor 4 give the same qualitative result), the exponent values decrease and at some point the dependence on t 0 ceases.This coincides with the beginning of the range of good values for t plat ∼ 1000.From this analysis we conclude for the depinning transition of the QEW elastic line, that Combining this with the value of ζ = 1.250 ± 0.005, obtained in Sec.IV A with a reliable exact steady-state method, this finally gives z = 1.433 ± 0.007; (22) and considering the statistical tilt symmetry of the model, which yields ν = 1/(2 − ζ), we obtain and F. Scaling relations around the critical point Now that we know the accurate critical exponents of the model, we test them in large scale simulations around the critical point, i.e., at force values around the critical force f c .
In Fig. 13(a) we present the behavior of v(t) with time, at different driving forces, for a string of size L = 2 22 = 4194304.In Fig. 13(b) we scale the data according to Eq. ( 15) with the exponents just obtained.Compared with Fig. 3 for a much smaller system, we now get a much better scaling, especially for large times, and since we are using now the new set of critical exponents.Consistently, Fig. 14(a) we show that the quality of the collapse in Fig. 13(b) is closely accompanied by the development of a plateau in γ(t) with time at f ≈ f c , for t 1000.For forces greater or smaller than f c , γ deviates from its flat behavior towards decreasing or increasing, respectively.Using the facts that in the scaling region we expect v(t, f with a new universal function F (y).This scaling is checked in Fig. 14(b) and, again, the quality of the collapse is closely accompanied by the development of a plateau in γ(t) at f ≈ f c , occurring for t 1000.This shows that our criterion for bounding the scaling region is consistent.

V. DISCUSSION
Power-law corrections to the non-equilibrium scaling have been observed in several other systems such as Ising, Potts, and XY models [57][58][59].Interestingly, it is found that these corrections are stronger for systems with disorder or frustration, or dilution [48].Corrections are usually used to improve the accuracy of the exponents obtained by the STD method.For the depinning transition we have shown that this practice does not much increase the accuracy of the STD method.We think that one reason is that we rely on an ad hoc model for the corrections, add extra parameters to the fit, and still have to decide where to start fitting the corrected formula (the non-universal microscopic time).We have shown that a better approach for depinning is to use scaling relations between the exponents that are violated in the presence of corrections.This allows us to find the truly assymptotic regime where we can expect to get very close to the true exponents.We think it may be useful to try this kind of approach in other problems.At large times, in the critical region, we can think of the relaxation of the velocity as "driven" by the small finite-size bias of the critical force fc−fc(ℓ) ∼ ℓ −1/ν at the scale ℓ = ℓ(t) ∼ t 1/z , where fc is the thermodynamic value.Although the geometry the interface beyond ℓ(t) still retains memory of the initial condition, these large-wavelength modes will not affect v(t).
A heuristic explanation for the power-law corrections can be proposed by relating Eq. ( 11) to the steady-state relation v ∼ ξ −β/ν .As we already pointed out, replacing ξ by the correlation length ℓ(t) ∼ t 1/z yields Eq. ( 11), obtained from general arguments.This suggests two simple pictures: (i) In the first picture we simply introduce phenomenological corrections to the dynamical exponent as ℓ ′ (t) ∼ t 1/z (1 + |c 3 |t −ǫ ) with c 3 and ǫ parameters of the correction.This type of correction has been introduced in studies of the off-equilibrium critical dynamics of the three-dimensional diluted Ising model, for instance, and attributed to the biggest irrelevant eigenvalue of the renormalization group in the dynamics [62].If we assume that v(t) ∼ ℓ ′ (t) −β/ν holds, it is easy to see that such a correction would lead to a corrected velocity v(t) ∼ t −β/zν (1 + |c 2 |t −σ ) with parameters c 2 and σ.Something similar would happen with w(t).Note, however, that if w(t) ∼ ℓ ′ (t) ζ also held, we would not explain the corrections observed, in the same time regime, in γ(t) = w(t)/v(t)t.In any case, it would be interesting to see whether such corrections for ℓ(t), which are certainly absent in the non-steady relaxation of a flat line in the EW model at T > 0 for instance (and more generally of non-Markovian Gaussian signals [63,64] related to manifolds with long-range elasticity), can be explained using the complicated coupling between Fourier modes induced by the non-linearity of the quenched disorder.It is known that "non-Gaussian" effects are indeed subtle [65,66].(ii) In the second picture we assume that the universal non-steady relaxation at the thermodynamic threshold f c is actually "quasistatically driven" by the finite-size bias of the critical force.In other words, the small velocity of the large string is roughly controlled by the steady-state dynamics of an effective string of increasing size ℓ(t) and width w we assume that this finite-size bias with respect to the thermodynamic value f c is positive, and scales with ℓ exactly as the finite-size critical force fluctuations around its mean value, i.e., f c − f c (L) ∼ L −1/ν , we get equivalent to Eq. ( 11).This picture is schematized in Fig. 15.The assumption of a positive bias f c − f c (L) ∼ L −1/ν can be justified from results using the quasistatically velocity-driven ensemble, where the interface is driven by a parabolic potential m 2 (vt − u).For small m, it was shown, via the functional renormalization group (FRG), that [56].Since the parabola imposes a characteristic length L m ∼ m −1 , beyond which the driven interface looks flat, we can write . Identifying L m with ℓ(t), we explain the assumed positive bias and again get Eq.( 11).This identification of lengths is supported by the fact that both L m and ℓ(t) (for the initially flat line) represent the same geometrical crossover.For lengths below ℓ(t)(or L m ), the string is steady-state quasi-equilibrated with the pinning landscape yielding the roughness exponent ζ, but it remains macroscopically flat because of the memory of the initial condition (due to the confinement of the parabola) at scales above ℓ(t) (or L m ).With the picture above we can now speculate about a possible mechanism for generating the scaling correction.Since the formula f c (ℓ) = f c −|c 1 |ℓ −1/ν is expected to hold only in the limit of large ℓ it is plausible to add corrections to it, which decay to zero in the large ℓ limit.If we add a correction of the type f c (ℓ) = f c − |c 1 |ℓ −1/ν (1 + |c 2 |ℓ −α ) with α > 0 and c 2 constants, and use v(ℓ) ∼ [f c − f c (ℓ)] β , we get v(t) ∼ t −β/zν (1 + |c 2 |t α/z ), using ℓ = ℓ(t) = t 1/z , the same type of correction we observe in our simulations.The correction just proposed is not completely ad hoc, but qualitatively consistent (the same sign and order of magnitude) with the corrections to the FRG formula f c (m) = f c − |c 1 |m 2−ζ observed at intermediate values of m in simulations in the quasistatically velocity-driven ensemble [56].It would thus be interesting, in the first place, to check quantitatively whether such corrections are directly related to the ones we report here, and in the second place to check whether one can explain them directly from the FRG-flow behavior.This picture for the corrections suggests that both finite-size steady-state simulations and finite-time non-steady simulations would yield, in particular, an effective exponent for β larger than the true one.This may explain the trend of discrepancies found in the literature.In any case, an analytical and more fundamental description of these effects would be welcomed.
The critical exponents associated with the depinning transition for the QEW model in one dimension as reported in previous works are shown in Table II cise than previous reported values and agrees well with the ones obtained in automaton models [8] and in molecular dynamics simulations [36], but disagrees appreciably with respect to random-field Ising model simulations [37] and other molecular dynamics simulations [52], which give about 30% larger values.All these works correspond to the steady-state dynamics.They thus rely on a proper steady-state equilibration and on a precise knowledge of the critical threshold, which is difficult to achieve, except in periodic samples [10,52].In steady-state simulations with dynamically generated disorder, such as the ones in Ref. [36], one should be cautious at long times or large center of mass displacements, because of the critical force extreme statistics, since, in this case, f c can be considered as the maximum among ∼ M/L ζ independent typical critical forces [53,54,69].Therefore, a finite-velocity steady-state might not exist at zero temperature if the critical force statistics tends to Gumbel's type for large M/L ζ for instance, as the interface will eventually be blocked (by virtue of the Middleton theorems [28]) at any finite-force.It is not clear how a finite temperature can diminish the dominant effect of these rare events at very large averaging times, or eliminate the sensibility to the tails of the distributions.At T = 0, one possibility is to perform a disorder average over finite time averages of the steady-state velocity, such that the distance covered by the center of mass is of order L ζ , and thus the critical force becomes typical again.The problem with this approach is that this critical force fluctuates and thus the control parameter f − f c also varies from sample to sample, complicating the direct estimation of β.In this respect finite samples with periodic disorder have the great advantage that we can calculate accurately the critical force for each sample using an exact algorithm [10].Steady-state simulations exploiting this algorithm were limited by the narrow scaling region for obtaining β [52] however, bounded by finite-size effects on one side, and by the effects of the crossover to the Edwards-Wilkinson equation at large velocities on the other side.Additionally, one should also be cautious with periodic boundary conditions: if M is small compared to the expected random-manifold width, i.e., M ≪ L ζ , we can have several crossovers from the random-manifold to the randomperiodic universality class.In this respect one should be warned that, for a fixed L, the crossovers to the randomperiodic class are not controlled only by M , but also by the velocity [51,61].
The STD method applied to the depinning transition has the advantage, on one hand, that it does not suffer from finite-size effects, as it assumes a growing correlation length which should lie well below the system size.On the other hand, it does not suffer from extreme statistics either, as the spatial region covered in the scaling region is at maximun of order L ζ and thus the critical force and configuration can not be rare.Moreover, we have shown that it actually detects the thermodynamical critical force corresponding to the random-manifold family M = kL ζ , separating the random-periodic case for k → 0, and the extreme Gumbel case k → ∞ [53].We have shown, however, that in order to have a long powerlaw decay in v(t) we must use a large system with the force precisely tuned to the thermodynamic critical force f c .Fortunately, as shown in Fig. 2, f c can be unambiguously and precisely determined from a finite-size analysis of the precisely known sample critical force f s c of periodic samples.It is also usually said that the STD method has the advantage that we do not need to (steady-state) equilibrate the system.The presence of power-law-like scaling corrections shows, however, that we do need to wait (or correct) until the truly asymptotic universal non-steady regime is reached, before determining the exponents.It is also worth noting that the STD method does not yield β or ζ directly, as the direct steady-state method does, but the combination of exponents β/νz.This is not very problematic however, since for the QEW model we accurately determine ζ from the analysis of large critical configurations [29], ν = 1/(2 − ζ), from the statistical tilt-symmetry exact relation, and ζ/z, from the power-law increase of the width w(t), as shown in Fig. 9.
It would thus be interesting to try the STD method in the computationally more convenient automaton models belonging to the QEW universality class.Such models are described in pioneering papers such as Ref. [8].Since by construction they involve a kind of coarse-grained dynamics, it would be instructive to study the non-steady dynamics of this model and see what happens in the mesoscopic time-regime reported here.

VI. CONCLUSIONS
We have studied the non-steady relaxation of a driven one-dimensional elastic interface at the depinning transition.Above a first, non-universal microscopic time regime, we have found a non-trivial long crossover towards the non-steady macroscopic critical regime, expected from general scaling arguments.We have shown that this mesoscopic time regime is robust under changes of the microscopic disorder, including its random-bond or random-field character, and can be fairly described as power-law corrections to the asymptotic scaling forms yielding the true critical exponents.These corrections may explain some numerical discrepancies found in the literature (as large as 30% for some exponents), for this universality class.In particular they explain the appearance of effective power laws in the nonsteady relaxation with exponents presenting a systematic bias with respect to the critical values.To improve the accuracy and consistency of the STD method for extracting critical exponents, we have implemented a practical criterion of consistency and tested it in large-scale (L = 2 25 ) simulations concurrently implemented on GPUs.In this way we obtained accurate exponents for the universality class of the paradigmatic continuum displacement quenched Edwards-Wilkinson equation.Accurate critical exponents are necessary tu succeed in classifying diverse experimental systems and theoretical models into universality classes.It is interesting to note that our estimates coincide, within error bars, with the critical exponents obtained for a stochastic sandpile model in the conserved directed-percolation class [70] The method applied in the present paper may be used to analyze exponents in different depinning universality classes, such as the long-range elasticity, quenched Kardar-Parisi-Zhang, or correlated disorder classes [71].
We believe that the features here observed might be experimentally relevant.In many experiments at low temperatures with magnetic or electric domain walls, or with contact lines of liquids, the dynamics of the interface is non-stationary, in the sense that the system keeps for a while a memory of the initial conditions (i.e. the memory of the preparation or "writing" of the domain wall is visible in the experimental time and length scales), and a growing dynamical correlation length is slowly developed in the presence of a driving field.The STD method might be thus applied experimentally to determine the critical field and exponents, and ultimately to determine the universality class of the system.In particular, we have shown here that scaling corrections are important in a mesoscopic time regime, between the so-called microscopic and macroscopic time-regimes.It would be important to determine whether this regime can be accessed experimentally.In this respect we note that in systems with weak pinning, where the microscopic timeregime t < t mic is determined by a large Larkin length ℓ c , we have ℓ(t mic ) ∼ t 1/2 mic ∼ ℓ c [43].If t mic is large enough, the mesoscopic time-regime might be observable.In such cases, our criterion of consistency may be useful, as we only require monitoring of the time dependence of the interface width and velocity simultaneously (by direct imaging for instance).This would allow corrections from the truly universal part of the relaxation to be disentangled and accurate critical steady-state exponents to be obtained experimentally.Finally, the inclusion of thermal fluctuations would be important in order to compare with some experiments.If the temperature is small enough, we expect a critical universal relaxation yielding the zero-temperature depinning exponents up to length-scales comparable with the thermal rounding length ξ T ∼ T −ψν/β or up to time-scales of the order of ξ z T .The STD method can thus also be used to determine ψ [39].

50 FIG. 1 .
FIG. 1. (Color online) Structure factor of the critical configuration averaged over N = 1000 samples of size L = 2048 and M = 13777.Using the self-affinity developed at long wavelengths, S(q) ∼ q −(1+2ζ) , we fit the roughness exponent ζ = 1.250 ± 0.005.Inset: Rescaled structure factor qualitatively displaying the accuracy of the fit, and the discretization effects at large q.

FIG. 2 .
FIG. 2. (Color online) Dependence of the finite-size critical force f s c on the longitudinal size L for periodic samples of size L × M with M = kL ζ .The dashed line corresponds to the thermodynamic limit fc = 1.916 ± 0.001.

FIG. 3 .
FIG. 3. (Color online) Non-stationary averaged string velocity as a function of time for different driving forces.(a) Raw data.The best power-law curve corresponds to f = f STD c = 1.915.(b) Rescaled data.The parameters f STD c = 1.915, β = 0.33, z = 1.5, and ν = 1.33 were used.The curves correspond to L = M = 8192, and averages were taken over N = 1000 samples.

FIG. 4 .FIG. 5 .
FIG. 4. (Color online) Structure factor of the string at different times during evolution at f = fc.(a) Raw data.(b) Rescaled data.The parameters z = 1.5 and ζ = 1.25 were used.The curves correspond to a system size L = M = 8192, and averages were taken over N = 10000 samples.

FIG. 6 .
FIG. 6. (Color online) Velocity at the corresponding critical forces of the four different cases defined in the text.(a) Raw data.Power-laws fits of the form v ∼ t β/νz for times t 30 give similar results for all curves, consistent with β/νz = 0.129 ± 0.002.(b) Velocity multiplied by t 0.129 .Dashed horizontal lines showing the apparent scaling regime are displayed to guide the eye.All curves correspond to L = 4096, M = 8192, and averages over N = 50000 samples.

FIG. 7 .
FIG. 7. (Color online) γ(t) as a function of time for Case 1 and different system sizes at f = fc = 1.915.The vertical dashed lines separate qualitatively the different dynamical regimes observed.The curves correspond to M = 32768, the averages were taken over N = 5000 samples for L = 512, 1024, 2048, N = 3000 samples for L = 4096, and N = 1700 samples for L = 8192.

FIG. 8 .
FIG.8.(Color online) γ(t) as a function of time for the RB case with Gaussian distributed disorder.The system size is L = 225 = 33554432.The continuous line is a fit using Eq.(18).

FIG. 9 .
FIG. 9. (Color online) String width w(t) as a function of time for the RB case with Gaussian-distributed disorder.The system size is L = 2 25 = 33554432.

FIG. 10
FIG. 10. (Color online) String velocity v(t) as a function of time for the RB case with Gaussian-distributed disorder.The system size is L = 2 25 = 33554432.

FIG. 11 .
FIG. 11. (Color online) String velocity v(t) as a function of time for the RF case with Gaussian-distributed disorder.The system size is L = 2 25 = 33554432.

FIG. 12 .
FIG. 12. (Color online) The combination of exponents 1−ζ/z and β/νz extracted from fittings to the data of Figs. 9 and 10, respectively.The horizontal axes shows the chosen values of three different arbitrary quantities for the fitting procedure.In black (full lines) we show the time tmic on and after which we fitted using the corrected expressions(19) and(20).In red (dotted lines) is shown the time tplat from which we can consider γ(t) to be constant, and we fitted w(t) and v(t) with pure power laws v ∼ t −β/νz and w ∼ t ζ/z , without corrections.In green (dot-dashed lines) we show the exponents as a function of t0, with t0 such that we fitted the data with pure power laws in the window [t0, 4t0].

FIG. 13 .
FIG.13.(Color online) String velocity v(t) as a function of time for the RB case with uniformly distributed disorder for which fc = 1.5652 using δt = 0.1.The system size is L = 4194304.In (a) we present the raw data, and in (b) v(t, f ) has been rescaled to vt β/νz and t to t 1/z |f − fc| ν .

FIG. 14 .
FIG.14.(Color online) γ(t) as a function of time for the RB case with uniformly distributed disorder for which fc = 1.5652 using δt = 0.1.The system size is L = 4194304.In (a) we present the raw data, and in (b) t has been rescaled.

FIG. 15 .
FIG.15.(Color online) Heuristic connection between steady state and the non-steady universal relaxation at the thermodynamical critical force fc.The circle represents the system and the short arrows indicate its non-steady time evolution.At large times, in the critical region, we can think of the relaxation of the velocity as "driven" by the small finite-size bias of the critical force fc−fc(ℓ) ∼ ℓ −1/ν at the scale ℓ = ℓ(t) ∼ t 1/z , where fc is the thermodynamic value.Although the geometry the interface beyond ℓ(t) still retains memory of the initial condition, these large-wavelength modes will not affect v(t).

TABLE II .
, together with the values obtained in the present work.Comparing previous estimates, the largest discrepancy is found for the β exponent.Our value β = 0.245 ± 0.006 is more pre-Representative values for the depinning exponent ζ, the dynamical exponent z, and the critical exponents ν and β, reported in the literature.The major dispersion is seen in the value of β.