Effectiveness of quantum annealing for continuous-variable optimization

The application of quantum annealing to the optimization of continuous-variable functions is a relatively unexplored area of research. We test the performance of quantum annealing applied to a one-dimensional continuous-variable function with a rugged energy landscape. After domain-wall encoding to map a continuous variable to discrete Ising variables, we first benchmark the performance of the real hardware, the D-Wave 2000Q, against several state-of-the-art classical optimization algorithms designed for continuous-variable problems to find that the D-Wave 2000Q matches the classical algorithms in a limited domain of computation time. Beyond this domain, classical global optimization algorithms outperform the quantum device. Next, we examine several optimization algorithms that are applicable to the Ising formulation of the problem, such as the TEBD (time-evolving block decimation) to simulate ideal coherent quantum annealing, simulated annealing, simulated quantum annealing, and spin-vector Monte Carlo. The data show that TEBD's coherent quantum annealing achieves far better results than the other approaches, demonstrating the effectiveness of coherent tunneling. From these two types of benchmarks, we conclude that the hardware realization of quantum annealing has the potential to significantly outperform the best classical algorithms if thermal noise and other imperfections are sufficiently suppressed and the device operates coherently, as demonstrated in recent short-time quantum simulations.

The application of quantum annealing to the optimization of continuous-variable functions is a relatively unexplored area of research.We test the performance of quantum annealing applied to a one-dimensional continuous-variable function with a rugged energy landscape.After domainwall encoding to map a continuous variable to discrete Ising variables, we first benchmark the performance of the real hardware, the D-Wave 2000Q, against several state-of-the-art classical optimization algorithms designed for continuous-variable problems to find that the D-Wave 2000Q matches the classical algorithms in a limited domain of computation time.Beyond this domain, classical global optimization algorithms outperform the quantum device.Next, we examine several optimization algorithms that are applicable to the Ising formulation of the problem, such as the TEBD (time-evolving block decimation) to simulate ideal coherent quantum annealing, simulated annealing, simulated quantum annealing, and spin-vector Monte Carlo.The data show that TEBD's coherent quantum annealing achieves far better results than the other approaches, demonstrating the effectiveness of coherent tunneling.From these two types of benchmarks, we conclude that the hardware realization of quantum annealing has the potential to significantly outperform the best classical algorithms if thermal noise and other imperfections are sufficiently suppressed and the device operates coherently, as demonstrated in recent short-time quantum simulations.
Combinatorial optimization problems with discrete variables are the main target of QA, but the basic idea of using quantum fluctuations to search for the lowestenergy state in the rugged energy landscape is also applicable to problems with continuous variables.An example is the early paper by Finnilla et al. [18] although they used an imaginary-time version of the Schrödinger equation without quantum coherence.Stella et al. [19] studied a few one-dimensional problems with and without rugged energy landscapes and compared QA with its classical counterpart of simulated annealing (SA) by directly solving the Schrödinger equation for QA and the Fokker-Planck equation for SA.
They found that the relative efficiency of QA and SA depends on the shape of the energy landscape.Simulated annealing is influenced by the height of the energy barrier, whereas QA reflects the structure of the energy spectrum, which affects the probability of quantum tunneling.An apparent difference in performance between SA and QA was clearly seen in a parabolic washboardlike potential with a single global minimum and many local minima as proposed by Shinomoto and Kabashima [20].These latter authors used a coarse-grained phe- * arai.s.ao@m.titech.ac.jp nomenological version of the problem to show that the best possible performance of SA can be achieved when the rate of temperature decrease is proportional to the inverse of logarithm, which results in a very slow logarithmic decrease of the residual energy as a function of the annealing time.Stella et al. discussed also QA phenomenologically for the same washboard-like potential by considering a coarse-grained tight-binding Hamiltonian under a semi-classical approximation, which consists of a harmonic potential and a Laplacian term with a coefficient representing a semi-classical WKB tunneling effect.They found that the residual energy shows a power law decrease, which is faster than the logarithmic decrease in SA as found by Shinomoto and Kabashima.Inack and Pilati studied the performance of QA by the projective quantum Monte Carlo method [21] and found that this imaginary-time formulation behaves better than SA in the long-time regime.Koh and Nishimori recently analyzed a problem with a Shinomoto-Kabashima-like potential with many local minima [22], which is known as the Rastrigin function in the field of nonconvex continuous optimization problems [23].They solved the timedependent Schrödinger equation directly numerically for QA and the master equation for SA and confirmed that the residual energy under QA decreases in a power law with an arbitrarily large power achievable by adjusting the annealing schedule, whereas SA has a smaller power.This is a simple but nontrivial example of a continuous optimization problem, in which QA has an advantage over SA.
The studies mentioned above were based on theoretical or numerical approaches.Several experimental investigations on real hardware have recently been reported using domain wall encoding [24] to represent continuous variables in terms of Ising spins.Examples include Abel et al. arXiv:2305.06631v2[quant-ph] 4 Oct 2023 [25,26] who observed quantum tunneling behavior from a local minimum to the global minimum in a double-well potential.In another paper, Abel et al. [27] have applied the same method to two-dimensional continuous optimization problems with many local minima using reverse annealing (RA) [28][29][30][31], which is a variant of QA starting from a classical spin configuration [32].They graphically compared the output of QA on the D-Wave hardware with data from classical continuous-variable optimization algorithms and found that the former reached the global minimum much more reliably than the latter.This is a very interesting development, but detailed analysis is still needed to quantitatively demonstrate a clear advantage of QA over state-of-the-art classical algorithms tailored for continuous-variable optimization.
In the present paper, we systematically investigate whether QA has an advantage over classical algorithms for optimization of the Rastrigin function in one dimension, adopting domain-wall encoding to map a continuous variable to Ising spins.The first part of our test compares the hardware D-Wave 2000Q with several wellestablished classical algorithms for continuous-variable optimization, Nelder-Mead (NM) [33], conjugate gradient descent (CGD) [34], basin-hopping (BH) [35], and differential evolution (DE) [36].We have used the D-Wave 2000Q, not the latest D-Wave Advantage, because the former produced better results than the latter as will be described later, in particular in Appendix B 1. The second part also runs the hardware, but we additionally test a classical algorithm to faithfully simulate quantum annealing, TEBD (time-evolving block decimation) [37,38] based on a matrix-product representation of the wave function, in order to clarify to what extent the hardware runs coherently.We also study several classical algorithms for combinatorial optimization, simulated annealing (SA) [39], simulated quantum annealing (SQA) [40], and spin-vector Monte Carlo (SVMC) [41], to investigate how classical algorithms are affected by the barrier height compared to the quantum case.
This paper is organized as follows.In Sec.II, we explain the Rastrigin function and how to implement continuous-variable optimization with domain wall encoding.In Sec.III, we perform benchmark tests for the Rastrigin function on the D-Wave 2000Q and several classical algorithms developed for continuous-variable optimization.In Sec.IV, we next compare D-Wave 2000Q, TEBD, SA, SQA, and SVMC, all of which are for discrete variable problems, to see how thermal effects influence the performance.Finally, Sec.V discusses the results and concludes the paper.A few technical details are described in the Appendixes.

II. THE PROBLEM AND METHODS
In this section, we explain the mathematical definition of the Rastrigin function and how to optimize the continuous function with the domain wall encoding method

A. Rastrigin function
The Rastrigin function is one of the standard benchmark functions for nonconvex optimization problems [23].We consider the one-dimensional Rastrigin function which consists of a harmonic potential term with a unique global minimum at x = 0 and a cosine term that gives many local minima as where k is the spring constant, and h 0 and w 0 are the control parameters of the height and width of each local minimum.The above mathematical formulation of the function was introduced by Koh and Nishimori [22] and they called it the Shinomoto-Kabashima-like potential.The original formulation of the Rastrigin function has no control parameters and corresponds to k = 2, w 0 = 1, and h 0 = 20 in Eq. (1).In this study, we do not adopt the original formulation, but Koh and Nishimori's formulation to examine the dependence of algorithm performance on the shape of the energy landscape.We illustrate the one-dimensional Rastrigin function with k = 0.5 and w 0 = 0.2 for different values of h 0 in Fig. 1.Due to the second cosine term, the simple gradient descent algorithm stops at a local minimum and fails to find the global minimum efficiently.

