Scaling Law for Irreversible Entropy Production in Critical Systems

We examine the Jarzynski equality for a quenching process across the critical point of second-order phase transitions, where absolute irreversibility and the effect of finite-sampling of the initial equilibrium distribution arise in a single setup with equal significance. We consider the Ising model as a prototypical example for spontaneous symmetry breaking and take into account the finite sampling issue by introducing a tolerance parameter. The initially ordered spins become disordered by quenching the ferromagnetic coupling constant. For a sudden quench, the deviation from the Jarzynski equality evaluated from the ideal ensemble average could, in principle, depend on the reduced coupling constant ε0 of the initial state and the system size L. We find that, instead of depending on ε0 and L separately, this deviation exhibits a scaling behavior through a universal combination of ε0 and L for a given tolerance parameter, inherited from the critical scaling laws of second-order phase transitions. A similar scaling law can be obtained for the finite-speed quench as well within the Kibble-Zurek mechanism.

Scientific RepoRts | 6:27603 | DOI: 10.1038/srep27603 irreversibility [23][24][25][26][27][28] . A process is called absolutely irreversible if there exists a path in phase space whose probability to occur in the forward direction is zero while that in the reverse direction is nonzero, or vice versa. A typical situation occurs when the accessible phase spaces for the system at the beginning and end of a protocol are not identical. This is indeed the case for the free expansion of initially confined particles whose accessible phase space is increased by removing the partitioning barrier.
In this work, we explore the fact that in systems driven through second-order phase transitions, both the absolute irreversibility and the convergence issue can take place in a single setup with equal significance. Using the scaling theory of phase transitions and numerical simulations, the deviation from the ideal ensemble average JE, i.e., the ensemble average of e −σ using infinite number of samples, is examined as a function of the system size and the reduced coupling constant. It exhibits a universal scaling behavior inherited from the critical scaling of the correlation length and the relaxation time in second-order phase transitions. This finding may provide a unique application of the FTs to study the dynamical properties of phase transitions. We note that the absolute irreversibility due to the spontaneous symmetry breaking is analogous to the ergodicity loss (within a finite observation time) in a finite system with metastable states. The latter has been explored theoretically by assuming partially equilibrated states 29,30 , and demonstrated in an experiment with a Brownian particle subject to a bistable potential of two optical traps 31 or similar experiments 32,33 in the context of Landauer's principle. In our work, we explore how the ergodicity breaking intensifies with the system size and reveal that the effects of spontaneous symmetry breaking are reflected in the scaling behavior of the deviation from the ideal ensemble average JE.
While the detailed arguments and analyses are discussed below, a summary of our discussion is as follows: On the one hand, a natural partitioning of the phase space emerges as a consequence of the ergodicity breaking in the ordered phase 34 in contrast to the partitioning externally imposed in the example of free expansion. The resultant absolute irreversibility is illustrated in Fig. 1 for the Ising model, which as the simplest model showing spontaneous symmetry breaking (SSB) will be used to anchor the rest of our discussions. It is expected that the SSB of  2 symmetry corresponds to halving the accessible phase space, resulting in 〈e −σ 〉 = 1/2 when the system is quenched from the equilibrium ordered to the disordered phase. On the other hand, in such a process the configurations with vanishing (spatial) mean order parameter give major contributions to 〈e −σ 〉, while such configurations are extremely rare in the initial equilibrium in the ordered phase. Here, such an insufficient sampling is accounted for by introducing a tolerance parameter to neglect some unlikely configurations and in general this leads to lower value of 〈e −σ 〉 than for the ideal case. We stress that an observation over a finite time in realistic experiments and numerical simulations inevitably leads to insufficient sampling. Our main result is that for a given tolerance the deviation from the ideal ensemble average JE, neither zero nor exactly one half, is still determined in terms of a universal combination of the reduced coupling constant and the system size. In the forward process, the system is initially at equilibrium with positive spontaneous magnetization, whereas in the backward process the initial equilibrium state has no magnetization. When the coupling J increases across the critical point, the system can have either positive or negative magnetization. The latter case has no corresponding forward path, which results in the absolute irreversibility.
j j and the set of all spin configurations by . The configuration space  consists of ±  and  0 ,  is naturally empty for odd N. For even N, 0  is negligible (probability measure is zero) in the thermodynamic limit and hereafter ignored.
We consider, for simplicity and to conform to the standard protocols of the JE, quenching processes where the coupling constant J(t) varies while temperature is kept constant. As usual, we define the reduced coupling constant by c where J c is the critical point. As the time t changes from To discuss absolute irreversibility, we will be mostly interested in quenching from the ordered (ε i < 0) to disordered (ε f > 0) phase. For simplicity, we consider symmetric quenching: Although the quenching process drives the system out of equilibrium, many physical effects are still described in terms of the initial and final equilibrium distribution functions where Z i/f are the respective partition functions. Since we start from the ordered phase at initial time, the allowed spin configurations are restricted either to  + or  − due to SSB. For keeping the discussion specific we take the spin configurations to be in  + giving the initial partition function f H S f  as usual. We will see that the restriction of the initial spin configurations has vital consequences. For large yet finite systems in the ordered phase, the free energy barrier is still high enough to make the transitions between +  and  − very rare, practically never happen at time scales sufficiently shorter than the ergodicity-breaking time (EBT) 34 . As long as the equilibrium state is concerned, ρ i can be safely confined within either +  or  − . Note that typical systems exhibiting phase transitions have the EBT much larger than the observation time. We also note that the equilibrium probabilities of magnetization per spin S are particularly useful.

