^†^†thanks: Present address: Arithmer Inc., R&D Headquarters, Terashimahonchonishi, Tokushima-shi, Tokushima 770-0831, Japan

Breakdown of the weak-coupling limit in quantum annealing

Yuki Bando Institute of Innovative Research, Tokyo Institute of Technology, Yokohama, Kanagawa 226-8503, Japan Ka-Wa Yip Center for Quantum Information Science & Technology, University of Southern California, Los Angeles, California 90089, USA Department of Physics & Astronomy, University of Southern California, Los Angeles, California 90089, USA Huo Chen Center for Quantum Information Science & Technology, University of Southern California, Los Angeles, California 90089, USA Department of Electrical & Computer Engineering, University of Southern California, Los Angeles, California 90089, USA Daniel A. Lidar Center for Quantum Information Science & Technology, University of Southern California, Los Angeles, California 90089, USA Department of Physics & Astronomy, University of Southern California, Los Angeles, California 90089, USA Department of Electrical & Computer Engineering, University of Southern California, Los Angeles, California 90089, USA Department of Chemistry, University of Southern California, Los Angeles, California 90089, USA Hidetoshi Nishimori International Research Frontiers Initiative, Tokyo Institute of Technology, Yokohama, Kanagawa 226-8503, Japan Graduate School of Information Sciences, Tohoku University, Sendai, Miyagi 980-8579, Japan RIKEN Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS), Wako, Saitama 351-0198, Japan

Abstract

Reverse annealing is a variant of quantum annealing, in which the system is prepared in a classical state, reverse-annealed to an inversion point, and then forward-annealed. We report on reverse annealing experiments using the D-Wave 2000Q device, with a focus on the $p=2$ $p$ -spin problem, which undergoes a second order quantum phase transition with a gap that closes polynomially in the number of spins. We concentrate on the total and partial success probabilities, the latter being the probabilities of finding each of two degenerate ground states of all spins up or all spins down, the former being their sum. The empirical partial success probabilities exhibit a strong asymmetry between the two degenerate ground states, depending on the initial state of the reverse anneal. To explain these results, we perform open-system simulations using master equations in the limits of weak and strong coupling to the bath. The former, known as the adiabatic master equation (AME), with decoherence in the instantaneous energy eigenbasis, predicts perfect symmetry between the two degenerate ground states, thus failing to agree with the experiment. In contrast, the latter, known as the polaron transformed Redfield equation (PTRE), is in close agreement with experiment. Thus our results present a challenge to the sufficiency of the weak system-bath coupling limit in describing the dynamics of current experimental quantum annealers, at least for reverse annealing on timescales of a $\mu$ sec or longer.

I Introduction

Quantum annealing is a quantum metaheuristic originally conceived as a method to obtain the global solution of an optimization problem by using quantum fluctuations to escape from local minima [1, 2, 3] (see Refs. [4, 5, 6, 7, 8] for reviews).

Commercial devices developed by D-Wave Systems that realize programmable quantum annealing in hardware have now been available for more than a decade [9]. Not only can these devices be used to solve optimization problems by physically realizing quantum annealing, but they are also able to perform physics experiments, for example, spin-glass phase transitions [10], the Kosterlitz-Thouless phase transition [11, 12], alternating-sector ferromagnetic chains [13], the Kibble-Zurek mechanism [14, 15, 16], the Griffiths-McCoy singularity [17], spin ice [18], the Shastry-Sutherland model [19], field theory [20], and spin liquids [21].

In this work we use a D-Wave 2000Q quantum annealer as a simulator of a simple spin system: the $p$ -spin model [22, 23, 24] with $p=2$ . But rather than using traditional “forward” annealing, wherein the Hamiltonian of a system initialized in the ground state of a transverse field evolves to a longitudinal target Hamiltonian [1], here we focus on reverse annealing, which is defined as the process of choosing an appropriate classical state as the initial state, starting from the target Hamiltonian and mixing it with the transverse field Hamiltonian to return to the original target Hamiltonian [25]. Reverse annealing is conceptually richer than forward annealing, since it introduces at least two additional parameters: the inversion point of the anneal and the pause duration. These parameters can be chosen to make forward annealing a special case of (the forward direction of) reverse annealing. Moreover, one can also iterate the process, leading to a strategy known as iterated reverse annealing [26, 12]. We note that there is some evidence that reverse annealing can outperform forward annealing in solving optimization problems [27, 28, 29, 26, 30, 31, 32, 33]. However, our focus in this work is not on algorithmic performance, but rather on using the rich playground provided by the $p=2$ $p$ -spin model under reverse annealing to answer the question of which of a variety of models best describes the results obtained from the D-Wave annealer. This question has been addressed before many times [34, 35, 36, 37, 38, 39, 40, 41, 13, 11, 12, 16], but as we shall show, the $p=2$ $p$ -spin model under reverse annealing allows us to rather clearly reveal a deficiency of one of the most popular and successful models, the weak-coupling adiabatic master equation (AME) [42]. A relatively recent model, with strong system-bath coupling, the polaron-transformed Redfield master equation (PTRE) [43, 44], provides a closer match to the empirical data we present.

We remark that the $p$ -spin model has been studied before in the context of reverse annealing, for $p=3$ [45]. The choice of $p=3$ was motivated by the fact that this model has a first order phase transition in the thermodynamic limit [22, 23, 24] (gap closing exponentially in the system size) and thus is a hard problem for conventional (forward) quantum annealing. However, this model does not have a direct physical embedding in current quantum annealing hardware, since it involves three-body interactions. While gadgets can be used to embed such a model in physical hardware supporting only two-body interactions, this comes at the cost of using up qubits and also introduces various errors [46]. Hence the study [45] was purely numerical. The $p$ -spin model with $p=2$ , which undergoes a second-order phase transition (polynomially closing gap), can be directly represented in the hardware graph of the D-Wave devices, and we do this here. The direct embedding allows us to avoid gadgetization errors and also gives us access to larger system sizes than the $p=3$ case: we experiment with up to $20$ fully connected (logical) qubits. Another reason for our choice of $p=2$ is the existence of ground-state degeneracy, which leads to additional information and the crucial ability to distinguish the AME from the other models, as will be described below.

The structure of this paper is as follows. We first define the $p$ -spin model with $p=2$ in Sec. II and describe the reverse annealing protocol we employed in this study. In Sec. III we present and discuss the empirical results from the D-Wave 2000Q device, for a variety of different parameter settings, including initial conditions, annealing time, pause duration, problem size, and the effect of iteration. This is followed in Sec. IV by an examination of numerical results based first on the closed-system Schrödinger equation, followed by quantum adiabatic master equation with both independent and collective system-bath coupling, and the PTRE. In all cases we compare and contrast the predictions of the various models to the empirical D-Wave data. Conclusions are presented in Sec. V, and appendixes provide further technical details, as well as additional results.

II Problem definition and reverse annealing protocol

In this section we describe the problem Hamiltonian and the reverse annealing protocol we used both in our experiments on the D-Wave device and in numerical simulations.

II.1 The $p$ -spin model with $p=2$

We consider a quantum annealing Hamiltonian comprising a driver Hamiltonian $H_{D}$ and a target Hamiltonian $H_{T}$ , the latter encoding the combinatorial optimization problem represented as an Ising model, as a function of dimensionless time $0\leq s(t)\leq 1$ :

H(s)=\frac{A(s)}{2}H_{D}+\frac{B(s)}{2}H_{T}\ .

(1)

Here, $A(s)$ and $B(s)$ are device-dependent annealing schedules. The D-Wave device used in the present experiment has $A(s)$ and $B(s)$ as shown in Fig. 1(a). We work in units where $\hbar=1$ and $k_{\mathrm{B}}=1$ except in Figs. 1(a), 6, 15, and 18, where we opted to use units where $h=1$ in accordance with the conventions of the D-Wave device documentation [47].

Refer to caption — Figure 1: (a) Annealing schedule of D-Wave 2000Q with the “DW_2000Q_6” solver. (b) Forward (FA) and reverse annealing (RA) protocols as defined in Eqs. (4) and (5), respectively. The latter incorporates a pause of duration $t_{p}$ . Here $t_{a}=1.2~{}\mathrm{\mu s}$ for forward annealing, and $\tau=1~{}\mathrm{\mu s}$ and $s_{\mathrm{{inv}}}=0.2$ for reverse annealing.

The final values of the schedule functions are $A(s=1)\approx 0$ and $B(s=1)>0$ so that ideally the ground state of $H_{T}$ is realized as the final state. The time dependence of $s(t)$ , in particular forward or reverse annealing as depicted in Fig. 1(b), corresponds to different variants of the general quantum annealing algorithm. The driver Hamiltonian $H_{D}$ is usually chosen as

H_{D}=-\sum_{i=1}^{N}\sigma_{i}^{x}

(2)

where $N$ is the number of qubits and $\sigma_{i}^{x}$ is the $x$ component of the Pauli matrix acting on the $i$ th qubit.

We study the $p$ -spin model,

H_{T}=-N\left(\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}^{z}\right)^{p}

(3)

with $p=2$ . This Hamiltonian couples every qubit to every other qubit, which exceeds the “Chimera” graph connectivity of the D-Wave 2000Q device. Thus the qubits in $H_{T}$ should be viewed as logical qubits, to be represented in the actual device by physical qubits. These physical qubits are coupled in ferromagnetic chains to form logical qubits, a procedure known as minor embedding [48] ¹¹1 Embedding a fully-connected graph was implemented in the standard way; see, e.g., Fig. 1 of Ref. [49].. We chose the ferromagnetic interaction between physical qubits in a logical qubit to be $J_{\text{F}}=-1.0$ , as discussed in App. A.

II.2 Reverse annealing protocol

The traditional protocol of forward annealing has

s(t)=\frac{t}{t_{a}}\ \ \ t\in[0,t_{a}],

(4)

where $t_{a}$ is the total annealing time; see the solid red line in Fig. 1(b). The initial state of forward annealing is the ground state of the driver Hamiltonian $H_{D}$ .

In reverse annealing, we start from $s=1$ and decrease $s$ to an intermediate value $s_{\mathrm{inv}}$ , pause, then increase $s$ to finish the process at $s=1$ . The explicit time dependence $s(t)$ realized on the D-Wave device is

s(t)=\begin{cases}1-\,\displaystyle\frac{\,t}{\tau}&0\leq t\leq t_{\rm inv}\\ 1-\displaystyle\frac{t_{\rm inv}}{\tau}=s_{\mathrm{inv}}&t_{\rm inv}\leq t\leq t_{\rm inv}+t_{p}\\ 2s_{\mathrm{inv}}-1-\displaystyle\frac{t_{p}}{\tau}+\displaystyle\frac{t}{\tau}&t_{\rm inv}+t_{p}<t\leq 2t_{\rm inv}+t_{p},\end{cases}

(5)

where $1/\tau$ is the annealing rate, $t_{\rm inv}$ is the inversion time defined as $t_{\rm inv}=\tau(1-s_{\mathrm{inv}})$ , and $t_{p}$ is the duration of the intermediate pause. Our reverse annealing begins from $s(t=0)=1$ and ends at $s(t=t_{a})=1$ , where

t_{a}=2t_{\rm{inv}}+t_{p}=2\tau(1-s_{\mathrm{inv}})+t_{p},

(6)

passing through the inversion point $s_{\mathrm{inv}}$ ; see the blue dashed line in Fig. 1(b). To distinguish between $\tau$ and $t_{a}$ , we henceforth refer to the former as the “annealing time” and the latter as the “total annealing time”. The initial state of reverse annealing is a classical state, usually a candidate solution to the combinatorial optimization problem, i.e, a state that is supposed to be close to the solution.

We can iteratively repeat the process of reverse annealing with the final state of a cycle as the initial condition of the next cycle. This protocol is called iterated reverse annealing [26] and we denote the iteration number by $r$ . Additional details are provided in App. B.

III Empirical results

In this section we report the results of performing reverse annealing experiments for the $p=2$ $p$ -spin model on the D-Wave 2000Q device.

III.1 Dependence on initial conditions

Figure 2 shows the empirical success probability, i.e., the probability that the final state observed is one of the two degenerate ground states. This is shown for $N=20$ and no pause. The initial condition is specified by the initial value of the normalized total magnetization,

\displaystyle m_{0}=\frac{1}{N}\sum_{i=1}^{N}\langle\psi_{0}|\sigma_{i}^{z}|\psi_{0}\rangle,

(7)

where $|\psi_{0}\rangle$ denotes the initial wave function (a classical state). Since the doubly degenerate ground state of the problem Hamiltonian $H_{T}$ has $\pm 1$ as the normalized magnetization, a value of $m_{0}$ closer to 1 (or $-1$ ) represents an initial condition exhibiting higher overlap with the ground state. In Fig. 2 we present results for a completely unbiased initial condition ( $m_{0}=0$ ) and initial conditions strongly biased toward the all-up state $\sigma_{i}^{z}=1~{}\forall i$ ( $m_{0}=0.8,0.9$ ). The state with $m_{0}=0$ is the highest excited state, and the state with $m_{0}=0.9$ is the first excited state.

III.1.1 $m_{0}=0$

Figure 2(a) shows the probability that the system reaches the all-up state (denoted “up”) and the all-down state ( $\langle\sigma_{i}^{z}\rangle=-1,~{}\forall i$ , denoted “down”). We call these the partial success probabilities.