B. Domain wall encoding for continuous function optimization
We briefly formulate QA and then describe the way to perform continuous function optimization by quantum annealing with domain wall encoding.
The generic Hamiltonian to be implemented by QA is where σz i and σx i are the Pauli operators acting at site i, J ij is the coupling constant on the edge (i, j) in the graph E describing the connectivity between qubits, h i is the local field at site i, and N represents the number of spins (qubits).The coefficients of the transverse field and the problem are denoted as A(s) and B(s), respectively, controlled by the time-dependent annealing parameter s = t/t a ∈ [0, 1].The initial state is the ground state of Eq. ( 4) with the coefficient −A(0)/2, and the system evolves from t = 0 to t = t a following the Schrödinger equation with the coefficients evolving as A(0) > 0 → A(1) = 0 and B(0) ≈ 0 → B(1) > 0. The precise forms of the scheduling functions of A(s) and B(s) depend on the physical devices [42].If the system evolves adiabatically, the final state of Eq. ( 2) approaches the ground state of Eq. ( 3).
Following Refs.[24][25][26], we use the domain-wall encoding method to rewrite a continuous variable optimization problem to a corresponding problem with Ising spins.In this method, we map the value of a continuous variable to the position of a domain wall within a spin chain.Concretely, we discretize a continuous variable x ∈ [x min , x max ] into N − 1 values as with the step size ∆x = and map each value x j to an N -spin state |x j ⟩ which has a domain wall at the jth bond: All spins at sites 1 to j are down, and the spins from j + 1 to N are up.
Using the fact that ⟨x j | (σ i+1 − σi )/2 |x j ⟩ = δ ij under this spin configuration, where δ ij is Kronecker delta, we can construct an N -spin Hamiltonian H RF such that its diagonal element in the basis |x j ⟩ corresponds to V RF (x j ) as To enforce that the ground state is a single domain wall state, we use a penalty Hamiltonian that punishes multi-kink and no-kink states as where the first term with J > 0 is the ferromagnetic coupling to align neighboring spins and the second term with h > J is the magnetic field at both ends of the chain to fix the spin at the first site in the down state and the spin at the N th site in the up state.The degenerate ground states of H DW are |x j ⟩ for j = 1, . . ., N − 1 and the penalty energy of each additional kink is 2J.Combining these two Hamiltonians, the domain wall encoding Hamiltonian for optimizing V RF (x) is given as where λ is a parameter that controls the relative strength of the penalty term H DW .For sufficiently small ∆x and λ, the ground state of H 0 corresponds to the global minimum of V RF (x).If λ is too large, the ground state of H 0 is no longer a single domain wall state.
Assuming that the potential function is differentiable and ∆x ≪ 1, Eq. ( 7) can be reduced to The above construction of the Ising Hamiltonian for continuous function optimization can be easily generalized to higher-dimensional problems [26].

III. COMPARISON OF QA WITH CLASSICAL ALGORITHMS FOR CONTINUOUS OPTIMIZATION
In this section, we show the results of benchmarks for QA and several classical algorithms designed for the optimization of continuous variable functions to investigate how QA applied to continuous variable optimization compares with dedicated classical algorithms.In the present section, QA means the hardware implementation, the D-Wave 2000Q.

A. Optimization algorithms
We run QA and the four classical algorithms for the optimization of continuous-variable functions: Nelder-Mead (NM) [33], conjugate gradient descent (CGD) [34], basin hopping (BH) [35], and differential evolution (DE) [36], which are easily performed on the classical computer and are commonly used for the optimization of continuous functions with a complex potential landscape.While the NM and CGD methods are for local optimization, BH and DE are for global optimization.The NM method searches for minimum by updating the polyhedron with the observed functional value.The CGD method utilizes the gradient of the potential function with the adjustable step size to find the minimum.Basin hopping performs two steps iteratively: local optimization around the proposed point and perturbation to search for the neighborhood domain of the current solution.Differential evolution is a kind of evolutionary computation and updates a proposed point by following the evolutionary procedure.Details of the classical algorithms are explained in Appendix A. Domain-wall encoding is applied only to QA on the D-Wave 2000Q, and classical algorithms use the continuous variable directly.
We apply these algorithms to Eq. ( 1) using the Scipy package [43].All numerical simulations for the classical optimization algorithms are performed on a laptop with Apple M1 Max chip with 10 CPU cores, which is one of the fastest chips [44].For QA, we optimize Eq. ( 9) with the D-Wave 2000Q, more precisely the "DW 2000Q 6" device.We have also run the same experiments on two types of the more recent D-Wave Advantage.Since "DW 2000Q 6" yielded the best results, we utilize it in the main text.The comparison of performance among a few models of D-Wave devices is described in the Appendix B 1.

B. Measured quantities
We compute two observables.The first is the ground state probability, where #GS is the number of cases where the global minimum was obtained by each algorithm, and n init stands for the number of trials (the number of different initial conditions).The second is the absolute error in the energy, where the brackets represent the average over the different initial conditions and V 0 = 0 is the minimum value of Eq. (1).

C. Parameter settings
For the classical optimization algorithms, we set n init = 10 3 and fix the search space as x ∈ [−3, 3] [45].For a fair comparison with QA, we regard the output x obtained by the classical optimization algorithms as the solution to Eq. (1) if the condition |x| < ∆x is satisfied.We control the execution time t ex by changing the maximum number of iterations of classical algorithms as t max = 2 i with the integer i running from 0 to 7. The execution time can be computed as t ex = tmax k=1 τ k where τ k represents the kth iteration time.
For the NM and CGD methods, the default parameter settings in Scipy are used.For BH, we apply the CGD method in the local optimization process.For DE, the number of populations is set to 10 and 120.Each algorithm stops on the way if it satisfies the stopping criterion as |x t − x t−1 | < 10 −7 (t = 1, . . ., t max ).
Our experimental settings are carefully designed for a fair comparison of QA and the classical algorithms in execution time and performance.Although there may be room for further improvements in execution time and performance, our studies are expected to lead to a general understanding of the features of QA for continuous function optimization against classical optimization algorithms.
For QA, we set the system size N = 211 and fix the coefficients as J = 1, h = 2, and λ = 1.The value N = 211, which corresponds to ∆x = 0.0286 with −x min = x max = 3, is large enough to represent the Rastrigin function with sufficient precision in its hills and valleys.The gradient of the potential in Eq. ( 1) is equal to 0 at the global minimum x = 0 and yields the vanishing of the magnetic field h j ′ = 0 where j ′ represents the place of the zero magnetic field.The vanishing of the magnetic field according to Eq. (11) gives two degenerate ground states of Eq. ( 9) as σz ± = (−1, . . ., −1, σz j ′ = ±1, 1, . . ., 1).The ground state σz − yields x j ′ = 0 with Eq. ( 5).On the other hand, the global minimum of Eq. (1) can not be recovered from the ground state σz + [46].Since the problem is negligible when N is large enough, we regard the two degenerate ground states of Eq. ( 9) as the equally qualified solution of Eq. (1).
We operate the execution time of QA by varying the annealing time 1(µs) ≤ t a ≤ 500(µs).The total execution time of QA is defined as t ex = t p + t s where t p is programming time and t s is sampling time which includes annealing time and readout time.The t ex can be obtained from "qpu access time" in the D-Wave Ocean SDK [42].To compute the observables, we repeat the annealing process 20 times independently and sample n reads = 10 3 spin configurations for each run.To reduce control errors, we apply flux bias calibration and gauge transformation [47,48].The flux bias is estimated before each run.A gauge is generated from a uniform distribution and fixed in each independent run.The ground-state probability is calculated from n reads = 10 3 spin configurations.The absolute error in the energy is computed from samples that meet the single-domain-wall condition.

D. Result for high energy barrier
We first consider the problem with a high energy barrier h 0 = 1.0.We show the dependence of the ground state probability and the absolute error in the energy on the execution time in Fig. 2. Data points denoted as "D-Wave (mean)" are the averaged result for 20 independent runs.Data with the maximum ground state probability among all runs are denoted as "D-Wave (best)".The error bars show the standard error computed by bootstrapping.Execution time t ex is always expressed in seconds.
Figure 2 (a) shows that around the region 0.7(sec) ≤ t ex ≤ 0.9(sec), the ground state probability obtained by "D-Wave (best)" can be larger than the other algorithms except for "DE (popsize=120)".The data for the NM and CGD methods are not shown in the left-side range beyond each left-most point because they failed to find the solution in this short-time region.Also, the data for the NM and CGD methods as well as for QA are not shown to the right of their respective rightmost points since they show saturation and no improvements in these long-time region [49].On the other hand, BH and DE find the solution in the longer time region 10(sec) ≤ t ex and their performance improves as the execution time increases.In the shorter time region, the minimum execution time for BH and DE cannot be reduced further than shown in Fig. 2 with our computational resources.
Figure 2 (b) shows that the absolute error in the energy obtained by QA is smaller than that of the NM and CGD methods and BH around the region 0.4(sec) ≤ t ex ≤ 0.9(sec).Therefore, QA can find the local minima closer to the global minima than the NM and CGD methods and BH in this region.The outputs of BH and DE with large populations concentrate around the global minima as we increase the execution time as we see in the right half of Fig. 2 (b).If the number of populations is small, DE falls into local minima even in the long execution time region.Similarly to the ground-state probability, the absolute error in the energy of QA saturates at a finite value.
For the higher energy barrier h 0 = 3.0, the qualitative behaviors of the classical algorithms are similar to those of h 0 = 1.0.On the other hand, the performance of QA deteriorates significantly, as will be described in Sec.IV A.