Absolute Irreversibility in the Ising Model
We first illustrate that for the Ising model 〈e −σ 〉 = 1/2, regardless of how the parameter is tuned in time. We remark that the argument here is exact. Following the previous work 28 , we characterize the break-down of the Jarzynski equality as S where λ S is the probability of the absolutely irreversible paths, i.e., the total probability of backward paths for which the corresponding forward paths are absent. Let ±  be the subset of  that cannot be reached by any forward path Γ fwd from  ± . Obviously, −  is the sym- A denote the set of spin configurations that can be reached by some Γ fwd starting from either +  or −  or both. Figure 2 summarizes the relation among these sets. Now we note that, by definition, any backward path Γ bwd starting from  ± will end up necessarily in   while Γ bwd from  may arrive at either +  or  − . This can be expressed as bwd the probability for a backward path Γ bwd starting from the set  of spin configurations to reach the set . From the symmetry we see that Scientific RepoRts | 6:27603 | DOI: 10.1038/srep27603 It implies that the probability of backward paths whose corresponding forward paths (starting from +  ) do not exist is given by We thus conclude that 〈e −σ 〉 = 1/2 for arbitrary quench protocols varying the parameter ε(t).

Tolerance Parameter
Here we introduce a tolerance parameter to account for the insufficient sampling, which is inevitable in experiments and numerical simulations. In addition, a properly defined tolerance parameter enables to study the dynamical properties of phase transitions as we shall see in the remaining part of this paper. In the so-called "sudden" (infinitely fast) quenching (more general cases are discussed below), the system does not have enough time to change its distribution over spin configurations, and hence the initial equilibrium distribution is preserved throughout the whole process. The work distribution is thus completely determined by ρ i , leading to is over the initial distribution ρ i (S). Recall that the initial spin configurations are restricted to +  due to the SSB. The exponential average of the entropy production σ = β(W − ΔF) follows easily from 〈e −βW 〉 by multiplying by the exponential of the free energy change given by βΔF = −log(Z f /Z i ). In realistic experiments and numerical simulations, spin configurations with exponentially small probability do not take actual effects. Therefore it is natural to ignore such spin configurations up to certain tolerance δ. (Also practically, configurations with very low probabilities that cannot be sampled sufficiently, have to be discarded in the actual analysis of experiments and numerical simulations.) Specifically, for a given probability distribution ρ and tolerance δ, we implicitly define the set of kept spin configurations  δ and the cutoff probability ρ cut by the following two conditions (see also Fig. 3): We introduce the short-hand notations δ i f /  for the initial/final configurations ρ i/f (S). The corresponding partition functions are given by The free energy change is not affected by tolerance, We thus obtain  (14) is one of our main results and manifests several features to be stressed 29,30 : (i) As illustrated schematically in Fig. 3, 〈e −σ 〉 δ depends crucially on the overlap of  δ i and  δ f . For finite δ, well separated initial and final distributions lead to vanishing 〈e −σ 〉 δ . For δ = 0, on the other hand,  , and hence 〈e −σ 〉 = 1/2 validating the heuristic analysis presented in Fig. 1 (14) also demonstrates how the convergence issue arises in quenching process of phase transitions. Namely, the dominant contributions to the ensemble average of e −σ comes from the spin configurations with larger ρ f (S) whereas the initial equilibrium is governed by those with larger ρ i (S). (iii) Equation (14) describes highly nonequilibrium processes merely in terms of equilibrium distributions, a remarkably simple way to study 〈e −σ 〉 δ .

