Simulated quantum annealing as a simulator of non-equilibrium quantum dynamics

Simulated quantum annealing based on the path-integral Monte Carlo is one of the most common tools to simulate quantum annealing on classical hardware. Nevertheless, it is in principle highly non-trivial whether or not this classical algorithm can correctly reproduce the quantum dynamics of quantum annealing, particularly in the diabatic regime. We study this problem numerically through the generalized Kibble-Zurek mechanism of defect distribution in the simplest ferromagnetic one-dimensional transverse-field Ising model with and without coupling to the environment. We find that,in the absence of coupling to the environment, simulated quantum annealing correctly describes the annealing-time dependence of the average number of defects, but a detailed analysis of the defect distribution shows clear deviations from the theoretical prediction. When the system is open (coupled to the environment), the average number of defects does not follow the theoretical prediction but is qualitatively compatible with the numerical result by the infinite time-evolving block decimation combined with the quasi-adiabatic propagator path integral, which is valid in a very short time region. The distribution of defects in the open system turns out to be not far from the theoretical prediction. It is surprising that the classical stochastic dynamics of simulated quantum annealing ostensibly reproduce some aspects of the quantum dynamics. However, a serious problem is that it is hard to predict for which physical quantities in which system it is reliable. Those results suggest the necessity to exert a good amount of caution in using simulated quantum annealing to study the detailed quantitative aspects of the dynamics of quantum annealing.