E. Result for low energy barrier
Next, we consider the problem with a lower energy barrier h 0 = 0.2. Figure 3 (a) shows an overall improvement of classical algorithms compared to the case of a higher energy barrier in Fig. 2 (a).Nevertheless, the performance of QA deteriorates slightly.The best result from QA does not reach the level obtained by BH and DE, even with a minimum execution time of about 0.8 seconds.The saturation of the ground-state probability of QA for a longer execution time is observed beyond about 0.9 seconds as in Fig. 2 (a) and is thus not drawn explicitly in Fig. 3 (a).Figure 3 (b) demonstrates that the absolute error in the energy generally decreases compared to the case of a high-energy barrier in Fig. 2 (b).The lowest value of the absolute error in the energy for the NM method is seen to be slightly lower than that of the CGD method in contrast to the case of Fig. 2 (b).This result implies that the one-dimensional structure with unique global minima fits the search strategy of the NM method.Except for this observation, the behaviors of the absolute error in the energy remain similar to those in Fig. 2 (b).
In summary, the performance of QA is slightly better than those of local optimization algorithms and is comparable to that from global algorithms in a limited region of execution time for the potential function with a higher energy barrier.On the other hand, if the energy barrier is low, we did not observe the clear benefits of QA over classical algorithms for the present problem of continuous function optimization.Whereas increasing the execution time enhances the performance of the global optimization algorithms, QA shows saturation beyond intermediate time scales.In the next section, we investigate whether this saturation originates in the problem in hardware implementation or in other reasons including the mapping of the continuous variable to Ising spins by the domain-wall encoding.

IV. COMPARISON OF QA WITH ALGORITHMS FOR DISCRETE-VARIABLE OPTIMIZATION
In this section, we compare the performance of the D-Wave 2000Q with several classical discrete-variable optimization algorithms, as well as with TEBD, a direct implementation of noise-free coherent QA on a classical computer, to examine how the D-Wave 2000Q compares with these other methods on the same footing of domainwall coding.The results may help us estimate how the hardware would perform if the noise level were significantly reduced.The classical algorithms used in this section are SA, SQA, and SVMC.We carefully distinguish between coherent QA simulated by TEBD and hardware implementation of QA, the latter being represented by the product name.For the sake of simplicity, the D-Wave 2000Q is referred to as the D-Wave hereafter.
We compute four quantities.The first one is the kink The second is the probability P const that the result satisfies the condition of single domain wall, and the third one is the residual energy Here, the brackets ⟨• • • ⟩ denote the sample average and E 0 is the ground state energy of Eq. ( 9).The last quantity is the ground state probability P GS that the ground state is successfully found.The kink density ρ represents the number of defects in spin chains.When a single domain wall exists, ρ = 1/N holds.Unlike the absolute error E abs used in Figs. 2 (b) and 3 (b), the residual energy E res is calculated from all samples, including the samples which do not satisfy the single domain wall condition.To understand the behavior of the entire system, we do not adopt E abs but use E res in this section.
Since the execution time of the D-Wave contains the measurement, initialization and other auxiliaryprocessing times, we mainly use the pure annealing time as the x axis of each figure for fair comparison.The experimental setting of the D-Wave is the same as those in Figs. 2 and 3.

A. D-Wave
Figure 4 shows the dependence of the observables obtained from the D-Wave on the annealing time for h 0 = 0.2, 1.0 and 3.0.The error bars represent the standard error obtained by bootstrapping and computed in 2 × 10 4 samples for ρ and E res and in 20 independent runs for P const and P GS .All error bars in figures of this section are computed in this way.The data for h 0 = 0.2 and h 0 = 1.0 are the same as those used in Figs. 2 and 3.

Kink density
Let us first look at the kink density ρ in Fig. 4 (a).We see that it needs some time that the data reaches ρ = 1/N , where the condition of single domain wall is satisfied.In other words, this is the minimum amount of time it takes to find the correct solution on average.If the kink density behaves as ρ ∝ t −b a below this threshold time, the annealing time required for the kink density to reach ρ = 1/N scales as t a ∝ N 1/b .Therefore, the value of the exponent b is important to understand how the performance of the domain wall encoding depends on N , and a larger b is preferable.Except for the case of h 0 = 3.0, Fig. 4 (a) demonstrates that the kink density gradually decreases as we increase the annealing time.I.A smaller value of the exponent reflects the effect of the higher energy barrier as seen in Table I (a).In addition, the higher energy barrier gives a larger absolute value of ρ, not just a smaller exponent.For h 0 = 3.0, the kink density behaves nonmonotonically.Such behavior has also been observed in several experiments on the D-Wave devices and is believed to be caused by the open-system nature of the hardware [48,50].
Kink densities for h 0 = 1.0 and 3.0 do not reach ρ = 1/N and saturate at finite values for long annealing time.Since the behaviors of all observables in Fig. 4 do not change even though we increase the annealing time beyond t a = 500(µs), we only show the data in 1(µs) ≤ t a ≤ 500(µs).

Single domain wall condition
Figure 4 (b) illustrates that almost all samples satisfy the single domain wall condition for h 0 = 0.2 beyond t a = 10(µs).On the other hand, the probability of satisfying the single domain wall condition saturates around P const ≈ 0.8 for h 0 = 1.0 and P const ≈ 0.1 for h 0 = 3.0.For those cases with higher energy barriers, the number of samples that satisfy the single domain wall condition deteriorates, and an excessive number of domain walls appear.region 1(µs) ≤ t a ≤ 10(µs) except for h 0 = 3.0.As with the kink density, the exponent of E res becomes smaller for higher energy barriers, see Table I (b).Higher energy barriers yield larger absolute values of the residual energy in addition to smaller exponents.

Residual energy
Next, we investigate the behaviors of E res in the large annealing time region t a ≥ 100(µs).The scaling analysis in this region gives us a hint on the performance of the D-Wave as an optimization solver.We omit the data for h 0 = 3.0 because few samples satisfy the single domain wall condition.We fit the polynomial function and show the exponents in Table II.For the problem with h 0 = 1.0, the slope is slightly larger than in the case of h 0 = 0.2.This result indicates interestingly that the process is more efficient in the problem with the higher energy barrier.
Although the residual energy has a polynomial decay as found in Koh and Nishimori [22], the exponent of E res is not close their result E res ∼ t −2 a .This discrepancy may come from the effects of noise on the hardware.In Fig. 4, we omit the data for t a ≥ 500(µs) because the behavior of E res hardly changes, and the effect of noise increases in the longer annealing time.This result corresponds to the fact in Figs. 2 and 3 that the observables hardly show improvements if we increase the execution time beyond a certain value and the data are omitted over t ex ≥ 0.9(sec).Therefore, the execution time region t ex ≥ 0.9(sec) corresponds to the annealing time region t a ≥ 500(µs).

Ground state probability
The ground state probability for h 0 = 1.0 is slightly larger than for h 0 = 0.2 but remains small for h 0 = 3.0 as seen in Fig. 4 (d).For h 0 = 0.2 and 1.0, the residual energy shows slightly decreasing behavior even in the large annealing time region t a ≥ 100(µs), but the ground state probability saturates around P GS ≈ 0.1.This result indicates that the transition to the ground states hardly appears, whereas the transition to the neighborhoods of the ground states occurs in this time scale.The ground state probability for h 0 = 3.0 is suppressed due to the existence of many defects with few samples satisfying the single domain wall condition as shown in Figs. 4 (a) and (b).