. (ii) Equation
The tolerance scheme (12) in terms of the microscopic spin configurations S is still difficult to implement in practice. For example, the tolerance parameter δ corresponding to the actual finite sampling is unknown or very difficult to estimate in most cases. For this reason, we introduce another operational tolerance scheme in terms of the macroscopic order parameter S: Given P μ (S) (μ = i, f), we implicitly define the interval of relevant magnetization  µ δ and the cutoff P cut by cut Note that the relation (14) does not depend on a particular tolerance scheme, and for the scheme (15) it reads as Below we will mainly use the tolerance scheme (15); the relation between the two tolerance scheme is discussed in the Supplementary Information.

Universal Scaling Behavior
We now examine 〈e −σ 〉 δ in Eq. (16) more closely with numerical simulations and analytical arguments, focusing on its scaling behavior inherited from the spontaneous symmetry breaking. According to Eq. (16), in the case of sudden quench it is enough to calculate the equilibrium distributions P i (S) and P f (S); no need to simulate the quenching dynamics.
To calculate the probability distribution of magnetization per spin, P(S), we perform a Monte Carlo simulation. Briefly, we (i) randomly generate initial spin configurations; (ii) choose one spin for flip at random and calculate the energies βH and βH′ in Eq.   Figure 4 displays the results of such a calculation and we can immediately see that for both 2 and 3 dimensional lattices 〈e −σ 〉 δ decreases as ε 0 is increased (at fixed L) or as L is increased (at fixed ε 0 ). Most remarkably, for a given tolerance δ, we find that 〈e −σ 〉 δ does not depend on ε 0 and L separately, but unexpectedly it is a universal function of ε ν L 0 with an exponent ν as shown in Fig. 5(a,b). Moreover we find that the exponent ν is nothing but the scaling exponent of the correlation length ξ ε  34 , Δ μ is related to the equilibrium susceptibility χ μ by χ ∆ = µ µ L / d . Note that beyond the above assumption to characterize the distributions using their peaks and widths, we do not assume any specific functional form of the distributions. Notably, the distribution can be highly non-Gaussian for small systems 37 and/or for systems with a continuous symmetry broken spontaneously 38 .
One can understand the arguments more rigorously along the lines of the large deviation theory 16,37,[39][40][41][42][43][44] : For sufficiently large systems (N = L d → ∞), the probability distributions behave as N is the so-called rate function. Thermodynamically, Φ μ (S) can be expressed in terms of the Helmholtz free energy density F μ (B) and the Gibbs free energy density G μ (S) by 37,34 where B is the external magnetic field.
2 which are consistent with the above arguments. This approximation becomes accurate in the thermodynamic limit and explored further below [Eqs (28)(29)(30)(31)(32)] to get an analytical universal expression of 〈e −σ 〉 δ in the limit. The universal tails of the work distribution near a quantum critical point has been revealed based on the large deviation theory 16 . Hence, investigating the overlap between P i (S) and P f (S) (and 〈e −σ 〉 δ in turn) more closely by means of the large deviation theory will be an interesting issue on its own, and we leave it open for future work. Instead, here we characterize the overlap by introducing the relative separation of the peaks of P i (S) and P f (S), The magnetization and susceptibility satisfy the standard scaling behaviors:  (21) and (22) together with the Rushbrooke scaling law 34 , α + 2β + γ = 2, one has the relative separation Using the Josephson hyperscaling law 34 does not depend on ε 0 and L separately but is a universal function of only the combination ε ν L 0 . This implies that 〈e −σ 〉 δ is also a universal function of ε ν L 0 alone, which is indeed confirmed by the numerical results shown in Fig. 5(a,b). Note that the hyperscaling law breaks down either in dimensions higher than the upper critical dimension d * = 4 or in the mean-field approximation. In such cases, where α = 0 and ν = 1/2, 〈e −σ 〉 δ is not necessarily a universal function of ε ν L 0 in general. For sufficiently large systems ε ν  L ( 1 ) 0 , one expects sharp distribution functions. Indeed, in this limit it follows that according to the hyperscaling law. It means that P i/f (S) have significant overlap with each other for a finite-size system. With the universal scaling behaviors of relative separation R at hand, let us now investigate 〈e −σ 〉 δ . For a given tolerance δ, the two asymptotic behaviors in Eqs (25) and (33) imply little and significant overlap between P i/f (S), respectively, and hence that 〈e −σ 〉 δ ~ 0 in the limit of ε ν  L 1 0 , while 〈e −σ 〉 δ ~ 1/2 in the opposite limit of ε ν  L 1 0 using the intuition provided by Eq. (16) and Fig. 3.

