Kibble-Zurek Mechanism for Nonequilibrium Phase Transitions in Driven Systems with Quenched Disorder

We numerically study the density of topological defects for a two-dimensional assembly of particles driven over quenched disorder as a function of quench rate through the nonequilibrium phase transition from a plastic disordered flowing state to a moving anisotropic crystal. A dynamical ordering transition of this type occurs for vortices in type-II superconductors, colloids, and other particle-like systems in the presence of random disorder. We find that on the ordered side of the transition, the density of topological defects $\rho_d$ scales as a power law, $\rho_d \propto 1/t_{q}^\beta$, where $t_{q}$ is the time duration of the quench across the transition. This type of scaling is predicted in the Kibble-Zurek mechanism for varied quench rates across a continuous phase transition. We show that scaling with the same exponent holds for varied strengths of quenched disorder. The value of the exponent can be connected to the directed percolation universality class. Our results suggest that the Kibble-Zurek mechanism can be applied to general nonequilibrium phase transitions.


INTRODUCTION
In systems that exhibit second-order phase transitions, a disordered phase on one side of the transition, such as a liquid or glass state, can be characterized by the presence of topological defects, while on the other side of the transition, there is an ordered or defect-free state such as a crystal.The defect-free ordered phase arises when the parameter controlling the transition is changed very slowly, so that the system remains in the adiabatic limit.If the rate of change through the transition increases, a portion of the topological defects do not have time to annihilate, and persist on the ordered side of the transition.A scenario for understanding the behavior for varied quench rates across a continuous phase transition is the Kibble-Zurek (KZ) mechanism, which predicts a universal power law for the defect density ρ d ∝ t −β q , where t q is the time duration of the quench through the transition, so that slower quenches produce fewer defects [1][2][3][4].The KZ mechanism has been studied in a variety of systems at the transition into an ordered phase, where the scaling exponents can be related to the universality class of the underlying phase transition [5][6][7][8][9][10][11].
The KZ scenario has generally been applied to systems that have equilibrium phase transitions; however, there have been recent proposals to use the KZ scenario to address transitions between different nonequilibrium steady states [12][13][14][15][16][17], such as those that occur in optical systems or for Rayleigh-Bènard convection, where defects can arise in otherwise hexagonal ordered lattices.Another class of nonequilibrium systems consist of assemblies of particles driven over quenched disorder that also exhibit behavior consistent with a continuous phase transition from a disordered state to an ordered state [18][19][20][21][22] or from a dynamical fluctuating state to a non-fluctuating state [23][24][25][26].An open question is whether the KZ scenario could also apply to nonequilibrium phase transitions for varied sweep rates through the transition.Nonequilibrium systems have several features that could make them ideal for studying the KZ scenario.They often contain very well defined topological defects, and the transition can occur at T = 0, so that thermal effects such as critical coarsening on the ordered side of the transition are absent.
One of the best examples of a system that shows evidence for a continuous nonequilibrium phase transition as a function of a continuously changing driving parameter from a disordered fluctuating state with a high density of topological defects to a dynamically ordered nonfluctuating state in which topological defects are scarce or absent is superconducting vortices driven over quenched disorder [18,19,[27][28][29][30][31][32][33][34][35][36][37][38].At T = 0 and in the absence of quenched disorder, a superconducting vortex system forms a triangular lattice free of defects; however, when quenched disorder is present, a disordered state containing numerous topological defects can appear that can be characterized as a vortex glass [39].A finite external drive in the form of an applied current causes vortices in the disordered state to depin and move [21,39].For drives above depinning, the system enters a strongly fluctuating or plastically flowing state in which there is a fluctuating number of topological defects, while at higher drives, there is a critical driving force above which the vortices dynamically order into a state with zero or a small number of topological defects [19,21,[29][30][31][32][33][34].This ordered state is not isotropic but takes the form of a moving anisotropic crystal [21,30,31,33,34] or a moving smectic state [29,[31][32][33][34][35].Here, a small number of topological defects in the form of dislocations are present which have their Burgers vectors aligned in the direction of the drive.
The critical drive amplitude at which the ordering transition occurs is a function of vortex density and disorder strength [21].The ordering transition has been observed experimentally by direct imaging [33,36], neutron scattering [27], changes in the noise [34,35,38], and changes in transport curve features [18, 21, 28-30, 34, 37].This same type of dynamical ordering transition is general to the broader class of particle-like assemblies driven over quenched disorder, and it has been studied for colloidal particles moving over disordered landscapes [20,21,[40][41][42], moving Wigner crystals [43], sliding charged systems [44], driven pattern forming states [45,46], driven dislocations [47], and magnetic skyrmions moving over quenched disorder [48].Since the transition from the disordered state to the ordered state is controlled by changing the driving amplitude, the number of topological defects on the ordered side of the transition can be measured for varied drive sweep rates across the transition to test the predictions of the KZ scenario.
Here, we numerically examine the density of defects in driven vortex and colloidal systems at T = 0 for varied driving sweep rates through the dynamical ordering transition.We find that the defect density obeys ρ d ∝ τ β q , where t q is the time required to go through the transition, consistent with the KZ scenario.The same behavior appears for various values of the quenched disorder and for both vortices and colloids moving over a random substrate.The exponent we find in both systems is β ≈ 0.39, which is consistent with an underlying transition that falls in a directed percolation (DP) universality class, which often describes nonequilibrium phase transitions from fluctuating states to non-fluctuating states [23,24,26].These results could be tested experimentally in superconducting vortex systems, colloids, or other driven systems with quenched disorder.Our results also imply that the KZ mechanism could be applied generally to nonequilibrium phase transitions in the same way that it has been applied to equilibrium phase transitions, and that it could be tested in many other types of driven systems that show dynamical continuous transitions from disordered to ordered states.

Modeling and characterization of the nonequilibrium phase transition
We use a well-established simulation model for superconducting vortices driven over random disorder.The vortices are represented as repulsively interacting point particles which undergo a dynamical ordering transition from a disordered plastic flow phase to a moving ordered crystal or smectic state [19,21,29,30,34,35].We consider a two-dimensional (2D) system of size L × L, with L = 36λ in units of the London penetration depth, containing a fixed number N of vortices produced by an external magnetic field B = Bẑ.The vortex density is n v = N/L 2 .The motion of vortex i at position R i is obtained by integrating the overdamped equation of motion where the damping constant η = 1.0.The vortex-vortex interaction force 0 /2πµ 0 λ 3 , µ 0 is the permeability of free space, and φ 0 = h/2e is the elementary flux quantum.The modified Bessel function K 1 falls off exponentially for R ij /λ > 1.0.The vortices are confined to the x, y plane and the magnetic field is applied in the z-direction.The quenched disorder is modeled by N p randomly placed harmonically attractive pinning sites of radius r p with k is the location of pinning site k, and R(p The pinning density is n p = N p /L 2 .The driving force F d is applied uniformly to all the vortices, and corresponds experimentally to a Lorentz force where J is the applied current.The initial vortex configurations are obtained using simulated annealing.We increase the magnitude of the dimensionless driving force from F D = 0 to a final maximum value over a total time of t q in increments of ∆F D = 0.002.Time is reported in dimensionless units achieved by scaling time by τ = η/f 0 .We characterize the system by measuring the average velocity per vortex in the driving direction, x and the fraction of sixfold coordinated vortices , where z i is the coordination number of vortex i obtained from a Voronoi construction.We also compute the fraction of topological defects P D = 1 − P 6 .

Numerical results
Figure 1 illustrates the effect of changing F D in the adiabatic limit for a system with n v = 1.0, n p = 0.5, F p = 0.4, and r p = 0.3.In Fig. 1(a) we plot the average vortex velocity V along with d V /dF D versus F D for t q = 7.5 × 10 6 .If t q is increased to a larger value, the curves remain nearly unchanged.Figure 1(b) shows the corresponding fraction of sixfold coordinated vortices P 6 versus F D , where we would have P 6 = 1 for a triangular lattice.A depinning transition occurs near F D = 0.06, and for F D < 0.225 we find P 6 ≈ 0.64, indicating a highly defected or disordered state in which the vortices are undergoing plastic flow.In Fig. 2(a) we show a Voronoi plot of a snapshot of the vortex positions from the system in Fig. 1   form an anisotropic triangular lattice of the type shown in Fig. 2(b) at F D = 1.5 in the ordered state.Here the lattice is aligned in the driving direction and contains a small number of aligned dislocations, forming a smectic state as studied previously [29-31, 33, 34].The critical force for the transition from the disordered state to the dynamically ordered phase is F c D ≈ 0.3.We do not observe any hysteresis across the transition, which has features consistent with a second-order transition.In Fig. 1(a), there is a peak in d V /dF D in the plastic flow phase near F D = 0.15 followed by a saturation near F D = 0.3 to d V /dF D ≈ 1.0, in agreement with previous studies of the dynamical ordering transition [18,28,30,35,37,38].The critical driving force for dynamic ordering also depends on the strength of the disorder, as shown in Fig. 1(c) where we plot a dynamical phase diagram for the system in Fig. 1(a,b) as a function of F p versus F D at t q = 7.5 × 10 6 .The disordered and ordered phases are highlighted, and the arrow illustrates the direction in which we sweep across the transition at different rates.Now that we have established the driving force which separates the disordered and ordered phases, we can cross this transition for varied t q .In Fig. 3(a) we plot P 6 versus F D for the system in Fig. 1 at t q = 7.5×10 6 , 3×10 6 , 7.5× 10 5 , 3×10 5 , 7.5×10 4 , 3×10 4 , 7.5×10 3 , 3×10 3 , and 750, showing that as the quench rate increases and t q becomes smaller, the fraction of six-fold coordinated particles on We illustrate Voronoi plots of the final F D = 1.5 state for the system in Fig. 3 at t q = 7.5 × 10 3 in Fig. 4(a), t q = 3 × 10 4 in Fig. 4(b), t q = 7.5 × 10 4 in Fig. 4(c), and 3 × 10 5 in Fig. 4(d), indicating that the number of remaining defects is larger for smaller t q .Figure 2(b) shows the same sample in the adiabatic limit with t q = 7.5 × 10 6 .The defects that appear on the ordered side of the transition generally take the form of dislocations composed of pairs of fivefold and sevenfold coordinated particles.At larger t q , the dislocations have their Burgers vectors oriented in the direction of drive, giving rise to a smectic ordering in which the system can be regarded as a set of one-dimensional (1D) moving channels.In Fig. 3(b) we plot P D = 1 − P 6 , the density of topological defects at the final F D = 1.5 state on the ordered side of the transition, versus t q .We find P D ∝ t β q , where the solid line indicates a fit with β = −0.385.The KZ scenario predicts a power law decay of P D with increasing t q .If we choose a final value of F D which is closer to but still above the critical transition drive F c , the scaling is not as good at smaller t q but we still find the same exponent for the larger values of t q .
The plots of P D versus t q in Fig 3 (c) indicate that the same scaling behavior remains robust over a range of different values of F p from F p = 0.2 to F p = 1.2.Since F c varies with F p , for each sample we select a final value of F D that is the same distance above F c for consistency.The straight line is a power-law fit with β = 0.39.For large F p , the scaling at larger values of t q is not as good because more defects become trapped by the pinning.
In addition to superconducting vortices, we have also considered colloidal particles driven over quenched disorder.Here we use the same type of overdamped simulations but replace the particle-particle interactions by the Yukawa or screened Coulomb potential U (r ij ) = A exp(−κr ij )/r ij [49], with κ = 1.0 and A = 0.1.The colloidal particles experience a stronger short-range repulsion than the vortices.In experiments, various types of random quenched disorder can be introduced along with an external driving force [42,50], and the number of topological defects can be measured with imaging techniques [50,51].In Fig. 5(a) we show P 6 versus F D in a colloidal system with colloidal density n c = 1.0, n p = 0.5, F p = 1.0, and r p = 0.35 at t q = 3 × 10 6 , 7.5 × 10 5 , 3.0 × 10 5 , 7.5 × 10 4 , 3.0 × 10 4 , 7.5 × 10 3 , 750, 300, and 75, where the final P 6 decreases for smaller t q .In Fig. 5(b) we plot the final value of P D at F D = 2.5 versus t q for the system in Fig. 5(a).The solid line is a power-law fit with β = −0.39,similar to the exponent obtained for the vortex case.

DISCUSSION
We can ask whether the scaling exponent we obtain can be related to possible universality classes of the dynamical transition.In the KZ scenario, the scaling exponent obeys β = (D − d)ν/(1 + zν), where ν and z are critical exponents associated with the universality class of the underlying phase transition, D is the spatial dimension of the system, and d is the dimension of the defect.For a 2D Ising model with D = 2, β = 2/3, which is not what we observe; however, since the defects in our sample are all aligned in the driving direction, our system is closer to coupled 1D channels, which would give β = 1/3.Recent simulations of certain quenched 2D spin ice models give β = 0.31 [52].A more likely candidate for the universality class of a nonequilibrium system is directed percolation (DP) [23], where the critical exponents depend on the effective dimension of the dynamics.If we assume dynamics of dimension 1+1, we would have z = 1.58 and ν = 1.097 [23], giving β = 0.401, which is in agreement with the values we observe.We argue that the 1 + 1 dynamics may be relevant since the topological defects are mostly aligned in a single direction on the ordered side of the transition, and the relevant length scale could be the distance between the defects in the direction of drive rather than perpendicular to the drive.
In the KZ scenario, the lag time between the nonequilibrium and equilibrium value is given by the freeze-out time t ∝ t zν q /(1 + zν) [4].Since the value of the driving force F D is a function of time, and since F D and t q are dimensionless, this implies that the defect fraction should scale as P D ∝ F D /t α q , where t q is the quench time.In Fig. 6(a) we illustrate this scaling for the vortex system from Fig. 3(a) and in Fig. 6(b) we show the scaling for the colloidal system from Fig. 5(a).In each case, the scaling exponent is α = 0.625.The KZ prediction gives zν/(1+zν) = α, where plugging in the scaling exponents for 1+1 directed percolation [23] leads to α = 0.632, consistent with the exponents we find.
To better understand why the dynamics is 1+1 dimensional, in Fig. 7(a) we illustrate the colloidal positions and trajectories on the disordered side of the transition for a subset of the system in Fig. 5 at F D = 0.25 and t Q = 7.5 × 10 5 , showing strongly disordered flow occurring in both the x and y directions.The corresponding Voronoi plot in Fig. 7(b) contains a high defect density, similar to what is observed in the vortex system in the disordered phase.Figure 7(c) shows the positions and trajectories of the colloids on the ordered side of the transition at F D = 1.5.Here the dynamics are strictly 1D in character, and the topological defects are all aligned in the direction of drive as shown in Fig. 7(d).In theoretical work for particles moving over random disorder in 2D, it is argued that the strongly driven case can be considered as a series of coupled 1D channels that slide past one another to form a moving smectic state [31,32].This 1D channeling could explain why the dynamics produce scaling exponents consistent with 1 + 1 DP rather than 2 + 1 DP.
In conclusion, we have investigated the evolution of the density of defects across a dynamical disorder to order nonequilibrium phase transition for driven particles moving over quenched disorder as we vary the quench rate.We find that the defect density scales as a power law with quench rate, in agreement with the predictions of the Kibble-Zurek scenario.For both superconducting vortices and colloidal assemblies, we find a scaling exponent of β ≈ −0.39, which is consistent with an underlying transition that falls in the 1+ 1 directed percolation universality class since the ordered system forms a moving smectic state in which the defects are aligned in 1D chains.Experimentally, our predictions could be tested in superconducting vortex systems using various imaging and transport measures that have previously been shown to be correlated with the number of defects in the vortex lattice [28].Direct imaging of the dynamics is feasible using colloidal systems, and it would also be possible to consider ac drives which would avoid edge effects.Our results suggest that the Kibble-Zurek scenario can be applied to more general non-equilibrium continuous phase transitions, particularly those between disordered and ordered states.Our results also imply that a system undergoing a directed percolation transition could exhibit features of the Kibble-Zurek scenario provided that some type of well-defined defect structure can be identified.

FIG. 1 .FIG. 2 .
FIG.1.Dynamical phase transition in a superconducting vortex system.(a) The velocity V versus applied drive FD (blue) and the corresponding d V /dFD vs FD (red) for a 2D superconducting vortex system with quenched disorder at vortex density nv = 1.0, pinning density np = 0.5, pinning force Fp = 0.4, and pinning radius rp = 0.3.The curves are obtained for a quench time of tq = 7.5 × 10 6 , which is considered the adiabatic limit.(b) The corresponding fraction of sixfold coordinated vortices P6 vs FD.A nonequilibrium transition from a disordered state to an ordered state occurs for FD ≈ 0.3.(c) Dynamical phase diagram as a function of Fp versus FD for the same system highlighting the disordered phase (green) and the dynamically ordered phase (blue).The arrow indicates the direction in which the transition is crossed at different sweep rates.

FIG. 6 .
FIG. 6.Defect density scaling in the vortex and colloid systems.Dimensionless scaling plots of PD versus FD/t α q , where tq is the quench time and where α = 0.625.(a) The vortex system from Fig.3(a).(b) The colloidal system from Fig. 5(a).

FIG. 7 .
FIG. 7. Disordered and 1+1 dimensional flow below and above the dynamical phase transition.(a) The colloid positions (red circles) and trajectories (lines) for a subset of the system in Fig. 5 with nc = 1.0, np = 0.5, Fp = 1.0, and rp = 0.35 at tq = 7.5 × 10 5 and FD = 0.25, showing 2D disordered flow.(b) The corresponding Voronoi plot showing a high defect density.(c) The colloid positions and trajectories in the same system at FD = 1.5 showing 1D channeling.(d) The corresponding Voronoi plot showing that all of the topological defects are aligned in the same direction.