Probing the Universality of Topological Defect Formation in a Quantum Annealer: Kibble-Zurek Mechanism and Beyond

The number of topological defects created in a system driven through a quantum phase transition exhibits a power-law scaling with the driving time. This universal scaling law is the key prediction of the Kibble-Zurek mechanism (KZM), and testing it using a hardware-based quantum simulator is a coveted goal of quantum information science. Here we provide such a test using quantum annealing. Specifically, we report on extensive experimental tests of topological defect formation via the one-dimensional transverse-field Ising model on two different D-Wave quantum annealing devices. We find that the results are in qualitative agreement with the theoretical predictions for a closed quantum system, but deviate from them quantitatively. Taking into account the effects of coupling to a bosonic environment reduces the quantitative deviations significantly but not completely. Other factors, such as random control errors and transient effects arising from the limited annealing time, are likely contributors to the deviations. In addition, we probe physics beyond the KZM by identifying signatures of universality in the distribution of the number of kinks, and again find good qualitative agreement with the quantum simulator results. To check whether an alternative, classical interpretation of these results is possible, we used the spin-vector Monte Carlo algorithm, a candidate classical description of the D-Wave device We find that it too provides a qualitative fit to both the predictions of the KZM and the quantum simulator (D-Wave device) results, but the degree of quantitative agreement with the data from D-Wave annealing devices is better for the KZM, a quantum theory, than for the classical spin-vector Monte Carlo. Our work provides an experimental test of quantum critical dynamics in an open quantum system, and paves the way to new directions in quantum simulation experiments.

The number of topological defects created in a system driven through a quantum phase transition exhibits a power-law scaling with the driving time. This universal scaling law is the key prediction of the Kibble-Zurek mechanism (KZM), and testing it using a hardware-based quantum simulator is a coveted goal of quantum information science. Here we provide such a test using quantum annealing. Specifically, we report on extensive experimental tests of topological defect formation via the one-dimensional transverse-field Ising model on two different D-Wave quantum annealing devices. We find that the results are in qualitative agreement with the theoretical predictions for a closed quantum system, but deviate from them quantitatively. Taking into account the effects of coupling to a bosonic environment reduces the quantitative deviations significantly but not completely. Other factors, such as random control errors and transient effects arising from the limited annealing time, are likely contributors to the deviations. In addition, we probe physics beyond the KZM by identifying signatures of universality in the distribution of the number of kinks, and again find good qualitative agreement with the quantum simulator results. To check whether an alternative, classical interpretation of these results is possible, we used the spin-vector Monte Carlo algorithm, a candidate classical description of the D-Wave device We find that it too provides a qualitative fit to both the predictions of the KZM and the quantum simulator (D-Wave device) results, but the degree of quantitative agreement with the data from D-Wave annealing devices is better for the KZM, a quantum theory, than for the classical spin-vector Monte Carlo. Our work provides an experimental test of quantum critical dynamics in an open quantum system, and paves the way to new directions in quantum simulation experiments.

I. INTRODUCTION
Quantum simulations are emerging to be one of the important applications of quantum annealing [1][2][3][4], quite different, and arguably more natural, than the original intent of using such devices for optimization, the subject of many recent studies [5][6][7][8][9][10][11][12][13][14][15]. Prominent examples include the simulation of the Kosterlitz-Thouless topological phase transition [16,17] and three-dimensional spin glasses [18] using the D-Wave quantum annealing devices, that have successfully reproduced the behavior of various physical quantities and the structure of the phase diagram, as predicted by classical simulations. Quantum simulation has also been pursued using other systems such as trapped ions [19][20][21].
Here we use D-Wave quantum annealers to perform quantum simulations of the Kibble-Zurek mechanism (KZM) [22,23], which predicts the kink (or defect 1 ) formation when a system crosses a phase transition point at a finite rate. While the theory of the KZM was originally formulated for classical phase transitions, it has been extended to describe quantum critical dynamics [24][25][26][27]. As the dominant paradigm to describe the universal dynamics of a quantum phase transition, it has motivated a wide variety of experimental and theoretical studies [28,29]. Laboratory tests of the KZM in quantum platforms have been carried out pursuing different quantum simulation approaches, e.g., using qubits to emulate free fermion models [30][31][32][33] and via a fully digital approach using Rydberg atoms [34].
Tests of the KZM can also be used to quantitatively assess the performance of a quantum device. Indeed, the KZM scaling is sensitive to, e.g., nonlinear driving schemes [35,36], disorder [37,38], inhomogeneities in the system [39][40][41][42] and decoherence [43][44][45][46]. The use of the KZM to assess the performance of a quantum annealer was studied by Gardas et al. [47], focusing on the onedimensional case on two previous generation D-Wave 2X devices, and by Weinberg et al. [48] on a current generation D-Wave 2000Q (DW2KQ) device, focusing on the Ising Hamiltonian on the two-dimensional square lattice.
Here we report on extensive DW2KQ experiments for the one-dimensional transverse-field Ising model, using two separate realizations of the device to perform quantum simulations of the predictions of the KZM for the kink density. We also test a recent theory for the kink density distribution developed by one of us [49], thus probing physics beyond the original KZM prediction of the average number of kinks [24][25][26][27]. Unlike Gardas et al. [47], our work finds a universal power-law scaling behavior of the average number of kinks. We choose the one-dimensional problem because departures from the ideal theoretical setting due to noise and other reasons would easily destroy ordering in one dimension, and therefore it is easy to detect the effects of imperfections in one dimension, implying that the data would clearly reveal open system effects. It is also an advantage of the one-dimensional problem that we can avoid the problem of embedding of the system on the Chimera graph of the D-Wave device [50]. Moreover, previous studies of both antiferromagnetic [51,52] and ferromagnetic chains [53] using previous generations of D-Wave devices obtained good agreement with open quantum systems theory [54].
Our work establishes the power-law scaling behavior of the average number of kinks, variance and third-order cumulant with the timescale in which the transition is crossed. In doing so, we provide a strategy to assess the behavior of quantum annealers, and find it to be well described within the framework of open system quantum dynamics. Our work thus provides an experimental test of quantum critical dynamics in an open quantum system. The universal power law scaling found in the cumulants of the kink-number distribution shows that signatures of universality beyond the KZM recently predicted in isolated quantum critical systems continue to hold in the presence of coupling to an environment. This paper is organized as follows. Background on the KZM and its generalization, the problem we study, and the experimental methods are described in Sec. II. The empirical results on the kink density are presented and compared with the generalized KZM theory in Sec. III, and Sec. IV similarly presents the kink distribution results. In Sec. V we address the question of whether classical models suffice to explain our empirical results. We do this by modeling the kink distribution using the classical Boltzmann distribution of the Ising spin chain, and by comparing the empirical results to the predictions of the classical spin-vector Monte Carlo model. We close the paper with a discussion in Sec. VI, including a comparison with Refs. [47,48], and conclude in Sec. VII. Additional materials are presented in the Appendixes.

