# Embedding quantum optimization problems using AC driven quantum ferromagnets

Gianni Mossi1,2, Vadim Oganesyan3,4, Eliot Kapit5 1 KBR, Inc., 601 Jefferson
St., Houston, TX 77002, USA.  2 Quantum Artificial Intelligence Lab. (QuAIL),
NASA Ames Research Center, Moffett Field, CA 94035, USA.  3 Department of
Physics and Astronomy, College of Staten Island, CUNY, Staten Island, NY
10314, USA  4 Physics program and Initiative for the Theoretical Sciences, The
Graduate Center, CUNY, New York, NY 10016, USA  5 Department of Physics and
Quantum Engineering Program, Colorado School of Mines, Golden, CO 80401

###### Abstract

Analog quantum optimization methods, such as quantum annealing, are promising
and at least partially noise tolerant ways to solve hard optimization and
sampling problems with quantum hardware. However, they have thus far failed to
demonstrate broadly applicable quantum speedups, and an important contributing
factor to this is slowdowns from embedding, the process of mapping logical
variables to long chains of physical qubits, enabling arbitrary connectivity
on the short-ranged 2d hardware grid. Beyond the spatial overhead in qubit
count, embedding can lead to severe time overhead, arising from processes
where individual chains “freeze” into ferromagnetic states at different times
during evolution, and once frozen the tunneling rate of this single logical
variable decays exponentially in chain length. We show that this effect can be
substantially mitigated by local AC variation of the qubit parameters as in
the RFQA protocol (Kapit and Oganesyan, Quant. Sci. Tech. 6, 025013 (2021)),
through a mechanism we call Symphonic Tunneling. We provide general arguments
and substantial numerical evidence to show that AC-driven multi-qubit
tunneling is dramatically faster than its DC counterpart, and since ST is not
a 1d-specific mechanism, this enhancement should extend to clusters of coupled
chains as well. And unlike a uniform transverse field, in higher dimensions
this method cannot be efficiently simulated classically. We explore schemes to
synchronize the AC tones within chains to further improve performance.
Implemented at scale, these methods could significantly improve the prospects
for achieving quantum scaling advantages in near-term hardware.

##  1 Introduction

While enormous progress has been made in gate model quantum computing in
recent years, quantum advantage for practical, real-world problems with such
devices has yet to be demonstrated. This is largely due to the extreme
precision requirements and noise sensitivity inherent to the gate model. In
contrast, _analog_ quantum optimizers, including quantum annealers (QA) [1, 2,
3, 4, 5, 6, 7, 8, 9, 10] and more recent Rydberg atom systems [11, 12, 13],
are much more resilient can routinely solve problems with hundreds or even
thousands of variables. While not computationally universal, a huge array of
NP-hard and NP-complete problems can be solved using these systems, with the
likely mechanism of quantum advantage lying in collective quantum tunneling as
a way to escape local minima once the system becomes glassy [14, 15, 16, 17].
Recent experimental results [10, 13, 18] on large problems in both neutral
atom and flux qubit analog devices suggest quantum scaling advantage over
classical heuristics, for problems which are short ranged and map natively to
the underlying hardware graph. However, demonstration of _broadly applicable_
quantum scaling advantage has remained elusive, particularly for problems
defined on abstract graphs with long ranged connectivity.

To attack these problems, encodings such as minor embedding are required [19,
20, 21], where single logical variables are mapped to long ferromagnetically
coupled chains (similar, antiferromagnetic schemes exist for Rydberg arrays
[22, 23]). These methods have an inescapable quadratic overhead in qubit
count, but more insidiously, they can introduce severe time penalties as well
[21, 24, 25, 26]. The essential mechanism for this in large problems is that
individual chains will “freeze” into ferromagnetic states at different times
during evolution (note of course that the whole system will eventually freeze
as the transverse field driving quantum dynamics is reduced to zero). In the
process, transverse field corrections and other effects can tip them into the
“wrong” state, and once frozen the tunneling rate of these single variables
decays exponentially in chain length. This leads to an ugly tradeoff in the
strength of the intra-chain coupling $J_{c}$. If it is too weak, the system
becomes more noise sensitive and can have (logically meaningless) broken
chains in its ground state. Larger $J_{c}$ values however, tend to freeze
earlier and have even worse scaling for multiqubit tunneling. In the worst
cases, the time to solution can scale exponentially in the total system size,
e.g. $\exp\left(N^{2}\right)$ for $N$ logical variables [25]. And unlike the
spatial overhead, this problem cannot be solved by simply making the system
larger; a more fundamental change to the quantum optimization protocol is
required.

We have demonstrated previously the use of randomised dynamical protocols to
ameliorate exponential scaling, both analytically and numerically in simple
0-dimensional problems of $N$ coupled qubits, where level structure exhibits
an avoided level crossing with a gap to other excited states. The Random Field
(or Radio Frequency) Quantum Annealing (RFQA) protocol[27] utilizes this
separation of scales to accelerate the mixing across the level crossing with
no appreciable heating, thus reducing the “difficulty exponent”

$\Upsilon\equiv\frac{1}{N}\log_{2}\Big{(}\frac{1}{\Gamma(N)}\Big{)},$ (1)

as extracted from the tunnelling rate $\Gamma(N)$ (more on this in the next
Section 2), which is constrained in those structureless models by the optimal
adiabatic value, e.g. $\Upsilon\geq\frac{1}{2}$ for the Grover problem [28].
While this optimal value is exponentially fragile to control precision, RFQA
is not. By contrast, quantum tunnelling of logical qubits based on the two
groundstates of the Ising chain of length $L$ is governed by a highly
structured Hamiltonian and we therefore expect to be able to accelerate the
dynamics more efficiently. To this end we explored a number of _structured_
protocols inspired both by basic understanding of correlated dynamics of small
spin clusters and practical considerations of hardware implementation.
Remarkably, we find significant reductions $\Upsilon$, seemingly to
arbitrarily small values (without appreciable heating or melting of
ferromagnetic order – see Section 4), as we vary drive strength and other
detials, see Fig. 1. We refer to the general scheme of structured/optimized
multifrequency protocols as Symphonic Tunnelling.

Figure 1: Difficulty exponent $\Upsilon$ obtained from the rate of complete
magnetization reversal (see Eq. 1 for definition and Sec. 4 for further
details) as a function of oscillations’ amplitude $\alpha$ for three different
driving protocols defined in Sec. 3, for three different values of transverse
field $\kappa$. A vanishing or negative $\Upsilon$ implies that an exponential
fit is no longer good (at least for the ranges simulated), indicating a
potential crossover to polynomial scaling.

To demonstrate the potential of RFQA/ST to mitigate the time overhead of
embedding, we report the results of an extensive series of numerical
simulations, focusing on the tunneling of a single one-dimensional ring in the
ferromagnetic (frozen) phase, where the base DC tunneling rate decays
exponentially in length $L$. We focus on single chains because, relying on AC
dynamics, RFQA cannot be simulated with quantum Monte Carlo [29, 30, 31, 32,
33], and thus embedded _problems_ are firmly out of reach of classical
simulation techniques at any meaningful scales. Likewise, while we expect RFQA
methods to accelerate problem solving for native problems, benchmarking that
in simulation is problematic because small instances (e.g. $N\leq 20$)
typically do not show the expected large-scale exponential difficulty scaling
and instances numerically mined for very small gaps are generally fragile to
perturbations [34]. The exponential scaling of collective tunneling in 1d
chains, in contrast, is well-controlled and obvious at small $L$.