I. INTRODUCTION
Quantum annealing was originally proposed as a metaheuristic to solve classical combinatorial optimization problems [1][2][3][4][5][6][7][8][9][10]. Recent years have seen its further development as a quantum simulator of materials [11][12][13][14][15][16][17][18][19][20]. In both of these applications, it is important to check the reliability of outputs from real quantum devices against the data obtained on classical computers by numerical solutions of the Schrödinger equation for closed systems (isolated from the environment) and those of quantum master equations for open systems (coupled to the environment). It is, however, difficult to carry out those numerical studies on classical hardware beyond moderate sizes because the dimension of the relevant Hilbert space increases exponentially with the number of qubits.
Classical simulations of quantum annealing using stochastic processes are often considered as an alternative tool to study the properties of quantum annealing for large systems. In the spin-vector Monte Carlo (SVMC) [21], one replaces spin- 1 2 Pauli operators in the Hamiltonian by classical rotors of unit length and stochastically updates the system state. This method is known to reproduce some features of the outputs from the quantum device, the D-Wave quantum annealer [22][23][24]. Simulated quantum annealing (SQA) is another powerful classical tool that uses the path-integral Monte Carlo [2,5,18,23,[25][26][27][28][29][30][31][32][33][34][35]. The latter method is in principle designed to simulate equilibrium properties of quantum systems without a sign problem [36]. However, in the context of quantum annealing, it is used to simulate dynamical behaviors of the transverse-field Ising model as initially attempted in Ref. [2]. Recent examples include Ref. [35], where the dynamics of a random Ising chain with transverse field was studied by SQA, and Refs. [37,38], where the performance of various embedding schemes was compared using SQA.
There are few a priori reasons to expect that the classical stochastic dynamics of SQA faithfully reproduces the quantum dynamics of quantum annealing, which often operates in the regime out of equilibrium or away from the adiabatic limit. It is nevertheless known that some aspects of the non-equilibrium dynamics of quantum annealing can be understood by SQA, notably in relation to incoherent quantum tunneling [34]. The present paper tries to shed light on this problem through the analyses of the simplest case of the non-random one-dimensional ferromagnetic Ising model in a transverse field. In particular, we study how far SQA reproduces the predictions of the generalized Kibble-Zurek mechanism on the distribution of defects after the system is driven across a critical point at a finite rate [39,40]. We also compare SQA with the direct numerical method for low-dimensional quantum systems, the infinite time-evolving block decimation (i-TEBD) combined with the quasi-adiabatic propagator path integral (QUAPI) [41,42], which is believed to pro-duce accurate results in a very short time region.
The problem of the generalized Kibble-Zurek mechanism in one dimension was studied extensively in Ref. [17], where it was seen that the behavior of the D-Wave device can be understood in terms of a quantum system under the influence of an environment. The present paper has a different point of view, i.e., not to examine the behavior of the D-Wave device but to study how far SQA is useful to explain the dynamical properties of the quantum system through comparison of the SQA data with the predictions of the generalized Kibble-Zurek mechanism as well as with the data from the i-TEBD-QUAPI.
The next section describes the problem and the methods of analyses. Section III reports the results. The final section is devoted to a summary and discussion.

II. PROBLEM AND METHODS
In the present section, we describe the problem to be studied, the methods of analysis, and the physical quantities to be observed.

A. Closed-and open-system Hamiltonians
We study the non-equilibrium dynamics, in particular the properties related to the original [43,44] and the generalized Kibble-Zurek mechanisms [39,40], of the ferromagnetic transverse-field Ising model in one dimension with a free boundary, Here L is the number of sites (qubits),σ µ i (µ = x, z) is the µ component of the Pauli matrix at site i, and J(t) and Γ(t) are the time-dependent coefficients for the target and driver terms, respectively. We use the linear time dependence of the coefficients for simplicity, where t a is the annealing time and the time t runs from 0 to t a . This system in the ground state undergoes an equilibrium second-order phase transition at a critical point J = Γ in the limit of infinite system size [45]. In addition to the above Hamiltonian representing a closed quantum system, we also study the open system described by the Hamiltonian, where V k is the coupling constant with the kth bosonic modeâ † i,k andâ i,k representing the environment, which is assumed to have the standard Ohmic spectral density J(ω) = 4π k V k 2 δ(ω − ω i,k ) = 2παω. Following Ref. [46], we assume J(ω) to be cutoff at some ω c , beyond which J(ω) is set to 0. The coefficient α describes the magnitude of dissipation caused by the coupling to the environment. The whole system of Eq. (4) is isolated from other degrees of freedom, kept at zero temperature, and in principle running under the unitary Schrödinger dynamics without the constraint of adiabaticity. SQA is supposed to simulate the dynamical behavior of this system.

B. Method
The problem defined above has been studied extensively for many years. Directly pertinent to the present paper are, first, the path-integral Monte Carlo simulation by Werner et al. [46], who studied the equilibrium phase diagram and critical exponents of the open system of Eq. (4), and found that critical exponents are different from those of the closed system, Eq. (1), but are independent of the coupling strength α(> 0). Also, Bando et al. [17] carried out experiments on the D-Wave quantum annealers and compared the data with the generalized Kibble-Zurek mechanism [39,40] to conclude that the data from the devices can be understood in terms of a system under the bosonic environment of Eq. (4), not of the closed system Eq. (1), and also that the classical SVMC is unable to explain the dynamical properties of the D-Wave devices.
We combine these two contributions and study the systems of Eqs. (1) and (4) by SQA, path-integral Monte Carlo simulations with time-dependent coefficients J(t) and Γ(t), and compare the data with those from the generalized Kibble-Zurek mechanism as well as from the more direct numerical method of the iTEBD-QUAPI [41,42], in which the time-dependent density operator is expressed by a matrix-product state after a Trotter decomposition and tracing out the bosonic degrees of freedom. Then a controlled truncation of basis states allows one to perform evaluation of physical quantities. See Refs. [41,42] for details.
To perform SQA, we first apply the Suzuki-Trotter decomposition [36] to the expression of the partition function for the Hamiltonian of Eq. (4). The partition function is then expressed as the corresponding partition function of a classical Ising model in two dimensions with long-range interactions in the Trotter (imaginary-time) direction [46], where Z 0 is the partition function of free bosons and .
Here, β eff = β/P with β the inverse temperature and P the number of Trotter slices along the imaginary-time axis. The coefficient in the second line is defined as . A periodic boundary condition S i (P + 1) = S i (1) is imposed along the Trotter axis. The closed system of Eq. (1) is recovered by setting α = 0. Zero-temperature properties of Eq. (4) can be simulated with sufficiently large values β and P with a fixed finite value of the ratio β eff = β/P [2,5,6]. We choose β eff = 1 and P = 4L to satisfy this condition as described in more detail below.
Although cluster updates have sometimes been used in related studies [32,35,46], we use the simple singleflip Metropolis-update dynamics because, first, it is nontrivial which stochastic dynamics better simulates quantum dynamics, and second, we do not expect to encounter the problem of slow relaxation, which hampers simulations of disordered systems based on simple single-spin updates [35], or we do not need to reach the equilibrium necessary for the study of critical exponents [46]. The annealing time t a is identified with the total number of Monte Carlo steps divided by the number of sites LP and will be denoted by τ MCS (= t a ). Values of time-dependent coefficients J(t) and Γ(t) are updated according to Eqs.

C. Physical quantities
Following Ref. [17], we measure several quantities by SQA to be compared with predictions of the generalized Kibble-Zurek mechanism and the iTEBD-QUAPI. All quantities are measured at the end of computation where Γ(t a ) = 0.
The main quantity of interest is the number of defects (spatially mis-aligned spin pairs) n, Notice that the true ground state of the effective Hamiltonian of Eq. (6) is perfectly ferromagnetic with n = 0. We measure the statistics of defects to be denoted by P SQA (n) and compare it with the prediction of the generalized Kibble-Zurek mechanism [39,40]. According to this theory, the defect distribution is bimodal and is well approximated by the Gaussian function Q(n) for large systems L ≫ 1, where κ 1 ≡ n and κ 2 ≡ (n − n ) 2 . Here angular brackets · · · denote the average. We quantify the difference between P SQA (n) and Q(n) by the L1 norm, We also follow Ref. [17] and check the proximity of P SQA (n) to the Boltzmann distribution of the classical Ising part of the original Hamiltonian in Eq. (1), where E(n) is the energy of the classical Ising model with n defects. We optimize the effective inverse temperature β BL in P BL (n) by minimizing the L1 norm of Eq. (9) (after replacing Q(n) with P BL (n)) such that P BL (n) approximates P SQA (n) as faithfully as possible. Notice that β BL is unrelated to any physical temperature but is a parameter to fit Eq. (10) to the data. The number of defects n in Eq. (7) is essentially equivalent to the residual energy per site E res , the difference between the achieved energy and the true ground-state energy, for which we use the expression of Ref. [35], Note that E res can also be regarded as the average number of defects, the average over the statistics of P SQA (n), which was denoted by κ 1 (the first cumulant or the average) above. The data for E res are to be compared with the predictions of the Kibble-Zurek mechanism, which is described qualitatively as follows. As the system approaches a critical point, the correlation length and relaxation time grow rapidly. When the system is not close enough to the critical point, the relaxation time is large but still shorter than the time scale of annealing, and the system has time to relax to equilibrium. As the system comes closer to the critical point, the relaxation time becomes long enough and the system has no time to relax and is effectively frozen with a non-vanishing residual energy. This picture leads to the asymptotic power law of the residual energy [47], Here, d = 1 is the spatial dimension and z and ν are the dynamical and correlation-length critical exponents, respectively.

III. RESULTS
We now present our results obtained by SQA.

A. Closed system
We first report the results for the closed system of Eq. (1). Figure 1 shows the residual energy E res defined in Eq. (11) as a function of the annealing time τ MCS for system sizes L = 64, 128, 256, and 512. We set β eff = β P = 1 and P = 4L. As shown in Appendix A, we have confirmed that the residual energy as a function of P converges for P = 4L. This condition is used throughout this paper. It is clearly shown in Fig. 1 that the residual energy follows a power law as predicted by the original Kibble-Zurek mechanism, The values of the exponent b are found as b = 0.543 ± 0.007 (L = 64), 0.519 ± 0.004 (L = 128), 0.506 ± 0.003 (L = 256), and 0.500 ± 0.003 (L = 512). These values, particularly the last one, are in close agreement with the prediction of the Kibble-Zurek mechanism Eq. (12) for this one-dimensional closed system, which has z = 1, ν = 1, and thus b = dν/(1 + zν) = 1 2 . This agrees with previous reports by the direct numerical method of the i-TEBD-QUAPI [41,42] as well as by SQA [32,35]. By contrast, as found in Ref. [17], the SVMC shows close, but slightly deviated, values b = 0.477 ± 0.005 and b = 0.482 ± 0.006, depending on the choice of the simulation temperature, 12.1 mK for the former value corresponding to the device at NASA Ames Research Center and 13.5 mK for the latter value for the device at Burnaby.
We proceed to the analysis of the distribution function of defects P SQA (n). Figure 2 shows the qth cumulants κ q (q = 1, 2, 3) of the defect distribution obtained from the final state at τ MCS of SQA for L = 64, 128, 256, and 512. The generalized Kibble-Zurek mechanism [39,40]  predicts that the second-and higher-order cumulants are proportional to the first-order cumulant (average) κ 1 . The data shows that this is indeed the case for τ MCS up to about 10 2 for any system size L but the tendency changes beyond τ MCS ≈ 10 2 . Comparison of data for different system sizes suggests that this deviation may possibly be a finite-size effect. Figure 3 shows the ratios of cumulants κ 2 /κ 1 and κ 3 /κ 1 for τ MCS up to about 10 2 for L = 512, which indicates the validity of proportionality of κ 2 and κ 3 to κ 1 up to τ MCS ≈ 10 2 . More quantitatively, we find κ 2 /κ 1 = 0.688 ± 0.005 and κ 3 /κ 1 = 0.394 ± 0.019. These values clearly deviate from the theoretical prediction for the closed system, κ 2 /κ 1 = 0.578 and κ 3 /κ 1 = 0.134 [39,40]. We therefore conclude that SQA does not faithfully reproduce the statistics of the defect distribution for the present closed system although the average or the residual energy is well described by SQA.
It may be useful to recall in passing that the data from the D-Wave devices showed similar proportionality but with values closer to the theory, κ 2 /κ 1 = 0.61 ± 0.03 and κ 3 /κ 1 = 0.23 ± 0.15 on the device at NASA and  4,118] for L = 512 for the closed system. The constants are evaluated as κ2/κ1 ≈ 0.688 ± 0.005 and κ3/κ1 ≈ 0.394 ± 0.019, respectively. κ 2 /κ 1 = 0.63±0.05 and κ 3 /κ 1 = 0.25±0.18 on the device at Burnaby [17]. In contrast, the SVMC data marginally indicate proportionality of κ 2 to κ 1 in the short time region but with significant deviations for longer annealing times [17]. The ratio κ 2 /κ 1 out of SVMC is close to 0.6 for the short-time region but the ratio κ 3 /κ 1 has large error bars and it is impossible to determine its value with reliability.
We further follow Ref. [17] and test if the defect distribution is closer to the Boltzmann or Gaussian distribution. Figure  appears to be close to the optimized Boltzmann distribution P BL (n), we see a slight deviation especially at τ MCS = 16, which is more clearly observed in the enhanced figure in Fig. 5, where we see a better fit by the Gaussian Eq. (8), similarly to the previous study by the D-Wave devices and SVMC [17].  Figure 6 shows the L1 norm of the difference between the defect distribution of SQA, P SQA (n), and the optimized Boltzmann distribution P BL (n) in Eq. (10) and the Gaussian distribution Q(n) in Eq. (8) as a function of the annealing time τ MCS for L = 512. In the time range up to τ MCS ≈ 10 2 , where the proportional relationship of cumulants is established, the L1 norm with the Gaussian is smaller, and thus SQA works as a Gaussian sampler rather than as a Boltzmann sampler, but for longer annealing times, e.g., τ MCS ≈ 10 3 , two functions show similar degrees of proximity to the data. Those results for the closed system indicate that the classical stochastic dynamics of SQA partly succeeds in reproducing the predictions on the quantum dynamics related to the original and generalized Kibble-Zurek mechanism. However, detailed quantitative analysis of the distribution function show deviations from the theory for the closed system.