II. THEORETICAL AND EXPERIMENTAL BACKGROUND
We first describe the problem to be studied, and then explain the predictions of the KZM, followed by our experimental methods for testing the theoretical predictions.
A. The problem studied The target system of our investigation is the onedimensional transverse-field Ising model defined and parameterized by the Hamiltonian, where L is the chain length (system size), s = t/t a , with the time t ∈ [0, t a ] and the final time t = t a being the annealing time. The absolute value of the interaction strength is chosen to be |J| = 1, and as explained below we consider both the ferromagnetic (J < 0) and antiferromagnetic (J > 0) cases. We adopt a free boundary condition as indicated by the upper bound L − 1 in Eq. (1). The functional forms of the annealing schedules A(s) and B(s) are shown in Fig. 1. This one-dimensional model has a second-order quantum phase transition at A(s) = B(s) for a timeindependent system, i.e., when s is regarded as a fixed parameter [55]. The system is in the ferromagnetic phase when A(s) < B(s) and is paramagnetic for A(s) > B(s). The system is initially in the paramagnetic phase since A(0) > 0 and B(0) ≈ 0 as seen in Fig. 1. The rate of change of the annealing schedules A(s) and B(s) is finite, and the system does not necessarily follow the instantaneous ground state even when the initial condition is chosen to be the ground state of the initial Hamiltonian H(0). Thus, at the end of the annealing schedule, when A(1) = 0 and B(1) > 0, the system is generally in an excited state with a number of kinks (defects) separating ferromagnetic domains (regions with aligned neighboring spins) when J < 0, or kinks separating antiferromagnetic domains (regions with anti-aligned neighboring spins) when J > 0.