When the initial condition is $m_{0}=0$ , the up and down probabilities are equal to within the error bars for all $s_{\rm{inv}}$ , as expected because the initial state is unbiased. Note that both success probabilities are close to $0.5$ for $s_{\mathrm{inv}}<0.4$ whereas they are zero for $s_{\mathrm{inv}}>0.5$ . The vertical dashed line at $s_{\mathrm{inv}}\approx 0.36$ indicates the minimum-gap point, i.e., the point where the energy gap between the ground and the second excited states becomes minimum, which we denote by $s_{\Delta}$ (the first excited state becomes degenerate with the ground state at the end of the anneal, hence it is the second excited state that is relevant). It is likely that the system remains close to its initial state for $s_{\mathrm{inv}}>0.5$ , where the transverse field and the associated quantum fluctuations are small, and the true final ground state with $m_{0}=\pm 1$ is hard to reach, since the initial state with $m_{0}=0$ has small overlap with the latter. The situation is different when $s_{\mathrm{inv}}<0.4$ , i.e., when the system traverses the minimum-gap region. In this case the system enters the paramagnetic (quantum-disordered) phase dominated by the driver $H_{D}$ , so that the initial condition is effectively erased. This renders the process similar to conventional forward quantum annealing, by which the two true ground states are reached with nearly equal probability $0.5$ . Recall that the phase transition around $s_{\Delta}$ is of second order in the present problem, and therefore the system finds (one of) the true final ground states relatively easily by forward annealing because of the mild, polynomial, closing of the energy gap, as we discuss in more detail later [see Fig. 6].

III.1.2 $m_{0}=0.8,0.9$ : up/down symmetry breaking for $0.3\lesssim s_{\mathrm{inv}}\lesssim 0.5$

For the initial states with $m_{0}=0.8$ and $m_{0}=0.9$ , the experimental results shown in Fig. 2(a) reveal significant differences between the probabilities of the final all-up and all-down states. Of course, the initial conditions $m_{0}=0.8,0.9$ have much larger overlap with the all-up state, and this fact alone suggests a mechanism for breaking the symmetry between the probabilities of the final state being the all-up or all-down state. However, it is remarkable that the success probability for the all-up state is almost 1 in a region $0.3\lesssim s_{\mathrm{inv}}\lesssim 0.5$ both to the left and to the right of $s_{\Delta}\approx 0.36$ (we discuss the additional asymmetry for $s_{\mathrm{inv}}\lesssim 0.3$ below). As we explain in Sec. IV.3, this behavior is inconsistent with a model of an open quantum system that is weakly coupled to its environment and is thus described by the adiabatic master equation [42]. However, it is consistent with both a quantum model of a system that is strongly coupled to its environment, as explained in Sec. IV.4, and with a simple semiclassical model captured by the spin-vector Monte Carlo algorithm [39, 50], as explained in App. F.

III.1.3 $s_{\mathrm{inv}}\gtrsim 0.5>s_{\Delta}$ : freezing

Note that for all three initial conditions the success probability eventually vanishes for $s_{\mathrm{inv}}\gtrsim 0.5>s_{\Delta}$ . The reason is that in all cases the system is initialized in an excited state, and remains in an excited state at the end of the anneal, since there is no mechanism for thermal relaxation to a lower energy state when $s_{\mathrm{inv}}$ is large. This is a manifestation of the phenomenon of freezing [42, 51, 37], i.e., the extreme slowdown of relaxation due to the fact that the system-bath interaction (nearly) commutes with the system Hamiltonian when the transverse field is very small. In addition, the annealing timescale is manifestly too slow for downward diabatic transitions. This is true even with the discontinuity in the derivative of the schedule depicted in Fig. 1(b). Despite this, the reversal of the anneal direction is apparently too slow in practice to have a non-adiabatic effect, or if diabatic transitions do occur, then they exclusively populate higher excited states.

III.1.4 $s_{\mathrm{inv}}\lesssim 0.3<s_{\Delta}$ : spin-bath polarization

It is also noteworthy from Fig. 2(a) that the initial condition results in different probabilities for the all-up and all-down states even in the paramagnetic region $s_{\mathrm{inv}}\lesssim 0.3<s_{\Delta}$ , where quantum fluctuations are large and the all-up and all-down states are expected to have the same probability in equilibrium. The system “remembers” the initial condition to a certain extent, an anomaly that may be attributable to spin-bath polarization [52, 53]. Namely, the persistent current flowing in the qubit body during the anneal produces a magnetic field that can partially align or polarize an ensemble of environmental spins local to the qubit wiring, with a much slower relaxation time than the anneal duration. Given the polarized initial condition $m_{0}=0.8$ or $0.9$ , this polarized spin-bath will be aligned with the all-up state even after the system crosses into the paramagnetic phase, thus preventing the system from equilibrating and explaining the observed memory effect. Spin-bath polarization is expected to be particularly pronounced under reverse annealing, since the strong polarization of the initial state will act to polarize the spin bath, more so than in forward annealing, where the initial state is unpolarized.

III.1.5 Total success probability

Figure 2(b) shows the total success probability, i.e., the sum of the final up and down probabilities. The resulting curve stays almost flat and close to $1$ for $s_{\mathrm{inv}}<0.4$ . This constant $1$ is a mixture of the two effects manifest in Fig. 2(a), the genuine effects of reverse annealing around $s_{\mathrm{inv}}\approx 0.4$ for $m_{0}=0.9$ and $0.8$ , and the effectively-forward-annealing-like behavior for $m_{0}=0$ . Indeed, the case with $s_{\mathrm{inv}}=0$ is essentially equivalent to standard, forward annealing: the inversion point is set at $s=0$ , so that the anneal restarts from a Hamiltonian that is completely dominated by the transverse field. That the empirical total success probability is $1$ in this case shows that the $p=2$ problem is easy also for forward annealing, in contrast to the $p=3$ case studied in Ref. [45], which numerically found very small values of success probability near $s_{\mathrm{inv}}=0$ . This may be explained by the very fast forward annealing in this parameter region in Ref. [45], which keeps the system almost unchanged from the quantum-disordered state at $s=s_{\mathrm{inv}}$ . Aside from this subtlety, our experimental data are consistent with the numerical results of Ref. [45], supporting the latter’s expectation that the system-environment interaction (through spin-boson coupling) prompts relaxation of the system toward the ground state around the minimum-gap region.

We shall see in Sec. IV.3.2 that all the salient features seen in Fig. 2(b), such as the shift to the left with decreasing $m_{0}$ , are captured by our open system simulations.

III.2 Dependence on annealing time

As seen from Eq. (6), the total annealing time $t_{a}$ is linearly dependent on the annealing time $\tau$ , and proportional to $\tau$ when $t_{p}=0$ . Figure 3 shows the success probability for different $\tau$ with $N=20$ , $m_{0}=0.9$ , and no pausing. The overall trend is similar to Fig. 2.

As shown in Fig. 3(a), an increase of $\tau$ leads to slightly higher all-up success probabilities for $s_{\mathrm{inv}}>s_{\Delta}$ but the other way around for $s_{\mathrm{inv}}<s_{\Delta}$ . This can be explained in terms of an increased relaxation to the ferromagnetic ground state near the minimum energy gap for larger $\tau$ and an enhanced relaxation to the paramagnetic ground state for $s_{\mathrm{inv}}<s_{\Delta}$ . As seen in Fig. 3(b), the total success probability stays close to $1$ for $s_{\mathrm{inv}}<s_{\Delta}$ and benefits slightly from increasing $\tau$ for $s_{\mathrm{inv}}>s_{\Delta}$ .

III.3 Effects of pausing

Figure 4 shows success probabilities for pausing time $t_{p}\in\{0,2,10\}~{}\mu s$ as a function of $s_{\mathrm{inv}}$ with $N=20$ , $m_{0}=0.9$ and $\tau=5~{}\mu s$ .

The overall trend is similar to Fig. 3, i.e., a longer pause leads to an increased relaxation, but the effects are more significant in the present case. Pausing at $0.5<s_{\mathrm{inv}}<0.6$ greatly improves the success probability from nearly 0 ( $t_{p}=0~{}\mu s$ ) to nearly 1 ( $t_{p}=2~{}\mu s$ and $10~{}\mu s$ ), implying that pausing in a relatively narrow region slightly past the minimum gap point is most effective. This is in line with previous findings [31, 54, 50].

Another significant difference between pausing and not pausing is that the partial success probabilities shown in Fig. 4(a) are nearly piecewise flat for $t_{p}>0$ , showing that the effect responsible for the anomaly disappears under pausing. This is consistent with the spin-bath polarization effect discussed in Sec. III.1, in that the pause provides the time needed for this polarization to relax to equilibrium, on a timescale of a few $\mu s$ .

III.4 Size dependence

Figures 5(a) and 5(b), respectively, show the partial and total success probability for different system sizes $N$ as a function of $s_{\mathrm{inv}}$ with $\tau=5~{}\mu s$ . The initial state for each $N$ is the first excited state with $m_{0}<1.0$ (closest to the ground state for which $m_{0}=1.0$ ), i.e., a state with one spin flipped. We observe a statistically significant slight non-monotonicity with $N$ in the region $s_{\mathrm{inv}}\geq s_{\Delta}$ of the up success probabilities [Fig. 5(a)] and the total success probabilities [Fig. 5(b)], namely, the ordering of the shift to the left with $N$ is: $\{4,8,16,12,20,24\}$ ; we do not have an explanation for this effect, though it could plausibly be an anomaly related to the minor embedding procedure. However, the down success probabilities are monotonic in $N$ , and the overall trend is clear, namely, the drop-off to zero success probability occurs at smaller $s_{\mathrm{inv}}$ as $N$ is increased.²²2 We have confirmed that this tendency persists for larger systems from preliminary data for up to $N=48$ . We did not collect data systematically for larger systems because of technical difficulties. This is consistent with the reduction in the position of the minimum gap, $s_{\Delta}$ , as a function of $N$ , shown in Fig. 6, which is tracked by the value of $s_{\mathrm{inv}}$ at which the total success probability equals $0.5$ . This suggests that the success probabilities (both partial and total) are sensitive to the location of the minimum quantum gap. See App. C for more details about the spectrum of the $p$ -spin problem.

Additionally, the asymmetry between the all-up and all-down states for $s_{\mathrm{inv}}<s_{\Delta}$ is enhanced with increasing $N$ . This is consistent with the formation of larger domains of spin-bath polarization, which would be expected to take longer to dissipate due to their size as $N$ increases.

Finally, the overlap of data for $N=20$ and $24$ suggests that large- $N$ effects have already converged at these sizes, i.e., that $N\sim 20$ is sufficiently large to infer the behavior at larger $N$ . Also, even the smallest systems with $N=4$ and $8$ already share qualitative features with larger systems.

III.5 Effects of iteration

We next study the effects of iteration on reverse annealing, i.e., how the success probability depends on the number of iterations $r$ . For this purpose, we used the reinitialize_state=False setting of the D-Wave device, meaning that the output state of the previous iteration was the initial state of the next iteration. Figure 7 shows the results for $r\in\{1,10,25,50\}$ as a function of $s_{\mathrm{inv}}$ , with $N=20$ , $m_{0}=0.9$ , $\tau=5~{}\mu s$ , and no pausing. The success probabilities improve in the region $s_{\mathrm{inv}}\geq s_{\rm\Delta}$ as the number of iterations $r$ increases, regardless of the initial state $m_{0}$ . We expect the relaxation to the low-energy state due to coupling to the environment to be induced by successive iterations, and the occupation of the ground states to increase correspondingly. The results in Fig. 7 confirm this expectation, in that the success probability is larger at given $s_{\mathrm{inv}}\geq s_{\Delta}$ as $r$ is increased.

Comparing Figs. 7(a) and (b) with $m_{0}=0.9$ and Figs. 7(c) and (d) with $m_{0}=0$ , we can see that the initial state with $m_{0}=0.9$ has a larger improvement in success probability with fewer iterations $r$ in the region where $s_{\mathrm{inv}}\geq s_{\Delta}$ . The initial state with $m_{0}=0$ deviates greatly from the ground states, and there are many excited states between it and the latter. Therefore, when the system transitions from the initial state of $m_{0}=0$ to other states by relaxation due to coupling to the environment, those excited states become populated. Many iterations are expected to be necessary to obtain a significant improvement in the success probability. On the other hand, there are only a few excited states between the initial state $m_{0}\lesssim 1.0$ and the all-up ground state. Therefore, it is expected that only a few iterations will be sufficient to make the transition to this ground state. Our results support this picture of improved performance under iterated reverse annealing.

IV Numerical Results

We next present closed and open-systems simulations and compare the results with the data presented in the previous section. We choose relatively small system sizes $N=4$ and $8$ to facilitate numerical computations. Because our experimental data do not show a strong size dependence, as discussed in Sec. III.4, the qualitative comparison our numerical results enable should suffice to draw relevant physical conclusions.

In all our numerical studies we simulate the logical problem directly, i.e., we do not simulate the embedded problem which replaces every logical spin by a ferromagnetic chain, as in our D-Wave experiments.

For simplicity, we focus on reverse annealing without pausing (i.e., $t_{p}=0$ ) in this section. We expect that in general pausing will lead to an overall success probability increase in open-system simulations.

IV.1 Closed system model

An analysis of the closed system case, while being unrealistic due to the absence of thermal effects, is instrumental in isolating the effect and importance of diabatic transitions in explaining our experimental results.

The time evolution of reverse annealing in a closed system is described as follows. The state after a single iteration (cycle) is

\ket{\psi(2t_{\text{inv}})}=U(2t_{\text{inv}},0)\ket{\psi(0)}\,,

(8)

where $\ket{\psi(0)}$ is the initial state and

U(2t_{\text{inv}},0)={\cal T}\exp\left[-i\int_{0}^{2t_{\text{inv}}}H(t^{\prime})dt^{\prime}\right]

(9)