B. Open system
We next discuss the results obtained for the open system defined in Eq. (4) and simulated by Eq. (6). Figure 7 shows the residual energy E res as a function of τ MCS for L = 64, 128, 256, and 512 under the same parameters as before, β eff = β P = 1 and P = 4L. Here, we fix the coupling strength with the boson field to α = 0.6 for all L. It is seen that the power decay pre-   [17,41,42]. We therefore checked the α-dependence of the residual energy by SQA, and the result is in Fig. 8 for L = 256. It is seen that the exponent b is a decreasing function of α (from b = 0.50 for α = 0 to b = 0.26 for α = 1.0), consistently with the results of i-TEBD-QUAPI reported in Refs. [17,41,42]. For the reader's convenience, we reproduce a figure from Ref. [17] as Fig. 15 of Appendix B of this paper, where the defect density (proportional to the residual energy) is drawn as a function of the annealing time. One observes there that b changes from b = 0.50 for α = 0 to b = 0.25 for α = 1.28, close to the values from SQA.
As suggested in Ref. [17], taking into account the fact that the i-TEBD-QUAPI reproduces quantum dynamics for very short times, this non-universality (dependence of the exponent on α) of the SQA data may imply that the system is still in a transient state and much longer annealing times for much larger systems may show a value of b independent α as expected from universality seen in the equilibrium Monte Carlo simulation [46]. Viewed differently, this consistency with the quantum-mechanical simulation by the i-TEBD-QUAPI may imply that SQA, a classical stochastic process, reproduces some aspects of the non-equilibrium quantum dynamics in the short time region for the present open system. We proceed to test if the distribution function of defects P SQA (n) for a typical fixed value of α = 0.6 is consistent with the prediction of the generalized Kibble-Zurek mechanism as we did for the closed system. Figure 9 shows the qth cumulants of the distribution as functions of the annealing time for L = 64, 128, 256, and 512. It is seen that the second cumulant is proportional to the first cumulant almost up to τ MCS ≈ 10 3 , a longer time range of proportionality than in the closed system depicted in Fig. 2. The behavior of the third cumulant is unstable due to insufficient statistics. Figure 10 shows the ratios of cumulants with the results κ 2 /κ 1 ≈ 0.598 ± 0.008 and κ 3 /κ 1 ≈ 0.185 ± 0.048. These values for the open system are closer to the theoretical prediction of the generalized Kibble-Zurek mechanism (κ 2 /κ 1 = 0.578 and κ 3 /κ 1 = 0.174 [40]) than in the case of SQA for the closed system discussed in the preceding section (κ 2 /κ 1 ≈ 0.688 ± 0.005 and κ 3 /κ 1 ≈ 0.394 ± 0.019). We have also calculated cumulant ratios with a different coupling strength α = 1.0 in L = 512 to see the α dependence. The result is κ 2 /κ 1 ≈ 0.569±0.015 and κ 3 /κ 1 ≈ 0.196 ± 0.135, indicating very weak depen- dence on α.
It was found in Ref. [17] that the numerical data by the i-TEBD-QUAPI for the open system indicate a similar value 0.6 for κ 2 /κ 1 independent of α but it is difficult to evaluate κ 3 /κ 1 . We have thus found that SQA for the open system successfully reproduces the prediction of the generalized Kibble-Zurek mechanism for the cumulant ratios, but the behavior of the average (residual energy) is non-universal, which is not compatible with the theory. Also, agreement with the data of the iTEBD-QUAPI in the short time region has been observed. As in the case of the closed system, we further check whether the distribution function is closer to the Gaussian or Boltzmann. Figure 11 shows the distribution function for τ MCS = 16, 256, and 4096 for L = 512 with the coupling constant α = 0.6 fitted to the Boltzmann distribution Eq. (10) with optimized effective temperature. Although the gross feature is captured by the Boltz- mann distribution, we find clear deviations. In fact, as seen in Fig. 12, the Gaussian distribution better matches the data, as anticipated from the smallness of the third cumulant.  Figure 13 shows the L1 norm between the defect distribution of SQA, P SQA (n), and the optimized Boltzmann distribution P BL (n) in Eq. (10) and the Gaussian distribution Q(n) of Eq. (8) as a function of the annealing time τ MCS for L = 512. It is seen that the data is closer to the Gaussian distribution than to the Boltzmann, similarly to the case of the closed system and in agreement with the generalized Kibble-Zurek mechanism which predicts small values of higher-order cumulants than the second order, meaning Gaussian approximately. Thus SQA is closer to a Gaussian sampler rather than a Boltzmann sampler in the present open system as well. A similar conclusion was drawn for the D-Wave device [17], and therefore we should be careful when we use quantum annealing (simulated or on the real device) for sampling purposes at least as long as the present one-dimensional system is concerned.  13. Distance, the L1 norm, between the defect distribution P SQA (n) of SQA and the optimized Boltzmann distribution P BL (n) (circle) and the Gaussian distribution Q(n) (square) for L = 512.