B. Kibble-Zurek mechanism and its extension
The Kibble-Zurek mechanism [27] describes the process of kink formation assuming that the ratio of the parameters A(s) and B(s) changes linearly as a function of time near the critical point A(s)/B(s) = 1. This assumption of linear time dependence in the vicinity of the critical point is reasonable because any analytical function can be expanded to linear order and we are interested in the system behavior near the critical point. Non-analytic driving schedules can also be accounted for within the KZM framework [35,36,42].
Let us state the main theoretical predictions of the KZM and its generalization to be tested on the D-Wave device.
1. The kink density ρ kink , the number of kinks divided by the system size, follows the formula where d is the spatial dimension, ν is the critical exponent for the correlation function, and z is the dynamical critical exponent. In our case d = ν = z = 1, and thus: 2. The qth cumulant κ q of the distribution function of the number of kinks P (n), divided by the average κ 1 = n , is independent of the annealing time. In particular, for the one-dimensional transverse-field Ising model the second and third cumulants satisfy 3. Since the third and higher-order cumulants are small relative to the first and second order ones, the distribution function can be well approximated by a Gaussian distribution Below we briefly describe how these formulas are derived based on Refs. [22,27] for item 1, and on Ref. [33] (its Supplementary Note 2 in particular) as well as on Ref. [49] for items 2 and 3, to provide pertinent physical background for our study. Readers interested only in the results can skip to Sec. II C.
Second-order continuous phase transitions are characterized by the divergence of the correlation length ξ and the relaxation time τ at the critical point. Specifically, as a function of the difference between the value of the control parameter λ and its critical value λ c , both quantities ξ and τ exhibit a power-law behavior where ε = (λ − λ c )/λ c , and ξ 0 and τ 0 are constants. The divergence of the relaxation time introduces a separation of time scales and allows one to describe the crossing of the phase transition as a sequence of stages: In the first stage, far from criticality where | | is not very small, the relaxation time is not large and system follows the instantaneous equilibrium state, the ground state in the context of quantum phase transition at zero temperature.
The system evolves adiabatically. Then, in the second stage, as | | becomes small, the relaxation time grows rapidly and the state of the system has no time to relax to the ground state, and the system becomes effectively frozen. As the parameter further changes, the system enters the final third stage, and | | again becomes large, thus the dynamics becomes adiabatic again. This is the so-called adiabatic-impulse approximation [27,56].
The key testable prediction of the KZM is that, after crossing the phase transition in the second stage, the average length scale in which the order-parameter is uniform is set by the equilibrium value of the correlation length when the system unfreezes at the point where the third stage is reached.
To formulate this idea quantitatively, consider a driving scheme such that the distance to the critical point varies linearly in time according to ε = (t − t c )/t a on a timescale t a , where t c denotes the instant when the system parameters cross the critical point. Equating the instantaneous equilibrium relaxation time τ (t) tô t − t c , the time elapsed after crossing the critical point, yields the freeze-out time scalet − t c = (τ 0 t zν a ) 1/(1+zν) , which yields the average correlation length asξ ≡ ξ(t) = ξ 0 (τ 0 /t a ) ν/(1+zν) . A kink may form at the interface between different domains of sizeξ. Then the kink density is given by the inverse of the volumeξ d and scales as a universal power-law [22,23] for a system in d spatial dimensions. This is Eq. (2). When the system size is L d , the average number of kinks is thus n = ρ kink L d . This picture applies both to classical and quantum phase transitions. In the quantum case, the relaxation time is identified with the inverse of the energy gap between the ground state and the first excited state that closes at the critical point, and the KZM describes the critical dynamics as well [24][25][26][27].
The above physical picture is quite generic and the result is valid independent of the details of the system Hamiltonian. If we restrict ourselves to the quantum phase transition of the one-dimensional transverse-field Ising model, more detailed information can be extracted on the distribution of kink numbers as follows [33,49].
The one-dimensional transverse-field Ising model with a periodic boundary can be solved (diagonalized) under periodic boundary condition by the Jordan-Wigner transformation [55], which rewrites the spin operators in terms of spinless fermion operators. Kinks appear in pairs under periodic boundary, and we therefore consider the number of kink pairs, which is described by the operator (for the ferromagnetic case) where L is the total number of sites and γ † k and γ k are creation and annihilation operators of fermions. The distribution function of kink pairs is defined by where ρ is the density matrix for the state after annealing. It helps to use the characteristic functionP (θ), the Fourier transform, Since kink pairs with different wave numbers are independent, the characteristic function is decomposed into a productP where We have used the Landau-Zener formula for the creation of a kink pair. Equation (11) indicates that the number of kink pairs follows the Poisson binomial distribution.
Then the cumulants are easily evaluated, the first three of which areκ In the long time scale limit t a /(2π 3 J), the first cumulant (the average), reduces tõ This is consistent with Eq. (3) of the KZM. Similarly, the second and third cumulants are evaluated to yield The cumulants κ q for the number of kinks can be derived from the above cumulants for the number of kink pairs as κ q = 2 qκ q , These give Eq. (4). The Gaussian distribution Eq. (5) follows from setting to zero all cumulants κ q with q ≥ 3, a reasonable approximation as their value is much smaller than κ 1 and κ 2 .

C. Experimental methods
We used two different DW2KQ devices, one located at the NASA Ames Research Center and the other at D-Wave Systems, Inc. in Burnaby. The latter is a lower-noise version of the former (for documentation see Ref. [57]). The D-Wave Chimera graph comprises × unit cells of sparsely connected K 4,4 bipartite graphs, for a total of 8 2 qubits, each coupled to up to 6 other qubits. We chose four chain lengths: L = 50, 200, 500, and 800. For each size we generated 200 instances of configurations of the one-dimensional chain with a free boundary by self-avoiding random walks starting from a randomly selected qubit on the Chimera graph. For each of these 200 instances, we carried out 1, 000 annealing cycles at a given annealing time t a . Thus, we generated 200, 000 samples for each t a and L, and recorded the distribution (histogram) of the kink density. The annealing time t a ranges from 1µs to 2ms, for a total of 33 values.
We tested three cases of the coupling parameter J: ferromagnetic (J = −1), antiferromagnetic (J = 1), and randomly chosen gauges. The latter starts from ferromagnetic interactions, then half of the qubits are chosen randomly and the signs of their interactions are flipped. This prescription is meant to cancel (unintended, devicespecific) local biases toward a specific direction at each qubit. As shown in Appendix A, the antiferromagnetic and random-gauge cases give almost identical results, while the ferromagnetic case tends to exhibit unstable behavior. We therefore show results for the antiferromagnetic case in the main text.