is the unitary time-evolution operator and ${\cal T}$ denotes forward time ordering. At the end of $r$ cycles, the final state $\ket{\psi(2rt_{\text{inv}})}$ can be expressed as:

\ket{\psi(2rt_{\text{inv}})}=U(2rt_{\text{inv}},0)\ket{\psi(0)}\,,

(10)

with

U(2rt_{\text{inv}},0)=\prod_{i=0}^{r-1}U(2(i+1)t_{\text{inv}},2it_{\text{inv}})\,.

(11)

Here, the final state of the $r$ th cycle is the initial state of the next cycle. This condition is shared by the experiments in the previous section. The solution states are doubly degenerate, and the total success probability at the end of $r$ cycles is

p(r)=|\braket{\psi(2rt_{\text{inv}})}{\text{up}}|^{2}+|\braket{\psi(2rt_{\text{inv}})}{\text{down}}|^{2}\,,

(12)

where $\ket{\text{up}}=\ket{\uparrow}^{\otimes N}$ and $\ket{\text{down}}=\ket{\downarrow}^{\otimes N}$ . For higher computational efficiency without loss of accuracy, we rotate the state vector into the instantaneous energy eigenbasis representation at each time step in our numerical simulations.

IV.1.1 Dependence on system size and annealing time

We initialize the state to have a single spin down, i.e., as the computational basis state $\ket{0001}$ [ $\ket{0}\equiv\ket{\uparrow}$ , $\ket{1}\equiv\ket{\downarrow},(m_{0}=0.5)$ ] for $N=4$ and $\ket{0000001}~{}(m_{0}=0.75)$ for $N=8$ , respectively. Note that in our simulations these are not exact eigenstates of $H(1)$ , due to a very small residual transverse field at $s=1$ , as in the D-Wave annealing schedule shown in Fig. 1.³³3 $A(1)=1.9\times 10^{-6}$ GHz, $B(1)=11.97718$ GHz, $A(1)/B(1)=1.58635004\times 10^{-8}$ . When the initial state $\ket{\psi(0)}$ is a computational basis state, there is a simple upper bound on the success probability achievable in closed system reverse annealing: the population of the initial state in the maximum-spin sector (see App. D). This upper bound explains why the following closed system results have relatively low success probabilities.

We plot in Fig. 8 the simulation results for the total success probability for various annealing rates, subject to the D-Wave annealing schedule shown in Fig. 1. We see that the success probability after a single cycle is non-negligible only when $\tau$ is small enough ( $\tau<1$ ns), in which case diabatic transitions to states with lower energies take place when $s_{\mathrm{inv}}<s_{\Delta}$ . However, the success probability remains small even in this case, and in any case much smaller than in our experimental results where thermal relaxation plays the dominant role.

IV.1.2 Dependence on initial state and number of iterations

We next choose the initial state $\ket{\psi(0)}$ with two spins down $\ket{0011}~{}(m_{0}=0)$ for $N=4$ and $\ket{00000011}~{}(m_{0}=0.5)$ for $N=8$ , respectively, i.e., the second excited states of $H_{T}$ for these respective system sizes.

Figure 9 (top row) shows that as expected, for $N=4$ and $r=1$ , the success probability is overall lower than the case when the initial state is the first excited state, Fig. 8. Only for the highest annealing rate (smallest $\tau$ , most diabatic) is the maximum success probability similar to the case of the first excited state. For $N=8$ , the success probability is much smaller, no matter how diabatic the process is. This indirectly confirms the dominant role played by thermal effects in our experimental results.

We plot in Fig. 9 (bottom row) the results for $r=2$ cycles, where we see that the success probability decreases compared to the $r=1$ case. This is consistent with the conclusions of Ref. [26].

Recalling the experimental results reported in the previous section, we may therefore safely conclude that, as expected, the closed-system picture is far from the experimental reality of the D-Wave device.

IV.2 Open system model: setup

For an open system, the state after the first cycle is

\rho(2t_{\text{inv}})=V(2t_{\text{inv}},0)\rho(0)\,,

(13)

where

V(t,0)={\cal T}\exp\left[\int_{0}^{t}\mathcal{L}(t^{\prime})dt^{\prime}\right]\,.

(14)

$\mathcal{L}(t)$ is the time-dependent Liouville superoperator, which generates a master equation of the form


$\displaystyle\frac{d\rho(t)}{dt}$	$\displaystyle=\mathcal{L}(t)\rho(t)$	(15a)
	$\displaystyle=i\bigl{[}\rho(t),H(t)+H_{\text{LS}}(t)\bigr{]}+\mathcal{D}\bigl{[}\rho(t)\bigr{]}\,.$	(15b)

Here, $H_{\text{LS}}(t)$ is the Lamb shift term and $\mathcal{D}$ is the dissipator. Again, the time dependence of the integrand in Eq. (14) is incorporated into $s(t)$ .

At the end of $r$ cyles, the final state $\rho(2rt_{\text{inv}})$ is expressed as:

\rho(2rt_{\text{inv}})=V(2rt_{\text{inv}},0)\rho(0)\,,

(16)

with

V(2rt_{\text{inv}},0)=\prod_{i=0}^{r-1}V(2(i+1)t_{\text{inv}},2it_{\text{inv}})\,.

(17)

We next consider both the weak and the strong coupling cases, which give rise to different Liouville superoperators.

IV.3 Weak coupling: adiabatic master equation

In this subsection we use the adiabatic master equation (AME), which holds under weak coupling to the environment [42]. In this limit $\mathcal{D}$ can be expressed in a diagonal form with Lindblad operators $L_{i,\omega}(t)$ :

	$\displaystyle{\mathcal{D}}[\rho(t)]$	$\displaystyle=$	$\displaystyle\sum_{i}\sum_{\omega}\gamma_{i}(\omega)\left(L_{i,\omega}(t)\rho(t)L^{\dagger}_{i,\omega}(t)\phantom{\frac{1}{2}}\right.$		(18)
			$\displaystyle\left.\qquad\qquad-\frac{1}{2}\left\{L_{i,\omega}^{\dagger}(t)L_{i,\omega}(t),\rho(t)\right\}\right)\,,$		(18)

where the summation runs over the qubit index $i$ and the Bohr frequencies $\omega$ [all possible differences of the time-dependent eigenvalues of $H(t)$ ]. This dissipator expresses decoherence in the energy eigenbasis [55]: quantum jumps occur only between the eigenstates of $H(t)$ [56].

We consider two different models of system-bath coupling: independent and collective dephasing [57]. In the first case, we assume that the qubit system is coupled to independent, identical bosonic baths, with the bath and interaction Hamiltonians being


$\displaystyle H_{B}$	$\displaystyle=\sum_{i=1}^{N}\sum_{k=1}^{\infty}\omega_{k}b_{k,i}^{\dagger}b_{k,i}\ ,$	(19a)
$\displaystyle H_{SB}^{\text{ind}}$	$\displaystyle=g\sum_{i=1}^{N}\sigma_{i}^{z}\otimes\sum_{k}\left(b_{k,i}^{\dagger}+b_{k,i}\right)\ ,$	(19b)

where $b_{k,i}^{\dagger}$ and $b_{k,i}$ are, respectively, the raising and lowering operators for the $k$ th oscillator mode with natural frequency $\omega_{k}$ . The rates appearing in Eq. (18) are given by

\gamma_{i}(\omega)=2\pi\eta g^{2}\frac{\omega e^{-|\omega|/\omega_{c}}}{1-e^{-\beta\omega}}\ ,

(20)

arising from an Ohmic spectral density, and satisfying the Kubo-Martin Schwinger (KMS) condition [58, 59] $\gamma_{i}(-\omega)=e^{-\beta\omega}\gamma_{i}(\omega)$ , with $\beta=1/T$ the inverse temperature. We use $\eta g^{2}=10^{-3}$ , the cutoff frequency $\omega_{c}=1$ THz, and the D-Wave device operating temperature $T=12.1\text{mK}=1.57$ GHz. With the independent dephasing assumption, the Lindblad operators are:

L_{i,\omega}(t)=\sum_{\epsilon_{b}-\epsilon_{a}=\omega}|\varepsilon_{a}(t)\rangle\langle\varepsilon_{a}(t)|\sigma^{z}_{i}|\varepsilon_{b}(t)\rangle\langle\varepsilon_{b}(t)|\,,

(21)

corresponding to dephasing in the instantaneous eigenbasis $\{|\varepsilon_{a}(t)\rangle\}$ of $H(t)$ . Similarly to the closed system case, we rotate the density matrix and Lindblad operators into the instantaneous energy eigenbasis at each time step in our numerical simulations. This keeps the matrices sparse, without loss of accuracy. For $N=8$ we truncate the system size to the lowest $n=18$ levels out of $256$ [ $18=2(1+8)$ : the total number of degenerate ground and first excited states at $s=1$ ].

We also consider the collective dephasing model, where all the qubits are coupled to a collective bath with the same coupling strength $g$ , which preserves the spin symmetry. In this case, the interaction Hamiltonian becomes

\displaystyle H_{SB}^{\text{col}}=gS^{z}\otimes B,

(22)

where

\displaystyle S^{z}=\sum_{i}\sigma_{i}^{z},\quad B=\sum_{k}\left(b_{k}^{\dagger}+b_{k}\right).

(23)

With this assumption, we can group together the Lindblad operators corresponding to different qubits $i$ into a single one:

L_{\omega}(t)=\sum_{\epsilon_{b}-\epsilon_{a}=\omega}|\varepsilon_{a}(t)\rangle\langle\varepsilon_{a}(t)|S^{z}|\varepsilon_{b}(t)\rangle\langle\varepsilon_{b}(t)|\,.

(24)

The resulting number of Lindblad operators is a factor of $N$ smaller than that of the independent system-bath coupling model.

The total success probability at the end of the $r$ cycles is

p(r)=\bra{\text{up}}\rho(2rt_{\text{inv}})\ket{\text{up}}+\bra{\text{down}}\rho(2rt_{\text{inv}})\ket{\text{down}}\,.

(25)

Any relaxation during the reverse annealing dynamics to the global instantaneous ground state or the instantaneous first excited state of $H(t)$ is beneficial, the latter since it becomes degenerate with the ground states $\{\ket{\text{up}},\ket{\text{down}}\}$ of $H_{T}$ at the end of the anneal (see App. C).

IV.3.1 Dependence on the annealing time.

Figure 10 (top row) shows the AME simulation results for the success probability as a function of $s_{\mathrm{inv}}$ with $N=4$ , various $\tau$ , and the initial state $\ket{0001}~{}(m_{0}=0.5)$ , using the independent and collective dephasing models.

The independent dephasing model leads to a maximum success probability of $1$ for large $\tau$ and small $s_{\mathrm{inv}}$ . However, the collective dephasing model leads to a maximum success probability of $1/4$ . The reason is the same as in the closed system simulations. For the initial state $\ket{0001}$ , only $1/4$ of its population is in the subspace of maximum quantum spin number. However, the global (instantaneous) ground state and the first excited state of the annealing Hamiltonian $H(s)$ both belong to the maximum-spin subspace. Since the collective dephasing model preserves the spin symmetry and the dynamics is restricted to each subspace, at most $1/4$ of the initial population can be relaxed to the instantaneous ground state and the first excited state, and reach the correct solutions at the end of the anneal (see App. D for more details). It is also noteworthy that coherent oscillations are visible for $\tau=1$ ns [compare with Fig. 8(a)], but not for the larger values of $\tau$ we have simulated. Recall that the experimental timescale in Sec. III is on the order of a $\mu$ s.

Comparing the simulation results in Fig. 10 and the experimental results in Fig. 3(b), we observe that they share the same main features. Namely, the success probability increases with $\tau$ ; and, as $\tau$ increases the maximum success probability can be achieved with a larger inversion point $s_{\text{inv}}$ . Most notably, the total success probability also drops to zero for sufficiently large $s_{\text{inv}}$ in our simulations. This is the freeze-out effect that is well captured by the adiabatic master equation (as reported in Ref. [42]), since thermal relaxation is suppressed when the transverse field magnitude becomes so small that the system and system-bath Hamiltonians effectively commute.

Note that the experimental results in Sec. III have a maximum success probability as high as $1$ , while our closed system simulations and open-system simulations with the collective dephasing model have success probabilities upper bounded by some constants $<1$ , as already discussed. The high total success probability observed in our experiments is evidence that, as expected, the dynamics in the D-Wave device do not preserve spin symmetry. Collective system-bath coupling cannot explain the experimental results, leaving independent system-bath coupling as the only candidate consistent with the experiments according to our simulations.

IV.3.2 Dependence on initial condition and the number of cycles.

For the initial state $\ket{0011}~{}(m_{0}=0)$ and with $r=1$ , the independent dephasing model still gives a maximum success probability of $1$ as seen in Fig. 10(c). The collective dephasing model (not shown) has a maximum success probability of $1/6$ , since the initial state has $1/{{4}\choose{2}}=1/6$ of its population in the maximum-spin subspace.

Comparing Fig. 10(a) and (c), the dependence of the success probability on $s_{\text{inv}}$ is similar to that of the initial state $\ket{0001}$ . The $\tau=1$ ns coherent oscillations visible for the latter are more attenuated for $\ket{0011}$ , but this was also the case for the closed system simulations [contrast Fig. 8(a) and Fig. 9(a) at $\tau=1$ ns]. The main feature distinguishing Fig. 10(a) and (c) is the shift to the left of the $\tau\geq 100$ ns curves for the initial state $\ket{0011}$ , i.e., as $m_{0}$ is reduced from $0.5$ to $0$ . This is consistent with our experimental results, as can be seen in Fig. 2(b).