Further, the key slowdown mechanism of freezing can be captured in a single
chain, and we average over random detuning as a proxy for the disordered
environment–and inability to find precise resonances–of real analog problems.
And since ST is fundamentally _not_ a 1d-specific mechanism, it will
generalize to clusters of coupled chains in larger systems. Finally, 1d chains
are amenable to matrix product simulation techniques (in this case, time-
evolving block decimation [35, 36, 37, 38, 39]), which allow us to extend our
simulations out to much larger system sizes than are possible in full-
wavefunction evolutions. We demonstrate parameter regimes where the chain
remains frozen–in that ferromagnetic correlations do not vanish at large
separation–but collective tunneling is accelerated to the point that
polynomial and exponential scaling cannot be distinguished from the data.

The rest of this paper is organized as follows. In the next section, we
introduce the basic 1d collective tunneling problem and the AC protocols we
will use to attack it. Following that, we present details of our numerical
simulation methods, including full wavefunction simulations and matrix product
methods. We then present numerical results for a range of system parameters
and protocols. While we do not have analytical arguments to predict
performance in arbitrarily large systems, the robustness of our observed
speedups (which persist out to the largest system sizes we can simulate)
strongly suggest significant performance improvements at the application
scale. We finally detail some of the considerations for a superconducting
hardware implementation, and offer concluding remarks.

##  2 Ising Model of reverse annealing – quasi-static DC protocol

We consider quantum dynamics of a single ferromagnetic ring, with total
Hamiltonian

$\displaystyle
H\left(t\right)=-\sum_{j}\left[\kappa_{j}\left(t\right)X_{j}+h_{j}\left(t\right)Z_{j}+J_{j}\left(t\right)Z_{j}Z_{j+1}\right],$
(2)

with time dependent couplings chosen to accelerate the mixing of two
ferromagnetic ground states without any appreciable excitation of excited
states. For ease of comparison we will benchmark all protocols using the
standard ”reverse annealing” quench, whereby we initialize the system in a
definite all-down classical groundstate with no transverse field, turn the
transverse field on and off slowly, and examine the probability of complete
magnetization reversal

$P(t_{f})\equiv\lvert\langle\uparrow\cdots\uparrow\rvert
U(t)\lvert\downarrow\cdots\downarrow\rangle\rvert^{2}.$ (3)

Importantly, we take great precautions to avoid exciting the problem out of
the low energy ferromagnetic doublet, i.e.

$P_{heat}\equiv 1-\lvert\langle\uparrow\cdots\uparrow\rvert
U(t)\lvert\downarrow\cdots\downarrow\rangle\rvert^{2}-\lvert\langle\downarrow\cdots\downarrow\rvert
U(t)\lvert\downarrow\cdots\downarrow\rangle\rvert^{2}\ll P(t_{f}).$ (4)

We do not attempt to measure heating during evolution, but since no cooling
mechanism is present to remove stray excitations, measuring it at the end is
sufficient.

In what we follows we focus on the statistically significant average probably
of success, whereby several instances of the problem are considered with
resultant $P(t_{f})$ averaged over the ensemble of static and dynamic
variations of potentials. We average $P\left(t_{f}\right)$ over fixed detuning
$h$ (corresponding to a uniform $Z$ bias $h/2L$) drawn from the uniform range
$\left\\{-2/\sqrt{L},2/\sqrt{L}\right\\}$ (this scaling choice is justified
below), for both the DC and RFQA cases. We calibrate the analysis of
complicated protocols defined below against the simplest uniform quasi-state
sweep characterized by the maximum value of transverse field $\kappa_{0}$ –
see Fig. 2 and, importantly, the value of detuning $W\equiv 2\sum_{i}h_{i}$.
In this simple case we can accurately estimate (see App. A of [27])

$P(t_{f})\approx\int
dW\mathcal{P}(W)\frac{\Omega_{0}^{2}}{W}\sin^{2}(Wt_{f})\to\pi\frac{\Omega_{0}^{2}}{W}t_{f}\equiv\Gamma_{0}t_{f},$
(5)

which ignores ramp time $t_{ramp}$ (in practice, we fix $t_{f}=6(2\pi
N)=6t_{ramp}$) and assumes are relatively smooth distribution of detuning
$\mathcal{P}(W)$ and sufficiently long $t_{f}$ to sample it efficiently. The
most important part of this expression is the many-body matrix element
connecting the two magnetization states which for a uniform field case is
well-captured by the scaling form

$\displaystyle\Omega_{0}\left(L\right)\propto\frac{\kappa_{0}}{L^{\kappa_{0}/J}}\left(\frac{\kappa_{0}}{J}\right)^{L-1},\;\;\Gamma_{0}\propto\frac{\Omega_{0}^{2}}{W}.$
(6)

This matrix element figures prominently in the standard textbook formulation
of symmetry restoration via path integral instantons, but can also be obtained
by direct high order perturbation theory; at large $L$ the polynomial
prefactor is largely irrelevant but included for completeness. Since
$\Omega_{0}(L)$ is normally exponentially small in $L$, the total number of
spins flipping, we can use the difficulty exponent
$\Upsilon=(1/L)\log_{2}(1/\Gamma(L))$ to help isolate the scaling advantage of
various protocols. Thus, the uniform field protocol discussed thus far
corresponds to $\Upsilon_{0}=2\log_{2}J/\kappa_{0}$, which only vanishes at
the quantum critical point, where the order parameter vanishes, and with it
the ability to use the chains for embedding logical qubits!

Before turning to the more powerful protocols that appear to deliver dramatic
reduction in $\Upsilon$ without destroying the logical qubits in the process
we close this section by defining two useful quantitative tools. One is the so
called ”time-to-solution” (TTS) that is a nice proxy[40] for $\Gamma$

$\displaystyle\mathrm{TTS}\equiv t_{f}\frac{\ln(1-0.99)}{\ln(1-\langle
P(t_{f})\rangle)}\sim\frac{t_{f}}{\langle P(t_{f})\rangle}\to 1/\Gamma.$ (7)

And finally, we also track the time-averaged order parameter correlation
function _during_ the waiting stage at largest separation available

$C_{ZZ}=\frac{1}{t_{wait}}\int_{t_{ramp}}^{t_{f}-t_{ramp}}\overline{\langle
Z_{0}(t)Z_{L/2}(t)\rangle}\,\mathrm{d}t.$ (8)

Tracking this quantity allows us to confirm that we do not leave the
ferromagnetic phase in AC evolution, as discussed below in the results
section.

Figure 2: The reverse annealing schedule explored in this work. The
ferromagnetic coupling $J$ is held fixed (red line), and the uniform
transverse field is ramped up and down to induce collective tunneling (blue
curve). In the RFQA cases, these parameters are locally modulated with
independent oscillating frequencies; the modulation of two (out of $L$)
transverse fields $\kappa_{j}\left(t\right)$ and $\kappa_{k}\left(t\right)$ is
shown in gold and green. Note in other protocols individual coupling strengths
$J_{j}$ can be oscillated, and/or these coupling strengths can be modulated by
the transverse field strengths of the coupled qubits (as in real quantum
annealers).

##  3 Dynamic protocols and other simulation details

This section details three novel modifications to the standard quasi-static
approach described above. All these approaches include modulations to the
three types of couplings already present in the Hamiltonian that are both weak
and slow so as to avoid exciting (heating) the system while dramatically
increasing the rate of mixing of the two ferromagnetic groundstates. The
seemingly inevitable success of this program was demonstrated by two of us
both analytically and numerically in the context of featuresless ”oracle”
problems[27] where the choice of driving terms was randomized and averaged
over, which each qubit assigned its own driver. By contrast, the protocols we
employ in this work are much less random, with total number of independent
frequencies scaling as $\sim\sqrt{L}$, i.e. with many qubits sharing the same
driver. This modification appears to dramatically improve the performance of
Symphonic Tunnelling, while at the same time precluding a straightforward
extension of prior analytic results for the fully random ensemble. Because of
its pivotal conceptual importance (certainly in our thinking) and clarity, we
first review the basic mechanism of the fully random protocol, RFQA, before
turning to the correlated variant employed here in the following paragraph.