III. AVERAGE KINK DENSITY
The average kink density as a function of the annealing time t a and for different sizes L is shown in Fig. 2(a) for the NASA device and Fig. 2(b) for the Burnaby device. We analyze the data for the time range t a ≤ 100 µs because the data beyond 100 µs show different, less stable, behavior. Likely reasons include the effect of 1/f noise, which becomes apparent at long annealing times, and a significant increase in the persistent current for t a > 100 µs [57], which reduces qubit coherence. See Appendix A for data beyond 100 µs and a more detailed discussion.
The KZM assumes that the number of kinks is at least 1 on average, which means that the inequality ρ kink > 1/L should hold. We therefore also exclude the data for L = 50 from the analysis because the kink density is too low: we find empirically (see Fig. 2) that ρ kink < 1/L = 0.02, which implies that the KZM does not apply.
A first qualitative observation from Fig. 2 is that the kink density obeys a power law. It is also clearly seen that the kink density is lower on the Burnaby device in Fig. 2(b) than on the NASA device in Fig. 2(a) for the same parameter values L and t a . This is in accordance with the 'low-noise' characteristics of the Burnaby device [57].
Delving deeper into quantitative aspects, we fit the data to the power law and evaluate the exponent α. The result is given in Ta  of α may not be attributed to the nonlinear functional form of A(s) and B(s). As the schedules can be effectively linearized, corrections to KZM behavior resulting from nonlinear passage across the critical point [35,36] may be ruled out. It is thus reasonable to suspect that the difference originates from deviations from unitary dynamics that are not accounted for in the theory. A natural first step is therefore to incorporate the coupling of qubits to the environment, for which we use a spin-boson model with the following Hamiltonian [58]: where H(s) is the original Hamiltonian of Eq. (1). Independent bosons (harmonic oscillators) with frequency ω i,k couple to the z component of the ith Pauli matrix. The coefficient C k is assumed to have an Ohmic spectrum, with the sharp cutoff frequency ω c and the coupling constant η.
The Ohmic spin-boson approach has been successfully used many times in modeling the dynamics of open system quantum annealing [59][60][61][62][63][64][65][66][67][68]. In particular, Ref. [53] reported a closely related open quantum systems study of transverse field Ising spin chains with alternating sectors of strong/weak ferromagnetic coupling, but this study did not include a comparison to KZM theory.
Ground-state (time-independent) properties of the above model have already been studied by quantum Monte Carlo simulations [69] and renormalization group methods [70,71]. The conclusion of these papers is that the quantum phase transition persists under a bosonic environment and the values of the critical exponents change from ν = z = 1 for the isolated system to ν = 0.64 and z = 1.99 for the system coupled to a zero temperature bosonic environment. Note the sharp contrast with other models of decoherence which do not alter the bare critical exponents and lead to environmentally-induced heating [43][44][45].
The modified values of the critical exponents ν and z are independent of the coupling constant η(> 0) in Eq. (19). We assume that the KZM applies to the present  open system case because KZM theory is developed based only on the divergence of the relaxation time near a critical point, without recourse to a microscopic Hamiltonian. We therefore apply the generic Eq. To achieve more precise quantitative agreement between theory and experiment, we would need to incorporate addi-tional elements that have not been taken into account so far. Such may include (i) finite temperature effects not considered in Refs. [69][70][71]; (ii) transient phenomena due to the finite annealing time; and (iii) control errors, i.e., imprecision in the parameter setting in the devices [15]. The impacts of the first two items (i) and (ii) listed above were studied by some of us in Refs. [72,73]. Extensive numerical computations using the time-evolving block decimation (TEBD) method, as well as infinite-TEBD (iTEBD) combined with the quasi-adiabatic propagator path integral (QUAPI) reveal that, as shown in Fig. 3, (a) the kink density approaches a temperaturedependent constant as t a becomes very large; (b) the kink density may behave non-monotonically as a function of t a in the transient time range if the temperature is finite; (c) the effective exponent α in ρ kink ∝ t a −α around a given time t a depends on the coupling strength even when the temperature is zero.
More precisely, Fig. 3(a) shows the temperature dependence of ρ kink for a fixed coupling strength η obtained by iTEBD with QUAPI. One can see that the curve of ρ kink for finite temperature deviates upwards from that of the zero temperature case with increasing t a and the deviation is more pronounced for higher temperatures. The results for the temperatures T = 1, 2, and 5 [in units of B(1)/2k B ] imply that ρ kink behaves nonmonotonically with t a and would approach the thermal average, [1 − tanh(B(1)/2k B T )]/2, as t a → ∞. Since our data in Fig. 2 do not show an approach to a constant, we may conclude that temperature effects in the form considered in Ref. [72] have not come into play in our data for the present range of annealing time. 2 Regarding the observation (c), Fig. 3(b) shows the dependence of the slope α on the coupling strength η at zero temperature. Note the transient effects, which increase in magnitude with η, and also extend to larger t a . We would expect the exponent α to approach a constant independent of the coupling constant for sufficiently large t a , if we assume consistency with the equilibrium computations in Refs. [69][70][71], which suggest a universal exponent α = 0.28 independent of η as mentioned above. We suspect that the deviations of our experimental result 0.20 and 0.34 for the exponent α from the theoretical equilibrium value of 0.28 are at least in part a result of transient effects. These effects are difficult to analyze in a precise way because the effective exponent changes as a function of the coupling constant η and the annealing time range, and is therefore non-universal as seen in Fig. 3(b) and Fig. 1 of Ref. [72].
Noise amplitude and control errors may qualitatively explain the difference between the NASA and Burnaby devices. The latter is a newer, low-noise model, with lower 1/f noise amplitude and more accurate control [57]. Better control in the specification of system parameters, the interaction strength J between neighboring qubits as well as the local longitudinal field (which is nominally zero in the present problem), results in less randomness in the problem parameters, which should lead to a more rapid decrease of the kink density as a function of annealing time on the Burnaby device, meaning a larger exponent α. Moreover, the less noisy Burnaby device has a value closer to that of the closed-system value, and so it stands to reason that the fact that the spin-boson value is intermediate between the noisier NASA device and the Burnaby device is a reflection of the fact that the former device is more closely described as an open system than the latter. The extracted exponents of 0.34 (Burnaby) and 0.20 (NASA) are at least consistent with this picture.
Note that randomness in J from location to location necessarily induces more kinks and eventually leads to a very slow inverse-logarithmic law, instead of a polynomial decay [37,38,74].