In Fig. 10(d), we consider the dependence on the number of cycles $r$ , using the independent system-bath model. For $r=2$ , we see that the results are similar to those of a single cycle. However, we do see a small improvement in the sense of a slight shift to the right of the curves with $\tau\geq 100$ ns compared with $r=1$ , which is consistent with the experimental result shown in Fig. 7(d). We note that this improvement was not observed in the closed system case, as can be seen by contrasting Figs. 9(a) and (c). Interestingly, there is also a small signature of coherent oscillations for $\tau\leq 100$ ns.

Finally, we note that compared to the closed systems case, the results depend much less on the initial condition.

IV.3.3 Size dependence and partial success probability.

In Fig. 11, we show the success probability for two different system sizes $N=4$ and $8$ . The initial state has a single spin down. We observe that the results do not depend much on the system size, with the total success probabilities slightly larger (shifted to the right) for $N=4$ . This is consistent with the experimental results for $N=4$ and $8$ shown in Fig. 5(b).

While the total success probabilities of our open-system simulations with independent dephasing are generally consistent with the experimental total success probabilities, our open-system simulations always produce symmetric partial success probabilities as shown in Fig. 11, in stark contrast with the experiments, where the all-up state is strongly favored over the all-down state in the region around the minimum gap (see all the top panels in the figures in Sec. III). Clearly, this reflects a significant failure of the AME model.

One reason for this discrepancy is that our simulations do not include a mechanism such as the polarized spin-bath which we believe explains the experimentally observed asymmetry for small $s_{\mathrm{inv}}$ (see Sec. III.1.4). Moreover, our adiabatic master equation simulations result in a symmetric final all-up and all-down population since any relaxation event during the anneal (at any $s$ ) to either the global instantaneous ground state $\ket{\epsilon_{0}(s)}$ or the instantaneous first excited state $\ket{\epsilon_{1}(s)}$ of $H(s)$ , which become degenerate at $s=1$ , eventually contributes amplitude to $\ket{\epsilon_{0}(s=1)}=({\ket{\uparrow}^{\otimes N}+\ket{\downarrow}^{\otimes N}})/{\sqrt{2}}$ or $\ket{\epsilon_{1}(s=1)}=({\ket{\uparrow}^{\otimes N}-\ket{\downarrow}^{\otimes N}})/{\sqrt{2}}$ . These two states have equal all-up and all-down populations, which is why the AME simulation results do not distinguish them. Notice the importance of the energy eigenbasis (i.e., the states $\ket{\epsilon_{0}(s)}$ and $\ket{\epsilon_{1}(s)}$ ) in this argument; this is special to the AME.

IV.4 Strong coupling: PTRE

Given the troubling discrepancy between the empirical partial success probability results and the AME simulations, in this subsection we use the Lindblad form polaron-transformed Redfield equation (PTRE) [44, 43]. The idea is to check whether a quantum model that is not subject to the weak coupling limit is capable of avoiding the discrepancy. This is motivated in part by previous studies, which showed that a hybrid AME/PTRE model works better than the pure AME in explaining the linewidth broadening phenomenon observed in a tunneling spectroscopy experiment [60, 44]. The PTRE holds under intermediate to strong coupling to the environment [44]. Unlike the AME, wherein decoherence is between energy eigenstates, decoherence in the PTRE is between computational basis states (eigenbasis of $\sigma^{z}$ ). But, unlike the singular coupling limit (SCL) [55], where decoherence is also between computational basis states, the PTRE is governed by a non-flat (and hence non-trivial) spectral density. This means that, unlike the SCL, the PTRE is sensitive to the particulars of the bath model, a fact we take advantage of below when we combine low and high frequency components into a single noise spectrum; see Eqs. (30) and (31) below. This type of hybrid spectrum approach has been used successfully [61] to explain a $16$ -qubit quantum annealing experiment with a small gap [62].

To be explicit, the PTRE is given by Eq. (15) with the polaron-frame system Hamiltonian

H(t)=\frac{1}{2}B(t)H_{T}\ ,

(26)

the polaron-frame dissipator

	$\displaystyle\mathcal{D}\left(\rho\right)$	$\displaystyle=\sum_{\alpha\in\left\{+,-\right\}}\sum_{\omega,i}\gamma_{\mathrm{p}}(\omega)\Big{(}L_{i,\omega}^{\alpha}(t)\rho L_{i,\omega}^{\alpha\dagger}(t)$
		$\displaystyle\quad-\frac{1}{2}\left\{L_{i,\omega}^{\alpha\dagger}(t)L_{i,\omega}^{\alpha}(t),\rho\right\}\Big{)}\ ,$		(27)

and the Lamb shift term

H_{\mathrm{LS}}(t)=\sum_{i,\alpha,\omega}L_{i,\omega}^{\alpha\dagger}(t)L_{i,\omega}^{\alpha}(t)S_{\mathrm{p}}(\omega)\ .

(28)

Here $\gamma_{\mathrm{p}}\left(\omega\right)$ and $S_{\mathrm{p}}(\omega)$ denote the polaron-frame noise spectrum and the corresponding principal value integral (discussed below), and the Lindblad operators are

L_{i,\omega}^{\alpha}(t)=\frac{A(t)}{2}\sum_{\varepsilon_{b}-\varepsilon_{a}=\omega}\bra{a}{\sigma^{\alpha}_{i}}\ket{b}|{a}\rangle\!\langle b|\ ,

(29)

where $\ket{a}$ is the eigenstate of $H(t)$ with eigenenergy $\varepsilon_{a}$ . Because $H(t)$ is diagonal in the $\sigma^{z}$ basis, the eigenstates $\ket{a}$ are classical spin states. This should be contrasted with the AME, for which the Lindblad operators involve the instantaneous eigenstates of the system Hamiltonian that includes the transverse field, as in Eq. (21).

Because the effective Hamiltonian $H(t)+H_{\mathrm{LS}}(t)$ [Eqs. (26) and (28)] and the Lindblad operators in Eq. (29) are diagonal in the $\sigma^{z}$ -basis, the PTRE can be simplified into an equation involving only the diagonal components of the density matrix, as shown in App. E. We use this form henceforth.

In contrast to the Ohmic noise spectrum we used for the AME [Eq. (20)], here we choose a hybrid model with a high-frequency Ohmic bath and a low-frequency component [63, 61]. Namely, $\gamma_{\mathrm{p}}(\omega)$ has a convolutional form

\gamma_{\mathrm{p}}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}G_{\mathrm{L}}(\omega-x)G_{\mathrm{H}}(x)d{x}\ ,

(30)

where


$\displaystyle G_{\mathrm{L}}(\omega)$	$\displaystyle=\sqrt{\frac{\pi}{2W^{2}}}\exp\left[-\frac{\left(\omega-4\epsilon_{L}\right)}{8W^{2}}\right]$	(31a)
$\displaystyle G_{\mathrm{H}}(\omega)$	$\displaystyle=\frac{4\gamma(\omega)}{\omega^{2}+4\gamma^{2}(0)}$	(31b)

which describe the low and high-frequency components, respectively. $W$ and $\varepsilon_{L}$ in Eq. (31a) are known as the macroscopic resonant tunneling (MRT) linewidth and reorganization energy. They are connected through the fluctuation-dissipation theorem: $W^{2}=2\varepsilon_{L}T$ . The quantity $\gamma(\omega)$ in Eq. (31b) is the standard spectral function for an Ohmic bath, given by Eq. (20). We assume that the low-frequency and high-frequency components share the same temperature.

We note that the PTRE can be approximately thought of as a multi-level extension of the noninteracting-blip approximation (NIBA) method [64, 61], which successfully modeled the open system dynamics in a multiqubit-cluster tunneling experiment [37]. Because in the low temperature limit NIBA only works for stronger coupling than PTRE [43, 65], the success of NIBA indicates the existence of a strong coupling region during the anneal.

IV.4.1 4-qubit case.

We first calibrate the simulation parameters using the empirical data from the $N=4$ case. We run the PTRE simulation with $\tau=5\;\mu$ s and $f_{c}=1$ THz, and vary $T\in\{6,30\}$ mK (step size of $1\;\mathrm{mK}$ ) and $W\in\{6,40\}$ mK (step size of $2$ mK), $\eta g^{2}\in\{2.5\times 10^{-i},5\times 10^{-i}\}_{i=1}^{5}$ . We pick the optimal parameters from this set such that the simulation results have the smallest distance from the D-Wave data in the following six cases: two cases where we start with $\ket{\uparrow\uparrow\uparrow\downarrow}$ and end up with the all-up and all-down states, and four cases where we start with $\ket{\uparrow\uparrow\downarrow\downarrow}$ and end up with the all-up, all-down, one-up and one-down states. Because all the numerical experiments are done on the same grid of inversion points $s_{\text{inv}}$ , we denote the measured population (corresponding to the grid-points) by a vector $P_{i}$ , where $i$ is the index for the six different cases. Then our parameter estimation procedure can be formally written as:

\min_{W,T,\eta g^{2}}\sum_{i}\|{\tilde{P}_{i}-P_{i}}\|\ ,

(32)

where $\tilde{P}_{i}$ is the simulation curve corresponding to the empirical data. The simulations results using the optimal parameters thus obtained are compared to the D-Wave data in Fig. 12. Panel (a) of this figure reproduces the empirical $N=4$ data from Fig. 5(a); panel (b) corresponds to the $m_{0}=0$ results shown in Figs. 2(a) and 7(c) (with $r=1$ ), though the latter are for $N=20$ . Crucially, unlike the AME the PTRE correctly describes the asymmetry in the populations of the all-up and all-down ground states when the initial state has one spin down, as is clearly visible in Fig. 12(a).

IV.4.2 8-qubit case.

Next, we simulate the 8-qubit case using the optimal parameters obtained above. Again, we consider the experiments with two different initial states and present the results in Fig. 13. Panel (a) of this figure reproduces the empirical $N=8$ data from Fig. 5(a); panel (b) corresponds to the $m_{0}=0$ results shown in Figs. 2(a) and 7(c) (with $r=1$ ), though the latter are for $N=20$ .

Similar trends are observed as in the 4-qubit case. The PTRE works best when the initial state has only a single spin flip from the ground state [Fig. 13(a)]. When an equal number of spins is up as down (the maximally excited state), the PTRE is less accurate but the qualitative trend is correct [Fig. 13(b)]. It is possible that fine-tuning of the PTRE parameters could further improve the quantitative agreement. However, we did not pursue this since the main goal has already been achieved: our simulation results demonstrate that the PTRE achieves a much better qualitative agreement with the experiments than the AME (Fig. 11). In particular, it correctly predicts the asymmetry in the populations of the two degenerate all-up and all-down ground states when the initial state is the first excited state, for both $N=4$ and $N=8$ . This is the case because in the PTRE decoherence happens between the classical spin states [Eq. (29)] while in the AME, decoherence happens between the eigenstates of the full Hamiltonian [Eq. (24)].

V Conclusions

In this work, we presented a comprehensive study of reverse annealing of the $p$ -spin model with $p=2$ , and compared empirical results from the D-Wave 2000Q used as a quantum simulator of this model, with two quantum master equations in the weak and strong coupling limits (the AME and PTRE, respectively). Results for one classical Monte Carlo based algorithm (SVMC-TF [50]) are presented in App. F; to keep the scope manageable we did not consider other popular models such as simulated quantum annealing (SQA) [66, 67, 68, 35, 69, 70] (also known as Path-Integral Monte Carlo Quantum Annealing in Refs. [71, 5] or Quantum Monte Carlo Annealing in Ref. [4]), or the recent Spin-Vector Langevin model [72]; these are interesting topics for a future study.

Our main observation is a stark failure of the AME, which successfully captured the D-Wave annealers’ behavior in a variety of previous studies [34, 40, 73, 60, 13, 50], in correctly describing the empirical “partial” success probabilities, i.e., the probability of finding either the all-up or all-down ground state, as opposed to the sum of the probabilities of these two degenerate ground states of the $p=2$ $p$ -spin model. As illustrated in Fig. 11, the AME predicts equal partial success probabilities, but the empirical data are strongly biased towards one or the other ground state, depending on the initial condition of the reverse anneal; see, e.g., Fig. 2 and all subsequent figures displaying empirical partial success probability data. The reason for this failure of the AME is that in the weak coupling limit that it describes, the two degenerate ground states are energy eigenstates that are equal superpositions of the all-up and all-down states with opposite signs, and there is nothing about the dynamics described by the AME that breaks the symmetry between these two states for the $p=2$ problem under consideration. The failure of the AME thus indicates that the weak coupling limit itself is the problem in the setting of the D-Wave 2000Q device subject to reverse annealing, and indeed, when we considered the PTRE (strong coupling) it exhibited the observed asymmetry; see, e.g., Fig. 12. This represents strong evidence of the breakdown of the weak coupling limit in the D-Wave 2000Q device, at least for reverse annealing on timescales of a $\mu$ sec or longer.

Our results demonstrate the importance of choosing the correct decoherence model when analyzing real-world devices, and that even when agreement is observed (as was the case for previous studies involving the AME), certain aspects that can reveal a way to distinguish between different decoherence models may remain hidden. The case in point is the difference between the partial probabilities of the two degenerate ground states, which forced us to conclude that strong coupling is the more appropriate decoherence model for the D-Wave 2000Q on the $\gtrsim 1\mu$ sec timescale. Had we considered only the total success probability, we would not have been able to distinguish between decoherence models with weak vs strong coupling.

