††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].

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] 111 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.

Figure 2: Empirical success probabilities for different initial conditions
$m_{0}$ (magnetization) in reverse annealing on the D-Wave 2000Q device as a
function of the inversion point $s_{\mathrm{inv}}$. The dashed line is the
minimum-gap point at $s_{\Delta}\approx 0.36$ for $N=20$. When
$s_{\mathrm{inv}}<s_{\Delta}$ the reverse direction of the anneal goes through
and past the minimum gap point, and then crosses it again during the forward
anneal. When $s_{\mathrm{inv}}>s_{\Delta}$ there is no crossing of the minimum
gap point. (a) Partial success probabilities for $m_{0}=0,0.8$, and $0.9$. (b)
Total success probabilities for the same set of $m_{0}$ as in (a). Here and in
all subsequent figures the labels “up” and “down” mean the populations of the
all-up state and the all-down state, respectively; $s_{\mathrm{inv}}$ was
incremented by steps of $0.02$, and error bars denote one standard deviation;
see App. B for details on the calculation of error bars.

###  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.

Figure 3: Empirical success probabilities for different annealing times $\tau$
in reverse annealing on the D-Wave 2000Q device as a function of
$s_{\mathrm{inv}}$. (a) Partial success probabilities for $\tau=5,10$, and
$20\mu s$. (b) Total success probabilities for the same set of $\tau$ as in
(a). The dashed line is the minimum-gap point at $N=20$, $s_{\Delta}\approx
0.36$. Figure 4: Empirical success probabilities for different pausing times
$t_{p}$ in reverse annealing as a function of $s_{\mathrm{inv}}$. (a) Partial
success probabilities for pause $t_{p}=0,2$, and $10~{}\mu s$. (b) Total
success probabilities for the same set of $t_{p}$ as in (a). The dashed line
is the minimum-gap point at $N=20$, $s_{\Delta}\approx 0.36$. Figure 5:
Empirical success probabilities for different system sizes $N$ in reverse
annealing as a function of $s_{\mathrm{inv}}$. The initial state for each $N$
is a state with one spin flipped. (a) Partial success probabilities for system
size $N\in\\{4,8,12,16,20,24\\}$. (b) Total success probabilities for the same
set of $N$ as in (a).  Figure 6: The position $s_{\Delta}$ of the minimum gap
$\Delta$ as a function of system size $N$, along with the value $s_{0.5}$ of
$s_{\mathrm{inv}}$ at which the total success probability $=0.5$ for each $N$.
The jump at $N=16$ is due to the non-monotonicity seen for the up
probabilities in Fig. 5(a). Inset: the minimum gap for each $N$. The gap obeys
a power law scaling $N^{-11/30}$.

####  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.222 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.

Figure 7: Empirical success probabilities for different number of iterations
$r$ in iterated reverse annealing as a function of $s_{\mathrm{inv}}$. Panels
(a) and (b) are the partial and total success probabilities with the initial
state $m_{0}=0.9$ (the first excited state), respectively. Panels (c) and (d)
are the partial and total success probabilities with the initial state
$m_{0}=0$ (the highest excited state). The dashed line is the minimum-gap
point $s_{\Delta}\approx 0.36$ for $N=20$.

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.333$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.

Figure 8: Total success probabilities as computed by solution of the
Schrödinger equation for a closed system with $r=1$ (single cycle) and a
single-spin down as the initial state. (a) $N=4$, (b) $N=8$. Here and in the
subsequent figures the nanoscecond timescale is set by the energy scale of the
D-Wave annealing schedule shown in Fig. 1. Note the different scale of the
vertical axes.

####  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: Total success probabilities as computed by solution of the
Schrödinger equation for a closed system with two spins down as the initial
state. (a) $N=4$ and $r=1$, (b) $N=8$ and $r=1$, (c) $N=4$ and $r=2$, (d)
$N=8$ and $r=2$. Note the different scale of the vertical axes.

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)\,,$

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.

Figure 10: Top row: Total success probabilities as computed by the AME as a
function of $s_{\mathrm{inv}}$ with (a) independent and (b) collective
dephasing (SB denotes system-bath coupling). Initial state:
$\ket{0001}~{}(m_{0}=0.5)$ and $r=1$, a single cycle. Bottom row: Total
success probabilities as computed by the AME with different $\tau$. Initial
state: $\ket{0011}~{}(m_{0}=0)$. (c) $r=1$, (d) $r=2$. Note the different
scale of the vertical axis of panel (b).

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).

Figure 11: Total and partial success probabilities as computed by the AME as a
function of $s_{\mathrm{inv}}$ for independent dephasing. The initial state
has one spin down and $\tau=1\;\mu$s. The main plot shows the total success
probabilities for $N=4$ and $8$. The inset shows the partial success
probabilities for $N=8$.

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).

Figure 12: PTRE simulation results versus the empirical results for the
4-qubit case. The initial states are (a)
$\ket{\uparrow\uparrow\uparrow\downarrow}$ and (b)
$\ket{\uparrow\uparrow\downarrow\downarrow}$. The best-fit simulation
parameters are: $\eta g^{2}=2.5\times 10^{-3}$, $T=25$ mK, $\tau=5\mu$ s,
$W=8$mK, and $f_{c}=1$THz. Here and in the subsequent figures the error bars
represent $1\sigma$ confidence intervals. Overall the PTRE results are in good
qualitative agreement with the empirical results for all cases shown. Note the
different scale of the vertical axes.

####  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$.

Figure 13: PTRE simulation results versus the experimental results for the
8-qubit case. The initial states are (a) the first excited state
$\ket{\uparrow\cdots\uparrow\downarrow}$ and (b) the maximally excited state
$\ket{\uparrow\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\downarrow}$.
The simulation parameters are identical to the optimal ones obtained from the
parameter estimation procedure in Fig. 12. Note the different scale of the
vertical axes.

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$.