IV. KINK DISTRIBUTION
We collected statistics of the kink density from 200, 000 samples for each t a at L = 800 as a test of the distribution function theory developed by one of us as described in Sec. II B and Ref. [49]. See also Ref. [33]. One of the important predictions of these references is that the ratio of the qth cumulant κ q (q ≥ 2) to the first cumulant κ 1 , the average, is independent of the annealing time t a [Eq. (4)]. Figure 4 shows the t a dependence of three cumulants [panels (a) and (b) for the NASA and Burnaby devices, respectively] and the ratios κ 2 /κ 1 and κ 3 /κ 1 [panels (c) and (d) for the NASA and Burnaby devices, respectively]. With the exception of κ 3 /κ 1 for the NASA device, these ratios indeed appear to be independent of t a , as predicted. The experimental values are extracted as κ 2 /κ 1 ≈ 0.61 − 0.63 and κ 3 /κ 1 ≈ 0.23 − 0.25. The theoretical predictions are 0.586 and 0.134, respectively, so the former ratio is closer to the theoretical prediction than the latter. A possible reason is the large uncertainty in statistics as reflected in the large error bars in Figure 4(c) and (d) for κ 3 /κ 1 . Indeed, the lower ends of the error bars of this ratio lie around 0.1, and the theoretically predicted value of 0.134 is within the error bars. Apart from this subtlety, the experimental data are consistent with the theory presented in Ref. [49] (see also Refs. [33,75]).
Since the third and higher order cumulants are much smaller than the second order cumulant, the distribution is well approximated by the Gaussian distribution function (5). Figure 5 shows the distributions at three values of t a . All three cases are very well approximated by this Gaussian, as drawn in solid, dashed, and dotted lines. Additional data are presented in Appendix B.
It is remarkable that we find such strong agreement between the closed-system quantum theory of Ref. [49] and the experimental results for the kink decay, cumulants, and distribution, given that the experiment is conducted at finite temperature and includes various other sources of errors, as discussed above. This suggests that these features are robust aspects of the kink statistics that lie beyond the KZM theory.

V. TESTS OF CLASSICALITY
In this section we address whether our empirical results and relatively good agreement with a quantum theory of the KZM can also be understood using a purely classical approach. We first consider a Boltzmann distribution of the kink density of the classical Ising chain, then the classical spin-vector Monte Carlo model [76], which has been successfully applied to at least partially describe the outcomes of experiments using the D-Wave devices in past studies [8,13,53,62,63,67,77], and also in recent theoretical studies of quantum annealing [78,79].