Our work shows that the PTRE – a first principles, fully quantum dynamical model with strong decoherence – achieves the best agreement overall with the empirical data among the various models we have tested. This does not necessarily imply that the D-Wave 2000Q behaves as a classical device on the $\gtrsim 1\mu$ sec timescale. Instead, our results should be interpreted to mean that models with strong decoherence can be successful in predicting the outcome of quantum annealing experiments when the chosen observable is not sensitive to quantum effects. An earlier study using a previous generation of the D-Wave devices already established evidence for entanglement [74], which was verified using AME simulations [60], and this is a clear-cut example of the measurement of an observable that is sensitive to quantum effects. We anticipate that new experiments based on much shorter anneal times than were available to us in this work will provide further evidence of quantum effects in experimental quantum annealing [75]. Our work highlights the importance of testing such new evidence using models that critically evaluate the strength of any claimed quantum effects.

Acknowledgements.

We thank Mohammad Amin and Masayuki Ohzeki for useful comments and Sigma-i Co., Ltd. for providing part of the D-Wave 2000Q machine time. The research is based partially upon work supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA) and the Defense Advanced Research Projects Agency (DARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, DARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Computation for some of the work described in this paper was supported by the University of Southern California Center for High-Performance Computing and Communications (hpcc.usc.edu).

Appendix A Optimal Value of Ferromagnetic Interactions within a Logical Qubit

The ferromagnetic coupling strength between physical qubits for a given logical qubit on the D-Wave device affects the final success probability. The appropriate value of $J_{F}$ depends on the system size $N$ [76]. We checked the success probability for various system sizes of the $p=2$ $p$ -spin model under traditional forward annealing. The result is shown in Fig. 14. It is clearly observed that the smallest value of coupling $J_{\text{F}}=-1.0$ we tried yields the best result for all $N$ values we tried, and we adopted this value in all our experiments. The smallest allowed value is $J_{\text{F}}=-2.0$ , but saturation is already observed for $J_{\text{F}}=-1$ .

Appendix B Reverse annealing details

In our experiments for single-iteration ( $r=1$ ) reverse annealing, we constructed 15 instances (different hardware embeddings of $H_{T}$ on the Chimera graph of the device) and generated 10 random gauges for each set of parameter values. A gauge (originally known as a spin inversion transformation [34]) is a given choice of $\{h_{i},\,J_{ij}\}$ in the general Ising problem Hamiltonian $H_{T}=\sum_{i=1}^{N}h_{i}\sigma^{z}_{i}+\sum_{i<j}^{N}J_{ij}\sigma^{z}_{i}\sigma^{z}_{j}$ ; a new gauge is realized by randomly selecting $a_{i}=\pm 1$ and performing the substitution $h_{i}\mapsto a_{i}h_{i}$ and $J_{ij}\mapsto a_{i}a_{j}J_{ij}$ . Provided we also perform the substitution $\sigma_{i}^{z}\mapsto a_{i}\sigma_{i}^{z}$ , we map the original Hamiltonian to a gauge-transformed Hamiltonian with the same energy spectrum but where the identity of each energy eigenstates is relabeled accordingly. In total, there are $2^{N}$ different gauges for an $N$ -spin problem. See Ref. [77] for a review and more details.

1000 annealing runs were performed for a given random gauge. Thus, the total number of runs for each set of annealing schedule parameter values is 150,000. In iterated reverse annealing with $r>1$ , we constructed 15 instances and generated 90 random gauges for each $r$ . Thus, the total number of samples at each $r$ and each $s_{\mathrm{inv}}$ is 1350.

We computed $1\sigma$ error bars by cluster sampling over instances. Namely, we regarded each of the $15$ instances as as a sample of size 10000 [the number of random gauges (10) $\times$ the number of iterations (1000)], and computed the standard error of the mean for each data point shown in our experimental figures.

We remark that the reverse annealing protocol we have adopted here for experiments on the D-Wave device is different from the one used in Ref. [45] for numerical computation. In particular, when $s_{\mathrm{inv}}\approx 0$ , the second part of Ref. [45]’s protocol has a very sharp increase of $s$ as a function of $t$ , which seems to have led to a significant drop of success probability. This is not the case in the present protocol, where the time derivative of $s(t)$ is a constant $\pm 1/\tau$ (or $0$ during pausing) irrespective of $s_{\mathrm{inv}}$ .

Appendix C Spectrum of the $p$ -spin problem with $p=2$ and scaling of the minimum gap with $N$

The spectrum of the $p$ -spin problem and the ordering of energies of different spin sectors can be found in Ref. [24]. For the particular case of $n=4,p=2$ and the D-Wave 2000Q annealing schedule, we plot the spectrum in Fig. 15.

We denote the energy gap between the instantaneous eigenenergies $\epsilon_{i}(s)$ and $\epsilon_{j}(s)$ by $\Delta_{ij}(s)=\epsilon_{i}(s)-\epsilon_{j}(s)$ , and the corresponding minimum energy gap by $\Delta_{ij}=\min_{s}\Delta_{ij}(s)$ . For $p=2$ , we are interested instead in the minimum energy gap between the ground state $\epsilon_{0}(s)$ and the second excited state $\epsilon_{2}(s)$ , denoted by

\Delta=\Delta_{20}=\min_{s}\Delta_{20}(s)=\min_{s}\epsilon_{2}(s)-\epsilon_{0}(s)\,,

(33)

since as can be seen from Fig. 15, the instantaneous first excited state $\epsilon_{1}(s)$ and the instantaneous ground state $\epsilon_{0}(s)$ converge at $s=1$ . This, of course, is true for every $N$ due to the double degeneracy of the ground state of the $p=2$ $p$ -spin model, which exhibits $\mathbb{Z}_{2}$ symmetry.

Figure 6 shows the value of $\Delta$ for $N\in\{4,\cdots,22\}$ , along with the position $s$ of the minimum gap, i.e.,

s_{\Delta}=\mathrm{argmin}\Delta_{20}(s)\,.

(34)

Appendix D Upper bound on the success probabilities for collective dephasing

If the dynamics preserve spin symmetry (for example in the closed system Schrödinger equation and in open-system simulations with collective system-bath coupling), then there exists a natural upper bound on the maximal success probabilities achievable in reverse annealing. The upper bound is the population of the initial state in the maximum-spin sector.

For the example of $N=4$ , there are two degenerate ground states ( $\ket{0000}$ and $\ket{1111}$ ) of the problem Hamiltonian $H_{T}$ and they both belong to the maximum-spin subspace of $S=S_{\max}={N}/{2}=2$ . While the computational basis state with one spin down, for example, $\ket{0001}$ ( $\ket{0}\equiv\ket{\uparrow}$ , $\ket{1}\equiv\ket{\downarrow}$ ) is a first excited state of $H_{T}$ , it does not belong to the maximum-spin subspace $S=2$ . A uniform superposition of the computational basis states with one spin down, i.e., $(\ket{0001}+\ket{0010}+\ket{0100}+\ket{1000})/2$ , however, does belong to the maximum-spin subspace $S=2$ . Unfortunately, the D-Wave device does not allow such an initialization.

In general, suppose that the initial computational basis state is $\ket{\psi(t=0)}_{\text{comp}}$ , with a particular magnetization $m_{0}$ [Eq. (7)]. It can be represented as a linear combination of states with a fixed value of total spin $S$ and magnetization $m_{0}$ :

\ket{\psi(t=0)}_{\text{comp}}=\sum_{S=S_{\min}}^{S_{\max}}a_{S}\ket{S,m_{S}=m_{0}}.

(35)

From angular momentum addition theory [78], we know that $S_{\min}=N/2-\lfloor N/2\rfloor$ and $S_{\max}=N/2$ . The total spin (integer or half-integer) $S\in\{S_{\min},S_{\min}+1,\dots,S_{\max}\}$ . The (unnormalized) magnetization is $M_{S}\in\{-S,-S+1,\dots,S\}$ , but we can also label the states using the normalized magnetization $m_{S}=\frac{2}{N}M_{S}$ ; then $\ket{S,m_{S}=m_{0}}$ is a simultaneous eigenstate of $\mathbf{S}^{2}$ and $S^{z}$ with eigenvalues

	$\displaystyle\mathbf{S}^{2}\ket{S,m_{S}=m_{0}}$	$\displaystyle=S(S+1)\ket{S,m_{S}=m_{0}}\,,$		(36)
	$\displaystyle S^{z}\ket{S,m_{S}=m_{0}}$	$\displaystyle=\Big{(}\frac{N}{2}m_{0}\Big{)}\ket{S,m_{S}=m_{0}}\,.$		(37)

In Eq. (35), $|a_{S}|^{2}$ is the initial state’s population in that spin subspace.

For the $p$ -spin Hamiltonian, in a closed system the state in each spin subspace evolves independently due to the preservation of spin symmetry of the Hamiltonian [26]. At the end of a single cycle of reverse annealing ( $r=1$ ), we have from Eq. (35)) that

U(2t_{\text{inv}},0)\ket{S,m_{S}=m_{0}}=\sum_{M_{S}=-S}^{S}c_{M}\ket{S,m_{S}=2M_{S}/N},

(38)

where $c_{M}$ is the amplitude developed by the other basis elements of the spin $S$ subspace at time $t=2t_{\text{inv}}$ .

Therefore, the final state of $1$ -cycle reverse annealing can be expressed as:

	$\displaystyle\ket{\psi}_{\text{fin}}$	$\displaystyle=U(2t_{\text{inv}},0)\ket{\psi(t=0)}_{\text{comp}}$
		$\displaystyle=\sum_{S=S_{\min}}^{S_{\max}}a_{S}\left(\sum_{M_{S}=-S}^{S}c_{M}\ket{S,m_{S}=2M_{S}/N}\right)\,.$		(39)

The all-up state ( $\ket{0}^{\otimes N}=\left|\uparrow\right>^{\otimes N}=\ket{\text{up}}$ ) and all-down state ( $\ket{1}^{\otimes N}=\left|\downarrow\right>^{\otimes N}=\ket{\text{down}}$ ) are both ground states of $H_{T}$ , and moreover they lie in the maximum-spin subspace $S=S_{\max}$ . In particular, $\ket{\text{up}}=\ket{S_{\max},1}$ and $\ket{\text{down}}=\ket{S_{\max},-1}$ . Therefore, projecting $\ket{\psi}_{\text{fin}}$ onto $\ket{\text{up}}$ gives:

	$\displaystyle\phantom{==}\braket{\text{up}}{\psi}_{\text{fin}}$
	$\displaystyle=\sum_{S=S_{\min}}^{S_{\max}}a_{S}\left(\sum_{M_{S}=-S}^{S}c_{M}\braket{\text{up}}{S,m_{S}=2M_{S}/N}\right)$		(40)
	$\displaystyle=a_{S_{\max}}c_{S_{\max}}\,.$		(41)

Similarly, $\braket{\text{down}}{\psi}_{\text{fin}}=a_{S_{\max}}c_{-S_{\max}}$ .

The total success probability is thus bounded by:

$\displaystyle p(r=1)$	$\displaystyle=\|\braket{\text{up}}{\psi}_{\text{fin}}\|^{2}+\|\braket{\text{down}}{\psi}_{\text{fin}}\|^{2}$
	$\displaystyle=\|a_{S_{\max}}\|^{2}(\|c_{S_{\max}}\|^{2}+\|c_{-S_{\max}}\|^{2})$
	$\displaystyle\leq\|a_{S_{\max}}\|^{2}\,,$	(42)

with the bound saturated when $(|c_{S_{\max}}|^{2}+|c_{-S_{\max}}|^{2})=1$ .

$\square$

The upper bound $|a_{S_{\max}}|^{2}$ is, as claimed above, the population of the initial state in the maximum-spin subspace. For the initial state $\ket{0001}$ we have $a_{S_{\max}}=1/2$ since $\ket{S=S_{\max}=2,m_{S}=0.5}=(\ket{0001}+\ket{0010}+\ket{0100}+\ket{1000})/2$ . Therefore, the total success probability is bounded by $|a_{S_{\max}}|^{2}=1/4$ . For the other examples in the main text, the initial state $\ket{0011}$ has $a_{S_{\max}}=1/{\sqrt{{{4}\choose{2}}}}=1/{\sqrt{6}}$ ; while the initial state of $\ket{00000001}$ has $a_{S_{\max}}=1/{\sqrt{8}}$ .

This conclusion is directly generalized to the open-system case under collective dephasing, where the population of each spin-sector is also preserved during the dynamics. This is illustrated in the main text for $N=4$ in Fig. 10(b), and also in Fig. 16, which shows the analogous result for $N=8$ , for which the initial state is $\ket{00000001}$ . In this case the maximum success probability of collective dephasing is $1/8$ , which is the population of the initial state in the maximal spin subspace. We remark that this case is different from Ref. [45], where the initial state (a superposition of computational basis state) lies completely inside the maximal spin subspace, in which case the collective dephasing model has a maximum success probability of $1$ .

Appendix E The PTRE equation involves only the diagonal components of the density matrix

As argued in Sec. IV.4, the PTRE equation involves only the diagonal components of the density matrix. To show this we first consider

M=[H_{\mathrm{e}},{\rho}]\ ,

(43)

where $H_{\mathrm{e}}\equiv H+H_{\mathrm{LS}}$ is a diagonal matrix, i.e., $H_{\mathrm{e}}^{ab}=\delta_{ab}H^{aa}_{\mathrm{e}}$ . The explicit form of $M$ can be written as $M_{ab}=\left(H_{\mathrm{e}}^{aa}-H_{\mathrm{e}}^{bb}\right){\rho}^{ab}$ , which implies that the diagonal elements of $M$ ( $M_{aa}$ ) vanish. Next, we consider