Kinks in low-energy states
To understand the reason for the significant performance deterioration observed for h 0 = 3.0, we examine the behavior of the low energy states of the classical Hamiltonian H 0 in Eq. ( 9).To evaluate how many singlekink states are properly encoded as low energy states of H 0 , we introduce n enc defined as where n kink (S) denotes the number of kinks in the spin state S, and S j is the spin state corresponding to the jth energy level.The value of n enc represents the index of the lowest energy level for which the number of kinks in the corresponding spin state is not equal to 1.If n enc > 0, the ground state is a single kink state.If n enc = 0 on the other hand, the ground state is a multi-kink state, indicating a failure of the encoding.The value of n enc depends on the relative magnitude of λH RF and H DW .
If we set λ to be sufficiently small relative to J, we can make n enc > 0.
Figure 5 shows the h 0 dependence of n enc obtained by dynamic programming [51].As h 0 increases, the values of the magnetic field h i determined by Eq. ( 17) become larger, resulting in a smaller n enc .For h 0 = 0.2 and 1.0, we have n enc = 196 and 160, respectively, indicating that most of the single kink states are encoded as low energy states of H 0 .For h 0 = 3.0, we have n enc = 26, indicating that the lowest energy states of H 0 are contaminated with multi-kink states.The reason for the extremely low constraint satisfaction probability P const for h 0 = 3.0 in Fig. 4 (b) is considered to be that the low energy multi-kink states trap the annealing process.Reducing λ increases n enc , which leads to an increase in P const .However, reducing λ also increases the thermal effect, and therefore, the optimization performance may not necessarily improve.The dependence of the observables on λ is shown in Appendix.B 2.

Summary
The values of observables obtained by the D-Wave depend on the height of the energy barrier.From the viewpoint as an optimization solver, achievement of the condition P const = 1 is important.For the problem with the higher energy barrier, the probability of satisfying the single-domain wall condition does not reach 1 in our experiments even if we increase the annealing time.The ground state probability also saturates at a finite value and does improve any further.Since the decoherence effects by thermal fluctuations may be dominant in the long annealing time region, these effects hamper the transition to the ground states and degrade the performance.Although the residual energy E res shows a polynomial decay in this region, the exponents are different from the direct simulation of the Schrödinger dynamics [22] due to the effects of thermal noise.

B. TEBD for coherent QA
In order to understand the origin of the suboptimal performance of the D-Wave described in the previous subsection, we now present the results of coherent QA based on TEBD to elucidate the differences between the D-Wave and the ideal QA that follows the Schrödinger dynamics without thermal noise.The TEBD algorithm is based on the Suzuki-Trotter decomposition [40] of the time evolution operator and the matrix product state (MPS) representation to approximately follow the time evolution of the wave function [52,53].Time and energy are measured in the unit of J in Eq. ( 8) and we adopt the natural unit where ℏ = 1.In this unit, if we choose B(1)/2 for the D-Wave as J, then t a = 1 corresponds to approximately 0.026 ns.The time step size of the Suzuki-Trotter decomposition is set to ∆t = 0.025/J.The maximum bond dimension of the MPS, which controls the accuracy of the approximation, is 32.We confirmed that the results converged well with respect to these parameters.The annealing schedules are the linear function, as A(s)/2 = 1 − s and B(s)/2 = s.
Figure 6 illustrates the results obtained by TEBD.A significant difference from the experimental results by the D-Wave in Fig. 4 is that the behavior of the kink density ρ exhibits only very slight dependence on the height of the energy barrier h 0 , which may originate in the weak dependence of coherent quantum tunneling effects on the barrier height.Moreover, in the case of coherent QA, for t a ≥ 10 3 , the single domain wall condition is almost completely satisfied, and the residual energy begins to decay quite rapidly, while in the case of the D-Wave in Fig. 4, the single domain wall condition is not satisfied even in the long time regime, and the residual energy does not decay.These results suggest that deviations from the ideal Schrödinger dynamics in the D-Wave lead to low performance, particularly for high h 0 .
In Figs. 6 (a) and (c), the kink density and the residual energy have the polynomial decay in the region t a ≤ 100, and the fitting exponents are summarized in Table I.In the limit where N is large and λ is small, the theory of the Kibble-Zurek mechanism [54][55][56]  In the simulation for Fig. 6, λ is set to 1, and therefore the magnitude of the longitudinal magnetic field is not negligible, see Eqs. ( 7) and ( 9), but the exponents show relatively mild dependence on h 0 .
As shown in Figs. 6 (b) and 6 (d), for t a > 100, the constraint satisfaction probability P const starts to increase from 0 and the ground state probability P GS also begins to increase.In this region, the residual energy deviates from the polynomial behavior and decreases more rapidly.
In a longer annealing time region (t a > 2 × 10 3 ), where the single domain wall condition is completely satisfied, P GS approaches 1.In this regime, the energy expectation value can be expressed as (17) and therefore the residual energy is proportional to 1 − P GS as We can confirm this behavior by plotting the relation between E res and 1 − P GS in Fig. 7 (a), as the slope is approximately 1 in the region where 1 − P GS < 0.1.
The annealing time dependence of 1 − P GS is shown in Fig. 7 (b), and in the region where 1 − P GS < 0.1, it can be seen that 1 − P GS follows an exponential decay as 1 − P GS ∝ exp(−Ct a ), which is reminiscent of the Landau-Zener formula.Consequently, the annealing time dependence of the residual energy can be expressed as E res = C ′ exp(−Ct a ), where C ′ and C are positive constants.The solid lines in Fig. 7 (c) represent the result of fitting this expression to the data in the range of 3 × 10 3 ≤ t a ≤ 10 4 .In Ref. [22], an analysis of QA by direct numerical solution of the time-dependent Schrödinger equation for a single particle system in a continuous space showed that for a linear annealing schedule, the residual energy decays as the inverse square of the annealing time beyond the Landau-Zener regime as suggested in Ref. [57].The deviation of the residual energy for h 0 = 0.2 from the exponential behavior in the longer-time region t a > 10 4 may be an indication of a transition to such a power law behavior, although the exponent is much larger than two.To confirm the existence of this transition and to determine the power law exponents, simulations for much longer annealing times are necessary, which is challenging due to the accumulation of numerical errors in TEBD.
To understand the effect of the height of the energy barrier on the coherent QA, we analyze the spatial distribution of kinks in the wave function after the process of QA has terminated.Figure 8 shows the local number of kinks (1−⟨σ z i σz i+1 ⟩)/2 at t a = 5×10 4 as a function of spatial coordinate, as well as the values of the ground state probability P GS and residual energy E res .For a small energy barrier h 0 = 0.2, kinks are distributed around the global minimum of the potential function at x = 0 without significantly avoiding the maxima.This result suggests that for the low energy barrier, the kink is more likely to move around, making it difficult to remain at the bottom of the potential minima.On the other hand, for a large energy barrier h 0 = 3.0, kinks localize in the vicinity of the minima and hardly exists in the vicinity of the maxima due to the high energy barrier.The increase in the population at the local (but not global) minima for h 0 = 3.0 leads to larger residual energy than that of h 0 = 0.2.However, the probability of obtaining the ground state P GS is higher when h 0 = 3.0.This behavior suggests that a lower energy barrier case is not necessarily easier to optimize, contrary to naive intuition.The characteristic of kink spatial distribution discussed here becomes indistinct when λ is larger, but it can still be observed.Indeed, for λ = 1, the optimization performance does not decrease even in the case of high energy barrier as shown in Fig. 6.
In summary, we have found that the optimization performance of coherent QA realized by TEBD is almost independent of the height of energy barrier h 0 , and the system converges to the ground state with high probability for sufficiently long annealing times.Therefore, the significant dependence of performance on h 0 and the premature saturation of the ground state probability in the data of the D-Wave can be attributed to imperfect realization of QA in the hardware, possibly mainly by thermal noise.As discussed in Appendix C, in order to fully understand the thermal effects on the hardware, it is necessary to analyze an open quantum system, not just an isolated system as done in this study.The performance of QA with a more coherent device would approach the results of TEBD simulations presented in this subsection and can be expected to be equivalent to or even better than that of global optimization algorithms such as BH and DE.Additional evidence will be provided in the following subsections to show that the h 0 dependence of data originates in thermal effects.

C. Simulated annealing
Next, we consider SA to investigate the effects of thermal fluctuations.In SA, we use a single-spin update with the Metropolis rule.The temperature is linearly decreased as with T 0 = 1 and T 1 = 10 −5 .The maximum number of Monte Carlo steps (MCS) is denoted as t MCS .The number of spin updates is t MCS N .We repeat SA 20 times independently and sample 10 3 spin configurations for each independent run.

Kink density
Figure 9 shows the results obtained by SA.As in Fig. 4 (a), Fig. 9 (a) indicates the power law behaviors of the kink density in the short annealing time region 2 ≤ t MCS ≤ 100.Unlike the results obtained from the D-Wave in Fig. 4 (a), the power law behaviors can be seen even in the case of high energy barrier, h 0 = 3.0.For h 0 = 0.2, the exponents are larger than those of the results obtained by the D-Wave and the coherent QA, see Table I.On the other hand, for larger h 0 with higher energy barrier, the exponents are larger than those of the D-Wave and smaller than those of the coherent QA for the problem.In this time region, the values of the kink density become larger than those of the D-Wave and take similar values to the case of coherent QA except for h 0 = 3.0.For h 0 = 3.0, the kink density does not reach ρ = 1/N even at t MCS = 10 5 .We need more annealing time to reduce the defects for the problem with the larger value of h 0 .