A. Boltzmann distribution and effective temperature of the kink distribution
A question of significant interest, e.g., due to the potential for quantum-assisted classical machine learning applications of quantum annealers as Boltzmann machines [80][81][82], is whether the kink distribution is thermal and well described by a Boltzmann distribution. Various previous studies have found mixed results in terms of trying to fit such thermal distributions to empirical quantum annealing data, an issue that is understood in terms of the fact that the distribution freezes once quantum fluctuations cease at low (but non-zero) values of the transverse field [65,[82][83][84]. The associated effective temperature is a relevant metric for quantifying how noisy one quantum annealer is relative to another. Given the result we found in Sec. III, that the Burnaby device generates fewer kinks than the NASA device, we expect the former to exhibit a lower effective temperature. In this section we address these issues, and provide an assessment of how well a trivial classical model matches the empirical kink distribution we have observed.
Let the empirical distribution be denoted by P (n; t a ), where n is the number of kinks, and let Q(n; β ) denote the Boltzmann distribution for the classical Ising model of a chain of length L. The latter is easily shown to be: where g(n) = L−1 n is the degeneracy of n kinks, Z = (e β +e −β ) L−1 is the partition function, E(n) = 2n+1−L is the dimensionless energy, and β is: We write β instead of β to indicate that this is a dimensionless effective inverse temperature reflecting all noise effects, not the inverse of the real physical temperature. To estimate β , we minimize (with respect to β ) the Kullback-Leibler (KL) divergence and the trace-norm distance between the experimental distribution and the Boltzmann distribution in Eq. (20). The KL divergence D KL and the trace-norm D TN are defined by using the empirical distribution P (n; t a ) and Q (n, β ) as follows: To obtain reliable estimates of the effective temperature, we first minimize the KL divergence to obtain the first approximation of the effective temperature 1/β , because the KL divergence turns out to be robust against data fluctuations. The β thus obtained is then used as the initial value in the effective temperature estimation based on the trace norm. Figure 6(a) shows the effective temperatures thus computed for L = 800. It is clear that, as expected, the NASA device has a larger effective temperature, 23% larger at t a = 100 µs, for example. This confirms the lower-noise aspect of the Burnaby device. The decrease of the effective temperature is consistent with more and more kinks being annihilated as the annealing time is increased (as seen in Fig. 2). Indeed, the kink density for the Boltzmann distribution is and given that we know that ρ kink ∝ t a −α , we expect 1/β to decrease as t a increases. Moreover, a larger α value then corresponds to a lower value of 1/β . As shown in Fig. 6(a), the effective temperature β obtained by equating Eq. (23) to the D-Wave data has almost the same value as the effective temperature obtained by Eq. (22b). Figure 6(b) compares the empirical data and the Boltzmann distribution with the optimized effective temperature at t a = 10 µs for L = 800. Although the optimized Boltzmann distribution captures the gross shape of the kink distribution, significant differences are apparent. This is consistent with the fact that, as already mentioned in Sec. IV, the actual distribution is close to Gaussian, as predicted by the quantum KZM theory. Thus, the deviation from the purely classical Boltzmann distribution of the kink density is to be expected. It is, however, possible that a closer agreement would be found with the quantum Boltzmann distribution obtained once quantum fluctuations freeze (at s < 1), but this distribution too would not be Gaussian [65]. We thus conclude that the kink distribution does not thermalize in accordance with equilibrium theory expectations, but is rather better described by the KZM and its generalization. Nevertheless, Fig. 6(b) indicates that the effective temperature obtained from our fitting procedure can serve as a reasonable proxy for quantifying the relative overall effect of noise for a comparison between different quantum annealing devices.

B. Test of a classical description by spin-vector Monte Carlo
A much more stringent model than a simple Boltzmann distribution is the standard classical model of the D-Wave devices, the spin-vector Monte Carlo (SVMC) model [76]. In our SVMC simulations, we replace the operators σ x i and σ z i by the components of a classical unit vector, sin θ i and cos θ i , respectively. The Hamiltonian is therefore written as: 3 where s n is a parameter representing time corresponding to the number of Monte Carlo steps, n. We choose the following parameterization of s n , s n = n t a N 0 (25) where N 0 is the number of Monte Carlo steps necessary to reproduce the kink density observed in the D-Wave device at 1 µs. The dimensionless parameter t a corresponds to the total annealing time, such that the total number of Monte Carlo steps is n = N 0 t a , and s n = 1 at the end of a simulation. In the present case, N 0 is 1000 and 1500 for the NASA and Burnaby devices, respectively. We use the actual NASA and Burnaby annealing schedules depicted in Fig. 1 for comparison of the DW2KQ data with SVMC results. We first set all local angles to θ i = π/2 and use a Metropolis move with the physical temperature of each device, T = 12.1 mK for NASA and T = 13.5 mK for Burnaby, and sequentially update each local state θ i . After the dynamical evolution from s = 0 to s = 1, we project the final state to +1 if 0 < θ i < π/2, and −1 if π/2 < θ i < π. We take 1, 000 samples for each t a for statistical analysis. Figure 7 shows the kink density ρ kink as a function of t a as obtained from the SVMC simulations. The power-law scaling seen for the D-Wave data in Fig. 2 and predicted from the KZM theory is observed here as well, but only for short annealing times 4 . For longer annealing times, the power law breaks down and a more rapid decay of   Fig. 7. Each exponent is obtained from a fit up to the L-dependent crossover point seen in Fig. 7. the kink-density sets in, with the crossover point increasing with chain length L. This is not the case for the DW2KQ results (see Appendix B). Moreover, the exponents extracted from the power-law regions, summarized in Table II, deviate substantially from the DW2KQ exponents summarized in Table I. It is also shown in Appendix C that the critical exponent ν assumes the value 1/2 in the SVMC model in contrast to the corresponding quantum value of ν = 1 for an isolated system and ν = 0.66 for a system coupled to a bosonic environment [69]. This implies that the closeness of the decay exponent α of the classical SVMC model to the quantum closed-system value of 0.5 is likely to be accidental. Another noticeable difference is that the kink density curves are all quite close, i.e. size-independent, for L ≥ 200 for the DW2KQ results, whereas the corresponding curves tend to differ for the SVMC simulations, with the kink density decaying more slowly for larger chain lengths.
A further test is provided by the cumulants, shown for SVMC in Fig. 8. The contrast with the DW2KQ   TABLE III. Results from SVMC model simulations at 50 mK for the exponent α of the power-law scaling describing the decay rate of the kink density as shown in Fig. 14. Each exponent is obtained from a fit up to the L-dependent crossover point seen in Fig. 14  data shown in Fig. 4 is clear, with the constancy of the ratio seen there, as predicted from generalized KZM theory, weaker in Fig. 8. We furthermore provide fits to the Gaussian distribution predicted by this theory to the SVMC simulation results in Fig. 9. Given the smallness of the third order cumulants, the Gaussian fits are unsurprisingly quite good, though not as good as to the DW2KQ data, shown in Fig. 5. Additional results for SVMC at the higher (and hence more classical) simulation temperature of 50 mK are provided in Appendix D. The overall trends are the same as those seen in Figs. 7-9, but the agreement with the predictions of generalized KZM theory is in fact closer than for the lower temperature simulations above. In particular, agreement with the power-law decay predictions for the first cumulant extend to larger t a values, as does the constancy of the cumulant ratios. The extracted decay exponent α is listed in Table III. The exponent α = 0.44 is closer to the experimental value 0.20 (NASA)/0.34(Burnaby) than the lowertemperature exponent α = 0.48 is, but the quantum- theoretical prediction 0.28 by the KZM is closer to the experimental result. One noteworthy qualitative difference is that for sufficiently large annealing times, the SVMC kink density deviates downward from the powerlaw (fewer kinks), while in the D-Wave case it deviates upward (more kinks). Finally, we also analyzed the D-Wave and SVMC data by computing their trace-norm distance from the Boltzmann distribution. Here the goal is to test the prediction of the adiabatic theorem for open quantum systems, that this distance decreases following a power law as a function of time t a for sufficiently large t a [85]. Figures 10  (a) and (b) show the result for L = 800. We computed the trace-norm distance according to Eq. (22b), which uses the optimized Boltzmann distribution Q (n; β ) in Eq. (20), and the kink distributions of D-Wave and the SVMC simulations as empirical distributions P (n). We fixed β of Q (n; β ) to that already obtained at t a = 100 µs (D-Wave) or t a = 100 (SVMC) because we are interested in how the trace-norm distance approaches the Boltzmann distribution at this annealing time, the largest reliable value available to us (see Appendix A). It is seen from Fig. 10 that the decrease of the tracenorm distance of SVMC fits an exponential (solid line), while the behavior of the D-Wave is closer to a power law for sufficiently large t a although the difference is not large. Thus it is reasonable to conclude that the D-Wave results are in closer agreement with the adiabatic theorem for open quantum systems than the classical SVMC simulation results.
Given all the discrepancies we have found, it is reasonable to conclude that the SVMC model does not explain the behavior of the DW2KQ devices reported here.  9. Histograms of the number of kinks computed using the SVMC model with the annealing schedule of (a) the DW2KQ at NASA at T = 12.1mK and (b) the DW2KQ in Burnaby at T = 13.5mK. The chain length is L = 800. From right to left, the annealing times t a are 1, 10, 100 µs, respectively. Solid, dashed, and dotted lines are the Gaussian distributions of Eq. (5).Contrast with the DW2KQ results shown in Fig. 5.

VI. DISCUSSION
We have reported on extensive one-dimensional transverse-field Ising model experiments performed using the NASA and low-noise Burnaby DW2KQ devices. We demonstrated that the kink density decays as a power law in the annealing time t a , in qualitative agreement with the theoretical prediction of the Kibble-Zurek mechanism (KZM). In more detail, we found that the exponent α describing the rate of power-law decay differs from the KZM prediction derived under the assumption of an isolated, closed quantum system. The difference between the theoretical value of α = 0.5 (which, coincidentally, is close to the outcome of the classical SVMC model as well), and the empirical values of α = 0.20 (NASA) and α = 0.34 (Burnaby), may be understood to a first approximation by modeling the coupling of the system to a bosonic environment with an Ohmic spectral density, which reduces the theoretical value to α = 0.28. Although it is difficult to quantitatively explain the remaining discrepancy, it is reasonable to suppose that it originates in other factors such as parameter control errors, and transient effects due to short annealing times. Indeed, the larger exponent (a faster decrease of the kink density) of the Burnaby device than the NASA device is consistent with the lower-noise characteristics of the former.
The lower noise aspect of the Burnaby device was further verified in our study by computing the effective tem-perature of the best-fit Boltzmann distribution at the end of the anneal, which shows that the effective temperature is about 23% higher on the NASA device. Note that this effective temperature reflects the combined effects of the dilution refrigerator temperature (which is in fact slightly higher for the Burnaby device) and a wide range of other noise sources including coupling to the environment and control errors.
Related work was reported by Gardas et al. [47] on the previous generation D-Wave 2X devices, at Los Alamos National Laboratory and in Burnaby. They also found a power-law decay of the kink density but the value of the exponent α depended strongly on experimental conditions such as the choice of the device and the sign of interactions (ferromagnetic or antiferromagnetic). The values reported range from α = 0.24 to α = 1.31, and they left the explanation for further work after listing several possible options. In contrast, our work quite definitively establishes that quantum simulation using the newer DW2KQ devices is capable of demonstrating and probing the KZM and its generalization, in particular using the lower-noise version in Burnaby with the DW 2000Q 5 solver.
Other closely related work is the recent experiment probing the two-dimensional transverse-field Ising model on the square lattice by Weinberg et al. [48], which demonstrated non-monotonicity in the kink density as a function of the annealing time t a . In the short annealing time regime, the kink density decreases as a function of t a , as in our case. The value of the exponent they found in this shorter time range, α = 0.74, is close to the theoretical value for an isolated system in two dimensions, α = 0.77. In contrast, for long annealing times, the kink density increases as t a increases. This kind of behavior is often referred to as anti-Kibble-Zurek scaling and can result from environmentally-induced heating [45,46]. Weinberg et al. also attribute this latter behavior to the effects of noise, including thermal fluctuations. Indeed, numerical calculations for the one-dimensional system shown in Fig. 3 and presented in Ref. [48] as well as in Refs. [43,72] show that the kink density can be nonmonotonic if the temperature is finite. A possible reason that our experiment did not find such non-monotonicity in the range t a ≤ 100 µs is that this annealing time is too short for the temperature effects to appear in the onedimensional problem. Our data in the range of very long annealing times up to 2000 µs show non-monotonicity in some cases and possibly reflect thermal and other deviations from the ideal quantum simulation, as discussed in Appendix A. We excluded this time range from our analysis since the data appears unstable with large uncertainties. In the short time region of the two-dimensional experiment of Ref. [48] where the KZM is likely to apply, the system seems to be much less susceptible to noise and the exponent α is close to the theoretical value of an ideal, isolated system as mentioned above. These observations suggest that how noise affects the system behavior strongly depends on the problem type as well as on the annealing time range.
We further investigated the distribution of the kink density at the end of the anneal, which encodes signatures of universality beyond the original predictions of KZM theory. We found agreement with the theoretical prediction that the ratio of the second and higher order cumulants κ q (q ≥ 2) to the first order cumulant κ 1 is independent of the annealing time t a [49]. We also found very good agreement with the prediction of a Gaussian distribution of the kink density. This agreement with a quantum theory constructed for a closed, isolated system suggests that these are robust features that remain largely intact even in the presence of coupling to an environment.
Given the history of challenges via classical models to experiments involving the D-Wave devices [76,86,87], and their rebuttals [61-64, 67, 88-90], we tested whether classical models alternatively explain the experimental data as well. We first tried a simple fit to a Boltzmann distribution but did not find very satisfactory agreement. We also tried the standard classical limit of the D-Wave device, the spin-vector Monte Carlo (SVMC) model [76], and found fair qualitative agreement. However, various important aspects of the data were not captured, and we found that overall the quantum theory of the generalized KZM provided better qualitative as well as quantitative agreement.
The two DW2KQ devices we tested therefore serve to a good approximation as quantum simulators of the onedimensional transverse-field Ising model under the influence of a dephasing Ohmic bosonic environment. These bosons do not represent thermal effects because we have observed neither an approach of the kink density to a constant, nor a non-monotonic behavior of the kink density as a function of annealing time (over the annealing time range where we have confidence in the reliability of the data). Instead, the bosons possibly correspond to dynamical fluctuations of the normal current flowing through Josephson junctions [91]. It is a difficult but interesting future direction of work to identify the nature of these fluctuations and to try to find a way to reduce them for better agreement with the closed quantum system limit.

VII. CONCLUSION
In the wake of a phase transition, topological defects form. The Kibble-Zurek mechanism predicts that the density of defects scales as a universal power-law with the time scale used to cross the transition. We have shown that this prediction can be tested on quantum annealers by using them for analog quantum simulation, i.e., as a test-bed for non-equilibrium statistical mechanics. Specifically, our work has tested the Kibble-Zurek mechanism in the one-dimensional transverse-field Ising model. Our analysis of the quantum annealer data shows that the behavior of the devices is consistent with the implementation of this model coupled to a bosonic environment. Our work thus provides experimental evidence of universal Kibble-Zurek scaling in an open quantum system, though a classical description by the spin-vector Monte Carlo model also provides fair qualitative agreement with the data.
By probing the full counting statistics of topological defects (kinks), we furthermore established signatures of universality beyond the original prediction of the Kibble-Zurek scaling, which is focused on the average kink number. In particular, we found that the power-law scaling with the annealing time of the average kink number is shared by its variance and the third cumulant. Our results thus indicate that the universal scaling recently predicted for the cumulants of the kink number distribution in an isolated quantum critical system also holds under open-system quantum dynamics.

ACKNOWLEDGMENTS
The research of YB, DL, and HN is based upon work supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.
Appendix A: Dependence of the kink density on the interaction type and annealing time Figure 11 shows the kink density obtained for three types of interactions up to t a = 2000 µs on the NASA DW2KQ device. Figure 11(a) displays the case of ferromagnetic interactions (J = −1), (b) antiferromagnetic (J = 1), and (c) random gauge. In the latter case we first set all the interactions to be ferromagnetic, randomly choose L/2 qubits, and change the sign of their interactions with their two nearest neighbors. In one dimension with free boundaries, these three cases are theoretically equivalent and the kink density should behave identically. Figures 11(a), (b), and (c) clearly indicate that the differences are small for the short time regime t a ≤ 100 µs. Beyond this regime marked deviations emerge in the ferromagnetic case, whereas the antiferromagnetic and random gauge data remain close even up to the longest annealing time t a = 2000 µs. This difference may imply the presence of a systematic bias toward ferromagnetic states in the D-Wave devices, which becomes prominent at larger annealing times. Antiferromagnetic and random-gauge interactions may be interpreted to have caused cancellations of such a bias. For these reasons we choose to use the data from antiferromagnetic interactions in the main text.
Additional data displayed in Fig. 12 show that the kink density obeys a power-law decay very accurately up to t a = 100 µs on both the NASA and the Burnaby DW2KQ devices with antiferromagnetic interactions and chain length L = 800. Beyond t a = 100 µs, deviations from a power law become apparent. This data set, along with Fig. 11, motivated us to use the empirical data only for t a ≤ 100 µs, in order to eliminate artifacts other than the coupling to bosonic environment.  We show in this Appendix that the critical exponent ν is 1/2 for the one-dimensional ferromagnetic SVMC model. The Hamiltonian is We have exchanged sin θ j and cos θ j from the conventional notation of Eq. (24) for later convenience. A periodic boundary is assumed. Since we are interested in how the system behaves as we decrease Γ from a very large value (where the system is in the paramagnetic phase with θ j ≈ 0 ∀j) toward a transition point, it is reasonable to expand the Hamiltonian to quadratic order as where we have ignored a constant term. Let us check if the paramagnetic state is stable by Fourier transformation,  It is observed that φ k = 0 ∀k is the stable state (ground state) configuration for Γ > 2J, and a second-order phase transition exists at Γ c = 2J. The correlation is given by where the angular brackets · · · stand for the statisticalmechanical average with respect to the Hamiltonian (C4). We have averaged over l using translation symmetry. Since the Hamiltonian (C4) represents independent Gaussian fields, we easily find Thus the correlation function (C5) becomes The behavior as r 1 is evaluated as Thus the exponent is ν = 1/2. Notice that the temperature is kept small but finite. If T = 0 exactly, no fluctuations exist classically and the spin configuration is fixed to θ j = 0 in the paramagnetic phase.  Table III. These values are somewhat closer to the DW2KQ values than those for the dilution refrigerator temperatures, and we also note that the qualitative behavior of the kink density curves seen in Fig. 14 is more like the DW2KQ results seen in Fig. 2, in that the curves for the two largest L values now nearly overlap. In these respects the warmer SVMC model is a closer match to the DW2KQ data than at the dilution refrigerator temperatures.