Consider an exponentially avoided level crossing or a first order phase
transition where $K$ spins must flip. Now imagine to every spin we add an
oscillating term $O_{j}\sin 2\pi f_{j}t$, where $O_{j}$ is a local operator.
At first order, when the energy difference of the two states approaches any
$\pm f_{j}$, the AC perturbation can resonantly mix the states, with a mean
Rabi frequency $c\Omega_{0}$, for some $O\left(1\right)$ constant $c$. This
accelerates the off-resonant tunneling rate by a factor of $K$, but amounts to
a dramatic underestimate, since there are also $4\binom{K}{2}$ two-frequency
processes at second order, $8\binom{K}{3}$ three-frequency terms at third
order, and so on; all are generically nonzero. If we assume that the _average_
$m$th order process scales as $\Lambda^{m}\Omega_{0}$ (for $\Lambda<1$, e.g.
exponential decay in $m$), then a simple incoherent sum of all contributing
processes produces a total energy-averaged tunneling rate

$\displaystyle\Gamma_{T}\propto\frac{\Omega_{0}^{2}}{W}\sum_{m=0}^{K}2^{m}\Lambda^{2m}\binom{K}{m}\propto\frac{\Omega_{0}^{2}}{W}\left(1+2\Lambda^{2}\right)^{K}.$
(9)

This can favorably shift the scaling exponent of the problem, and of
particular interest to us is the case where $\left(1+2\Lambda^{2}\right)^{K}$
grows faster than $\Omega_{0}^{2}$ decays and this simple prediction breaks
down. For structureless problems examined previously[27] adiabatic scaling
acts as a speedlimit, hence breakdown of perturbation theory is unphysical.
However, typical problems of interest are structured and one might expect
physical consequence of such a breakdown, e.g. a breakdown of exponential
scaling of mixing in favor of a faster process. We _numerically_ demonstrate
such crossovers for multiple protocols below, where exponential and polynomial
scaling cannot be distinguished from the data, though we emphasize that these
are numerical results and we do not have any analytical guarantees this
scaling would persist as $L\to\infty$. However, the existence of this
dynamical regime for $L\approx 40$ is of keen experimental and practical
interest to ongoing efforts to engineer quantum platforms.

We now compare the standard uniform-transverse-field reverse annealing
protocol (the “DC protocol”) with three different implementations of the RFQA
protocol. These are differentiated by (i) the operators that are being
oscillated, and (ii) the distribution of $L$ random frequencies $\\{f_{j}\\}$
and initial phases $\\{\phi_{j}\\}$ that parametrize these oscillations. We
define a parameter $\alpha$, common to all protocols, that fixes the amplitude
of the oscillations. Note that this parameter determines the _relative_
magnitude of the oscillations and it multiplies the base parameter (e.g.
coupling $J_{j}$ or transverse field $\kappa_{j}$) in the full time-dependent
$H\left(t\right)$. In all cases the longitudinal bias field
$h_{j}\left(t\right)=-h$ is kept fixed throughout evolution. The three RFQA
protocols we simulate here are:

  * •

X-RFQA with randomly-paired oscillations: transverse field strengths are
locally modulated, so and $J_{j}(t)=1$ while
$\kappa_{j}(t)=\kappa(t)(1+\alpha\sin(f_{j}t+\phi_{j}))$ where $\kappa(t)$ is
the ramp schedule described in the DC protocol. $O(\sqrt{L})$ random
frequencies and phases are generated independently and these are randomly
assigned to the spins in the chain (each tone is thus repeated
$O\left(\sqrt{L}\right)$ times). This protocol could be implemented in neutral
atoms by locally modulating laser intensities.

  * •

ZZ-RFQA with randomly-paired oscillations: local ferromagnetic coupling
strengths are modulated. Here $\kappa_{j}(t)=\kappa(t)$ is the schedule
described in the DC protocol and interactions are oscillated like
$J_{j}(t)=1+\alpha\sin(f_{j}t+\phi_{j})$. We sample the random frequencies and
phases in the same way as the X-RFQA case. This protocol can be implemented in
quantum annealers through local AC flux control of the coupling terms.

  * •

XZZ-RFQA with randomly-paired oscillations: the transverse fields follow the
same schedule $\kappa_{j}(t)=\kappa(t)(1+\alpha\sin(f_{j}t+\phi_{j}))$ of the
X-RFQA protocol, while ferromagnetic interactions are modulated based on the
transverse field strengths of the coupled qubits, as
$J_{j}(t)\equiv(1-\alpha\sin(f_{j}t+\phi_{j}))(1-\alpha\sin(f_{j+1}t+\phi_{j+1}))$
This reflects real flux qubit physics [6], where the stronger the transverse
field is, the lower the susceptibility to longitudinal flux, so increasing a
local transverse field decreases couplings to that qubit. This protocol could
be implemented in flux qubit hardware through time-dependent local flux
control of transverse field strengths.

Notice that in all cases we repeat each tone $O\left(\sqrt{L}\right)$ times at
random spatial locations–$\sqrt{L}/3$, to be precise (when this quantity is
not an integer, each tone is repeated either ${\rm
floor}\left(\sqrt{L}/3\right)$ or ${\rm ceiling}\left(\sqrt{L}/3\right)$ with
appropriate probability, and the process halts when all $L$ sites have
frequency/phase pairs assigned). We choose this pairing structure for a few
reasons. First, for all three protocols, synchronizing tones causes the
corresponding AC terms to interfere constructively, increasing the boost to
collective tunneling. However, if the same tone is repeated too many times, it
can cause large mean value swings in $\kappa$ and/or $J$, leading to phase
transitions into the paramagnetic state. This shows up as more significant
heating and the vanishing of the ferromagnetic order parameter (defined
above), which are both very undesirable in the context of embedded chains for
problem solving. We found $O\left(\sqrt{L}\right)$ pairing to be a kind of
“sweet spot” between reaping the benefits of synchronization without
significant heating or phase transitions out of the ground state doublet (and
further, relevant for experimental implementation, smaller modulation
amplitudes $\alpha$ are required to obtain the same change in $\Upsilon$ with
tone repetition). Finally, in a real implementation of $N$ qubits in a 2d
lattice, a requirement of only $O\left(\sqrt{N}\right)$ unique frequencies
significantly reduces the signal generation and control complexity.

Each random frequency is drawn from the range
$f_{j}\in\left\\{1/2L,1/L\right\\}$ with random phase $\phi_{j}$. This choice
significantly mitigates heating from the AC drives. The $1/\sqrt{L}$ scaling
of the detuning range ensures that an $O\left(1\right)$ fraction of the RFQA
trials are driving tunneling between states separated by any energy difference
less than $O\left(\sqrt{L}\left<f_{j}\right>\right)$; as discussed extensively
in [27], multiqubit transitions are only accelerated by RFQA if the energy
difference between the competing states falls in the window in which the
frequency combinations are dense, consistent with the core mechanism being an
exponential proliferation of weak resonances. This slowly decaying “resonance”
condition is in stark contrast to DC protocols (the “population transfer” or
“reverse annealing” simulated here) where given an exponentially decaying
minimum gap $\Omega_{0}$, $P\left(t_{f}\right)$ can only reach appreciable
values if the detuning $h$ is $O\left(\Omega_{0}\right)$ or less. Finally, to
regularize behavior at small $L$–and consequently provide a larger range of
reliable data to fit–we also perform our simulations with a single transverse
field, at site 0, weakened by a constant prefactor. This reduces the
degeneracy splitting $\Omega_{0}$ by a prefactor without changing the scaling
with $L$; absent this step, at larger $\kappa$ values and small $L$, tunneling
occurs too quickly to allow us to employ long enough runtimes to use small
frequencies and thus, mitigate heating from the AC drives. For fair
comparisons the same single field weakening is used in all simulations.