N^{\omega,\alpha}_{i}=L_{i,\omega}^{\alpha}{\rho}L_{i,\omega}^{\alpha\dagger}-\frac{1}{2}\left\{L_{i,\omega}^{\alpha\dagger}L_{i,\omega}^{\alpha},{\rho}\right\}\ ,

(44)

where the Lindblad operators $L_{i,\omega}^{\alpha}(t)$ were defined in Eq. (29). To simplify the notation, we omit the indices $\omega$ , $\alpha$ and $i$ in the following discussion. Let us denote $A_{ab}\equiv\frac{A\left(t\right)}{2}\bra{a}{\sigma_{i}^{\alpha}}\ket{b}=A^{*}_{ab}$ . The first term in Eq. (44) can be written as

\displaystyle L{\rho}L

\displaystyle=\sum_{a,b,a^{\prime},b^{\prime}}A_{ab}A_{a^{\prime}b^{\prime}}{\rho}^{bb^{\prime}}|{a}\rangle\!\langle a^{\prime}|\ ,

(45)

where ${\rho}^{bb^{\prime}}=\bra{b}{\rho}\ket{b^{\prime}}$ . The second term in Eq. (44) is


$\displaystyle-\frac{1}{2}\left\{L^{\dagger}L,{\rho}\right\}$	$\displaystyle=-\sum_{a,b,b^{\prime}}\frac{1}{2}A_{ab^{\prime}}A_{ab}\{\|{b^{\prime}}\rangle\!\langle b\|,{\rho}\}$	(46a)
	$\displaystyle=-\sum_{a,b,b^{\prime}}\frac{1}{2}A_{ab}A_{ab^{\prime}}\left(\|{b^{\prime}}\rangle\!\langle b\|{\rho}+{\rho}\|{b^{\prime}}\rangle\!\langle b\|\right)\ .$	(46b)

If we assume ${\rho}$ is diagonal, then all the diagonal elements of $\dot{\rho}$ depend only on the diagonal elements of ${\rho}$ , with an explicit form

\dot{\rho}^{aa}=\sum_{b\neq a}\gamma_{\mathrm{p}}(\omega_{ba})Z_{ab}\rho^{bb}-\sum_{b\neq a}\gamma_{\mathrm{p}}(\omega_{ab})Z_{ba}\rho^{aa}\ ,

(47)

where $Z_{ab}=A^{2}(t)\sum_{\alpha,i}|\bra{a}{\sigma_{i}^{\alpha}}\ket{b}|^{2}/4$ , and $\omega_{ab}=\omega_{a}-\omega_{b}$ . Thus, if the initial state is diagonal (a classical spin state), we can consider the Pauli master equation [79] for the diagonal elements only [Eq. (47)] to speed up the computation. The matrix form of this equation is


	$\displaystyle\begin{pmatrix}\dot{\rho}_{00}\\ \dot{\rho}_{11}\\ \vdots\end{pmatrix}=T\begin{pmatrix}\rho_{00}\\ \rho_{11}\\ \vdots\end{pmatrix}$		(48a)
	$\displaystyle T\equiv\begin{pmatrix}-\sum_{b\neq 0}\gamma_{\mathrm{p}}(\omega_{0b})Z_{b0}&\gamma_{\mathrm{p}}(\omega_{10})Z_{01}&\cdots\\ \gamma_{\mathrm{p}}(\omega_{01})Z_{10}&-\sum_{b\neq 1}\gamma_{\mathrm{p}}(\omega_{1b})Z_{b1}&\cdots\\ \vdots&\vdots&\ddots\\ \end{pmatrix}\ .$		(48b)

Before proceeding, we note that the sparsity of $T$ is determined by the sparsity of $Z_{ab}$ , whose full size is $2^{N}\times 2^{N}$ . To calculate the number non-zero elements in $Z_{ab}$ , we first notice that for each $\sigma_{i}^{\alpha}$ operator, $2^{N-1}$ elements are non-zero (recall that $\alpha\in\{+,-\}$ ). The number of $\sigma_{i}^{\alpha}$ operators is $2N$ . So there are $N2^{N}$ non-zero elements in $Z_{ab}$ . If we add the number of diagonal elements in $T$ , the total number of non-zero elements in $T$ is $(N+1)2^{N}$ . As a result, the sparsity of the transfer matrix $T$ is $(N+1)/2^{N}$ .

The transfer matrix $T$ in Eq. (48) provides the incoherent “tunneling” rate between classical spin states. Because the total Hamiltonian is symmetric under permutations of qubits, we may also want to calculate the tunneling rate between the spin coherent states [80]

|\theta,\phi\rangle=\bigotimes_{i=1}^{n}\left[\cos\left(\frac{\theta}{2}\right)|0\rangle_{i}+\sin\left(\frac{\theta}{2}\right)e^{i\varphi}|1\rangle_{i}\right]\ .

(49)

However, because in the current form of the PTRE the coherence between computational basis decays exponentially, the spin coherent state cannot survive.

Complementing the results shown in the main text, the differences (in absolute value) between the simulation and experimental results are shown in Fig. 17.

Appendix F Fully classical simulations using the spin-vector Monte Carlo algorithm

In an effort to better understand the asymmetric partial success probability observed in our experiments, we also performed fully classical simulations of the same problem using the spin-vector Monte Carlo (SVMC) [39] algorithm and a new variant with transverse-field-dependent updates (SVMC-TF) [50], where this variant was successfully used to explain empirical D-Wave results for a particular $12$ -qubit instance. For the $p$ -spin problem, we replace the Hamiltonian of Eq. (1) by a classical Hamiltonian:

\mathcal{H}(s)=-\frac{A(s)}{2}\left(\sum_{i}^{N}\sin\theta_{i}\right)-\frac{B(s)N}{2}\left(\frac{1}{N}\sum_{i}^{N}\cos\theta_{i}\right)^{p}\,.

(50)

Each qubit $i$ is replaced by a classical $O(2)$ spin $\vec{M}_{i}=(\sin\theta_{i},0,\cos\theta_{i})$ , $\theta_{i}\in[0,\pi]$ . For the purpose of reverse annealing, we also need to specify the $t$ dependence of $s(t)$ . The concept of time is here replaced by the number of Monte Carlo sweeps: we replace $\tau$ by a specified number of total sweeps. The total number of sweeps is then $2\tau(1-s_{\text{inv}})$ , in analogy to the total annealing time in Eq. (6).

To simulate the effect of thermal hopping through this semi-classical landscape with inverse temperature $\beta$ , we perform at each time step a spin update according to the Metropolis rule. In SVMC, a random angle $\theta_{i}^{{}^{\prime}}\in[0,\pi]$ is picked for each spin $i$ . Updates of the spin angles $\theta_{i}$ to $\theta_{i}^{{}^{\prime}}$ are accepted according to the standard Metropolis rule associated with the change in energy ( $\Delta E$ ) of the classical Hamiltonian. For the $p$ -spin problem, $\Delta E$ cannot be expressed in a simple form as in case of the Ising problem Hamiltonian [39].

In SVMC-TF, the random angle $\theta_{i}^{{}^{\prime}}=\theta_{i}+\epsilon_{i}(s)$ is picked in a restricted range

\displaystyle\epsilon_{i}(s)\in\left[-\min\left(1,\frac{A(s)}{B(s)}\right)\pi,\min\left(1,\frac{A(s)}{B(s)}\right)\pi\right].

(51)

The goal of SVMC-TF is to restrict the angle update for $A(s)<B(s)$ , and imitate the freeze-out effect discussed in Sec. III.1.3. The full SVMC and SVMC-TF algorithms for reverse annealing are summarized in App. I, including the expression for $\Delta E$ .

A simple intuitive way to visualize the semiclassical dynamics described by the SMVC and SVMC-TF algorithms is to consider the energy landscape defined by $\mathcal{H}(s)$ [Eq. (50)] when equating all angles ( $\theta_{i}\equiv\theta$ ). We plot the resulting surface in Fig. 18 for the simple case of $N=1$ . For $s_{\mathrm{inv}}<0.3948$ , Metropolis updates to either spin direction happen frequently and are equally likely since there is no potential barrier separating them. For $1\gg s_{\mathrm{inv}}>0.3948$ , under Metropolis updates the system prefers staying in the original well and escaping to the opposite well is unlikely, due to the potential barrier. This preference explains the asymmetry part of partial success probabilities in the experimental results. For $1>s_{\mathrm{inv}}\gg 0.3948$ , Metropolis updates are very rare due to the rate suppression in SVMC-TF.

F.1 Total success probability

We report on SVMC and SVMC-TF simulations for $N=4$ and $8$ . We choose the initial condition with a single spin down. In terms of angles, for $N=4$ the initial angles are $\{0,0,0,\pi\}$ . We again use the temperature $T=12.1$ mK. The classical analogue of the annealing time is chosen to be $\tau=10^{3}$ and $10^{4}$ sweeps.

In Fig. 19 we display the simulation results for the total success probability using SVMC and SVMC-TF. The number of samples is $10^{4}$ .

SVMC gives high total success probabilities even with large inversion points $s_{\text{inv}}$ . A large $s_{\text{inv}}$ value means that during the whole reverse annealing process the ratio $A(s)/B(s)$ is small. In the D-Wave device, it is expected that for small $A(s)/B(s)$ the dynamics freezes, which makes it difficult to reach the ground state(s) when the initial state is excited. With the number of sweeps increased from $\tau=10^{3}$ to $10^{4}$ , we observe that even for very large inversion points $s_{\text{inv}}$ the total success probability of SVMC can be as high as $1$ . This is because the angle updates in SVMC are completely random and thus with a sufficient number of sweeps, it is possible for the state to flip to the correct solutions.

However, in SVMC-TF, the range of angle updates is restricted for $A(s)/B(s)<1$ . The restricted angle updates (freeze-out effect) prevent the state from flipping to the correct solutions. Therefore the total success probability for large inversion points $s_{\text{inv}}$ is basically zero in SVMC-TF, regardless of the number of sweeps. This is also what we observe in the empirical data and in the adiabatic master equation simulations, as discussed in Sec. IV.3.1.

F.2 Partial success probability

We compare the partial success probability obtained from SVMC and SVMC-TF in Figs. 20(a) and (b), respectively, for $N=4$ and $\tau=10^{3}$ . It is again seen that SVMC-TF more accurately captures the early freezing than SVMC. The empirical data from Fig. 5(a) is reproduced in Fig. 20(c) and (d), where is it compared to SVMC-TF at an optimized number of sweeps (explained in App. H). This yields a semi-quantitative agreement, in particular the correct trend and transition locations for the unequal up and down partial success probabilities. However, the agreement is clearly better for $N=4$ than for $N=8$ despite the optimization, which suggests the interesting possibility that SVMC-TF becomes less accurate at higher numbers of spins. Also noteworthy is that for $N=8$ we observe a deviation from the empirical data for $s_{\text{inv}}\lesssim 0.3$ , where there exists a small but clear difference in the probabilities of all-up and all-down states, whereas the SVMC-TF data do not show such a trend. As discussed in Sec. III.1.4, we attribute the difference for $s_{\text{inv}}\lesssim 0.3$ to the spin-bath polarization effect, which is not modeled in our SVMC-TF simulations.

In Fig. 21 we display SVMC-TF reverse annealing simulation results of partial success probabilities for $N=16,32$ with $\tau=10^{3},10^{4}$ sweeps. For both sizes shown, the regime of high partial success probability for the all-up state is shifted slightly to higher $s_{\textrm{inv}}$ for $\tau=10^{4}$ sweeps than for $\tau=10^{3}$ . This is consistent with the trend in the experimental results observed in Fig. 3(a). However, the trend with system size is inconsistent with the empirical partial success probabilities for the up case shown in Fig. 5(a): the numerical results for $N=16$ and $N=32$ are virtually indistinguishable (apart from statistical fluctuations), while the empirical data shows that $N=16$ is not yet large enough for convergence. This (small) failure of the SVMC-TF model may hint at an interesting way in which to identify a “quantum signature” in experimental quantum annealing [34, 40]. However, we did not pursue this direction further since we cannot rule out that further fine-tuning of the SVMC parameters will result in a closer match with the empirical data. To further explore ways in which a classical model such as SVMC-TF, or a model with strong coupling to the bath such as PTRE, might fail in describing the empirical data, we focus on excited state populations in the next Appendix.

Appendix G Excited states

In this Appendix we present results comparing D-Wave data to simulations for the population in low-lying excited states. Our goal is not to be comprehensive, but rather to highlight agreements and discrepancies between the empirical and the numerical results.

The overall conclusion of this Appendix is that because of a persistence of small discrepancies even for the PTRE and SVMC-TF, in particular their failure to accurately predict the excited-state populations as shown below, further work is needed in order to improve open-system models.

In both the empirical and the simulation data, the population of the $i$ ’th excited up (down) state is obtained after summation over all permutations of computational basis states with $N-i$ and $i$ up, or $i$ and $N-i$ down spins, where $N$ is the number of spins. We chose $\tau=5\mu s$ for the D-Wave experiments.

G.1 $N=4$ with a maximally excited initial state

We compare in Fig. 22 the AME, PTRE, and SVMC-TF simulation results to the empirical data for the initial state $\ket{0011}$ . The open-system parameter settings are the same as in the ground state simulations reported above, for all three simulation methods, but the number of sweeps used in SVMC-TF is optimized for the closest agreement with the D-Wave data, as explained in App. H.

Since the initial state has no bias toward spin up or down, it is surprising that the D-Wave data exhibits an asymmetry between probability of ending in states with one up or one down spin. Since the same anomaly is not observed for $N=8$ (see Fig. 24), we attribute it to an unexplained peculiarity associated with the embedding of the $N=4$ problem. All three simulation methods correctly predict no distinction between the up and down states. The predicted position of the peak in the first excited state population is different for the three simulation methods; the AME [panel (a)] predicts a position that is roughly the average of the empirically observed peaks for the case where the final state has a down or an up spin, while the PTRE [panel (b)] and SVMC-TF [panel (c)] are shifted to the right and left, respectively. It is not possible to say which prediction is correct due to the aforementioned anomaly.