Single domain wall condition
Figure 9 (b) illustrates that few samples satisfy the single domain wall condition in the short-time region where the kink density has the polynomial decay.While this result is similar to the result obtained by the coherent QA in the short-time region as seen in Fig. 6 (b), a clear difference from the D-Wave in Fig. 4 (b) is observed.Since the system is not sufficiently equilibrated in the short annealing time region, many defects remain.For the D-Wave result in Fig. 4 (b), P const has a finite value except for h 0 = 3.0 in this region.Although this D-Wave result is similar to the behavior of SA in the intermediate annealing time region 10 2 ≤ t MCS ≤ 10 3 in Fig. 9 (b), it is impossible to reproduce all D-Wave data quantitatively by SA, for example, the residual energy as shown below.
In the long annealing time region t MCS ≥ 10 4 , almost all samples satisfy the single domain wall condition except for h 0 = 3.0.For h 0 = 3.0, P const is suppressed.Similarly to the D-Wave data in Fig. 4 (b), SA is also affected by the height of the energy barrier.

Residual energy
As with the kink density in Fig. 9 (a), Fig. 9 (c) shows that the residual energy scales polynomially in the short annealing time region 2 ≤ t MCS ≤ 100.The exponents are smaller than the D-Wave results as seen in Table I (b).The exponents are larger than those by the coherent QA for h 0 = 0.2 and smaller for h 0 = 1.0 and 3.0.These results imply that the dynamics based on thermal fluctuations are hampered by the higher energy barrier, and quantum fluctuations are more efficient than thermal fluctuations for higher barriers.
The residual energy of SA decreases polynomially in the long annealing time regions 2 × 10 3 ≤ t MCS ≤ 10 5 for h 0 = 0.2 and 10 4 ≤ t MCS ≤ 10 5 for h 0 = 1.0 where P const ≈ 1 holds.According to Table II, the slope of the fitting line becomes smaller for the problem with the higher energy barrier.The higher energy barrier impedes the transition to the ground states by thermal fluctuations.On the other hand, the exponents obtained by the D-Wave show the opposite trend.Note that P const = 1 does not hold for h 0 = 1.0 in Fig. 4 (b) for the D-Wave.
Whether or not this is related to the quantumness of the D-Wave is not clear at this point.Although the exponents of SA are better than those of the D-Wave, the inverse square scaling of E res is not observed.

Ground state probability
Figure 9 (d) shows that the ground state probability starts to increase monotonically from t MCS = 10 2 .The performance depends on the height of the energy barrier.The higher the energy barrier is, the smaller the success probability becomes.The behavior is different from that in Figs. 4 (d) and 6 (d).Since the ground state probability is not close to 1 even in the long annealing time region, the exponential decay of E res observed in the coherent QA in Fig. 7 (c) does not appear in SA as seen in Fig. 9 (c).

Summary
The observables of SA clearly depend on h 0 while the observables of the coherent QA hardly depend on h 0 .This is a significant difference between SA and coherent QA.In coherent QA, quantum tunneling effects play an important role in the state transition between local minima.On the other hand, SA shows deteriorated performance under higher energy barriers.This behavior is consistent with the physical picture that SA goes over the energy barrier with thermal fluctuations.Also in the case of the D-Wave, dependence of the observables on h 0 can be clearly seen.Premature saturation of P const and P GS at lower values appear in the long annealing time region for the D-Wave.Such saturation of these observables does not exist in SA.This suggests that the degraded performance of the D-Wave may not be due to thermal noise alone, but to other factors as well [58,59].

D. Simulated quantum annealing
Next, we study SQA and see its difference from SA and other protocols in some detail.Simulated quantum annealing is the algorithm to simulate some aspects QA on a classical computer.In SQA, we apply the Suzuki-Trotter decomposition [40] to the partition function of Eq. ( 2).The quantum system at a finite temperature is then mapped to the effective classical Ising model as where M is the Trotter number, β = 1/T is the inverse temperature, and σ ik ∈ {±1}(i = 1, . . ., N, k = 1, . . ., M ) is the classical Ising configuration.The coefficients of the problem and transverse field are linearly controlled as where t = 0, . . ., t MCS − 1.We set the effective temperature β/M = 1 and M = 10 3 .We utilize a single-spin update with the Metropolis rule.We sweep the spin configurations M N times for each MCS.We perform SQA 20 times independently.As shown in Ref. [60], the final spin configurations in Trotter replicas are utilized for computing the observables.

Kink density
Figure 10 shows the results obtained from SQA.In Fig. 10 (a), the kink density does not depend on the value of h 0 in the shortest annealing time region 2 ≤ t MCS ≤ 3.This behavior is similar to that of the coherent QA in Fig. 6 (a).In the short annealing time region 4 ≤ t MCS ≤ 100, the polynomial decay of ρ can be seen.Except for h 0 = 3.0, SQA has the largest exponent in annealing protocols for discrete-value optimization, as shown in Table I.The exponents of ρ become smaller as the height of the energy barrier increases.The tendency is shared by other annealing protocols except for the coherent QA.For h 0 = 1.0 and 3.0, SQA requires less annealing time to satisfy the single domain wall condition as compared to SA, and thus classically-simulated quantum fluctuations are more effective than thermal fluctuations in the higher barrier case.This result reflects the difference in the dynamics of both methods and can also be understood from the behaviors of P const in Fig. 10 (b).For h 0 = 0.2, we cannot observe a clear difference between SA and SQA.Compared with the coherent QA, the dependence of ρ on h 0 appears.The dynamics of SQA is fundamentally different from the dynamics of coherent QA [61] and may be regarded to be rather closer to the D-Wave behavior and SA.(b).The conditionP const ≈ 1 holds in the long annealing time region t MCS ≥ 10 4 including the case of h 0 = 3.0.This behavior is different from the D-Wave and SA in Figs. 4 and 9.No suppression of P const for h 0 = 3.0 is observed for SQA, which is similar to coherent QA.

Residual energy
Figure 10 (c) illustrates that the residual energy decreases polynomially in the short annealing time region.The range of fitting is the same as that in Fig. 10 (a) to extract the values in Table I.As seen in this Table, as the height of the energy barrier increases, the exponents of SQA decrease.The exponents of SQA are larger than those of SA.The result indicates that SQA can solve the problem more efficiently than SA.For h 0 = 3.0, the exponent of SQA is smaller than that of the coherent QA.Although SQA can simulate some aspects of quantum effects, the higher energy barrier degrades the performance, unlike the coherent QA.In the long annealing time region 10 3 ≤ t MCS ≤ 10 4 , the polynomial decay of E res can be observed.As seen in Table II, SQA can solve the problem more efficiently than the D-Wave and SA since SQA gives larger values of the exponent than the D-Wave and SA.In our experiments, the exponential decay of E res as in coherent quantum annealing is not observed in SQA.To see the behavior, we may need to investigate longer annealing time regions where P GS ≈ 1 holds.

Ground state probability
Figure 10 (d) shows that the success probability becomes smaller as the height of the energy barrier increases in the large annealing time region t MCS ≥ 10 4 .Since the values of P GS obtained by SQA are larger than those by SA in the region t MCS ≥ 2 × 10 3 , the effect of the energy barrier is weakened by classically-simulated quantum fluctuations.

Summary
The observables obtained by SQA depend on the value of h 0 , which is different from the case of coherent QA.This attributes to the fact that the dynamics of SQA does not follow the Schrödinger dynamics but is based on classical Monte Carlo updates under thermal fluctuations.Compared with SA, SQA suffers less performance degradation under the higher energy barrier, possibly by partly simulating tunneling effects [62].

E. Spin-vector Monte Carlo
Spin-vector Monte Carlo is a classical algorithm proposed to explain the outputs from the D-Wave device [41].In SVMC, a quantum spin is replaced by a classical rotor of unit length with a continuous angle variable The total Hamiltonian can then be described as : The annealing schedules of A(s)/2 and B(s)/2 are the same as those in SQA.The temperature in Monte Carlo process is set in our simulation to T = 10 −5 and the number of all spin updates is written as N t MCS .Each continuous angle is initialized to θ i = π/2, and spins are randomly chosen and updated with the Metropolis rule as In addition, we apply the transversefield-dependent (TFD) update [63,64], which is proposed to capture the freezing phenomena in the D-Wave device.In the TFD update, the current angle is updated as θ i → θ i + ϵ i (s) where ϵ i (s) is the random variable dependent on the annealing parameter and the energy coefficients, The range of updating the angle gradually decreases in the late annealing process A(s) < B(s).We execute SVMC 10 times independently and sample 10 3 spin configurations for each run.