##  4 Results

###  4.1 Acceleration of Tunnelling

We study the transverse-field values of $\kappa=0.7,0.8,0.9$, with AC
modulation amplitudes $\alpha$ between 0 and 0.5 for X-RFQA and ZZ-RFQA, and
0.3 for XZZ-RFQA. These peak amplitudes are chosen based on the observed
response of the system’s $L$-scaling to the various protocol choices. We
observe that in all the cases we considered, RFQA oscillations significantly
accelerate the global-spin-flip tunneling process at fixed scale $L$, as
confirmed by the increase of the average tunnelling probability compared to
the analogous DC protocol shown in Fig. 3 for $L=18$, and by the exponential
fits shown in Fig. 1 for a wide range of parameters.

In order to capture the $L$-dependence of this improvement we study the
approximate TTS ratio $t_{f}/\langle P(t_{f})\rangle$ vs $L$ for a fixed
choice of $\kappa,\alpha$ and RFQA variant (see 7). For all the DC protocols
and at least the smaller-$\kappa$, smaller-$\alpha$ RFQA protocols we observe
what appears _prima facie_ to be an exponential scaling of the TTS with $L$.
We fit the approximate TTS ratio $t_{f}/\langle P(t_{f})\rangle$ obtained from
the numerical data with a two-parameter fitting function $f(L)=aL2^{\Upsilon
L}$. At small $L$ there is of course some inherent ambiguity in the choice of
polynomial prefactor in such fits. We chose this form because in fits where
the polynomial prefactor was allowed to vary (e.g. $L^{c}$ instead of $L$) the
best fit exponent $c$ was always close to 1, and because the empirical scaling
form in Eq. 6 combined with $W\propto L^{-1/2}$ predicts prefactor exponents
close to 1 as well. We chose to fit all data using the same functional form
for consistency and to make the easiest comparisons. From this fit we obtain
the difficulty exponent $\Upsilon$. For the DC protocol, we used standard
exact diagonalization methods to check that the observed TTS scaling exponent
closely tracks the prediction $TTS\propto W/\Delta_{\min}^{2}$ obtained from
the empirical gap $\Delta_{\min}$ between the two dressed quasi-degenerate
ground states of the TFIM in zero longitudinal field, at fixed $\kappa$ (see
Appendix B for a demonstration of this).

We observe that for all RFQA protocols, the difficulty exponent $\Upsilon$
decreases as $\alpha$ is increased from zero to positive values (Fig. 1), with
XZZ-RFQA showing the most dramatic improvements and X-RFQA the most modest.
For large enough values of $\kappa$ and $\alpha$, the X- and XZZ-RFQA
protocols exhibit a crossover from positive to negative values of the
exponent, suggesting the end of the exponentially-diverging regime of the TTS.
We expect the X-RFQA variant, being the least performant of the three, will
exhibit an analogous crossover at larger $\kappa,\alpha$ values of the ones we
studied here. For all the three RFQA protocols we employes Time-Evolving Block
Decimation (TEBD) [35] methods in order to extend the calculations of the
average tunnelling probability to larger system sizes $L\geq 20$, for
$\kappa=0.9$ and $\alpha$ values close to the crossover point, obtaining
results largely consistent with the ones at smaller $L$ (see Fig. 3). The
shifts in the fitted difficulty exponents produced by aggregating the
additional TEBD data are very small: $0.048\rightarrow 0.062$ for the X- ,
$-0.026\rightarrow-0.018$ for the ZZ-, and $-0.0009\rightarrow-0.00005$ for
the XZZ-RFQA variants respectively.

Figure 3: Average tunnelling probability $\langle P(t_{f})\rangle$ vs system size $L$ at $\kappa=0.9$. In the main plot: average tunnelling probability curves for the DC protocol ($\alpha=0$) and RFQA protocols with values of $\alpha$ closest to where the the difficulty exponent $\Upsilon$ is changing sign. For X- (blue dots), ZZ- (orange crosses) and XZZ-RFQA (green diamonds) these are respectively $\alpha=0.5,0.4$ and $0.15$). Note that the average tunnelling probability for the RFQA protocols shown here is approximately constant in (or very weakly dependent on) the system size $L$, while for the DC protocol it is clearly exponentially decaying. The data for the larger system sizes $L\geq 20$ (dashed segments) were obtained using TEBD. Inset: $\langle P(t_{f})\rangle$ vs $L$ plot in the positive $\Upsilon$ regime. Shown here is the case with $\kappa=0.7$, with a choice of $\alpha=0.4$ for the X- and ZZ-RFQA protocols, and $\alpha=0.2$ for the XZZ-RFQA protocol. In this regime all probability curves decay exponentially with $L$, albeit with different rates which favour the RFQA protocols over the DC-driven one. Protocol | $\kappa$ | $\alpha$ | $\Upsilon_{ST}$ | $2\gamma_{dis}$ | $\left<Z_{0}Z_{L/2}\right>_{0}$ | $\left<Z_{0}Z_{L/2}\right>_{dis}$
---|---|---|---|---|---|---
X | 0.8 | 0.5 | 0.28 | 0.39 | 0.88 | 0.87
| 0.9 | 0.5 | 0.05 | 0.15 | 0.81 | 0.8
ZZ | 0.8 | 0.4 | 0.14 | 0.19 | 0.88 | 0.84
| 0.8 | 0.5 | 0.01 | 0.12 | 0.88 | 0.81
| 0.9 | 0.3 | 0.05 | 0.18 | 0.81 | 0.78
| 0.9 | 0.4 | -0.03 | 0.08 | 0.81 | 0.75
XZZ | 0.8 | 0.2 | 0.05 | 0.07 | 0.88 | 0.84
| 0.8 | 0.25 | -0.03 | 0.06 | 0.88 | 0.81
| 0.9 | 0.1 | 0.11 | 0.19 | 0.81 | 0.79
| 0.9 | 0.15 | 0.0 | 0.06 | 0.81 | 0.77
Table 1: Evidence for an irreducibly AC nature of tunneling acceleration in
this system, by comparison of the fitted TTS exponent to the exponential decay
of the DC average $\left<\Omega_{0}^{2}\right>$ over the instantaneous
modulation of Hamiltonian parameters appropriate to that protocol. Two point
correlations are reported for $L=20$; see text for more details.

###  4.2 DC parameter fluctuations and ferromagnet-paramagnet phase
transitions do not explain the observed speedup

In all three of our simulated protocols, the transverse field and coupling
terms are modulated fairly substantially. Consequently, in rare events, one
will find appreciable shifts in the mean values of $\kappa$ or $J$, and in
those cases the instantaneous collective tunnel splitting $\Omega_{0}$ will be
much larger than in the corresponding uniform field case, potentially even
large enough to induce a transition to the paramagnetic phase. A skeptical
reader could very reasonably ask if rare large fluctuations are sufficient to
explain the speedups we observe, and since the challenge of implementing
simultaneous AC modulation of all terms is significant, it is important to
rule out more prosaic explanations for the fast tunneling we report.