The partial success probability results are shown in Fig. 23 for the second excited state as well as the ground state, for the AME and SVMC-TF. From (a) we see that the AME is qualitatively but not quantitatively in agreement with the empirical results for the second excited state, where it predicts an $s_{\text{inv}}$ value that is too large for the onset of the rise in the population of this state, and this rise is also somewhat too steep. The same is true for SVMC-TF [panel (b)] but it is in slightly closer agreement than AME for the second excited state. We also show results for $N=8$ [panel (c)], where SVMC-TF continues to exhibit good qualitative agreement with the empirical results.

G.2 $N=8$ with a maximally excited initial state

In Fig. 24 we display the results for $N=8$ with $\ket{\uparrow\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\downarrow}$ as the initial state, for the probability of ending in the first, second, or third excited state, using the PTRE and SVMC-TF.⁴⁴4We do not present AME results since the computational cost of $N=8$ highly excited states using the AME is prohibitive: all $256$ eigenstates need to be taken into account. This time the empirically observed population in the states with $i$ spins up or $i$ spins down ( $1\leq i\leq 3$ ) is identical, as expected, i.e., we do not observe the anomaly mentioned above for $N=4$ .

The top row shows the results for the PTRE, with the same set of optimal parameters as explained in Sec. IV.4. The agreement is relatively poor, in that both the magnitude and the position of the population peak is missed, both being systematically overestimated.

The bottom row shows the results for SVMC-TF, with the optimal number of sweeps as determined in App. H. The agreement is somewhat better than for the PTRE, in that both the peak’s position and magnitude are closer to the empirical data, though the agreement in the peak’s position deteriorates with increasing excitation level.

G.3 First excited state as initial state

As a final test, we checked the PTRE and SVMC-TF for initial states with one spin down, i.e., the first excited state. The results for $N=4$ and $N=8$ are shown in Fig. 25, and are in reasonable qualitative agreement. The agreement is overall somewhat better for the PTRE, especially for $N=8$ . For $N=4$ the PTRE predicts a small non-zero probability for the first excited down state near the minimum gap point, which is absent in the empirical data and in the SVMC-TF results. This suggests that the PTRE slightly overestimates the incoherent tunneling rates for small $N$ .

Appendix H SVMC-TF’s dependence on the number of sweeps

For the $i$ ’th eigenstate population, the $\ell_{2}$ norm between the experiment and simulation data series (over a range of $s_{\text{inv}}$ ) is:

\ell^{i}_{2}=\sqrt{\sum_{k=1}^{44}(x^{i,\text{exp}}_{k}-x^{i,\text{sim}}_{k})^{2}}\ ,

(52)

where for the $k$ ’th data point, $x^{i}_{k}$ is the population of the $i$ ’th eigenstate at $s_{\text{inv}}=0.02k$ at the end of the reverse anneal (recall that the experimental data is sampled at every point from $\{0.02,0.04,\dots,0.88\}$ ).

To rigorously evaluate the differences between the experiment data and simulation data, we consider the sum of such $\ell_{2}$ norms over all eigenstates of the Hamiltonian, i.e., $\ell_{2}=\sum_{i}\ell^{i}_{2}$ . We plot $\ell_{2}$ vs sweeps in Fig. 26. We observe that for the initial states shown, with the exception of $\ket{\uparrow\cdots\uparrow\downarrow}$ , we can find an optimal sweep number. Considering both initial states at the respective system size, the optimal number of sweeps is $450$ for $N=4$ and $150$ for $N=8$ .

Appendix I Pseudocode for SVMC and SVMC-TF

Algorithm 1 SVMC and SVMC-TF (

p

-spin reverse annealing)

\mathcal{H}(s)=-\frac{A(s)}{2}(\sum_{i=1}^{N}\sin\theta_{i})-\frac{B(s)N}{2}\left(\frac{1}{N}\sum_{i=1}^{N}\cos\theta_{i}\right)^{p}

, with a specified reverse annealing dependence

s(t)

procedure SVMC

for

k=1

K

(number of samples) do

Initial computational basis state:

\ket{0}_{i}\rightarrow\theta^{k}_{i,t}:0

\ket{1}_{i}\rightarrow\theta^{k}_{i,t}:\pi

for

t=1

T

(total number of sweeps) do

for

i=1

n

(number of qubits) do

Randomly choose a new angle

\theta_{i,t}^{{}^{\prime}k}\in[0,\pi]

Calculate

\Delta E

\Delta E\leq 0

then

\theta_{i,t}^{k}\rightarrow\theta_{i,t}^{{}^{\prime}k}

else if

p<\exp(-\beta\Delta E)

where

p\in[0,1]

is drawn with uniform probability. then

\theta_{i,t}^{k}\rightarrow\theta_{i,t}^{{}^{\prime}k}

end if

end for

Take the mean of

K

samples:

\bar{\theta}_{i,t}=(\sum_{k=1}^{K}\theta_{i,t}^{{}^{\prime}k})/K

end for

end procedure

return

\bar{\theta}_{i,t}

procedure SVMC-TF

for

k=1

K

(number of samples) do

Initial computational basis state:

\ket{0}_{i}\rightarrow\theta^{k}_{i,t}:0

\ket{1}_{i}\rightarrow\theta^{k}_{i,t}:\pi

for

t=1

T

(total number of sweeps) do

for

i=1

n

(number of qubits) do

\theta_{i,t}^{{}^{\prime}k}=\theta_{i,t}^{k}+\epsilon_{i}(s(t/T))

a random