Effect of Tolerance.
Ideal sampling (δ = 0) always gives 〈e −σ 〉 δ = 1/2. However, for ε ν L 0 fixed, large tolerance (δ ~ 1) naturally leads to 〈e −σ 〉 δ ~ 0, because it decreases the overlap between the initial and final distributions P i/f (S) in Fig. 3. In the parameter space of ε ν L 0 and δ, 〈e −σ 〉 δ ~ 1/2 in the limit ε δ → ν L ( , ) (0, 0) 0 , while 〈e −σ 〉 δ ~ 0 in the opposite limit of ε δ → ∞ ν L ( , ) ( , 1) 0 as depicted in Fig. 5(c,d). Evidently, a crossover of 〈e −σ 〉 occurs as a result of combined effects of finite size and tolerance. One can locate the crossover boundary δ δ ε = ν ⁎ The crossover boundaries are illustrated by the thick red lines in Fig. 5(c,d). Here we note that the universality of 〈e −σ 〉 δ as a function of ε ν L 0 is valid up to the critical scaling of phase transitions. In the critical region, any physical quantity has both regular and singular parts, and only the singular part exhibits the scaling behavior. Therefore, 〈e −σ 〉 δ can also in general have a regular part, which does not obey the universal dependence on ε ν L 0 . Indeed, for fixed δ, Fig. 5(a,b) show three slightly different curves for different values of L. In Fig. 5(c,d) the contours plot 〈e −σ 〉 δ averaged over different values of L.
Thermodynamic Limit. Figure 5 shows that 〈e −σ 〉 δ is suppressed exponentially for sufficiently large systems Therefore, in practice δ δ  ⁎ always and 〈e −σ 〉 δ tends to vanish in the thermodynamic limit. Using Eq. (16), one can examine the tendency more closely: Thus, the exponential of entropy production at finite tolerance exponentially vanishes in the thermodynamic limit.