We first consider DC fluctuations in $\Omega_{0}$. As argued above, in the DC
case the detuning-averaged tunneling rate $\Gamma\propto\Omega_{0}^{2}/W$;
this remains true with disorder if we replace
$\Omega_{0}^{2}\to\left<\Omega_{0}^{2}\right>_{dis}$. Here, the disorder
average is over random modulations of the appropriate terms equivalent to
taking instantaneous time slices of $H\left(t\right)$ in Eq. 2 for the
appropriate protocol. To check whether this effect is able to explain the
speedup, we used exact diagonalization with $h=0$ to compute
$\left<\Omega_{0}^{2}\right>$ for eight protocol/parameter choices (all at or
close to the potential scaling crossover regime), with 2000 random disorder
realizations for each datapoint and $L$ running from 8 to 20. To compare to
the AC-driven tunneling rate $\Upsilon$, we numerically fit

$\displaystyle\left<\Omega_{0}^{2}\left(L\right)\right>=A\frac{2^{-2\gamma_{dis}L}}{L}.$
(10)

This is the same scaling form used to extract $\Upsilon$; as shown in Table 1,
the resulting exponents are all larger by a shift of 0.05-0.13, well outside
any fitting or sampling uncertainty here. DC fluctuations in $\Omega_{0}$ thus
cannot explain the tunneling rates we observe. Further, since $\Omega_{0}$ is
exponentially sensitive to mean value shifts, $\left<\Omega_{0}^{2}\right>$
returns significantly larger shifts (compared to the uniform field) than
$\left<\Omega_{0}\right>$, and the _median_ value of $\Omega_{0}$ at each $L$
shows very little change compared to the uniform case.

We now address the second possibility, transitions into the paramagnetic
phase. In the clean limit this transition has a gap that scales as $1/L$; with
disorder this becomes a stretched exponential [41, 42, 43], which is still
substantially larger than the exponentially decaying gaps of transitions in
the ferromagnetic phase. Fast tunneling driven by such transitions would not
suggest utility for embedded problems, since in the paramagnetic phase long-
ranged correlation (and thus, energetic awareness of the logical problem) is
lost. Such transitions would also likely lead to signficant heating through
the Kibble-Zurek mechanism. To ensure that this is not happening, as mentioned
earlier we tracked the ferromagnetic order parameter
$\left<Z_{0}Z_{L/2}\right>$, averaged over the “waiting” phase of AC-driven
evolution (and over detuning and frequency distributions). As shown in FIG. 4,
and for DC eigenstates in Table 1, this order parameter is modestly reduced by
the modulation terms but we see no evidence that it will vanish at large $L$.
Interestingly, the AC averages are slightly _larger_ than the DC averages;
this is because the DC average is taken at degeneracy but the AC average
includes small detuning terms $h$ that bias the system further toward
ferromagnetism. We thus conclude that global mixing with paramagnetic states
does not explain the speedup we report.

Figure 4: Time-averaged two-point correlation function $C_{ZZ}$ in Eq. 8 vs
system size $L$, at fixed $\kappa=0.9$. The time-averaging is taken over the
“plateau” part of the protocol in between the two ramps, when the global
transverse field value $\kappa$ does not depend on time. The the main plot
shows the curves for the DC protocol ($\alpha=0$, solid black line) and for
each of the RFQA protocols, for the value of $\alpha$ closest to the where the
TTS scaling exponent crosses from positive to negative values. For X- (blue
dots), ZZ- (orange crosses) and XZZ-RFQA (green diamonds) these values are
respectively $\alpha=0.5$, $0.4$ and $0.15$. The data for the larger system
sizes $L\geq 20$ (dashed segments) were obtained using TEBD. In the inset, the
curves for the largest values of $\alpha$ we studied. For X-, ZZ- and XZZ-RFQA
respectively these are $\alpha=0.5,0.5$ and $0.3$.

##  5 Implementation Prospects and Discussion

Through extensive numerical simulations, we have demonstrated significant
acceleration of multiqubit tunneling in transverse field Ising chains through
variations of RFQA, assuming $O\left(1/\sqrt{L}\right)$ energy uncertainty (in
comparison to an exponentially decaying minimum gap). In many of our
simulations, exponential and polynomial scaling of the average tunneling time
could not be distinguished from the data, in contrast to the clear exponential
decay of all uniform field cases. While we reported only 1d simulations here,
preliminary full state evolution studies of weakly coupled rings showed
similar speedups when compared to single rings with the same total number of
qubits; we will present those results in future work. As ST is not a
1d-specific mechanism, this is not surprising, and suggests these results
should generalize to embedded _problems_ , which are far out of reach of
simulation.

Of course, the ultimate goal of these AC methods is to address two of the core
physics problems in quantum annealing (embedding overhead and computational
weakness of a uniform transverse field), in a simple and scalable manner.
While any RFQA implementation is undoubtedly more complex than current uniform
transverse field hardware, as we are only varying flux-tunable
$X_{j}/Z_{j}Z_{k}$ terms already present in the system, no direct changes to
the qubit or coupler hardware are required and the complexity lies entirely on
the control side. Further, we have shown that reusing the same frequencies
across the lattice has performance benefits, and $O\left(\sqrt{N}\right)$
unique tones reused across $O\left(N\right)$ qubits is sufficient to rapidly
accelerate tunneling while avoiding local transitions out of the ferromagnetic
phase for individual chains. This favorable condition reduces our methods’
implementation complexity to a local addressing and/or signal routing problem,
which while difficult is still substantially easier than the expected
requirements for topological error correction codes [44].

With regard to real implementations, it is likewise important to note that all
of these simulations have been noise-free (though they do incorporate random
detuning as a proxy for the many energy uncertainty sources in real
experiments). This may strike the reader as a significant oversight in a paper
making claims about improving near-term hardware. We chose to forgo noise
simulation for a few reasons. First, we expect our results to apply (with
suitable tuning and modifications) to both quantum annealers and neutral
atoms, which have extremely different noise models that are both difficult to
simulate in their own ways.

In neutral atoms, the error model (see the supplemental information of [13])
is dominated by laser noise (amplitude and phase fluctuations of the drive
beams), decay from the Rydberg state (a leakage error) and atom loss.
Amplitude fluctuations111Here, we assume these fluctuations can be local or
global; in previous experiments global beams were used, but obviously local
control would be involved in any implementation of the protocols we discuss in
this work. in the drive lasers–weak, random modulation of the transverse
fields–are going to be negligible compared to the larger modulations used to
implement RFQA; phase variations correspond to slowly fluctuating $Z$ biases
that stymie attempts to hit exponentially narrow resonances (at sufficiently
large system size) but are less deleterious to transitions active over wider
detuning ranges. Leakage and atom loss are more serious, in that they change
the problem graph, and they are a common challenge to all neutral atom
protocols. Of course, an RFQA implementation cannot solve these problems,
though we do not expect it would meaningfully increase their rates, and since
they are comparatively slow error sources, _any_ methods capable of finding
the solution in reduced quantum evolution time can reduce their impact.

