Distribution of Kinks in an Ising Ferromagnet After Annealing and the Generalized Kibble-Zurek Mechanism

We consider the annealing dynamics of a one-dimensional Ising ferromagnet induced by a temperature quench in finite time. In the limit of slow cooling, the asymptotic two-point correlator is analytically found under Glauber dynamics, and the distribution of the number of kinks in the final state is shown to be consistent with a Poissonian distribution. The mean kink number, the variance, and the third centered moment take the same value and obey a universal power-law scaling with the quench time in which the temperature is varied. The universal power-law scaling of cumulants is corroborated by numerical simulations based on Glauber dynamics for moderate cooling times away from the asymptotic limit, when the kink-number distribution takes a binomial form. We analyze the relation of these results to physics beyond the Kibble-Zurek mechanism for critical dynamics, using the kink number distribution to assess adiabaticity and its breakdown. We consider linear, nonlinear, and exponential cooling schedules, among which the latter provides the most efficient shortcuts to cooling in a given quench time. The non-thermal behavior of the final state is established by considering the trace norm distance to a canonical Gibbs state.


I. INTRODUCTION
Nonequilibrium phenomena occupy a prominent role at the frontiers of physics, where few and highly-valuable paradigms are able to provide a description making use of equilibrium properties. Notable examples include linear response theory and the fluctuation-dissipation theorem [1], fluctuation theorems and work relations valid far from equilibrium [2,3], and the Kibble-Zurek mechanism [4][5][6][7]. We focus on the latter, as it provides a framework to analyze the course of a phase transition and the breakdown of adiabatic dynamics leading to the formation of topological defects. In this context, the system of interest exhibits different collective phases as a control parameter is varied across a critical value. This parameter is the temperature in thermal phase transitions but can be identified by other quantities such as a magnetic field, or the density of particles in the system. The crossing of a continuous phase transition is characterized by the divergence of the (equilibrium) relaxation time in the neighborhood of the critical point, known as critical slowing down. As a result, whenever the phase transition is driven in a finite quench time τ Q , adiabaticity is broken [8].
A scenario of spontaneous symmetry breaking is characterized by the presence of a manifold of degenerate ground states in the low-symmetry phase of the system. During the course of the phase transition, causally disconnected regions of the system may single out different ground states, leading to the formation of domains and the creation of topological defects at the resulting interfaces. A familiar example in this context is the cooling of a paramagnet below the Curie temperature, resulting in domains with a homogenous local magnetization and separated by domain walls. According to the Kibble-Zurek mechanism (KZM) the mean number of defects decays as a function of the time scale in which the transition is crossed. Specifically, a universal power-law scaling is predicted when the control parameter is driven linearly in time. The finite-time cooling of an Ising ferromagnet has been used as a paradigmatic testbed to explore KZM physics [9][10][11][12][13], which provides useful heuristics in adiabatic quantum optimization and quantum annealing [14][15][16][17][18]. Generalizations of KZM have been established that account for disorder [19,20], nonlinear driving protocols [21][22][23] as well as inhomogeneous systems [23][24][25][26][27][28][29][30][31][32][33], see Refs. [8,34] for a review. These developments have inspired novel protocols in adiabatic quantum computation [25,[35][36][37][38][39]. While the early formulation of KZM was focused on classical systems, following decades of research, the applicability of the KZM in the quantum domain has been established by a combination of analytical, numerical, and experimental studies [8,27,34].
Beyond the mean number of kinks, one may wonder whether the full kink number distribution exhibits universal behavior. The latter is directly accessible in many experiments and can be as well probed via single-qubit interferometry [40]. The kink number distribution has recently been shown to exhibit signatures of universality beyond the KZM in a family of models known as quasi-free fermion systems, that include paradigmatic in-stances such as the one-dimensional transverse-field Ising and XY models, and the Kitaev chain [41]. In particular, not only the mean number of defects but as well the variance, third centered-moment, and any cumulant of the kink number distribution of higher order have been shown to scale following a universal power-law with the quench time [41][42][43]. This prediction has been experimentally explored using a trapped-ion for the quantum simulation of critical dynamics in momentum space [43]. Universal features of kink number statistics in the onedimensional transverse-field quantum Ising model have also been reported using D-Wave quantum annealers as quantum simulators [44]. It is thus natural to wonder whether the distribution of topological defects in classical systems is as well universal. Indeed, a framework to account for the distribution of topological defects generated across a classical continuous phase transition has been put forward and predicts a binomial distribution, in agreement with numerical simulations for the timedependent Ginzburg-Landau theory [45]. In higher dimensional systems, further evidence for the presence of universality in the full counting statistics of topological defects has been provided by the study of the vortex number distribution in a newborn holographic superconductor, which was predicted to be Poissonian [46]. At the time of writing, experimental evidence of Poissonian vortex statistics has been reported after cooling of an atomic Bose gas into a superfluid in finite time [47].
In this work, we characterize the exact kink-number distribution of a one-dimensional classical ferromagnet cooled in finite time. Specifically, we consider the onedimensional Ising model with no magnetic field and evolving under Glauber dynamics. The mean number of defects in this setting has been studied by Krapivsky [10], see as well [12,13] for related work.
Here, using a ring topology endowed with translational invariance, the kink number distribution is studied. We calculate explicit expressions for the general two-point correlator for spins separated by a lattice-point distance n in the same limit. The first three cumulants of the kink number distribution are explicitly shown to be equal and described by a universal power-law with the quench time, indicating that the slow cooling of an Ising ferromagnet under the Glauber dynamics yields a Poissonian kink-number distribution. The relevance of these findings to finite-annealing times is verified by numerical simulations, in which we consider three different families of cooling schedules: linear, nonlinear, and exponential quenches.
From the outset, we note that the only critical behavior of a one-dimensional Ising ferromagnet is exhibited at zero temperature. A cooling schedule cannot possibly involve the crossing of the critical point considered in the original studies of KZM. However, the KZM prediction can be extended to account for "half-quenches" ending at the critical point [48][49][50][51][52]. At the same time, the finitetemperature treatment endows the dynamic with coarsening. The nonequilibrium dynamics is thus governed by a coexistence of KZM universality and coarsening. Their contribution can generally be discriminated by considering the time scales involved. In some instances, such as the artificial spin ice [53], the discrimination is not possible. However, for the one-dimensional Ising ferromagnet with non-conserved order-parameter (magnetization), domain growth due to coarsening scales with the square root of the time of evolution. In our study, we can uniquely identify signatures of critical scaling for different kinds of quenches (linear, nonlinear, exponential, etc.), ruling out the effects of coarsening.
The equilibrium critical behavior of the one dimensional Ising ferromagnet is peculiar in that it does not exhibit a power-law divergence of the correlation length as a function of the proximity to the critical point, known in higher dimensional continuous phase transitions. As a result, the correlation length critical exponent ν is not well-defined. However, we shall see that critical scaling with the finite driving time governs the cumulants of the kink distribution as it does in the generalized KZM.

II. KINKS DISTRIBUTION IN AN ISING FERROMAGNET
A one dimensional Ising ferromagnet in the absence of an external magnetic field is described by the Hamiltonian H = − j J ij σ i σ j where the spin at a site j can take any of the two values σ j = ±1. The ferromagnetic character stems from J ij ≥ 0. We first discuss the distribution function of the number of kinks and its characteristic function. Given a one-dimensional spin chain, the number of kinks in a given configuration can be studied via the number operator which can take integer value k ∈ [0, N ]. We assume periodic boundary conditions (in the case of an open chain, the upper limit of the sum is N − 1 instead of N and k ∈ [0, N − 1]). We shall be interested in the distribution of the number of kinks where denotes the state of the system. For its characterization, we shall resort to the characteristic function As the kink number takes integer values, using the Fourier transform yields The kink distribution is accessible in experiments and can be measured, e.g., via single-qubit interferometry [40].
We shall focus on the distribution of kinks in the nonequilibrium state resulting from cooling the Ising ferromagnet in a finite time. For its analysis, it will prove useful to use the cumulants κ j (j ∈ N) of the distribution. The cumulant generating function of the kink number distribution is the logarithm of the characteristic function and admits the expansion In particular, we shall focus on the mean given by κ 1 = N , the variance κ 2 = N 2 − N 2 and the thirdcentered moment κ 3 = (N − κ 1 ) 3 .

III. COOLING BY GLAUBER DYNAMICS OF AN ISING FERROMAGNET: EXACT SOLUTION
To describe the finite-time cooling of the Ising ferromagnet we shall consider its evolution under Glauber dynamics [54,55]. Specifically, we consider the nonequilibrium quenching process in which the evolution of the Ising chain is described as a reversible Markov process obeying the detailed balance condition P eq ( σ)w i ( σ) = P eq ( σ (i) )w i ( σ (i) ), where σ = (σ 1 , ..., σ i , ..., σ N ) is the current system state and σ (i) = (σ 1 , . . . , −σ i , . . . , σ N ) denotes the same state with the ith spin flipped. The probabilities P eq given by the Boltzmann distribution where the partition function is given by Z = {σi=±1} e −βH( σ) . The flipping rate w i ( σ) of spin i is obtained by direct substitution of Eq. (6) into the detailed balance condition: = e −βσi j Jij σj e βσi j Jij σj (8) where the final equality is obtained by substitution of the hyperbolic identity for the exponential. We shall focus on the uniform coupling case with nearest neighbor interactions, i.e., J ij = Jδ i,i+1 . Equation (7) then implies the most general flipping rate in this case to be where we have normalized the rates by setting the limit w i → α/2 as T → ∞ and defined γ = tanh(2βJ) for convenience. The Glauber dynamics of Ising ferromagnets have been the subject of an extensive literature. We focus on the non-equilibrium case, in which both γ = γ(t) and α = α(t) act as control parameters of the temperature and local flipping barrier, respectively. Making use of translational invariance, consider the correlator between nearest-neighbor spins W 1 = σ i σ i+1 . Previous work by Krapivsky indicates that in the slow-cooling regime this correlator takes the form [10].
(up to a logarithmic factor in τ Q ) where C is a constant dependent on the cooling schedule specifics, τ Q is the time taken in total for the temperature to pass from an effectively infinite value to T = 0 and the power-law exponent δ is set by the dynamic critical exponent z, the cooling schedule, and system dimensionality; see as well [12,13].
The study of kink statistics requires the calculation of moments N m of the kink operatorN , each of which involves the evaluation of expressions proportional to correlators up to and including 2m individual spins. In the classical Ising model, it is known that even correlators may also be decomposed into an alternating sum of twopoint correlators, parameterized solely by their separation n and their time dependence under translational invariance [55]. In the limit of slow cooling as W 1 → 1, it is expected that the correlator W n of particles separated by a distance n also converges with the same power law in τ Q . In this limit, the dependence of correlators W n should linearize their functional dependence on the distance n. This motivates us to start considering the behavior of two-point functions of the form In the sections to come, we derive explicit asymptotic expressions for W n in the case of linear, algebraic, and exponential quenches from the generating function. We shall use these results to establish the universal form of the kink number distribution and the scaling of its cumulants with the quench time.
We also perform extensive dynamical simulations of the ferromagnetic Ising-Glauber model with various cooling schedules. The Glauber dynamics, which is equivalent to the so-called heat-bath method in Monte Carlo simulations, can be easily implemented numerically. To take into account the stochastic and local nature of the spin dynamics, at each fundamental step, a spin σ i that is randomly chosen from the system is to be updated according to the Glauber transition dynamics. Specifically, a random number r uniformly distributed in the interval [0, 1] is generated from a pseudo-random number generator. The chosen spin is flipped, i.e. σ i → −σ i if this random number satisfies r < w i ( σ), where w i ( σ) is given by the Glauber acceptance rate Eq. (10) with α set to 1. To properly compare simulation results from different system sizes N , we define a time-step in our simulations as consisting of N single spin-updates described above. The system is initialized in a random spin configuration and then cooled down by tuning the control parameter γ at each time step according to the cooling schedule. For each cooling speed τ Q , a large number of independent cooling simulations are performed and observables are computed from instantaneous snapshots of the spin configurations.
Periodic boundary conditions are used in all our numerical simulations presented below. When a final configuration is generated, the kink numberN is measured in different configurations, obtaining the the moments N , N 2 , N 3 from the Monte Carlo average. The cumulants are then calculated using the identities where · · · denotes average over independent annealing simulations.

A. Cumulant Generating Function
Consider the probability distribution for a given spin configuration P ( σ, t). Under Glauber dynamics, this distribution evolves according to the master equation where · σ denotes an expectation over spin realizations σ, and {η i } is a set of Grassmann variables satisfying and we assume an infinite chain for simplicity. Explicit correlation functions can be derived from the generating function via the identity, applicable for an even number n of indices i j : After the propagation of the generating function by way of (15) in the manner described by Aliev [58], the form of Ψ({η i }; t) induced by the Grassmann variables {η i } allows the explicit expression by differentiation of correlators in the form of an alternating sum of products of two-point functions. More concisely, this may be encoded as where W i1,i2,...,in is an antisymmetric 2n × 2n matrix whose elements are defined in terms of the two-point cor- In Eq.
The Pfaffian equivalence induces an equivalent structure to the Wick contraction for fermionic field operators. A power-series expansion of the kink number characteristic function involves these correlators Thus, we can formally write the characteristic function in terms of the generating function Ψ = Ψ({η i }; t) (23) Its logarithm, lnP (θ), is the cumulant generating function. Specifically, the j-th cumulant κ j of the kink number distribution can be found as Let us consider three first terms with j = 1, 2, 3. One readily finds the mean number of kinks as the first cumulant, i.e, Similarly, the correlator between two spins that are n sites apart in the presence of translational invariance will be denoted by W n = 1 N i σ i σ i+n . The explicit computation of higher-order cumulants is somewhat laborious. Here, we simply quote the result for the second and third cumulant derived in the appendix A. The second cumulant equals the variance of the number of kinks and reads The third cumulant equals the third centered moment and its explicit computation yields At this stage, we can analyze the general features of the kink distribution in the binomial model. The latter is associated with N Bernoulli trials describing the presence of a kink at the interface between different spins with a success probability p. The first three cumulants of the binomial distribution are given by N p, N p(1 − p) and N p(1 − p)(1 − 2p). From the expression of the mean, one can identify the probability p in terms of the spin-spin correlator as p = 1 2 (1 − W 1 ). An analogous identification holds for κ 2 , with (26) having first term κ 2 = 1 4 (1 − W 2 1 ). Regarding the κ 3 the first two terms in (27) are consistent with the binomial expression as . We shall revisit the connection with the binomial distribution in a different framework, that of the generalized KZM, in Section X.
To summarize this section, we have obtained exact expressions -Eqs. (25), (26) and (27) -for the first three cumulants of the kink number distribution in terms of the n-site correlator W n . A crucial observation is that these equations, together with the ansatz for the leading power-law behavior of W n in Eq. (13), yield, to leading order in 1/τ Q suggesting that the kink-statistics becomes Poissonian in this limit. We next turn our attention to the explicit analysis for specific cooling protocols.

IV. FINITE TIME COOLING
We consider an infinite ring, with an uncorrelated initial state corresponding to the high-symmetry phase satisfying σ i 0 = σ i σ j 0 = · · · = 0 and no local flipping barrier (thus α(t) = α 0 = 1). In this case, the generating function reads [56,58] (29) in terms of [58] and I ν (x) denotes the νth modified Bessel function of the first kind. The indices f i denote the (ordered) index of each spin. In the case of an infinite chain, the correlators depend only on the distance n between successive spins. In addition, the successive modified Bessel functions may reduce the expression via the identity Thus, evaluating at the instant where T = 0, we may write which provides the exact integral representation of the n-site correlator. Using this result, together with those derived in the preceding section, we next describe the kink statistics resulting from different cooling schedules. Specifically, under Glauber dynamics the flipping rate is dictated by the parameter γ = tanh(2βJ) and we shall consider different functional forms for the variation of this parameter in time [10][11][12][13].

A. Linear Quench
For a linear cooling schedule, we consider where τ Q denotes the total time taken to cross from the initial condition to the T = 0 state. After evaluation at where we have defined η = (τ 2 Q − t 2 )/τ Q for convenience, and used the approximation to the exponential factor justified by the exponential rolloff of the contribution for η = O( √ τ Q ). Asymptotic solution of the integral (34) can be exactly computed as shown in Appendix B and yields The above result in conjunction with our results for κ 1 in Eq. (25), yields By comparing the amplitude of the leading and subleading terms, one concludes that the power-law behavior sets in for quench times where the superindex indicates that this time scale characterizes the first cumulant.
In agreement with Eq. (28), to leading order in a 1/τ Q expansion, we further find which suggests that the distribution becomes Poissonian distribution in the limit of arbitrarily slow cooling. Fig. 1 shows the probability distribution functions of kink number P (n) obtained from the Glauber dynamics simulations for finite quench times. It is worth noting that the distribution of kinks is well described by binomial distribution B(n, p) with parameters n = κ 1 p, p = 1 − κ 2 /κ 1 . Fig. 2 shows the scaling of the first three cumulants as a function of the annealing time. Numerically, the cumulants are calculated by averaging the recorded kink numbers using Eqs. (14). Our numerical results clearly show that the first two cumulants follow a power-law dependence with the annealing rate τ Q . A relatively larger fluctuation can be seen in the data points of the third cumulant κ 3 , which is expected due to the enhanced statistical error in the numerical calculation of higher-order moments. Nonetheless, the trend of κ 3 still roughly follows the power law. The power-law exponents obtained from nonlinear least-squares fitting are 0.239 ± 0.001, 0.229 ± 0.001, and 0.189 ± 0.007 for the first three cumulants. These values are close to the theoretically predicted value 1/4 in Eq. (39) but exhibit some deviations from it. Thus, only the first cumulant is governed by the leading 1/τ 1/4 Q term in this range, while the subleading corrections are important for κ 2 and κ 3 .

B. Nonlinear Algebraic Quench
The nonlinear passage across a critical point has been proposed to suppress the mean number of defects generated in a phase transition, as it yields a power-law dependence on the quench time with a tunable exponent [21][22][23]. This feature is also found in the finite-time cooling of an Ising ferromagnet under Glauber dynamics [10,12,13] and we next study its effect on the distribution of kinks and the cumulant scaling. To this end, we consider the algebraic cooling schedule parameterized by α. The analogous form of equation (34) is thus where now Again asymptotically approximating (41) in the limit of large τ Q we find the general expression where 1+α . The expression for κ 1 is thus Comparing again leading and subleading amplitudes, we find A straightforward exercise verifies that (43), (44) and (45) coincide with the respective linear schedule expressions (36), (37) and (38) in the special case α = A = 1.
To leading order in 1/τ Q , where 1+α . Thus, the expressions analogous to (39) are, keeping only the leading order in 1/τ Q , Cumulants of the kink distribution thus exhibit a powerlaw scaling with the quench time. In particular, the power-law exponent α 2(1+α) increases within its range [0, 1/2] as the parameter α of the cooling protocol is increased.
Minimizing κ j (j = 1, 2, 3) with respect to α we find that the optimal value of α which yields the minimum value of the cumulants (49) Comparing the mean number of kinks resulting from this optimized nonlinear schedule and the linear case in Eq. (39), we find that the latter leads to an enhanced suppression by a factor Said differently, for a given cooling time τ Q it is possible to reduce the mean number of kinks with respect to the linear schedule by using an algebraic schedule. This finding is reminiscent of the suppression of the mean number of excitations in the quantum dynamics of isolated critical systems, in which the dynamics is unitary and thus preserves entropy along the evolution [21][22][23]. Here, we further note that the same conclusion applies to higherorder cumulants. However, as discussed in Sec. VI, a sudden quench outperforms these schedules. Fig. 3 shows the distribution function of kink-number P (n) obtained from Glauber dynamics simulations for three different annealing rates. Again, the kink-statistics is well described by binomial distribution parametrized

C. Exponential Quench
Both linear and algebraic cooling schedules lead to a power-law scaling of the cumulants of the kink distribution. We next consider an exponential quench in the form suggested by Krapivsky [10] 1 Here, b, β > 0 are positive real coefficients and B = exp(b) is a normalization factor ensuring γ(0) = 0 and γ(τ Q ) = 1. Making the substitution η = 2h(τ Q , τ ), we find the rather cumbersome integral expression where Γ(a, b) denotes the (upper) incomplete gamma function. Asymptotic solution of the integral (52) then leads to the result for W n giving   The kink-statistics is also well described by binomial distribution parametrized by n = κ 1 p, p = 1 − κ 2 /κ 1 . Fig. 6 shows the first three cumulants as functions of the annealing rate using parameters β = 2, b = 1 and B = 1 for the exponential cooling schedule. Consistent with the analytical prediction Eq. (54), we find all three cumulants exhibit a power-law dependence as a function of log(τ Q ). However, contrary to a constant exponent 1 2β = 0.25, our best nonlinear least-square fit gives three different values 0.245 ± 0.001, 0.241 ± 0.002, 0.114 ± 0.040 for the exponents. In particular, the third cumulant is substantially different from the theoretical prediction. For β = 2, the log(τ Q ) 1/2β term is changing very slowly in the τ Q range of the simulation, making the cumulants only weakly dependent on log(τ Q ); thus, the τ Q dependence of cumulants is dominated by the 1/ √ x term in this τ Q range.
This makes the higher-order cumulant fitting in this cooling schedule more sensitive to numerical uncertainties. A remarkable feature of the exponential schedule is that its cooling efficiency surpasses that of the optimized nonlinear schedule. Indeed, taking the ratio of (49) over (54) we find that is, the exponential schedule leads to a logarithmic suppression of the mean kink density with the quench time over the optimized nonlinear schedule.

V. THERMAL EQUILIBRIUM AT ARBITRARY TEMPERATURE
We briefly consider the equilibrium kink number distribution that will be relevant to the following sections devoted to sudden quenches and non-thermal behavior. We consider an arbitrary inverse temperature β ≥ 0. In this case, it is known that the k-point correlator takes the form [56,59] with It follows that at equilibrium two-point correlator at distance n equals Using the expressions (25), (26) and (27) for κ j (j = 1, 2, 3), one obtains We note that these expressions are equivalent to those of the binomial distribution B(N, p) with the kink formation probability In the infinite temperature case, the distribution describes as well that of the quantum Ising chain [60].

VI. FAST AND SUDDEN QUENCHES
We next consider the behavior of the system under a rapid quench and note that each cooling schedule yields in this limit Taking the series expansion of I n at η → 0, we find that The integral in (67) is equal to the lower incomplete gamma function γ(a, b), which gives a Taylor series, the leading factor of which is of the form b a Γ(a)e −b . Taking the leading factors of both summations, observing that e −τ Q ≈ 1 − τ Q for fast quenches and taking the minimal power in τ Q yields Substituting (68) into the expressions for κ 1 , κ 2 and κ 3 , while taking the leading powers gives The kink distribution upon completion of the quench in the limit of vanishing τ Q is that of a ferromagnet at infinite temperature. According to Eq. (65), for β = 0 the kink formation probability is p = 1/2, as expected. As a result, the cumulant values of a binomial distribution B(N, 1/2) in Eq. (62) are recovered, i.e., κ 1 = N 2 , κ 2 = N 4 , κ 3 = 0. We further notice that the values in Eqs. (69) also agree with those of the binomial distribution in (62) when the kink formation probability reads which captures the leading correction away from the sudden limit due to the finite value of the quench time τ Q . We note in passing that another interesting dynamical phenomenon related to sudden quench is domain coarsening [61]. It is generally believed that coarsening systems exhibit dynamical scaling, i.e., the typical domain size grows algebraically with time L ∼ t 1/z , where z is a dynamical exponent that is independent of microscopic details of the system. In the 1D Ising chain, the typical domain size is simply related to the average distance between kinks, hence L ∼ 1/κ 1 . The scaling hypothesis thus implies a power-law behavior for the first cumulant. Interestingly, our extensive Glauber dynamics simulations show that all three cumulants follow a diffusive scaling law: κ j ∼ t −1/2 , corresponding to an exponent z = 2; see Appendix C for more details. As a caveat, it should be noted that the phenomenology of sudden thermal quenches differs from that of sudden quenches in quantum phase transitions in isolated spin chains. Indeed, the sudden quench followed by an evolution time leads to a lower density of defects than the linear, nonlinear, and exponential schedules for a given total duration of the process.

VII. NON-THERMAL BEHAVIOR
One may wonder whether the non-equilibrium state resulting from the finite-time cooling of a ferromagnet is effectively thermal. To that end, one can compute the distance between an equilibrium thermal distribution of kinks P β (n) with inverse temperature β as a free parameter and the numerically obtained distribution P (n) for given P (n) = P τ Q (n). The proximity between the two distributions can be quantified by a distance. We consider the trace-norm distance Minimizing it with respect to the free parameter β, one can identify the effective temperature β * that best approximates the non-equilibrium state with distance D * TN . The equilibrium distribution P β (n) is obtained from standard Monte Carlo simulation using Glauber dynamic spin update and the trace-norm distance is minimized using golden search method. For concreteness, we focus on the case of a linear quench protocol. The numerical simulations for a chain of L = 500 spins indicate non-thermal behavior in the final state for quench times τ Q ∈ [10 2 , 10 4 ], see Fig. 7. For these parameters, the minimum trace norm distance remains in the interval D * TN ∈ [0.07, 0.09]. The canonical Gibbs state that best approximates the final state is characterized by an inverse temperature that scales as a power law of the quench time with exponent −0.096±0.001. We note that in a Gibbs state, the mean kink number equals [40] Assuming D TN = 0, and comparing this expression with the result for a linear quench (39) suggests that the power-law scaling is effective and results from linearizing the logarithmic dependence (taking the Boltzmann constant k B = 1) over the studied range of quench times. Naturally, the explicit dependence of T * on τ Q varies with the cooling schedule.

VIII. LIMIT OF NUMERICAL SIMULATION
The deviations between the analytical results and the numerical data observed in the histograms and cumulant scaling behavior come from the fact that the slow cooling limit τ Q → ∞ and the thermodynamic limit of infinite system size, both considered in the analytical approach, are not accessible in the numerical simulation. In the low temperature, long-time limit, the topological defects are exponentially scarce, making finite-size effects significant in the slow cooling regime. An arbitrarily long cooling time will simply bring a finite system close to equilibrium and the non-equilibrium physics can not be fully captured. Therefore, to explore the non-equilibrium physics in the slowing cooling limit, the infinite size limit is also required, which is beyond reach in the numerical simulation. The asymptotic behavior of the system in slow cooling limit can still be analyzed by observing the scaling of the correlator W n , which has a stronger dependence on the cooling rate. It can be seen that W n converges to the predicted expression 1 − nC τ −δ Q in the long time limit. This dependence will further lead to the cumulant behavior predicted by the analytical calculation. Note that in the range of quench times τ Q ∈ [10 2 , 10 4 ], the system size does not have a significant impact on the results.

IX. CONNECTION TO THE KIBBLE-ZUREK MECHANISM: MEAN NUMBER OF KINKS
KZM predicts the mean density of defects upon completion of a cooling schedule making use of equilibrium properties [8]. We first recall the equilibrium correlation length of the one-dimensional Ising model is [62] where ξ 0 is the lattice spacing. The one-dimensional ferromagnet thus differs from the standard setting in higher dimensions, where correlation length exhibits a powerlaw scaling as a function of the proximity to the critical point. By contrast, the relaxation time under Glauber dynamics exhibits the conventional power-law divergence [10,12] where τ 0 is a microscopic constant. As the critical point at T = 0, the system exhibits critical slowing down and the dynamics can be expected to be nonadiabatic for any finite quench time. For the sake of illustration, we focus on the linear cooling schedule KZM invokes the adiabatic impulse approximation, according to which the relaxation time is in an early stage small enough so that the system quickly adjusts to the instantaneous equilibrium configuration with γ = γ(t).
The growth of the relaxation time close to the critical point gives rise to the effective freezing of the order parameter of the system. KZM estimates the mean size of the domains (out-of-equilibrium correlation length) after cooling in finite time by the equilibrium value of the instantaneous correlation length at freezing, the so-called freeze-out timet. To estimate the freeze-out timet, we match the instantaneous equilibrium relaxation time to the time left until reaching the critical point τ Q − t, that is, For the linear schedule, the solution is given bŷ By an analogous procedure, one can estimate the freezeout-time for other schedules such as the algebraic and the exponential one. Using the relation between the correlation length and the relaxation time KZM predicts the mean domain size after cooling to be given byξ which for z = 2 yields the power-law scaling The accuracy of the KZM in accounting for the finitetime cooling of the Glauber dynamics has been discussed in [10,12]. We next focus on physics beyond KZM associated with the kink number statistics.

X. BEYOND THE KIBBLE-ZUREK MECHANISM: KINK NUMBER STATISTICS
A growing body of results [41][42][43][44][45][46] suggests that the signatures of universality govern the kink number distribution and not only its mean value. To appreciate this, in a classical setting, it suffices to assume that the formation of kinks at different locations is described by independent stochastic events [45].
The key tenet of KZM is that the cooling dynamics sets the average length scale of domains to be given by the equilibrium correlation length evaluated at the freezeout timeξ. By contrast, to generalize KZM we consider that the effect of the cooling is to partition a system of size L = N ξ 0 into "proto-domains" of the same length scaleξ over which the order parameter stabilizes. At the boundary between adjacent domains, kinks form with a given probability p. Conversely, with probability (1 − p) no kink is formed and the two adjacent proto-domains coalesce to form a larger domain. The number of boundaries between proto-domains determines the number of stochastic events for kink formation set by (the floor of) where the second-equality holds for the Ising ferromagnet. Assuming kink formation events at different locations to be uncorrelated leads to a kink number distribution associated with N b independent and discrete random Bernoulli variables. Upon assuming the success probability p to be the same at different locations, the distribution takes the binomial form see Figure 9. In one spatial dimension, Similarly, higher order cumulants of the binomial distribution read As a result, all cumulants are predicted to follow the same universal-power-law scaling predicted by the mean number of kinks, in agreement with the numerical simulations and analytical calculations we have reported. In the binomial distribution, cumulant ratios are determined by the kink formation probability p, independent of the quench time. Accordingly, the numerical values of the cumulant ratios found in Fig. 2 are expected to be consistent with a well-defined probability for kink formation, and thus, with the binomial distribution. According to the generalized KZM this probability is independent of the quench time, as shown in Figure 10. From κ 2 /κ 1 = 1 − p, one finds p = 0.32. Using any of the ratios involving the third cumulant, κ 3 /κ 2 , one finds the close value p = 0.27, though we recall that the power-law scaling of κ 3 deviates from the KZM prediction. These values are comparable to those observed in other one-dimensional systems, such as the overdamped Ginzburg-Landau model (p ≈ 0.42) and the transverse-field quantum Ising model, well in isolation (p = 0.41) [41][42][43], or coupled to a bath (p ≈ 0.37 − 0.39) [44],

XI. DISCUSSION AND CONCLUSION
We have analyzed the kink number distribution of an Ising chain thermally annealed under Glauber dynamics. While generally it is well described by a binomial distribution, in the limit of slow annealing kink statistics becomes Poissonian. We have explicitly computed the two-point function and used it to derive the low-order cumulants of the distribution. Specifically, the mean number of kinks, the variance, and the third centered moment are identical and given by a power-law with the quench time in the limit of slow cooling.
The one-dimensional Ising model does not exhibit a phase transition and in the absence of a magnetic field becomes degenerate only at zero temperature. The annealing schedules we have reported involve positive temperatures, approaching only degeneracy at infinite time. As a result, the annealing of the ferromagnet does not involve the crossing of the critical point. The situation is similar to that in recent experiments with colloidal monolayers that probe only "half of the transition" [50,51]. In principle, such a scenario does not preclude the appearance of KZM scaling [48,49,52], although under slow cooling, the dynamics is inextricably woven with coarsening [10,12,13,63,64].
The Ising ferromagnet in one spatial dimension does not exhibit a power-law divergence of the equilibrium correlation length. As a result, the correlation length critical exponent ν is not defined. This precludes the application of the KZM in its original form [8]. However, the appearance of power-law behavior can be established using the adiabatic-impulse approximation [10,12,13], a core tenet of KZM. Focusing on kink number fluctuations, we have characterized the full kink number distribution that exhibits signatures of universality as predicted by the generalized KZM [45]. The dependence of the cumulants of the kink number distribution on the annealing time varies with the schedule. When the temperature is a linear function of time, all cumulants are shown to scale with a power-law of the quench time. When the annealing schedule involves a polynomial variation of the temperature with time, a modified power law is observed. This generalizes to arbitrary cumulants the scaling prediction for the mean number of defects resulting from the nonlinear passage across a critical point [21,22,45,65]. We have found corrections to the power-law behavior of cumulants when the temperature decays exponentially as a function of time, see as well [10,12,13] for the mean number. At variance with previous studies exploring Berezinskii-Kosterlitz-Thouless phase transition [66,67] and holographic systems [68], the logarithmic corrections to KZM scaling that we have reported stem directly from the annealing schedule.
We have further analyzed the dependence of the cooling efficiency for a given quench time as a function of the choice of the cooling schedule. Nonlinear quenches are shown to reduce the residual density of kinks below the value obtained under a linear quench. This result is consistent with previous findings on nonlinear quenches [10,12,13,21,22,45,65]. The nonlinearity can be optimized to maximize this suppression as suggested in [22]. An exponential cooling protocol proves even more efficient than the optimal nonlinear quenches in suppressing kink formation. Yet, a sudden quench to zero temperature with subsequent evolution for the same total time surpasses all to this end. The cooling scenario thus exhibits a different phenomenology from that observed in quantum phase transitions, where sudden quenches enhance defect formation over finite-time protocols.
In the opposite limit of sudden and nearly sudden quenches the cumulant values are those of a binomial distribution. As a result, the shape of the distribution varies from binomial to Poissonian as the cooling rate is decreased. Yet, even in the scaling regime for slow quenches, we have shown that in the regime of generalized scaling behavior the final state is nonthermal, by establishing the trace distance between the resulting kinknumber distribution upon completion of the quench and the corresponding one for a canonical Gibbs state. The thermal state that best approximates the final state exhibits an effective temperature that scales as an inverse power-law of the quench time.
We hope that the current findings motivate new stud-ies of the finite-time annealing dynamics of a ferromagnet beyond the KZM. A natural generalization involves the inclusion of disorder, which is known to turn the powerlaw scaling on the quench time under a linear schedule into a logarithmic dependence [9,11,19,20]. An analogous description may be invoked in one-spatial dimension [57] and the full counting statistics, as well as the role of the schedule, remain unexplored in this context. Similarly, one can envision studies in higher spatial dimensions, as well as with continuum and gauge symmetries [69].
We close by pointing out the relevance of the cooling dynamics of classical spin models in the benchmarking of quantum annealers and quantum simulators. By embedding Ising models in quantum annealing devices, tests of the KZM have been used to benchmark their performance [44,70,71] and a study of the kink statistics can provide a stringent test, helping to elucidate the kind of dynamics emulated in these devices [44]. In addition, we note that in other setups the kink number distribution can be directly measured making use of single-qubit interferometry, whereby an auxiliary qubit is used to probe the state of the Ising ferromagnet [40].
Appendix A: Computation of the second and third cumulants For convenience we introduce the notation ∂ n i1,...,in ≡ ∂ n ∂η i1 , ..., ∂η in . (A1) Using the cumulant generating function, the variance of the kink number is given by The third cumulant equals the third centered-moment and is given by Equations (25), (A2) and (A3) thus lead to the well-known results for the cumulants in terms of moments of the distribution and show the consistency of using the logarithm of (23) as the cumulant generating function.

Explicit Calculation of κ2
We begin with the expression (A2), and substitute in the explicit correlators W n as derived previously from Ψ. Lastly, we use translational invariance of the system to dispense with the index m to find a condensed form of the expression.
Note that the expression on the right hand side of the final equality of (A4) is also valid for the zero-order term, i.e. when n = m and so W 0 = 1. Taking this out of the summation and using the symmetry of the ring, (taking N even for convenience) we have In the case of W n being described by an expression of the form (13), the resulting limiting expression for κ 2 is then

Explicit Calculation of κ3
Returning to the expression for the third cumulant Making use of translational invariance in the system allows us to dispense with the first indices, defining the n'th and m'th spins relative to their distance from the first: A lengthy differentiation process of the terms containing Ψ, or equivalently collecting all index pairs in a Fermionic Wick contraction yields: Terms proportional to W 3 1 immediately cancel, as do those symmetric under the transformations n → −n and n ↔ m, with an analogous counterpart with opposite sign (i.e. the second, third and fourth brackets within the first summation), leading to Collecting alike terms in the second summation yields Making use of the periodic boundary conditions to reduce the sum, we find For a further simplification, it proves convenient to redefine indices as l = |n − m|, to find that there are (N + 1)/2 terms in which l = 0, 2((N + 1)/2 − 1) = (N − 1) terms such that l = 1, 2((N + 1)/2 − 2) = (N − 3) with l = 2, etc.
With this redefinition, we can write Evaluating the first term in the second summation, and collecting alike terms gives where we recognize the first two terms as . Once again, if in the τ Q → ∞ limit Eq. (13) holds, we have

Appendix B: Calculation of the 2-Point Correlator
In this appendix, we detail the computation of the two-point correlator W n (τ Q ) in Eq. (32) for the Ising model under Glauber dynamics and slow quenches. Before dwelling on specific cases, we note that by applying the recursion formula for modified Bessel functions the integral W n (τ Q ) to be written as

Linear Cooling
In the case of α = 1, the integral reduces to Split the integral into two parts, and define f = f (τ Q ) a function to be optimized later The Taylor and asymptotic expansions of the modified Bessel function I n (x) are given by where a k (ν) denotes a member of the class of polynomials defined by the general formula a n (ν) = (4ν 2 − 1)(4ν 2 − 9)...(4ν 2 − (2n − 1) 2 ) 8 n Γ(n + 1) .
Plugging the upper expression (B6) into the lower integral, the first term becomes where the approximation is justified pre-emptively by the choice of f (τ ) , namely that it should go to zero in the limit of large τ . In such a limit, we have that lim x→0 Solving the integral exactly, we find it to be equal to where Γ(a, b) denotes the upper incomplete gamma function. Turning our attention to the second integral, and making use of the asymptotic form (B6) we find it to be given by The integral part of the expression (B10) is exactly solvable, and the resulting form is Examining the first sum in (B9), we find: leading to the full expression for W n Now we go about optimizing f (τ Q ). We know that lim τ Q →∞ W n (τ Q ) = 1, since an infinite quench has thermal motion allowing spins to align, and continued coarsening dynamics to ensure that any excitation in the system is eventually removed. This requires that f (τ Q ) √ τ Q → ∞ as τ Q → ∞. Furthermore, we wish to have an expression with complete Γ functions, which do not depend on τ Q . Thus, in the limit of τ Q → ∞, this leads to the condition that f (τ Q ) → 0 as τ Q → ∞. Thus, we pick as a suitable choice f (τ Q ) = τ − 1 4 Q . We see upon expansion in inverse powers of τ Q that the second and third Γ functions are exponentially suppressed in powers of τ Q . Applying these conditions while taking the leading power in τ Q , and thereafter simplifying the Γ function using standard identities leaves us with the final expression: In the limit of slow cooling, T final density of defects is governed by the power-law behavior (B15)

Algebraic Cooling
We turn our attention now to the general case of algebraic cooling. Restating the integral (41) We proceed by making the substitution x = cη(2τ Q ) − α 1+α , where c = ( A 1+α ) 1 1+α is defined for convenience. Applying this substitution, and defining c = ( A 1+α ) 1 1+α (1 − A 1+α ), we find: Splitting the integral in the same fashion as before, we find that the lower contribution becomes approximately equal to n ∞ k=0 1 Γ(k + n + 1)k!
while the upper contribution is n 2π/c(2τ Q ) x k+ 3 2 . (B19) As in the previous section, we find conditions on f (τ ), which turn out to be lim τ →∞ τ , and taking the limit as before, we find that: giving the final density of defects equal to to leading order.

Exponential Cooling
In the case of exponential cooling, we begin with equation (37) W n = n where, in this case, Using the substitution u = b 1/β (1−t/τ Q ) , the integral becomes.
Inverting as in [10], we find that: and so the integral becomes Splitting the integral (B30) into an upper and lower part, while defining η = 2τ Q (ξb/ ln(τ Q )) 1/β to substitute in the upper contribution, while keeping the first contribution in terms of η for convenience gives n ∞ k=0 1 Γ(k + n + 1)k! where we have defined c = (1 − Bb 1/β β Γ(−1/β, b)) β for convenience. Taking a cue from [10] and noting that the upper part of the integral converges to zero for all ξ > 1, we find that its contribution can be replaced by n β Solving the integral (B32) leads to the expression The condition implied by both (B33) to remove the divergence and the expansion of the lower expression (B31) is that while still lim τ Q →∞ f (τ Q ) = 0. Thus, we may pick f (τ Q ) = ln(τ Q )τ 1−β Q , and find the final expression for exponential cooling to be The density of defects is, therefore