Critical dynamics of classical systems under slow quench

We study the slow quench dynamics of a one-dimensional nonequilibrium lattice gas model which exhibits a phase transition in the stationary state between a fluid phase with homogeneously distributed particles and a jammed phase with a macroscopic hole cluster. Our main result is that in the critical region ({\it i.e.}, at the critical point and in its vicinity) where the dynamics are assumed to be frozen in the standard Kibble-Zurek argument, the defect density exhibits an algebraic decay in the inverse annealing rate with an exponent that can be understood using critical coarsening dynamics. However, in a part of the critical region in the fluid phase, the standard Kibble-Zurek scaling holds. We also find that when the slow quench occurs deep into the jammed phase, the defect density behavior is explained by the rapid quench dynamics in this phase.

Introduction. -In the past few decades, extensive studies have been carried out to understand the phase ordering dynamics of classical systems with equilibrium [1] and nonequilibrium steady states [2][3][4] when the system is quenched infinitely fast from a disordered state to an ordered one. Slow quench (or annealing) dynamics have also attracted some attention in the recent years, and have been invoked to understand the defect structures in the early universe [5,6] and more generally, systems exhibiting second-order phase transitions [7,8] in both classical [9][10][11][12][13][14][15] and quantum [16][17][18][19][20] settings. Moreover, a number of experiments investigating the relationship between defect density and quench rates have also been performed in a variety of systems such as the liquid crystals, superfluid 3 He, superconductors, Bose-Einstein condensates and colloidal systems, see [21] for a recent review. The defect density at the end of the quench is generally found to decay as a power law in the inverse quench rate (although some systems such as the 2D XY model exhibit non-algebraic decay [22,23]).
Much of the body of work on slow quench dynamics appeals to the Kibble-Zurek argument [5,7] which states that if the control parameter is varied slowly across the critical point, the system stays close to the steady state (adiabatic regime) until its relaxation time becomes longer than the quench time after which the dynamics are hypothesised to remain "frozen" (impulse regime) until the critical region is crossed. Thus, the Kibble-Zurek argument describes the quench dynamics before the critical point is crossed. Recent works have elucidated the slow annealing dynamics when the system is quenched far from the critical point to an ordered phase and argued that the defect density is determined by coarsening dynamics [13]. The dynamics in the 2D XY model in which the system is quenched slowly from the high-temperature disordered phase to the low-temperature critical phase have also been studied [22]. However, a complete study of a system with isolated critical point in which one can examine the change in the dynamical behavior when a slow quench occurs below, at and above the critical point has not been carried out. Moreover, all the studies mentioned above deal with systems with equilibrium steady state (but, see [24] that considers a nonequilibrium quantum system).
Here we consider the slow quench dynamics of the jamming transition that occurs in diverse settings such as vehicular traffic [25], cellular traffic [26] and granular media [27]. We study a classical nonequilibrium system in one dimension which shows a jamming transition in the stationary state [28]. The steady state of this model is known exactly [29], and some results for the coarsening dynamics [4,30] and stationary-state dynamics [31] have also been obtained. However, the dynamics of this model under slow annealing have not been studied and here we address this problem using numerical simulations. We find that the standard Kibble-Zurek scaling explains our results in a part of the critical region but close to the critical point and for quenches deep in the jammed phase, the defect density decay can be understood using the corresponding results for rapid quench dynamics [13].
Models. -We consider a one-dimensional lattice with periodic boundary having L sites and N particles. The total number of particles in the system is conserved. Due to hard-core interactions between the particles, each site can have at most one particle. As shown in fig. 1, a particle hops to its right empty neighbor with a rate u(n) where n is the number of holes in front of it (unidirectional model ).
It is useful to map this model to a zero range process (ZRP) with L = N sites and N = L − N mass units by considering particles in the lattice gas model as sites and holes as mass clusters [29]. Thus, in ZRP, a site can contain any number of particles and a particle hops out of a site containing n particles to its left neighbor with a rate u(n), see fig. 1 for an illustration. For an infinitely large system in the stationary state, the probability that a site contains n particles is given by [29] where In the above equation, g(ω) = ∞ n=0 ω n f (n) is the generating function of f (n) and the fugacity ω is calculated from the fugacity-density relation, From the above equation, it can be seen that the fugacity is an increasing function of the density . For certain choices of u(n), the fugacity reaches its maximum value (given by the radius of convergence of the generating function g(ω)) at a finite density c . Then for > c , the excess density − c resides at a single site while for < c , each site supports a density [29].
For the hop rate the ZRP shows the phase transition described above in the -b plane for b > 2. Correspondingly, the lattice gas model exhibits a phase transition at the critical point given by  where ρ = N/L is the total particle density. For b < b c , the typical hole cluster length is of order unity (fluid phase) while for b > b c , a macroscopically long hole cluster coexists with gaps that are power law distributed as n −b (jammed phase) [29]. It has been shown that close to the critical point, the static correlation The stationary-state dynamics have also been studied and it has been found that at the critical point, the steadystate density fluctuations decay on a time scale that grows as L z where the dynamic exponent z = 3/2 [31].
In the following, we also consider a bidirectional model in which the particle first chooses either the left or the right neighbor with equal probability and then hops with a rate that depends on the vacancies in the chosen direction provided the target site is empty. The steady state obeys detailed balance and is the same as that in the unidirectional model [29]. As a result, the correlation length exponent ν is given by (6); however, the dynamic exponent z = 2 in this case [31].
To study the slow quench dynamics, we introduce time dependence in the hop rates and write where, for simplicity, we work with linearly varying b given by The quench protocol was carried out by changing the parameter b from zero (in the fluid phase) to a final value b τ = b c (critical point) and 2b c (jammed phase) keeping 26003-p2 0.4 0 .6 0 .8 the density fixed at ρ c given by (5). Our quantity of interest is the domain wall density (which is the interface between the particle and hole) in the lattice gas model.
From the inset of fig. 1, we see that the number of domain walls is equal to two times the number of occupied sites in the ZRP. For finite inverse quench rate τ , as the system is far from its steady state and has more domain walls than in the stationary state, we consider the excess defect density given by where p(0, t) is the probability that a site is empty at time t in ZRP when the parameter b is time-dependent and p(0) is given by (1). In Monte Carlo simulations of the models described above, we measured δρ d (t) for system sizes in the range 15000-20000 and averaged the data over 2000-4000 independent initial conditions.
Results. -When the quench rate τ −1 is very small, the parameter b changes very slowly allowing the system to relax to the stationary state. But for faster quench, the system is farther from the steady state. Indeed, as fig. 1 shows, the excess defect density δρ d (t) decreases with increasing τ . Our objective here is to understand how δρ d (τ ) decays with τ when the system is quenched slowly to b τ .
Dynamics in the fluid phase close to the critical point: When the system is far from the critical point, it relaxes quickly. But as the critical point is approached, the relaxation time increases and at time t * , the remaining time t r = τ − t * to reach the critical point becomes comparable to the relaxation time in the stationary state which can be expressed as [5,7]  so that t r ∼ τ zν 1+zν . As the system falls out of equilibrium at t * , in numerical simulations, we picked the time t * to be the one where the excess domain wall density is 10 −3 and found the time t r . Using the exponents ν and z quoted in the last section, we find that for the unidirectional model, the time t r ∼ τ 3/(2bc−1) for 2 < b c < 3 and as τ 3/5 for b c ≥ 3 which is in good agreement with the numerical data in the main panel of fig. 2. (We have checked that our scaling results are not affected if t * is determined by the criterion that the excess defect density 10 −3 .) Assuming that the system does not evolve after time t * , the defect density after crossing the critical point is posited to be δρ d (τ ) ∼ ξ −1 * ∼ τ −ν/(1+νz) [5,7]. However, we find that the Kibble-Zurek scaling works well for times smaller than τ but not at or after the critical point is crossed (see below). Our simulation results for the unidirectional model shown in the inset of fig. 2 for a quench to the critical point demonstrate that the excess defect density is of the form where f 1 (x) and f 2 (x) are the scaling functions, and t * < t < τ. We also find that the above scaling form breaks down at times of order t * . Dynamics at the critical point : Close to the critical point, the system undergoes critical coarsening during the time interval t r [13]. We therefore turn to a discussion of the fast quench dynamics of the ZRP in which the system initially in the fluid phase is quenched instantaneously to the critical point. To distinguish between the quantities obtained using slow and fast quench, in the following, we useˆto refer to quantities obtained using the fast quench protocol.
In [30], critical coarsening dynamics of the ZRP have been investigated in mean-field geometry and in one 26003-p3 dimension. In the latter case, numerical simulations indicate that a measure of the domain length grows with time as t 1/ẑ with the coarsening exponentẑ = 3 (5) for unidirectional (bidirectional) model and a scaling argument suggests that the probability δp(0, t) ∼ t −α with the exponentα = (b − 2)/ẑ for b > 3 (see (35) of [30]). While our numerical results forẑ are in agreement with those of [30], the exponentα for the distribution of empty sites is not consistent with that in [30]. Instead, our numerical results shown in figs. 3 and 4 suggest that the exponentα varies with b c for 2 < b c < 3 but it is a constant otherwise, and the best regression fit giveŝ Using the above results and recalling that the coarsening process initiates at time of order t * , we obtain δρ d (τ ) ∼ t −α r ∼ τ −zνα/(1+zν) when the system is quenched to the critical point. Here we have ignored the dependence on the quench depth (i.e., b c −b(t * )) since our numerical results in fig. 3 suggest that the long-time dynamics are independent of it. More explicitly, for the unidirectional model, we have while, for the bidirectional model, we get The above predictions for the excess defect density are consistent with the numerical results shown in figs. 5 and 6 when b τ = b c . We have also checked that the above scaling predictions hold when the system is quenched in the vicinity of b c . Dynamics deep in the jammed phase: For b τ b c , the system undergoes coarsening in the jammed phase after the time t c + t r where b(t c ) = b c [13]. As the time scale t c ∼ τ but t r is sublinear in τ , the time left until the quench is of order τ . Thus, we expect the defect density to simply decay as where the coarsening exponent in the ordered phase,ẑ O = 2 (ẑ O = 3) for unidirectional (bidirectional) model [4,30]. Figure 7 shows that our numerical results for b τ = 2b c are consistent with the above scaling prediction. Moreover, the excess defect density δρ d (t) is of a scaling form similar to (11) as attested by the data collapse shown in the inset of fig. 7.
Conclusions. -Slow quench dynamics have been studied extensively when the control parameter is changed across the critical point of a second-order phase transition in equilibrium systems [7,8], and recent works have considered slow quenches deep into the ordered phase [13] or the critical phase [22]. Here we have performed, to our knowledge, the first quantitative study of the slow quench dynamics when a nonequilibrium system with isolated critical point is quenched in the critical region.
The Kibble-Zurek argument assumes that in the critical region, the dynamics are frozen since the relaxation times are much longer than the time available. As shown in the inset of fig. 2, this argument indeed holds when the system is close to the critical point; more precisely, we find that the Kibble-Zurek scaling is obeyed over a time window [t * , ft * ], f < 1, when the system is quenched to the critical point. However, it does not apply during f t * < t < τ (and similarly, in the symmetric time window beyond the critical point in the jammed phase). Instead, we find that   the defect density decays as a power law although with an exponent smaller than or equal to the Kibble-Zurek prediction, see figs. 5 and 6. While the dynamics outside the critical region involve only the stationary-state dynamic exponent (z) and the deep quench in the ordered phase is determined by the coarsening exponent (ẑ O ), the critical point quench dynamics are more complex involving both static fluctuations and critical coarsening. For the class of models considered here, a comparison of the decay exponents shows that the excess defect density decays faster in the unidirectional model which has a nonequilibrium steady state than in the bidirectional model with equilibrium steady state. We also obtain continuously varying exponents for 2 < b c < 3 where the mass fluctuations are anomalous but constant otherwise [29].
A detailed exploration of other nonequilibrium and equilibrium systems with critical point annealing should inform us better about the dynamics in the impulse regime and is a goal for the future. * * *