Finite-Speed Quenching
So far we have examined the nonequilibrium dynamics of the entropy production for the sudden quenching of the reduced coupling. Here we discuss the more general case of a finite-speed quench. As we already illustrated at the Scientific RepoRts | 6:27603 | DOI: 10.1038/srep27603 beginning, for zero tolerance one still has 〈e −σ 〉 = 1/2 exactly. Now we examine the case with a finite tolerance. We will employ the Kibble-Zurek mechanism (KZM) and demonstrate that the nonequilibrium dynamics of 〈e −σ 〉 δ reflects the equilibrium scaling properties of the system even in this case. Recall that due to the convergence issue in large systems, the numerical simulation for the finite-speed quench is very demanding computationally. The analysis based on the KZM may provide an outlook.
We consider a linear quenching process where the coupling constant varies in time as ε(t) = ct for t ∈ (−∞, ∞) with the finite quenching speed c. We assume a finite tolerance δ, up to which relatively improbable spin configurations are ignored.
In general, the nonequilibrium dynamics caused by such a quench can be very complicated. Here we adopt the spirit of the KZM and make the adiabatic-impulse approximation. The KZM was originally put forward to study the cosmological phase transition of the early Universe 45,46 and later extended to study classical phase transitions in condensed matter systems 47,48 . Recently, it was also found to apply to the Landau-Zener transitions in two-level quantum systems 49,50 and the dynamics of second-order quantum phase transitions 51,52 . As a theory of the formation of topological defects in second-order phase transitions, the KZM establishes accurate connections between the equilibrium critical scalings and the nonequilibrium dynamics of symmetry breaking.
The key idea of the KZM is to classify the dynamics into two distinct regions, the adiabatic and impulse regimes. Recall that, in a second-order phase transition, the relaxation time τ provides the single time scale analogous to the correlation length ξ which is the only length scale of the system. The former is naturally related to the latter by where z is the dynamic critical exponent, and diverges at the critical point (ε = 0). In the early stage of the linear quenching process (t → −∞), the system is far from the critical point, and its relaxation time τ(t) ~ |ε(t)| −zν → 0 is short enough so that it can adjust itself to the change in the reduced coupling actuated by the driving. In this sense, the dynamics is adiabatic. On the contrary, as the system approaches the critical point (ε(t) → 0), the relaxation time τ diverges and the relaxation of the system is too slow to follow the external change. The system essentially remains the same and the external driving can be regarded as impulsive. In the far future (t → ∞), getting away from the critical point, the relaxation time decreases back and the dynamics becomes adiabatic again. The crossover between the two dynamical regimes occurs at t = ±t c when the relaxation time becomes comparable with the time scale, ε ε  / , for ε(t) to develop, namely, τ ε ε ± = t t ( ) / c c . It follows from (33) that t c ~ c −zν/(1+zν) and ε(t c ) ~ c 1/(1+zν) . The slower the quenching process, the closer the crossover point ε(t c ) comes to the critical point.
The adiabatic-impulse approximation drastically simplifies the picture of the nonequilibrium dynamics of the system. In the early stage from t = −∞ to −t c and the later stage from t c to ∞ of the quenching process, the dynamics is adiabatic (or quasi-static), and the work is equal to the free-energy change. It means that there is no entropy production in these periods. In contrast, the system is completely frozen during the interval −t c ≤ t ≤ t c . The dynamics during this interval is equivalent to the case of the sudden quenching discussed above. Therefore, 〈e −σ 〉 δ for the whole quenching process is determined solely by the dynamics during the impulse interval from −t c to t c . As a consequence, the same analysis for sudden quenching processes given above can be applied for the finite-speed quenching processes, but with ε 0 replaced by ε(t c ) ~ c 1/(1+zν) .

Conclusion
We have found that, near a critical point, for a given tolerance parameter δ the ensemble average of e −σ follows a scaling law in terms of a universal combination of the reduced coupling constant ε 0 of the initial state and the system size L: ε ξ = ν L L/ 0 0 with ν being the critical exponent of the correlation length ξ 0 . As noted previously [18][19][20][21][22][23][24][25][26][27][28] , the Jarzynski equality may break down for many practical and intrinsic reasons. Its breakdown in our case is peculiar as the deviation is determined by an universal combination of L and ε 0 , which is inherited from the equilibrium scaling behavior of second-order phase transitions. It is stressed that such a universal scaling behavior is not limited to the sudden quenching but holds in general due to the critical slowing down. Our findings may provide a unique application of the Jarzynski equality to study the dynamical properties of phase transitions.