For quantum annealers, the picture is more complex, and one which was reviewed
in [27]. The error model in quantum annealers primarily consists of
quasistatic control error modeled by small random fluctuations in the problem
Hamiltonian parameters and transverse field strengths, $1/f$-like noise along
$Z$ for each qubit, and comparatively strong interaction with a cold (but not
zero temperature, e.g. $k_{B}T\ll J$ but $k_{B}T\gg\Omega_{0}$ at small
avoided crossings) bath. RFQA does not _directly_ mitigate control error or
$1/f$ longitudinal noise, the effects of which are at least partially captured
by our detuning averaging. It is possible that RFQA could _indirectly_ reduce
the impact of both error sources by allowing the user to work with larger
intra-chain ferromagnetic coupling $J$ (which can lead to rapid freezing in
the uniform field case), increasing the local gap and overall energy scale of
the problem. Interaction with the cold bath is more complex, and likely
prohibitively difficult to simulate in the AC-driven regime [45, 46]; in a
work by two of us [47], a Hilbert space of five hundred states was used to
simulate the bath interaction with two qubits, and this method does not scale
to long evolution times. The original RFQA work established that AC driving
accelerates bath-assisted phase transitions as well, in line with the benefits
seen in uniform field annealing with a pause in the schedule [17], but it is
possible that the AC driving could also amplify harmful bath effects though an
as-yet unknown mechanism. Since we cannot rigorously simulate the most
important and interesting open system effects in a flux qubit implementation,
we found it unnecessary to simulate other forms of noise for this work.

We want to close with an interesting question that this work raises, but
certainly does not settle. Namely, the generic expectation (with some notable
counterexamples [48]) of first order quantum phase transitions between ground
states is that the collective spin rearrangements defining the transition are
exponentially slow in the number of participating degrees of freedom, e.g.
$\Omega_{0}\propto 2^{-cK}$. What we ask is the following: are there
physically realistic–and hopefully, application relevant–first order phase
transitions where, under the influence of AC driving, the nature of the
transition does not change (e.g. all order parameters that are finite across
the transition remain so, and the system does not meaningfully heat), but the
collective tunneling rate crosses over to polynomial scaling even for
arbitrarily large system sizes? Our results suggest that this is _possible_
but we make no claims based on numerical evidence alone that such scaling must
persist at $L\to\infty$. A conclusive answer to this question would be of
significant interest for quantum optimization, and dynamical many-body physics
more generally.

##  6 Acknowledgements

We would like to thank Steven Dissler, Andrew King, Glen Mbeng, Eleanor
Rieffel, Paul Varosy and Steven Weber for useful discussions around this
project. We also would like to thank Zhijie Tang for assisting in early
simulations of RFQA in 1d chains. EK and VO were jointly supported by DARPA
under the Reversible Quantum Machine Learning and Simulations program,
contract HR00112190068. EK’s research in this area was also funded by NSF
grant PHY-1653820. GM would like to acknowledge support from the NASA Ames
Research Center and from DARPA under IAA 8839, Annex 128. Resources supporting
this work were also provided by the NASA High-End Computing (HEC) Program
through the NASA Advanced Supercomputing (NAS) Division at Ames Research
Center.

##  Appendix A Linear regime

The main text includes an explanation as to why one should see a linear growth
regime $\langle P\left(t_{f}\right)\rangle=\Gamma\,t_{f}$ for the average
tunnelling probability as a function of the final time of the protocols. By
analogy, the prefactor $\Gamma$ can informally be interpreted as an “average
tunnelling rate”. This linear-growth prediction is confirmed in Fig. 5 by
means of example, where we compare $\langle P\left(t_{f}\right)\rangle$ vs the
rescaled final time $t_{f}/L$ in the case of X-RFQA for
$\kappa=0.8,\alpha=0.3$. By using this linear-growth expression for $\langle
P\left(t_{f}\right)\rangle$ in Eq. (7) one can see that in the large-$L$ limit
we have that $TTS\sim 1/\Gamma$ and it would in principle be possible to
extract the difficulty exponent $\Upsilon$ by studying the average rate
$\Gamma$ of the linear regime. In practice, however, this becomes quickly very
expensive since $\Gamma$ is typically exponentially decaying with $L$, and
extracting its value requires fitting a linear function with a slope that even
at moderately large values of $L$ becomes almost indistinguishable from zero.
In order to do so reliably, then one would have to perform simulations for
times $t_{f}$ that grow exponentially with $L$. For this reason we decide to
extract the difficulty exponent $\Upsilon$ by fitting the scaling of the TTS
directly.

Figure 5: Emergence of the linear regime of the average tunnelling probability
$\langle P(t_{f})\rangle=\Gamma\,t_{f}$, for the X-RFQA protocol at
$\kappa=0.8$ and $\alpha=0.3$. Dashed lines represent linear fits. The
difficulty exponent $\Upsilon$ and the average rate $\Gamma$ (_i.e._ the slope
of the linear growth) are related through Eq. (1).

##  Appendix B DC protocol: minimum-gap scaling exponent and difficulty
exponent

We studied the scaling of the empirical gap between the two quasi-degenerate
dressed ferromagnetic states at given $\kappa$ using exact diagonalization
methods. We extracted the scaling exponent $\hat{\Upsilon}$ of the quantity
$1/\Delta_{\min}^{2}$ by fitting it with the function
$f(L)=a2^{\hat{\Upsilon}L}$, with fit parameters $a,\hat{\Upsilon}$. The TTS
datapoints obtained from the dynamical simulations of the DC protocol for the
same given $\kappa$ were then fitted with the single-parameter function
$TTS(L)=b\,2^{\hat{\Upsilon}L}$ by using the $\hat{\Upsilon}$ value extracted
before. Fig. 6 shows that the scaling exponent $\hat{\Upsilon}$ of the inverse
square gap $1/\Delta_{\min}^{2}$ corresponds to the difficulty exponent
$\Upsilon$ defined by the TTS’s exponential scaling with $L$. Formally,
$\lim_{L\rightarrow\infty}\frac{1}{L}\log_{2}(1/\Delta^{2})=\lim_{L\rightarrow\infty}\frac{1}{L}\log_{2}TTS=\Upsilon$.

Figure 6: Comparison between the scaling of the TTS for the DC protocol and
the inverse-square-gap prediction obtained from the ($\kappa$-dressed)
empirical gap $\Delta_{min}$ between the two ferromagnetic states of the
($h=0$) transverse-field Ising chain. Datapoints correspond to dynamical data
from the numerics, while dashed lines are the single-parameter fit
$f(L)=b2^{\hat{\Upsilon}L}$, where $\hat{\Upsilon}$ is the scaling exponent of
the quantity $1/\Delta_{min}^{2}$, as explained in the main text.

##  Appendix C TTS and difficulty exponent

As explained in the main text, in order to calculate the difficulty exponent
$\Upsilon$ (shown in Fig. 1) for a particular protocol and fixed values of
$\kappa,\alpha$, we plot the large-$L$ approximation to the TTS,
$t_{f}/\langle P\left(t_{f}\right)\rangle$ vs the system size $L$ and fit the
data with the two-parameter fitting function $f(L)=aL2^{\Upsilon L}$,
obtaining $\Upsilon$. Fig. 7 (left) shows this for the ZZ-RFQA protocol at
$\kappa=0.8$ and $\alpha$ between zero (the DC protocol) and $0.5$. The fitted
$\Upsilon$ decreases monotonically with $\alpha$. Analogous behaviours are
observed for all the other RFQA protocols in all parameter ranges we have
studied.

The distribution of the tunneling probabilities $P(t_{f})$ depends on the
random choices longitudinal fields $h$, and (for the AC-driven protocols) on
the random frequencies and phases. It is undesirable that this distributions –
and as a consequence, the value of the TTS – should be dominated by rare
events. If this were the case, then any observed acceleration of tunneling
could be indicative of an atypical behaviour of the RFQA protocols. In order
to rule out this possibility we computed the median of the tunnelling
probability $P_{med}\left(t_{f}\right)$ and fitted the ratio
$t_{f}/P_{med}\left(t_{f}\right)$ vs $L$ using the same fitting function as
before. The results for the ZZ-RFQA protocol at $\kappa=0.8$ and
$\alpha=0-0.5$ are shown in Fig. 7 (right). Even though we obtain slightly
different scaling exponents, the overall picture is unchanged: the “difficulty
exponent” $\Upsilon_{med}$ obtained in this way connects to the DC value at
$\alpha=0$ and crosses into negative values at for large enough $\alpha$.