Kink density
Figure 11 shows the results of SVMC.In Fig. 11 (a), nonlinear behaviors of ρ are observed.In the region of very short annealing time 2 ≤ t MCS ≤ 3, the kink density does not depend on the values of h 0 .Dependence of ρ on the values of h 0 appears for t MCS ≥ 3.For h 0 = 3.0, the kink density saturates at a larger value than that of SA in the region t MCS ≥ 50.This result represents that SVMC is more affected by the high energy barrier than SA.1.0.Compared with other annealing protocols, SVMC needs more annealing time to satisfy P const ≈ 1.The result indicates that the dynamics of SVMC is slower than that of other annealing protocols.

Residual energy
As in ρ in Fig. 11 (a), Fig. 11 (c) illustrates that nonlinear behavior of E res emerges and the saturation of E res for h 0 = 3.0 appears.In the long annealing time region, the polynomial decay of E res can be observed.The fitting regions are 5 × 10 4 ≤ t MCS ≤ 5 × 10 5 for h 0 = 0.2 and 10 5 ≤ t MCS ≤ 5 × 10 5 for h 0 = 1.0.In these regions, P const ≈ 1 holds as seen in Fig. 11 (b).From Table II, the exponents of SVMC are the largest among all annealing protocols.Therefore, SVMC can reach the ground states faster than other annealing protocols in the long time region.Compared with SQA in the region where the single domain wall condition holds, the values of E res obtained by SVMC are smaller for h 0 = 0.2 and 1.0.

Ground state probability
Figure 11 (d) demonstrates that the ground state probability starts to increase more slowly than SA and SQA.The higher energy barrier suppresses the increase of P GS .At a later stage, SVMC produces a larger value of P GS than SQA does for h 0 = 0.2 and 1.0 in the region where P const ≈ 1 holds.

Summary
Values of observables in SVMC depend on the value of h 0 as in SA and SQA.The dynamics of SVMC is more affected by the higher energy barrier.This is a consequence of the fact that the dynamics of SVMC is activated by thermal fluctuations in the semi-classical potential function.Although the behaviors of the observables are unique to SVMC and are different from those obtained by the D-Wave, the performance of SVMC is the best in the region where P const ≈ 1 holds except for the coherent QA.This result suggests that the dynamics of SVMC is more effective than other classical annealing protocols for the Rastrigin function in the space where only a single domain wall exists.
Considering SVMC as a practical solver for continuous function optimization, a smaller absolute value of the execution time is required, not just larger scaling exponents.In the next section, we analyze the dependence of observables on the absolute execution time.

F. Actual execution time
Finally, we show the dependence of observables on the actual execution time in seconds for different annealing protocols in Fig. 12.We use the same data as those presented in this Sec.IV for k = 0.5, w 0 = 0.2, and h 0 = 1.0.The definition of t ex is the same as that defined in Sec.III B. Since the TEBD is the deterministic algorithm and different from the other annealing protocols based on the Monte Carlo method, we avoid to show the TEBD 10 3 10 2 10 1 10 0 10 1 10 2 10 3 10 4 t ex (s) data for direct comparison in this section.Figure 12 (a) illustrates that SQA reaches ρ = 1/N most quickly among different annealing protocols [65].Although the D-Wave reaches the lowest value of ρ around 0.4(sec) ≤ t ex ≤ 0.9(sec), its value of ρ remains suboptimal in the longer time region due to the noise and other imperfections as observed more clearly in Fig. 4 (a).
Around the moderately-short time region 0.4(sec) ≤ t ex ≤ 0.9(sec), P const reaches close to 1 by the D-Wave as shown in Fig. 12 (b).If the goal is to satisfy the condition of single domain wall as quickly as possible in the actual time of execution, the D-Wave is the choice.
Figure 12 (c) shows that the D-Wave has the smallest value of E res around 0.4(sec) ≤ t ex ≤ 0.9(sec).Around 2(sec) ≤ t ex ≤ 370(sec), SQA performs best.From the analysis of the exponents in Table II, SVMC surpasses SQA in the long time region, and the difference of SQA and SVMC widens with time.To satisfy the single domain wall condition, SVMC needs more time, see Figs. 12 (a) and (b).However, once the condition is satisfied, the residual energy decreases most rapidly by SVMC.
Figure 12 (d) indicates that the success probability of the D-Wave has the largest value around 0.4(sec) ≤ t ex ≤ 0.9(sec), P GS ≈ 0.124 as indicated by the horizontal line.The D-Wave reaches this value most quickly among all annealing protocols.However, it shows saturation at around this finite value, and the performance hardly improves with time.In Fig. 2, we see that the ground state probability obtained by BH and DE are large P GS > 0.9 even in the longest time region 10(sec) ≤ t ex ≤ 100(sec).As far as our implementation is concerned, other annealing protocols, SA, SQA, SVMC, are similarly less effective than the global opti-mization algorithms for continuous-variable optimization since the former do not achieve such high values in this long-time region.
Overall, in the short time window of 0.4(sec) ≤ t ex ≤ 0.9(sec), the D-Wave is the best choice.At around the intermediate time t ex ≈ 2(sec), SQA is competitive.Beyond this time region, dedicated global optimization algorithms for continuous variables, DE and BH, keep improving with time, whereas classical metaheuristics for discrete optimization remain relatively mediocre.The D-Wave stays stagnant.It is nevertheless expected from the data of TEBD that the ideal QA without noise is likely to be most competitive, if realized at the hardware level, also in the intermediate and long-time regions.

V. CONCLUSION
We have performed extensive benchmarks of QA and classical optimization algorithms for the continuousvariable Rastrigin function in one dimension with a rugged energy landscape.By using the domain wall encoding, the continuous variable problem was represented by a one-dimensional Ising chain.
We have first compared the D-Wave with classical algorithms for continuous-variable optimization with local or global updates and have demonstrated quantitatively that, in the case of higher energy barrier, the performance of the D-Wave is better or comparable to that of the local optimization algorithms in a limited region of execution time and is comparable to the global optimization algorithms in the same time region.For longer execution times, the global algorithms show outstanding performance whereas the D-Wave remains stagnated.For the potential function with lower energy barrier, no clear advantage of the D-Wave has been observed.While the performance of the global optimization algorithms keeps improving with increasing execution time, no improvements were identified in the D-Wave.
For better understanding of the performance of QA, in its hardware realization of the D-Wave as well as in an ideal setting without noise, we next compared QA with protocols for discrete-variable optimization, TEBD for coherent QA, SA, SQA, and SVMC.The observed values of physical quantities show a clear dependence on the height of the energy barrier except for TEBD for coherent QA.This is natural for SA, SQA, and SVMC because these classical algorithms make explicit use of thermal noise to go over the barrier.In the case of the D-Wave, its barrier-height dependence implies that the process realized in the D-Wave is explicitly affected by thermal noise.Particularly in the case of the highest barrier with h 0 = 3.0, we found no sign of the successful solution by the D-Wave, whereas coherent QA by TEBD reached the ground state with high probability for relatively long annealing time with little-to-no dependence on the barrier height.This result suggests that the performance degradation of the D-Wave at long annealing time may be significantly reduced by improvements in the hardware to reduce noise and possibly other imperfections as demonstrated in recent experiments of quantum simulation [48,66].Stated otherwise, quantum annealing has the potential to become quite competitive even for continuous-variable problems under the domain wall encoding.This is achieved by quantum tunneling effects, as implied by the independence of the barrier height observed in the TEBD data.Classical algorithms, including those for continuous-variable optimization, would struggle to find a solution if the barrier become even higher than we tested, while QA, if realized coherently, would be much less affected by the height.This conclusion has become possible by the analysis of the continuous-variable problem with an explicit and easy control of the barrier height, which is generally difficult to achieve in problems expressed directly by the Ising model.In other words, we have demonstrated the effectiveness of coherent quantum tunneling effects explicitly and quantitatively in an optimization problem.Simulated quantum annealing may be regarded as a quantum-inspired classical algorithm that partly takes into account quantum effects.Nevertheless, comparison of data for TEBD in Fig. 6 and SQA in Fig. 10 suggests that a complete realization of quantum effects in QA is necessary for better results.
This paper represents a first step toward systematic and quantitative studies of continuous-variable optimization by QA in comparison with a series of well-established classical algorithms.In particular, the present onedimensional problem with a periodic structure of the barrier is one of the simplest examples of continuousvariable optimization.The global optimization algorithm DE based on the genetic algorithm may find it more diffi-cult to reach the correct solution for aperiodic potentials than for the periodic Rastrigin function.Also, higherdimensional problems with the complex energy landscape as studied qualitatively in Refs.[25][26][27], which are generally difficult to analyze by TEBD, may be the interesting next target of systematic studies.