Figure 14:  Success probability of forward annealing with annealing time
$\tau=1.0~{}\mu s$ as a function of the coupling between physical qubits in a
logical qubit for different system sizes. The solid lines are Bezier curves.
The penalty value we used is $-J_{\text{F}}$, i.e., ferromagnetic.

##  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)

Figure 15: Full spectrum of four qubits for the $p=2$ $p$-spin model, subject
to the annealing schedules shown in Fig. 1.

##  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$.

Figure 16: Success probability as a function of $s_{\mathrm{inv}}$ for
independent and collective dephasing, for $N=8$ spins. The initial state has a
single spin flipped and $\tau=1\;\mu$s.

##  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.

Figure 17: (a) and (b): The population differences between the experimental
and PTRE simulation results shown in Fig. 12 for $N=4$ qubits. (c) and (d):
The population differences between the experimental and simulation results
shown in Fig. 13 for $N=8$ qubits. Note the different scale of the vertical
axes.

##  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.

Figure 18: The semiclassical potential landscape corresponding to the
Hamiltonian of Eq. (50) for $N=1$. Ground states correspond to $\theta=0,\pi$.
The two saddle points are at $s=0.3948$ and $\theta=\pi/2\pm 0.31518\pi$.

###  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}$.

Figure 19: Total success probabilities as computed by SVMC and SVMC-TF. (a)
$\tau=10^{3}$ sweeps, (b) $\tau=10^{4}$ sweeps. Error bars are 2$\sigma$ over
$10^{4}$ samples.

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.

Figure 20: Partial success probabilities for the one-down initial state as
computed by (a) SVMC for $N=4$, (b) SVMC-TF for $N=4$. Panels (c) and (d) show
the partial ground state (GS) success probabilities for SVMC-TF (“Sim.”) at
optimized sweeps numbers of (c) $\tau=450$ for $N=4$ and (d) $\tau=150$ for
$N=8$, which yields a qualitatively good agreement with the experimental data
(“Exp.”). In panels (a) and (b) $\tau=10^{3}$ sweeps and error bars are
$2\sigma$ over $10^{4}$ samples. In panels (c) and (d) error bars are
$1\sigma$ over $10^{4}$ samples.

Figure 21: Partial success probabilities as computed by SVMC-TF. (a) $N=16$,
$\tau=10^{3}$ sweeps, (b) $N=16$, $\tau=10^{4}$ sweeps, (c) $N=32$,
$\tau=10^{3}$ sweeps, (d) $N=32$, $\tau=10^{4}$ sweeps. Error bars are
2$\sigma$ over $10^{4}$ samples.

Figure 22: Population in the first excited state as a function of
$s_{\text{inv}}$ for D-Wave (Exp) vs simulations (Sim): (a) the AME, (b) the
PTRE, and (c) SVMC-TF with $\tau=450$ sweeps. The initial state is
$\ket{\uparrow\uparrow\downarrow\downarrow}$. Error bars denote $1\sigma$ over
$10^{4}$ samples. The slight difference in SVMC-TF peak heights in panel (c)
is a numerical artifact. Note the different scale of the vertical axis of
panel (c).

Figure 23: Population in the ground and second excited states as a function of
$s_{\text{inv}}$ for D-Wave (Exp) vs simulations (Sim): (a) the AME and (b)
SVMC-TF with $\tau=450$ sweeps, for the initial state
$\ket{\uparrow\uparrow\downarrow\downarrow}$. (c) Ground and fourth excited
states for D-Wave vs SVMC-TF with $\tau=150$ sweeps, for the initial state
$\ket{\uparrow\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\downarrow}$.
Error bars denote $1\sigma$ over $10^{4}$ samples.

Figure 24: Population in the first, second, and third excited states vs
$s_{\text{inv}}$ for D-Wave vs the PTRE [(a)-(c)] and SVMC-TF with $\tau=150$
sweeps [(d)-(f)]. The initial state is
$\ket{\uparrow\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\downarrow}$.
Error bars denote $1\sigma$ over $10^{4}$ samples. For the PTRE the simulation
parameters are identical to the optimal ones obtained from the parameter
estimation procedure in Fig. 12. Note the different scale of the vertical
axes.

###  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.444We 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$.

Figure 25: Comparison of experimental data (Exp) and simulation data (Sim) for
the initial state $\ket{\uparrow\uparrow\uparrow\downarrow}$. (a) PTRE, (b)
SVMC-TF with $450$ sweeps. (c) PTRE, (d) SVMC-TF with $150$ sweeps, for the
initial state
$\ket{\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}$ .
The population shown is for the state with all but one spin up (“1st excited
up”) and all but one spin down (“1st excited down”). Error bars denote
$1\sigma$ over $10^{4}$ samples.

##  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$.

Figure 26: Sum of $\ell_{2}$-norms (between simulation and experiment data
series) over all the eigenstates vs the number of sweeps used in SVMC-TF
simulations. (a) Subplot legend shows the initial states. (b) For $N=4$, SVMC-
TF produces simulation data best matched with the empirical date at $450$
sweeps, while for $N=8$ the optimal number is $150$.

##  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$ to $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$ to $T$ (total number of sweeps) do

for $i=1$ to $n$ (number of qubits) do

Randomly choose a new angle $\theta_{i,t}^{{}^{\prime}k}\in[0,\pi]$.

Calculate $\Delta E$.

if $\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

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$ to $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$ to $T$ (total number of sweeps) do

for $i=1$ to $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$.

if $\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

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

if $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).