Figure 7: Left: The $L$-dependence of the approximate TTS ratio $t_{f}/\langle
P\left(t_{f}\right)\rangle$ obtained from the numerical simuations is fitted
with the two-parameter function $f(L)=aN2^{\Upsilon L}$ (dashed lines) in
order to extract the difficulty exponent $\Upsilon$, in the exemplar case of
the ZZ-RFQA protocol for $\kappa=0.8$ and $\alpha$ values between 0 and 0.5.
Right: the ratio $t_{f}/P_{med}(t_{f})$ defined through the _median_
tunnelling probability $P_{med}(t_{f})$ is fitted with the same fitting
function, yielding $\Upsilon_{med}$.

##  Appendix D Heating

In this work we use RFQA to accelerate tunnelling from one ferromagnetic state
of a quantum Ising chain to the other. Ideally, this should happen without
exciting the system out of the the quasi-degenerate ferromagnetic doublet. We
use the _heating_ quantity $P_{heat}$ defined through Eq. (4) in order to
assess the amount of excitations produced by the various protocols. In
particular, we compare $P_{heat}$ with the tunnelling probability
$P\left(t_{f}\right)$ in order to rule out the possibility that $P_{heat}\gg
P\left(t_{f}\right)$, which was observed in [27] to obscure the possible
advantage achieved by RFQA. For all DC protocols, and all RFQA protocols at
large enough $\alpha$, we have that $\langle P_{heat}\rangle<\langle
P(t_{f})\rangle$ for all the system sizes studied. Some RFQA protocols go
through an intermediate regime at small $\alpha$ where $\langle
P\left(t_{f}\right)\rangle\lesssim\langle P_{heat}\rangle$, for the larger
values of $L$. This is due to a comparatively larger amount of heat than what
is observed in the DC case for the same $\kappa$, combined with too modest an
increase of the tunnelling probability (See Fig. 8). A qualitatively analogous
behaviour is observed by considering the median heat and the median tunnelling
probability.

Given the empirically-observed behaviour of $P(t_{f})$ and $P_{heat}$, one
could reasonably imagine that $P_{heat}\gg P\left(t_{f}\right)$ would
eventually hold in the $L\rightarrow\infty$ limit. For the DC case, this issue
can arguably be solved by using better adiabatic ramps, but for the AC-case
this prediction seems hard to confirm or deny conclusively without a
theoretical description of the AC-induced excitation processes that holds for
non-perturbative values of $\kappa$ and $\alpha$, which we currently lack.
Nevertheless, we do not expect the heating we observe _at given size_ to
significantly affect the results presented in this work.

Figure 8: Comparisons between average $P(t_{f})$ (solid black line), median
$P(t_{f})$ (black crosses), average $P_{heat}$ (dashed red line) and median
$P_{heat}$ (red squares) in DC- and AC-driven protocols, as functions of the
system size $L$. Leftmost column: X-RFQA protocol for $\kappa=0.7$ and
$\alpha=0,0.2,0.4$. Center column: ZZ-RFQA protocol for $\kappa=0.8$ and
$\alpha=0,0.1,0.4$. Rightmost column: ZZX-RFQA protocol for $\kappa=0.9$ and
$\alpha=0,0.1,0.2$.

##  Appendix E Considerations for TEBD in this system

For the one-dimensional model studied in this work we use Time-Evolving Block
Decimation [35], as implemented by the ITensor library [38, 39] to access
larger system sizes. Unlike the exact-state simulation code used for the
smaller $L$, TEBD uses matrix-product states (MPS) with finite bond dimension
$\chi$ in order to approximate the time-evolving wavefunction. Physical
results are accessed by increasing $\chi$ until convergence is reached (up to
the desired numerical accuracy). Unlike other MPS-based methods such as the
Density Matrix Renormalization Group (DMRG), where the one-dimensional area
law allows for an efficient representation of the ground state of a local
Hamiltonian by MPSs and _e.g._ for the calculation of intensive quantities (of
order $O(1)$), the calculation of quantities that exhibit an exponential decay
with the system size (such as the tunnelling probability $P(t_{f})$ of our
model) require increasingly better approximations of the time-evolving
wavefunction, even in cases where the entanglement is not diverging with $L$.