\epsilon_{i}(s(t/T)\in[-\min\left(1,\frac{A(s(t/T)}{B(s(t/T)}\right)\pi,\min\left(1,\frac{A(s(t/T)}{B(s(t/T)}\right)\pi]

Calculate

\Delta E

\Delta E\leq 0

then

\theta_{i,t}^{k}\rightarrow\theta_{i,t}^{{}^{\prime}k}

else if

p<\exp(-\beta\Delta E)

where

p\in[0,1]

is drawn with uniform probability. then

\theta_{i,t}^{k}\rightarrow\theta_{i,t}^{{}^{\prime}k}

end if

end for

Take the mean of

K

samples:

\bar{\theta}_{i,t}=(\sum_{k=1}^{K}\theta_{i,t}^{{}^{\prime}k})/K

end for

end procedure

return

\bar{\theta}_{i,t}

procedure Projection onto computational basis

0\leq\theta_{i,t}^{k}\leq\pi/2

then

\ket{\psi^{k}(t)}_{i}=\ket{0}

else if

\pi/2\leq\theta_{i,t}^{k}\leq\pi

then

\ket{\psi^{k}(t)}_{i}=\ket{1}

end if

end procedure

Remark:

\Delta E

due to the update of qubit

k

, can be expressed in the following form:

\displaystyle\Delta E_{k}

\displaystyle=\left(-\frac{A(s)}{2}\left(\sum_{\begin{subarray}{c}i=1\\ i\neq k\end{subarray}}^{N}\sin\theta_{i}+\sin\theta^{{}^{\prime}}_{k}\right)-\frac{B(s)N}{2}\left(\frac{1}{N}\left(\sum_{\begin{subarray}{c}i=1\\ i\neq k\end{subarray}}^{N}\cos\theta_{i}+\cos\theta^{{}^{\prime}}_{k}\right)\right)^{2}\right)-\mathcal{H}(s)\,.

References

Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Phys. Rev. E 58, 5355 (1998).
Farhi et al. [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem., Science 292, 472 (2001).
Santoro et al. [2002] G. E. Santoro, R. Martoňák, E. Tosatti, and R. Car, Theory of quantum annealing of an Ising spin glass, Science 295, 2427 (2002).
Das and Chakrabarti [2008] A. Das and B. K. Chakrabarti, Colloquium: Quantum annealing and analog quantum computation, Rev. Mod. Phys. 80, 1061 (2008).
Morita and Nishimori [2008] S. Morita and H. Nishimori, Mathematical foundation of quantum annealing, J. Math. Phys. 49, 125210 (2008).
Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
Hauke et al. [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Perspectives of quantum annealing: Methods and implementations, Rep. Prog. Phys. 83, 054401 (2020).
Crosson and Lidar [2021] E. J. Crosson and D. A. Lidar, Prospects for quantum enhancement with diabatic quantum annealing, Nature Reviews Physics 3, 466 (2021).
Johnson et al. [2011] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Quantum annealing with manufactured spins, Nature 473, 194 (2011).
Harris et al. [2018] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Phase transitions in a programmable quantum spin glass simulator, Science 361, 162 (2018).
King et al. [2018] A. D. King, J. Carrasquilla, J. Raymond, I. Ozfidan, E. Andriyash, A. Berkley, M. Reis, T. Lanting, R. Harris, F. Altomare, K. Boothby, P. I. Bunyk, C. Enderud, A. Fréchette, E. Hoskinson, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich, Y. Sato, A. Y. Smirnov, L. J. Swenson, M. H. Volkmann, J. Whittaker, J. Yao, E. Ladizinsky, M. W. Johnson, J. Hilton, and M. H. Amin, Observation of topological phenomena in a programmable lattice of 1,800 qubits, Nature 560, 456 (2018).
King et al. [2021a] A. D. King, J. Raymond, T. Lanting, S. V. Isakov, M. Mohseni, G. Poulin-Lamarre, S. Ejtemaee, W. Bernoudy, I. Ozfidan, A. Y. Smirnov, M. Reis, F. Altomare, M. Babcock, C. Baron, A. J. Berkley, K. Boothby, P. I. Bunyk, H. Christiani, C. Enderud, B. Evert, R. Harris, E. Hoskinson, S. Huang, K. Jooya, A. Khodabandelou, N. Ladizinsky, R. Li, P. A. Lott, A. J. R. MacDonald, D. Marsden, G. Marsden, T. Medina, R. Molavi, R. Neufeld, M. Norouzpour, T. Oh, I. Pavlov, I. Perminov, T. Prescott, C. Rich, Y. Sato, B. Sheldan, G. Sterling, L. J. Swenson, N. Tsai, M. H. Volkmann, J. D. Whittaker, W. Wilkinson, J. Yao, H. Neven, J. P. Hilton, E. Ladizinsky, M. W. Johnson, and M. H. Amin, Scaling advantage over path-integral monte carlo in quantum simulation of geometrically frustrated magnets, Nature Communications 12, 1113 (2021a).
Mishra et al. [2018] A. Mishra, T. Albash, and D. A. Lidar, Finite temperature quantum annealing solving exponentially small gap problem with non-monotonic success probability, Nature Communications 9, 2917 (2018).
Gardas et al. [2018] B. Gardas, J. Dziarmaga, W. H. Zurek, and M. Zwolak, Defects in Quantum Computers, Sci. Rep. 8, 4539 (2018).
Weinberg et al. [2020] P. Weinberg, M. Tylutki, J. M. Rönkkö, J. Westerholm, J. A. Åström, P. Manninen, P. Törmä, and A. W. Sandvik, Scaling and Diabatic Effects in Quantum Annealing with a D-Wave Device, Phys. Rev. Lett. 124, 090502 (2020).
Bando et al. [2020] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, and H. Nishimori, Probing the universality of topological defect formation in a quantum annealer: Kibble-zurek mechanism and beyond, Phys. Rev. Research 2, 033369 (2020).
Nishimura et al. [2020] K. Nishimura, H. Nishimori, and H. G. Katzgraber, Griffiths-McCoy singularity on the diluted Chimera graph: Monte Carlo simulations and experiments on quantum hardware, Phys. Rev. A 102, 042403 (2020).
King et al. [2021b] A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, and A. Lopez-Bezanilla, Qubit Spin Ice, Science 373, 576 (2021b).
Kairys et al. [2020] P. Kairys, A. D. King, I. Ozfidan, K. Boothby, J. Raymond, A. Banerjee, and T. S. Humble, Simulating the shastry-sutherland ising model using quantum annealing, PRX Quantum 1, 020320 (2020).
Abel and Spannowsky [2021] S. Abel and M. Spannowsky, Quantum-field-theoretic simulation platform for observing the fate of the false vacuum, PRX Quantum 2, 010349 (2021).
Zhou et al. [2021] S. Zhou, D. Green, E. D. Dahl, and C. Chamon, Experimental realization of classical ${\mathbb{z}}_{2}$ spin liquids in a programmable quantum device, Phys. Rev. B 104, L081107 (2021).
Derrida [1981] B. Derrida, Random-energy model: An exactly solvable model of disordered systems, Phys. Rev. B 24, 2613 (1981).
Gross and Mezard [1984] D. J. Gross and M. Mezard, The simplest spin glass, Nucl. Phys. B 240, 431 (1984).
Bapst and Semerjian [2012] V. Bapst and G. Semerjian, On quantum mean-field models and their quantum annealing, J. Stat. Mech.: Theory Exp. 2012 (6), P06007.
Perdomo-Ortiz et al. [2011] A. Perdomo-Ortiz, S. E. Venegas-Andraca, and A. Aspuru-Guzik, A study of heuristic guesses for adiabatic quantum computation, Quantum Inf. Process. 10, 33 (2011).
Yamashiro et al. [2019] Y. Yamashiro, M. Ohkuwa, H. Nishimori, and D. A. Lidar, Dynamics of reverse annealing for the fully connected $p$ -spin model, Physical Review A 100, 052321 (2019).
Chancellor [2017] N. Chancellor, Modernizing quantum annealing using local searches, New J. Phys. 19, 023024 (2017), 1606.06833 .
Ohkuwa et al. [2018] M. Ohkuwa, H. Nishimori, and D. A. Lidar, Reverse annealing for the fully connected $p$ -spin model, Physical Review A 98, 022314 (2018).
Ottaviani and Amendola [2018] D. Ottaviani and A. Amendola, Low Rank Non-Negative Matrix Factorization with D-Wave 2000Q, arXiv:1808.08721 (2018), arXiv:1808.08721 .
Venturelli and Kondratyev [2019] D. Venturelli and A. Kondratyev, Reverse quantum annealing approach to portfolio optimization problems, Quantum Mach. Intell. 1, 17 (2019).
Marshall et al. [2019] J. Marshall, D. Venturelli, I. Hen, and E. G. Rieffel, Power of Pausing: Advancing Understanding of Thermalization in Experimental Quantum Annealers, Phys. Rev. App. 11, 044083 (2019).
Pelofske et al. [2020] E. Pelofske, G. Hahn, and H. N. Djidjev, Advanced anneal paths for improved quantum annealing, in 2020 IEEE International Conference on Quantum Computing and Engineering (QCE) (IEEE Computer Society, Los Alamitos, CA, USA, 2020) pp. 256–266.
Rocutto et al. [2021] L. Rocutto, C. Destri, and E. Prati, Quantum semantic learning by reverse annealing of an adiabatic quantum computer, Adv. Quantum Technol. 4, 2000133 (2021).
Boixo et al. [2013] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Experimental signature of programmable quantum annealing, Nat. Commun. 4, 2067 (2013).
Boixo et al. [2014] S. Boixo, T. F. Ronnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nat. Phys. 10, 218 (2014).
Smolin and Smith [2014] J. A. Smolin and G. Smith, Classical signature of quantum annealing, Frontiers in Physics 2, 52 (2014).
Boixo et al. [2016] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Computational multiqubit tunnelling in programmable quantum annealers, Nat Commun 7 (2016).
Albash et al. [2015a] T. Albash, T. F. Rønnow, M. Troyer, and D. A. Lidar, Reexamining classical and quantum models for the D-Wave One processor, Eur. Phys. J. Spec. Top. 224, 111 (2015a).
Shin et al. [2014] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani, How “quantum” is the D-Wave machine?, arXiv:1401.7087 (2014).
Albash et al. [2015b] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, Consistency tests of classical and quantum models for a quantum annealer, Phys. Rev. A 91, 042314 (2015b).
Denchev et al. [2016] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, What is the computational value of finite-range tunneling?, Phys. Rev. X 6, 031015 (2016).
Albash et al. [2012] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, Quantum adiabatic Markovian master equations, New J. Phys. 14, 123016 (2012).
Xu and Cao [2016] D. Xu and J. Cao, Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: A polaron transformation approach, Frontiers of Physics 11, 110308 (2016).
Chen and Lidar [2020a] H. Chen and D. A. Lidar, HOQST: Hamiltonian Open Quantum System Toolkit, arXiv:2011.14046 (2020a).
Passarelli et al. [2020] G. Passarelli, K.-W. Yip, D. A. Lidar, H. Nishimori, and P. Lucignano, Reverse quantum annealing of the $p$ -spin model with relaxation, Phys. Rev. A 101, 022331 (2020).
Babbush et al. [2013] R. Babbush, B. O’Gorman, and A. Aspuru-Guzik, Resource efficient gadgets for compiling adiabatic quantum optimization problems, Ann. Physik 525, 877 (2013).
[47] D-Wave System Documentation – QPU-Specific Characteristics.
Choi [2008] V. Choi, Minor-embedding in adiabatic quantum computation: I. the parameter setting problem, Quantum Inf. Process. 7, 193 (2008).
Albash et al. [2016] T. Albash, W. Vinci, and D. A. Lidar, Simulated-quantum-annealing comparison between all-to-all connectivity schemes, Phys. Rev. A 94, 022327 (2016).
Albash and Marshall [2021] T. Albash and J. Marshall, Comparing relaxation mechanisms in quantum and classical transverse-field annealing, Phys. Rev. Appl. 15, 014029 (2021).
Amin [2015] M. H. Amin, Searching for quantum speedup in quasistatic quantum annealers, Phys. Rev. A 92, 052323 (2015).
[52] D-Wave System Documentation – Spin-Bath Polarization Effect.
[53] T. Lanting, M. H. Amin, C. Baron, M. Babcock, J. Boschee, S. Boixo, V. N. Smelyanskiy, M. Foygel, and A. G. Petukhov, Probing environmental spin polarization with superconducting flux qubits, arXiv:2003.14244 .
Chen and Lidar [2020b] H. Chen and D. A. Lidar, Why and when pausing is beneficial in quantum annealing, Phys. Rev. Appl. 14, 014100 (2020b).
Albash and Lidar [2015] T. Albash and D. A. Lidar, Decoherence in adiabatic quantum computation, Phys. Rev. A 91, 062320 (2015).
Yip et al. [2018] K. W. Yip, T. Albash, and D. A. Lidar, Quantum trajectories for time-dependent adiabatic master equations, Phys. Rev. A 97, 022116 (2018).
Lidar and Whaley [2003] D. A. Lidar and K. B. Whaley, in Irreversible Quantum Dynamics, Lecture Notes in Physics, Vol. 622 (Springer, Berlin, 2003) p. 83.
Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, J. Phys. Soc. Jpn. 12, 570 (1957).
Martin and Schwinger [1959] P. C. Martin and J. Schwinger, Theory of many-particle systems. I, Phys. Rev. 115, 1342 (1959).
Albash et al. [2015c] T. Albash, I. Hen, F. M. Spedalieri, and D. A. Lidar, Reexamination of the evidence for entanglement in a quantum annealer, Phys. Rev. A 92, 062328 (2015c).
Smirnov and Amin [2018] A. Y. Smirnov and M. H. Amin, Theory of open quantum dynamics with hybrid noise, New J. Phys. 20, 103037 (2018).
Dickson et al. [2013] N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris, F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple, P. Chavez, F. Cioata, T. Cirip, P. deBuen, M. Drew-Brook, C. Enderud, S. Gildert, F. Hamze, J. P. Hilton, E. Hoskinson, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Lanting, T. Mahon, R. Neufeld, T. Oh, I. Perminov, C. Petroff, A. Przybysz, C. Rich, P. Spear, A. Tcaciuc, M. C. Thom, E. Tolkacheva, S. Uchaikin, J. Wang, A. B. Wilson, Z. Merali, and G. Rose, Thermally assisted quantum annealing of a 16-qubit problem, Nat. Commun. 4, 1903 (2013).
Amin and Averin [2008] M. H. S. Amin and D. V. Averin, Macroscopic Resonant Tunneling in the Presence of Low Frequency Noise, Phys. Rev. Lett. 100, 197001 (2008).
Dekker [1987] H. Dekker, Noninteracting-blip approximation for a two-level system coupled to a heat bath, Phys. Rev. A 35, 1436 (1987).
Grifoni and Hänggi [1998] M. Grifoni and P. Hänggi, Driven quantum tunneling, Phys. Rep. 304, 229 (1998).
Rieger and Kawashima [1999] H. Rieger and N. Kawashima, Application of a continuous time cluster algorithm to the two-dimensional random quantum Ising ferromagnet, Eur. Phys. J. B 9, 233 (1999).
Martoňák et al. [2002] R. Martoňák, G. E. Santoro, and E. Tosatti, Quantum annealing by the path-integral Monte Carlo method: The two-dimensional random Ising model, Phys. Rev. B 66, 094203 (2002).
Bapst and Semerjian [2013] V. Bapst and G. Semerjian, Thermal, quantum and simulated quantum annealing: analytical comparisons for simple models, J. Phys.: Conf. Ser. 473, 012011 (2013).
Rønnow et al. [2014] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup, Science 345, 420 (2014).
Bando and Nishimori [2021] Y. Bando and H. Nishimori, Simulated quantum annealing as a simulator of nonequilibrium quantum dynamics, Phys. Rev. A 104, 022607 (2021).
Santoro and Tosatti [2006] G. E. Santoro and E. Tosatti, Optimization using quantum mechanics: quantum annealing through adiabatic evolution, J. Phys. A 39, R393 (2006).
Subires et al. [2021] D. Subires, F. J. Gómez-Ruiz, A. Ruiz-García, D. Alonso, and A. del Campo, Benchmarking quantum annealing dynamics: the spin-vector langevin model (2021), arXiv:2109.09750 .
Albash et al. [2013] T. Albash, D. A. Lidar, M. Marvian, and P. Zanardi, Fluctuation theorems for quantum processes, Phys. Rev. E 88, 032146 (2013).
Lanting et al. [2014] T. Lanting, A. J. Przybysz, A. Y. Smirnov, F. M. Spedalieri, M. H. Amin, A. J. Berkley, R. Harris, F. Altomare, S. Boixo, P. Bunyk, N. Dickson, C. Enderud, J. P. Hilton, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, R. Neufeld, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, S. Uchaikin, A. B. Wilson, and G. Rose, Entanglement in a quantum annealing processor, Phys. Rev. X 4, 021041 (2014).
King et al. [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Coherent quantum annealing in a programmable 2000-qubit Ising chain, arXiv:2202.05847 (2022).
Venturelli et al. [2015] D. Venturelli, S. Mandrà, S. Knysh, B. O’Gorman, R. Biswas, and V. Smelyanskiy, Quantum optimization of fully connected spin glasses, Phys. Rev. X 5, 031040 (2015).
Job and Lidar [2018] J. Job and D. Lidar, Test-driving 1000 qubits, Quantum Sci. Tech. 3, 030501 (2018).
M.E. Rose [1995] M.E. Rose, Elementary Theory of Angular Momentum (Dover, New York, 1995).
Lidar [2019] D. A. Lidar, Lecture notes on the theory of open quantum systems, arXiv preprint arXiv:1902.00967 (2019).
Muthukrishnan et al. [2016] S. Muthukrishnan, T. Albash, and D. A. Lidar, Tunneling and Speedup in Quantum Optimization for Permutation-Symmetric Problems, Phys. Rev. X 6, 031010 (2016).

$\displaystyle p(r=1)$	$\displaystyle=\|\braket{\text{up}}{\psi}_{\text{fin}}\|^{2}+\|\braket{\text{down}}{\psi}_{\text{fin}}\|^{2}$
	$\displaystyle=\|a_{S_{\max}}\|^{2}(\|c_{S_{\max}}\|^{2}+\|c_{-S_{\max}}\|^{2})$
	$\displaystyle\leq\|a_{S_{\max}}\|^{2}\,,$	(42)

Breakdown of the weak-coupling limit in quantum annealing

Abstract

I Introduction

II Problem definition and reverse annealing protocol

II.1 The pp-spin model with p=2p=2

II.2 Reverse annealing protocol

III Empirical results

III.1 Dependence on initial conditions

III.1.1 m0=0m_{0}=0

III.1.2 m0=0.8,0.9m_{0}=0.8,0.9: up/down symmetry breaking for 0.3≲sinv≲0.50.3\lesssim s_{\mathrm{inv}}\lesssim 0.5

III.1.3 sinv≳0.5>sΔs_{\mathrm{inv}}\gtrsim 0.5>s_{\Delta}: freezing

III.1.4 sinv≲0.3<sΔs_{\mathrm{inv}}\lesssim 0.3<s_{\Delta}: spin-bath polarization

III.1.5 Total success probability

III.2 Dependence on annealing time

III.3 Effects of pausing

III.4 Size dependence

III.5 Effects of iteration

IV Numerical Results

IV.1 Closed system model

IV.1.1 Dependence on system size and annealing time

IV.1.2 Dependence on initial state and number of iterations

IV.2 Open system model: setup

IV.3 Weak coupling: adiabatic master equation

IV.3.1 Dependence on the annealing time.

IV.3.2 Dependence on initial condition and the number of cycles.

IV.3.3 Size dependence and partial success probability.

IV.4 Strong coupling: PTRE

IV.4.1 4-qubit case.

IV.4.2 8-qubit case.

V Conclusions

Acknowledgements.

Appendix A Optimal Value of Ferromagnetic Interactions within a Logical Qubit

Appendix B Reverse annealing details

Appendix C Spectrum of the pp-spin problem with p=2p=2 and scaling of the minimum gap with NN

Appendix D Upper bound on the success probabilities for collective dephasing

Appendix E The PTRE equation involves only the diagonal components of the density matrix

Appendix F Fully classical simulations using the spin-vector Monte Carlo algorithm

F.1 Total success probability

F.2 Partial success probability

Appendix G Excited states

G.1 N=4N=4 with a maximally excited initial state

G.2 N=8N=8 with a maximally excited initial state

G.3 First excited state as initial state

Appendix H SVMC-TF’s dependence on the number of sweeps

Appendix I Pseudocode for SVMC and SVMC-TF

References

II.1 The $p$ -spin model with $p=2$

III.1.1 $m_{0}=0$

III.1.2 $m_{0}=0.8,0.9$ : up/down symmetry breaking for $0.3\lesssim s_{\mathrm{inv}}\lesssim 0.5$

III.1.3 $s_{\mathrm{inv}}\gtrsim 0.5>s_{\Delta}$ : freezing

III.1.4 $s_{\mathrm{inv}}\lesssim 0.3<s_{\Delta}$ : spin-bath polarization

Appendix C Spectrum of the $p$ -spin problem with $p=2$ and scaling of the minimum gap with $N$

G.1 $N=4$ with a maximally excited initial state

G.2 $N=8$ with a maximally excited initial state