IV. DISCUSSION
We have performed numerical tests to check if simulated quantum annealing (SQA) is able to describe the non-equilibrium quantum dynamics of the simple onedimensional ferromagnetic transverse-field Ising model across its critical point using the generalized Kibble-Zurek mechanism and the direct quantum numerics of i-TEBD-QUAPI as tools to measure the degree of success of SQA. This is a highly non-trivial problem because the classical algorithm SQA has no a priori reason to reproduce quantum dynamics, although SQA has often been used to simulate quantum annealing, in particular for large systems. It is nevertheless known that the phenomenon of incoherent quantum tunneling is well simulated by SQA if there exists a finite number of energy barriers between local minima [34], a prototypical energy landscape of a simple first-order phase transition. We have studied the problem from the point of view of non-equilibrium dynamics across a critical point (secondorder transition point) in the open system keeping in mind that the data from SQA is known to agree with the prediction of the original Kibble-Zurek mechanism in the closed one-dimensional system without disorder [35].
We presented detailed quantitative comparison of defect distribution with the asymptotically exact theory of the generalized Kibble-Zurek mechanism and the i-TEBD-QUAPI, the latter of which is expected to faithfully reproduce quantum dynamics in the very short time region in the present one-dimensional system. Our results indicate that some, but not all, of the dynamical properties of the quantum system can be reproduced by SQA for the one-dimensional system: In the absence of coupling to the environment, SQA correctly describes the annealing-time dependence of the residual energy, but the ratios of cumulants of defect distribution clearly deviate from theoretical values. When the system is open, the residual energy does not follow the prediction of the Kibble-Zurek mechanism but shows compatibility with the numerical result by the iTEBD-QUAPI. The ratios of cumulants for the open system turn out to be closer to the theoretical values of the generalized Kibble-Zurek mechanism than in the case of the closed system. The distribution of defects is closer to the Gaussian function than to the Boltzmann distribution in both closed and open systems, which is also the case in the D-Wave quantum annealers Ref. [17], meaning that SQA and the D-Wave device serve as Gaussian samplers, not as Boltzmann samplers.
It is remarkable that the classical stochastic dynamics of SQA ostensibly reproduces some of the properties of quantum dynamics for the present system. Nevertheless, it is hard to predict when it is reliable for what physical quantities, which is a serious problem in practical applications. It is therefore necessary to stay very cautious when one uses SQA to clarify the dynamical properties of quantum annealing in a detailed quantitative way. It is also desirable to establish a theoretical framework to explain the present numerical observations. Though there exists a formal mapping between classical stochastic dynamics and a quantum-mechanical system [48,49], this is still far from sufficient to understand the results found in the present paper. It is also known that the convergence condition of SQA in the long-time limit has a very similar expression to that of quantum annealing [8,50]. Whether or not this is a coincidence is an interesting topic closely related to the present work.
is enough to reach convergence of the residual energy for those system sizes and annealing time. In Fig. 15, we reproduce Fig. 4(b) of Ref. [17], where the density of defects (kinks), proportional to the resid-ual energy, is plotted as a function of the annealing time obtained from the extensive numerical computation of i-TEBD-QUAPI. The exponent b for the decay of the residual energy, Eq. (13), is seen to depend on the coupling constant α with the bosonic environment in reasonable agreement with the SQA result described in Sec III B. Taken from Ref. [17].