## References

  * [1] AB Finnila, MA Gomez, C Sebenik, C Stenson, and JD Doll.  Quantum annealing: a new method for minimizing multidimensional functions.  Chemical physics letters, 219(5-6):343–348, 1994.
  * [2] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  Physical Review E, 58(5):5355, 1998.
  * [3] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser.  Quantum computation by adiabatic evolution.  arXiv:quant-ph/0001106, 2000.
  * [4] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  Quantum adiabatic evolution algorithms versus simulated annealing.  arXiv preprint quant-ph/0201031, 2002.
  * [5] Arnab Das and Bikas K Chakrabarti.  Colloquium: Quantum annealing and analog quantum computation.  Reviews of Modern Physics, 80(3):1061, 2008.
  * [6] Mark W Johnson, Mohammad HS Amin, Suzanne Gildert, Trevor Lanting, Firas Hamze, Neil Dickson, R Harris, Andrew J Berkley, Jan Johansson, Paul Bunyk, et al.  Quantum annealing with manufactured spins.  Nature, 473(7346):194, 2011.
  * [7] Sergio Boixo, Troels F Rønnow, Sergei V Isakov, Zhihui Wang, David Wecker, Daniel A Lidar, John M Martinis, and Matthias Troyer.  Evidence for quantum annealing with more than one hundred qubits.  Nature Physics, 10(3):218, 2014.
  * [8] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * [9] Tameem Albash and Daniel A. Lidar.  Adiabatic quantum computing.  arXiv:1611.04471, 2017.
  * [10] Tameem Albash and Daniel A Lidar.  Demonstration of a scaling advantage for a quantum annealer over simulated annealing.  Physical Review X, 8(3):031016, 2018.
  * [11] Pascal Scholl, Michael Schuler, Hannah J Williams, Alexander A Eberharter, Daniel Barredo, Kai-Niklas Schymik, Vincent Lienhard, Louis-Paul Henry, Thomas C Lang, Thierry Lahaye, et al.  Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms.  Nature, 595(7866):233–238, 2021.
  * [12] Pascal Scholl, Hannah J Williams, Guillaume Bornet, Florian Wallner, Daniel Barredo, L Henriet, Adrien Signoles, Clément Hainaut, Titus Franz, S Geier, et al.  Microwave engineering of programmable x x z hamiltonians in arrays of rydberg atoms.  PRX Quantum, 3(2):020303, 2022.
  * [13] Sepehr Ebadi, Alexander Keesling, Madelyn Cain, Tout T Wang, Harry Levine, Dolev Bluvstein, Giulia Semeghini, Ahmed Omran, J-G Liu, Rhine Samajdar, et al.  Quantum optimization of maximum independent set using rydberg atom arrays.  Science, 376(6598):1209–1215, 2022.
  * [14] Alejandro Perdomo-Ortiz, Salvador E Venegas-Andraca, and Alán Aspuru-Guzik.  A study of heuristic guesses for adiabatic quantum computation.  Quantum Information Processing, 10(1):33–52, 2011.
  * [15] Andrew D King, Juan Carrasquilla, Jack Raymond, Isil Ozfidan, Evgeny Andriyash, Andrew Berkley, Mauricio Reis, Trevor Lanting, Richard Harris, Fabio Altomare, et al.  Observation of topological phenomena in a programmable lattice of 1,800 qubits.  Nature, 560(7719):456–460, 2018.
  * [16] Nicholas Chancellor.  Modernizing quantum annealing using local searches.  New Journal of Physics, 19(2):023024, 2017.
  * [17] Jeffrey Marshall, Davide Venturelli, Itay Hen, and Eleanor G Rieffel.  Power of pausing: Advancing understanding of thermalization in experimental quantum annealers.  Physical Review Applied, 11(4):044083, 2019.
  * [18] Andrew D King, Jack Raymond, Trevor Lanting, Richard Harris, Alex Zucca, Fabio Altomare, Andrew J Berkley, Kelly Boothby, Sara Ejtemaee, Colin Enderud, et al.  Quantum critical dynamics in a 5,000-qubit programmable spin glass.  Nature, pages 1–6, 2023.
  * [19] Vicky Choi.  Minor-embedding in adiabatic quantum computation: I. the parameter setting problem.  Quantum Information Processing, 7(5):193–209, 2008.
  * [20] Vicky Choi.  Minor-embedding in adiabatic quantum computation: Ii. minor-universal graph design.  Quantum Information Processing, 10(3):343–353, 2011.
  * [21] Mario S Könz, Wolfgang Lechner, Helmut G Katzgraber, and Matthias Troyer.  Embedding overhead scaling of optimization problems in quantum annealing.  PRX Quantum, 2(4):040322, 2021.
  * [22] Minhyuk Kim, Kangheun Kim, Jaeyong Hwang, Eun-Gook Moon, and Jaewook Ahn.  Rydberg quantum wires for maximum independent set problems with nonplanar and high-degree graphs.  arXiv preprint arXiv:2109.03517, 2021.
  * [23] Minh-Thi Nguyen, Jin-Guo Liu, Jonathan Wurtz, Mikhail D Lukin, Sheng-Tao Wang, and Hannes Pichler.  Quantum optimization with arbitrary connectivity using rydberg atom arrays.  PRX Quantum, 4(1):010316, 2023.
  * [24] Jeffrey Marshall, Eleanor G Rieffel, and Itay Hen.  Thermalization, freeze-out, and noise: Deciphering experimental quantum annealers.  Physical Review Applied, 8(6):064025, 2017.
  * [25] Ryan Hamerly, Takahiro Inagaki, Peter L McMahon, Davide Venturelli, Alireza Marandi, Tatsuhiro Onodera, Edwin Ng, Carsten Langrock, Kensuke Inaba, Toshimori Honjo, et al.  Experimental investigation of performance differences between coherent ising machines and a quantum annealer.  Science advances, 5(5):eaau0823, 2019.
  * [26] Matthew Kowalsky, Tameem Albash, Itay Hen, and Daniel A Lidar.  3-regular 3-xorsat planted solutions benchmark of classical and quantum heuristic optimizers.  arXiv preprint arXiv:2103.08464, 2021.
  * [27] Eliot Kapit and Vadim Oganesyan.  Noise-tolerant quantum speedups in quantum annealing without fine tuning.  Quantum Science and Technology, 6(2):025013, 2021.
  * [28] Christof Zalka.  Grover’s quantum searching algorithm is optimal.  Physical Review A, 60(4):2746, 1999.
  * [29] Sergei V Isakov, Guglielmo Mazzola, Vadim N Smelyanskiy, Zhang Jiang, Sergio Boixo, Hartmut Neven, and Matthias Troyer.  Understanding quantum tunneling through quantum monte carlo simulations.  Physical review letters, 117(18):180402, 2016.
  * [30] Evgeny Andriyash and Mohammad H Amin.  Can quantum monte carlo simulate quantum annealing?  arXiv preprint arXiv:1703.09277, 2017.
  * [31] Zhang Jiang, Vadim N Smelyanskiy, Sergei V Isakov, Sergio Boixo, Guglielmo Mazzola, Matthias Troyer, and Hartmut Neven.  Scaling analysis and instantons for thermally assisted tunneling and quantum monte carlo simulations.  Physical Review A, 95(1):012322, 2017.
  * [32] Zhang Jiang, Vadim N Smelyanskiy, Sergio Boixo, and Hartmut Neven.  Path-integral quantum monte carlo simulation with open-boundary conditions.  Physical Review A, 96(4):042330, 2017.
  * [33] Andrew D King, Jack Raymond, Trevor Lanting, Sergei V Isakov, Masoud Mohseni, Gabriel Poulin-Lamarre, Sara Ejtemaee, William Bernoudy, Isil Ozfidan, Anatoly Yu Smirnov, et al.  Scaling advantage in quantum simulation of geometrically frustrated magnets.  arXiv preprint arXiv:1911.03446, 2019.
  * [34] Elizabeth Crosson, Edward Farhi, Cedric Yen-Yu Lin, Han-Hsuan Lin, and Peter Shor.  Different strategies for optimization using the quantum adiabatic algorithm.  arXiv preprint arXiv:1401.7320, 2014.
  * [35] Guifré Vidal.  Efficient simulation of one-dimensional quantum many-body systems.  Phys. Rev. Lett., 93:040502, Jul 2004.
  * [36] Guifré Vidal.  Efficient classical simulation of slightly entangled quantum computations.  Phys. Rev. Lett., 91:147902, Oct 2003.
  * [37] Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, and Claudius Hubig.  Time-evolution methods for matrix-product states.  Annals of Physics, 411:167998, 2019.
  * [38] Matthew Fishman, Steven R. White, and E. Miles Stoudenmire.  The ITensor Software Library for Tensor Network Calculations.  SciPost Phys. Codebases, page 4, 2022.
  * [39] Matthew Fishman, Steven R. White, and E. Miles Stoudenmire.  Codebase release 0.3 for ITensor.  SciPost Phys. Codebases, pages 4–r0.3, 2022.
  * [40] Itay Hen, Joshua Job, Tameem Albash, Troels F. Rønnow, Matthias Troyer, and Daniel A. Lidar.  Probing for quantum speedup in spin-glass problems with planted solutions.  Phys. Rev. A, 92:042325, Oct 2015.
  * [41] Daniel S Fisher.  Critical behavior of random transverse-field ising spin chains.  Physical review b, 51(10):6411, 1995.
  * [42] Jacek Dziarmaga.  Dynamics of a quantum phase transition in the random ising model: Logarithmic dependence of the defect density on the transition rate.  Physical Review B, 74(6):064416, 2006.
  * [43] Tommaso Caneva, Rosario Fazio, and Giuseppe E Santoro.  Adiabatic quantum dynamics of a random ising chain across its quantum critical point.  Physical Review B, 76(14):144427, 2007.
  * [44] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland.  Surface codes: Towards practical large-scale quantum computation.  Phys. Rev. A 86, 032324, 2012.
  * [45] Daniel Jaschke, Lincoln D Carr, and Inés de Vega.  Thermalization in the quantum ising model?approximations, limits, and beyond.  Quantum Science and Technology, 4(3):034002, 2019.
  * [46] Eliot Kapit, Pedram Roushan, Charles Neill, Sergio Boixo, and Vadim Smelyanskiy.  Entanglement and complexity of interacting qubits subject to asymmetric noise.  Physical Review Research, 2(4):043042, 2020.
  * [47] Eliot Kapit and Vadim Oganesyan.  Small logical qubit architecture based on strong interactions and many-body dynamical decoupling.  arXiv preprint arXiv:2212.04588, 2022.
  * [48] C. Laumann, R. Moessner, A. Scardicchio, and S. L. Sondhi.  The quantum adiabatic algorithm and scaling of gaps at first order quantum phase transitions.  Phys. Rev. Lett. 109, 030502, 2012.