Basin-Hopping
The Basin-hopping (BH) method is a widely used global optimization algorithm for nonlinear multimodal optimization problems, where the local optimization algorithm often fails to find the global minimum due to the rugged energy landscape [35].The BH method is combined with the gradient-based local optimization and the perturbation based on Monte Carlo random walks with the Metropolis rule [70,71].Starting from an initial point x (1) , we apply the local optimization to map a current point x (t) to the nearest minima y (t) in the local optimization process.Next, we add a perturbation to the current solution as x new = y (t) + ϵ where ϵ is a random coordinate generated from the uniform distribution U (a, b).We perform the local optimization once again and attain another local minimum y new .Following the Metropolis rule, we accept the proposed point as y (t+1) = y new or reject it as y (t+1) = y (t) with the probability, min 1, exp(−(f (y new ) − f (y (t) ))/T ) where T is the temperature to adjust the scale of the energy.The above iteration continues until t = t max .In our experiments, we utilize the default setting a = −0.5, b = 0.5, and T = 1, and adopt the CGD method in the local optimization process.

Differential evolution
The differential evolution (DE) method is a stochastic search algorithm for minimizing real-valued global optimization problems based on the evolutionary strategy [36,72].In DE, we initialize a population vector for each generation G as x G i ∈ R n (i = 1, . . ., p, G = 0, . . ., G max ), where p is the number of populations.The population vector is generated from the uniform distribution to cover the entire parameter space.Next, we create a new population vector with different ones as , where r 1 , r 2 , and r 3 are indices of the population determined by the "crossover" rule and F > 0 is the amplification factor.The new trial vector u G+1 i can be obtained by mixing the target vector x i and the new population vector xG+1 i with the "crossover" rule [73].If the new trial vector has a lower objective function value than the target vector, the new trial vector supersedes the target vector and is inherited by the next generation.The above process is called the "selection".The algorithm stops at G = G max or the mean objective function value over all populations becomes lower than the threshold ϵ.In our experiments, we adopt the "best/2/bin" strategy [36,72] and set ϵ = 10 −7 .The number of populations G max corresponds to the maximum number of iterations t max in the main text.The other hyperparameters are used in the default settings.
Figure 13 illustrates the dependence of observables on the annealing time for a few different quantum annealers.In addition to the "DW 2000Q 6" device, we utilize the "Advantage system 4.1", "Advantage system 6.1", and "Advantage2 prototype 1.1" devices and compare their performance.
Regardless of the type of annealers, the polynomial decay of ρ can be observed in the short annealing time region 1(µs) ≤ t a ≤ 10(µs) as shown in Fig. 13 (a).The fitting exponents are summarized in Table III.The largest exponent of ρ is obtained from the D-Wave 2000Q.Except for the D-Wave Advantage 6.1, the kink density gradually decreases as we increase the annealing time.The kink density obtained by the D-Wave Advantage 6.1 shows similar behavior to the D-Wave Advantage 4.1 in the region 1(µs) ≤ t a ≤ 10(µs) and starts to increase at t a = 30 (µs).This behavior may be due to decoherence.In the long annealing time region t a ≥ 100(µs), no device reaches the target ρ = 1/N .In this region, the D-Wave 2000Q takes similar values as the D-Wave Advantage2 Prototype 1.1 does.
Figure 13 (b) illustrates that the D-Wave 2000Q has the largest P const in the region 1(µs) ≤ t a ≤ 100(µs).The D-Wave Advantage2 Prototype 1.1 takes similar values of P const to those obtained by the D-Wave 2000Q in the region 100(µs) ≤ t a ≤ 500(µs).In the long annealing time region, most samples obtained by the above two devices satisfy the single domain wall condition.Both devices yield the larger value of P const than the D-Wave Advantage 4.1 and 6.1.
Figure 13 (c) shows that the residual energy decays polynomially irrespective of the types of the quantum annealer in the short annealing time region as also seen in Fig. 13 (a III, the performance gap appears in P GS .This result implies that the D-Wave Advantage2 Prototype 1.1 may reach low energy states, but not necessarily the ground states. The observable values obtained by different quantum annealers show different behaviors.We generally conclude that the D-Wave 2000Q is suitable for optimizing the Rastrigin function under the present problem setting.IV.In this region, the coefficient λ hardly affects the behavior of the system.On the other hand, in the long annealing time region t a ≥ 100(µs), the saturation seen in ρ for λ = 1.0 does not occur in λ = 0.1 and 0.5.The small λ weakens the effect of the problem term Eq. ( 10) and enhances the effect of the constraint term Eq. (8).
Figure 14 (b) demonstrates the value of λ hardly affects P const in the short annealing time region.Similarly to the kink density in Fig. 14 (a), the effect of λ can be seen in the long annealing time region t a ≥ 100(µs).Premature finite saturation of P const appears for λ = 1.0.The number of samples that satisfy the single domain wall Both axes are the same as those in Fig. 4. condition increases by adjusting λ.In the long annealing time region, the effect of quantum fluctuation becomes weaker due to the decoherence.In this region, thermal fluctuations affect the system.Since small values of λ largely reduce the height of the energy barrier and the lower energy barrier is favorable for thermal fluctuations as shown in Fig. 9, the number of samples that satisfy the single domain wall condition increases with decreasing λ in this region.Figure 14 (c) shows that the polynomial decay of E res can be observed in the short annealing time region 1(µs) ≤ t a ≤ 10(µs).From Table IV, the values of λ hardly affects the exponents of E res as in the kink density.In the long annealing time region t a ≥ 100(µs), the residual energy decreases polynomially and the fitting exponents are written in Table IV (b).A smaller λ yields a larger exponent.This result implies that a smaller λ leads to lower energy states and reduces the hardness of optimization due to the higher energy barrier.Figure 14 (d) illustrates that larger values of λ give larger P GS .This result is in contrast to the behavior of E res in Fig. 14 (c).Smaller values of λ decrease the values of the magnetic field physically encoded into the D-Wave.Since the bit width to represent the problem Hamiltonian is limited in the D-Wave, the actual implemented value in the D-Wave may be different.Due to this discretization error, the actual energy landscape slightly changes.Because the success probability is more affected by variation of the energy landscape than the residual energy, the success probability for the smaller λ is degraded.The D-Wave is a device to reproduce the system behavior under the time dependent Hamiltonian H(s) in Eq. ( 2).However, as the hardware is an open system, its dynamics is affected by thermal noise.To describe such dynamics, the hypothesis of quasi-static time evolution and freeze-out phenomenon has been proposed [74].In this hypothesis, during the annealing process, the density matrix of the system at a given time is approximately equal to the density matrix of the thermal equilibrium state of the Hamiltonian at that time.The dynamics freezes at a time denoted as s * , when the relaxation time to equilibrium becomes much larger than the time scale of annealing, and the density matrix of the final result becomes exp(−H(s * )/T phys ), where T phys denotes the physical device temperature that is about 13.5 mK.
This can be interpreted as the density matrix at the thermal equilibrium state of H 0 with a temperature of 2T phys /B(s * ).By determining this rescaled temperature from the output samples of the D-Wave, it is possible to estimate the freeze-out point s * .It should be noted that if the obtained effective temperature value is too high, the estimation of s * may become inaccurate and it would be appropriate to estimate the temperature using H(s * ) instead of H 0 .We now discuss the consistency between the experimental results of the D-Wave in this study and the hypothesis of quasi-static time evolution and freeze-out.The effective temperature T eff is defined as the value obtained by solving the equation E(T ) = ⟨H 0 ⟩ for T , where E(T ) is the internal energy of the thermal equilibrium state of H 0 at temperature T and ⟨H 0 ⟩ is the energy expectation value of the data of the D-Wave.The annealing time dependence of the effective temperature from the experimental results of the D-Wave is shown in Fig. 15.For h 0 = 0.2, the effective temperature decreases toward T phys while for h 0 = 1.0 and 3.0, the effective temperature saturates at a constant value several times higher than 2T phys /B(1).
Regarding T eff as 2T phys /B(s * ), the freeze-out point s * is estimated as the value obtained by solving the equation 2T phys /B(s) = T eff for s.The relation between s * and 1/T eff is shown in Fig. 16 and the annealing schedule function B(s) used to determine s * and the schedule function of the transverse field A(s) are also shown.For h 0 = 0.2, at t a = 500(µs) the freeze-out point s * reaches approximately 0.6, where B(s * ) ≫ A(s * ), indicating that the classical term is dominant in H(s * ).For To examine how closely the D-Wave data matches a thermal equilibrium value, we compare the data of the D-Wave and the observables of the thermal equilibrium of H 0 at the obtained effective temperature as shown in Fig. 17.For h 0 = 0.2 and for t a > 10, where s * > 0.5, the ground state probability P GS and the constraint satisfaction probability P const are close to the value of the thermal equilibrium state, implying that the hypothesis of quasi-static time evolution is valid.On the other hand, in the region where s * < 0.5 for all h 0 values, P GS of the D-Wave is larger than that of the thermal equilibrium.Furthermore, P const of the D-Wave is smaller than that of thermal equilibrium.This result suggests that in the case of s * < 0.5, the hypothesis of quasi-static time evolution may not hold, or that the final state of the annealing is affected by the transverse field.
A similar analysis of the effective temperature is performed for the results of TEBD, as shown in Fig. 18.As in the case of the D-Wave, P GS is larger, and P const is smaller than the value of the corresponding thermal equilibrium state.This result suggests that it may be generally observed that the single kink state probability is larger, and the ground state probability is smaller, in non-equilibrium states compared to the thermal equilibrium states at the corresponding effective temperature.
In summary, the annealing time dependence of the effective temperature varies significantly depending on h 0 .For h 0 = 0.2, the effective temperature decreases toward the hardware temperature while it saturates at higher temperatures for h 0 = 1.0 and h 0 = 3.0.Moreover, the value of effective temperature changes non-monotonically with h 0 .The mechanism of this nontrivial dependence of the effective temperature on h 0 remains to be clarified.In the region where the freeze-out point is small and not in thermal equilibrium of H 0 , the behaviors of P GS and P const exhibit similar characteristics to coherent QA, al-though the underlying mechanism remains unclear.For a deeper understanding of the results, analysis taking the effects of thermal noise into account is necessary.

FIG. 1 :
FIG. 1: Illustration of the Rastrigin function with k = 0.5 and w 0 = 0.2 for different values of h 0 .

FIG. 2 :
FIG. 2: Dependence of observables on the execution time for k = 0.5, w 0 = 0.2, and h 0 = 1.0.The vertical axes denote the observables: (a) the probability of ground state and (b) the absolute error in energy.The insets enlarge the region 0.25(sec) ≤ t ex ≤ 1(sec).

FIG. 3 :
FIG.3: Dependence of observables on the execution time for k = 0.5 and w 0 = h 0 = 0.2.Both axes are the same as those in Fig.2.The insets enlarge the region 0.25(sec) ≤ t ex ≤ 1(sec).

FIG. 4 :
FIG. 4: Dependence of the observables obtained by the D-Wave on the annealing time for three values of h 0 .The vertical axes are the observables: (a) the kink density, (b) the probability satisfying the single domain wall condition, (c) the residual energy, and (d) the ground state probability.The horizontal axes represent the annealing time.

Figure 4 (
Figure 4 (c) indicates that the residual energy polynomially decreases as E res ∝ t −c a in the short annealing time

FIG. 5 :
FIG. 5: Number of single kink states encoded in the lowest energy states of the classical Hamiltonian H 0 as a function of the height of energy barrier h 0 .The parameters of the Rastrigin function are set to k = 0.5 and w 0 = 0.2, and the weight of H RF is set to λ = 1.0.

FIG. 6 :
FIG.6: Results of TEBD simulations for different values of h 0 as functions of annealing time.The vertical axes are the same as those in Fig.4.The annealing time t a is in the unit of J and does not correspond to that of Fig.4.

FIG. 7 : 2 FIG. 8 :
FIG. 7: Results of TEBD simulations.(a) Relation between the residual energy and the deviation of the ground state probability from unity.The dashed line denotes (1 − P GS )(E 1 − E 0 ) where E 1 is the first excited energy of H 0 .(b) Deviation of the ground state probability from unity as a function of the annealing time.(c) Residual energy and its fitting curve as a function of the annealing time.The solid lines represent the curves of E res = C ′ exp(−Ct a ) fitted to the data in the range of 3 × 10 3 ≤ t a ≤ 10 4 .

FIG. 9 :
FIG. 9: Dependence of the observables obtained from SA on the annealing time for different values of h 0 .The vertical axes are the same as those in Fig. 4. The horizontal axes represent the Monte Carlo step.

FIG. 10 :
FIG. 10: Dependence of the observables obtained by SQA on the annealing time for different values of h 0 .Both axes are the same as those in Fig. 4.

Figure 10 (
Figure 10 (b) shows that few samples satisfy the single domain wall condition in the short-time region similarly to the cases of SA and coherent QA in Figs. 9 (b) and 6

Figure 11 (FIG. 11 :
Figure 11 (b)  shows that most samples satisfy the single domain wall condition t MCS ≥ 10 5 for h 0 = 0.2 and

FIG. 12 :
FIG.12: Dependence of the observables on the execution time for different annealing protocols with k = 0.5, w 0 = 0.2, and h 0 = 1.0.Both axes are the same as those in Fig.4.

Figure 13 (
Figure 13 (d) demonstrates that the D-Wave 2000Q gives the largest value of P GS .Although the D-Wave 2000Q and the D-Wave Advantage2 Prototype 1.1 give similar exponents in E res in the long annealing time region from TableIII, the performance gap appears in P GS .This result implies that the D-Wave Advantage2 Prototype 1.1 may reach low energy states, but not necessarily the ground states.

2 .
Figure14shows the dependence of the observables on the annealing time for different values of λ.All data are obtained from the D-Wave 2000Q.As in the main text, we refer to the D-Wave 2000Q as the D-Wave.The data for λ = 1.0 is the same as used in Fig.4.In Fig.14(a), the polynomial decay of ρ can be seen in the short annealing time region 1(µs) ≤ t a ≤ 10(µs) for different values of λ.The fitting exponents are summarized in TableIV.In this region, the coefficient λ hardly affects the behavior of the system.On the other hand, in the long annealing time region t a ≥ 100(µs), the saturation seen in ρ for λ = 1.0 does not occur in λ = 0.1 and 0.5.The small λ weakens the effect of the problem term Eq. (10) and enhances the effect of the constraint term Eq.(8).Figure14(b) demonstrates the value of λ hardly affects P const in the short annealing time region.Similarly to the kink density in Fig.14 (a), the effect of λ can be seen in the long annealing time region t a ≥ 100(µs).Premature finite saturation of P const appears for λ = 1.0.The number of samples that satisfy the single domain wall

FIG. 14 :
FIG. 14: Dependence of the observables attained by the D-Wave on the annealing time for different values of λ.Both axes are the same as those in Fig.4.
Appendix C: Comparison of the data of the D-Wave2000Q and thermal equilibrium state

FIG. 15 :
FIG. 15: Effective temperature obtained from the data of the D-Wave as a function of annealing time.The dashed horizontal line represents the physical temperature of the D-Wave.
FIG.16: Relation between the inverse effective temperature and freeze-out point.The lines represent the annealing schedule functions used to determine s * from T eff .

FIG. 17 :FIG. 18 :
FIG. 17: Comparison between data of the D-Wave (shown by symbols with lines) and the thermal equilibrium of H 0 at the corresponding effective temperature (shown by filled triangles).(a) Constraint satisfaction probability and (b) ground state probability as functions of annealing time.

TABLE II :
The fitting exponents of the residual energy in the long annealing time region.
).From TableIII, the largest exponent of E res can be obtained by the D-Wave 2000Q.In our experiments, the D-Wave 2000Q gives the lowest value of FIG.13: Dependence of the observables on the annealing time for the different quantum annealers.Both axes are the same as those in Fig.4.

TABLE III :
The fitting exponents of the observables for the different quantum annealers.The "E res (short)" represents the residual energy in the short annealing time region 1 ≤ t a ≤ 10.The "E res (long)" is the residual energy in the long annealing time region 100 ≤ t a ≤ 500.resamongall devices.In the long annealing time region t a ≥ 100(µs), we can see a polynomial decay of E res .The fitting exponents obtained by the D-Wave 2000Q and the D-Wave Advantage2 prototype 1.1 take similar values to those in TableIII.The result indicates that qualitative behaviors of the two devices show little differences.The D-Wave Advantage 4.1 has the largest exponent among all devices, from which one may expect Advantage 4.1 to yield smaller values of E res than other models if the device would operate with minimal noise in the longer annealing time region. E

TABLE IV :
The fitting exponents of the observables for different values of λ.The fitting region is the same as used in TableIII.