# The Advantage of pausing: parameter setting for quantum annealers

Zoe Gonzalez Izquierdo  QuAIL, NASA Ames Research Center, Moffett Field,
California 94035, USA  USRA Research Institute for Advanced Computer Science,
Mountain View, California 94043, USA  Shon Grabbe  QuAIL, NASA Ames Research
Center, Moffett Field, California 94035, USA  Husni Idris  Aviation Systems
Division NASA Ames Research Center, Moffett Field, California 94035, USA
Zhihui Wang  QuAIL, NASA Ames Research Center, Moffett Field, California
94035, USA  USRA Research Institute for Advanced Computer Science, Mountain
View, California 94043, USA  Jeffrey Marshall  QuAIL, NASA Ames Research
Center, Moffett Field, California 94035, USA  USRA Research Institute for
Advanced Computer Science, Mountain View, California 94043, USA  Eleanor
Rieffel  QuAIL, NASA Ames Research Center, Moffett Field, California 94035,
USA

###### Abstract

Prior work showed the efficacy of pausing mid-anneal: such a pause improved
the probability of success by orders of magnitude in a class of native problem
instances and improved the time-to-solution in a class of embedded problem
instances. A physics-based picture provides qualitative suggestions for where
pausing mid-anneal is effective, for the interplay between annealing schedule
parameters and other annealing properties and parameters such as embedding
size and strength of the ferromagnetic coupling $|J_{F}|$, and for the
conditions under which pausing can improve the time-to-solution. Here, through
studies on an updated annealing architecture that has higher connectivity than
previous annealers, and on multiple embedded problem classes, we are able to
confirm various aspects of this picture. We demonstrate the robustness of the
optimal pause parameters across platforms and problem classes, explore how to
set $|J_{F}|$ to optimize performance in different scenarios, and provide
experimental evidence that short pauses trump longer overall annealing times
in time-to-solution. We also identify the number of different coefficients in
a problem as a predictor of problem hardness, and explore its interplay with
the optimal $|J_{F}|$ and embedding size. Based on these results we are able
to present qualitative guidelines for parameter setting in quantum annealers.

##  I Introduction

The theory behind quantum annealing has now existed for over two decades
Finnila  _et al._ (1994); Farhi  _et al._ (2000); Kadowaki and Nishimori
(1998); Farhi  _et al._ (2001); Aharonov  _et al._ (2007); Santoro  _et al._
(2002), and actual physical quantum annealers for more than one Harris  _et
al._ (2010, 2010); Johnson  _et al._ (2011); Boixo  _et al._ (2013); Bunyk
_et al._ (2014); Jiang  _et al._ (2018); Hauke  _et al._ (2020). Yet, the much
sought-after speedup for solving optimization problems with these devices has
remained elusive Rønnow  _et al._ (2014); Boixo  _et al._ (2014); Katzgraber
_et al._ (2015); Mandrà and Katzgraber (2018); Pearson  _et al._ (2019). One
way of improving performance for optimization problems is pausing mid-anneal
Marshall  _et al._ (2019); Venturelli and Kondratyev (2018); Passarelli  _et
al._ (2019); Winci  _et al._ (2020); Chen and Lidar (2020); Gonzalez Izquierdo
_et al._ (2021). By making clever use of the thermalization that takes place
in finite temperature annealers and taking into consideration the dynamics of
the system throughout the annealing process, this technique was first shown to
increase the ground state (GS) success probability, $p_{success}$, for native
problems Marshall  _et al._ (2019), and later confirmed to work for embedded
problems and adjusted to also improve the time-to-solution Gonzalez Izquierdo
_et al._ (2021). This makes pausing a necessary ingredient of benchmarking
studies going forward.

A quantum annealing solver has multiple parameters or variables whose values
can be set by the user, including choice of annealing schedule, total anneal
time and embedding parameters. Adding a pause further complicates this
landscape as two more parameters (location and duration of the pause) come
into play. Guidance in what regime to look for the optimal or good parameters
is therefore valuable for efficiently boosting success probability.
Benchmarking studies that help inform parameter setting for families of
problems have been a very useful tool to navigate these questions Rieffel  _et
al._ (2015); Vinci and Lidar (2016); Albash and Lidar (2018); Quiroz (2019);
Grant  _et al._ (2021), and they become crucial as quantum annealers offer
further opportunities to be controlled by user input.

As the landscape becomes more complicated, exploring a large range of values
and combinations becomes prohibitively costly. The aspiration to make
annealers useful for a multitude of applications adds to the issue, with
understanding how certain characteristics of the problems interact with
annealing parameters coming into play. Moreover, as newer and more diverse
quantum annealing devices become available Weber  _et al._ (2017), we need to
consider how their potential differences affect performance and parameter
setting.

We explore these questions in the present work, by testing several predictions
on optimal parameter regions across two annealing devices and three
optimization problems, to help provide parameter setting guidance in a variety
of scenarios. In particular, in this paper we:

  * •

Experimentally show that a short pause mid-anneal, in the region of pause
location 111For definitions of the pause location and other concepts related
to quantum annealing, see Sec. II. $s_{p}=0.3$-0.5, can improve time-to-
solution for ensembles of three different non-native problems as well as for
two quantum annealing devices with different architectures, further confirming
that the general location at which a pause will help remains consistent with
the physical picture describing the dynamics at different times in the anneal
Gonzalez Izquierdo  _et al._ (2021).

  * •

Confirm that within that region, larger embedding sizes and larger
ferromagnetic couplings will lead to an earlier optimal location, and vice
versa.

  * •

For a problem class with optimal annealing time higher than the hardware
minimum, experimentally show that introducing a short pause brings the optimal
down to at most the hardware minimum. This provides evidence for the fact that
extra time is only helpful in some regions of the anneal, and adds support to
the idea that, even with access to shorter annealing times, pausing will still
be beneficial.

  * •

Provide evidence for coefficient heterogeneity (i.e. the number of different
$h_{i}$ and $J_{ij}$ values in the logical problem) being a strong predictor
of problem hardness (partly due to experimental precision limitations). We do
this by considering three problems which QUBOs have different degrees of
coefficient heterogeneity, which gets amplified through the potential
coefficient splits over multiple qubits and couplers that occur during
embedding.

  * •

Show the more secondary role that problem size plays in problem hardness,
finding its effect to become apparent once coefficient heterogeneity stays
relatively constant.

  * •

Qualify the correlation between vertex model222We use vertex model to refer to
the set of physical qubits representing a single logical variable after minor-
embedding. The word “chain” is sometimes used in the literature with the same
meaning, but we prefer vertex model given that their structure is often not
chain-like. size and optimal $|J_{F}|$, showing that it gets complicated by
coefficient heterogeneity, rather than simply larger vertex models benefiting
from stronger couplings.

  * •

Through the above points, provide insights for parameter setting across
diverse classes of problems and hardware architectures, both general and more
specific based on knowledge of the logical and embedded problem, thus removing
the need to explore all regions when trying to optimize performance for
optimization problems.

  * •

Define an information sharing problem within the context of Collective
Autonomous Mobility Idris  _et al._ as one of the three problems analyzed in
this paper, and formulate it as a QUBO. In this problem, information sharing
among vehicles is optimized to mitigate communication bandwidth constraints in
high density traffic operations such as Urban Air Mobility Thipphavong  _et
al._ .

The rest of the paper is organized as follows. In Sec. II we give a brief
overview of quantum annealing, including a description of the two devices we
will be using. Sec. III introduces the three optimization problems (for which
QUBO mappings can be found in App. B, C, and D), as well as the metrics and
parameters used in the experimental runs. We present results for the same set
of instances solved on two different devices in Sec. IV.1, while Sec. IV.2
contains results for the three different problems solved on the newer device.
The findings from these results, as well as future research avenues are
discussed in Sec. V. A summary of our main points and closing remarks can be
found in Sec. VI.

##  II Background: Quantum Annealing

Quantum annealing is a quantum metaheuristic for optimization, and quantum
annealers are quantum hardware designed to run this metaheuristic. Any
classical cost function $C(x)$ that is a polynomial over binary variables
$x\in\\{0,1\\}^{n}$ can, with the addition of auxiliary variables, be turned
into a quadratic cost function. Problems with quadratic cost functions over
binary variables without additional constraints are called quadratic
unconstrained binary optimization (QUBO) problems, and are solvable by quantum
annealers (subject to size constraints imposed by each particular device).

Quantum annealing relies on the fact that the solution to such an optimization
problem can also be understood as the GS of an Ising problem Hamiltonian
Morita and Nishimori (2008)
$H_{p}=\sum_{<ij>}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}+\sum_{i}h_{i}\sigma_{i}^{x}$,
where $\sigma_{i}^{x}$ and $\sigma_{i}^{z}$ are individual Pauli matrices
acting on spin $i$, $J_{ij}$ represents the strength of the coupling between
spins $i$ and $j$, and $h_{i}$ that of the bias on spin $i$. Combinatorial
optimization problems can be expressed in Ising form by a straightforward
mapping between QUBO and Ising (mapping binary variables 0 and 1 to spin
variables $\pm 1$).

Then, to find the GS of $H_{p}$ (and thus the solution to the corresponding
QUBO) quantum annealing is carried out by evolving the system under a time-
dependent Hamiltonian $H(s)=A(s)H_{d}+B(s)H_{p}$, where $H_{d}$ is a driver
Hamiltonian (most commonly a transverse field
$H_{X}=-\sum_{i}\sigma_{i}^{x}$), $s$ is a dimensionless time parameter that
ranges from $0$ to $1$ over the course of a single anneal, and $A(s)$ and
$B(s)$ are device-dependent functions determining the strength of the driver
and problem Hamiltonians, respectively. These functions are such that $A(0)\gg
B(0)\approx 0$ and $B(1)\gg A(1)\approx 0$, so that the time-dependent
Hamiltonian evolves from $H_{d}$ to $H_{p}$ throughout the anneal. The system
is then initialized in the easy to prepare GS of $H_{d}$, and expected to
remain near the GS of $H(s)$ throughout, to finally yield the GS of $H_{p}$
(or equivalently, the solution to the optimization problem) upon final
measurement. More information about quantum annealing generally, including
mappings of optimization problems to QUBO can be found in Choi (2020); Rieffel
_et al._ (2015); Lucas (2014).

While the functions $A(s)$ and $B(s)$ cannot be changed by the user in
currently available devices, the dependence $s(t)$ can be modified under
certain constraints: $s$ must be a linear, non-decreasing function of $t$,
starting at 0 and ending at 1, with a maximum slope of $1\mu s^{-1}$, and a
limited number of slope changes allowed. This leads to different so-called
schedules $s(t)$. In particular, we will use schedules that include a pause,
that is, a period during which $ds/dt=0$, meaning that $H(s)$ is held constant
and the anneal is effectively paused.

For most application problems, since the hardware has restricted qubit
connectivity, the resulting QUBO is unlikely to conform to this hardware
connectivity, and a one-to-one correspondence between logical variables and
physical qubits does not exist. Instead, we must perform minor embedding Choi
(2008), which enables coupling between logical variables in the QUBO graph by
representing each of them by an appropriate set of physical qubits that will
allow the required connectivity. Following standard terminology in graph
theory, each such set of physical qubits is called a vertex model for its
corresponding logical variable.

Because all qubits in the vertex model should act as a single variable (i.e.
be aligned, otherwise the vertex model is considered to be broken, and the
solution containing it will be discarded), they are ferromagnetically coupled
to promote this collective behavior. We use the same coupling strength
$|J_{F}|$ for all the couplings within a vertex model ($J_{F}$ is always
negative, so we typically refer to its magnitude $|J_{F}|$). Problems that do
not require embedding because their structure matches that of the hardware are
called native problems for that hardware.

While $|J_{F}|$ can be set to a large value such that the embedded problem
preserves the GS of the logical problem, and analytical bounds on this value
can be obtained Choi (2008), too large a $|J_{F}|$ can reduce quantum
annealing performance. Physically there is an energy limit on the Hamiltonian
as a whole, and too large a $|J_{F}|$ relative to other parameters would mean
that all of the problem parameters could reduce performance due to precision
issues and noise in implementation. Furthermore, the energy spectrum
throughout the anneal varies with the value of $|J_{F}|$, and its effect on
the annealing often requires careful case-by-case consideration Choi (2020);
Rieffel  _et al._ (2015); Fang and Warburton (2019); Venturelli  _et al._
(2015); Marshall  _et al._ (2020, 2022). Thus, optimally setting the
ferromagnetic coupling $|J_{F}|$ is a challenging task. Prior work has shown
there is a sweet spot for this value. Physically this makes sense because a
stronger $|J_{F}|$ makes it less likely for individual qubits within a vertex
model to flip, which helps to avoid breaking the vertex model, but too large a
$|J_{F}|$ makes it increasingly costly for the vertex model qubit values to
flip together, potentially preventing the system from leaving a non-optimal
configuration. To boost the probability of success, $|J_{F}|$ must strike the
right balance, leading to better chances of arriving at—and staying in—the
correct configuration.

The schedule $s(t)$ can also significantly affect performance. Of particular
interest to us are schedules that include a pause where for some sub-interval
$s(t)$ is constant (i.e., $H(s)$ is constant for a specified time). It was
first observed for an ensemble of native problems that a pause at a location
(generally) insensitive to the instance specifics boosts the success
probability by orders of magnitude Marshall  _et al._ (2019), and a physical
picture was presented explaining this effect. It was later shown that it also
applies to embedded problems, and that with the right pause location and
duration Albash and Marshall (2021a) the time to solution can also be improved
Gonzalez Izquierdo  _et al._ (2021).

###  II.1 Quantum annealing devices

We use two different quantum annealing devices for our experimental runs. The
first part of our results will focus on a comparison across these devices, and
the same set of instances will be run in both. In the second part, where we
shift the focus to a comparison across problems, only one will be used.

The two devices we use are D-Wave 2000Q (DW2K), and D-Wave Advantage (DWA).
DWA is the newer of the two. While we discuss some of their features relevant
to our study below, more detailed technical information can be found in DWa
(2019) for DW2K, and in McGeoch and Farré (2020) for DWA.

Let us first talk about the similarities and differences between the two
quantum annealers that we use for these experimental runs. Both devices
perform qualitatively the same process, implementing the time-dependent
Hamiltonian

$H(s)=A(s)\sum_{i}\sigma_{i}^{x}+B(s)\left(\sum_{<i,j>}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}+\sum_{i}h_{i}\sigma_{i}^{z}\right),$
(1)

with $s$ the dimensionless time parameter going from 0 to 1 to perform one
anneal. Their driver Hamiltonian is the same, a transverse field over all
qubits that produces quantum fluctuations. Other quantitative details,
however, differ.

The annealing functions $A(s)$ and $B(s)$ present somewhat different shapes,
as shown in Fig. 12 (top) in Appendix A. The $B(s)$ is fairly similar in both
devices until the later part of the anneal, where DW2K ends up higher, and
$A(s)$ decays earlier and more rapidly for DWA than for DW2K. Their relative
strength is what determines the different regions in the dynamics of the
system, so it is more convenient to look at two dimensionless scales as
defined in Marshall  _et al._ (2017); $Q(s)=A(s)/B(s)$ corresponding to
quantum fluctuations and $C(s)=k_{B}T/B(s)$ to thermal ones. Note that the
temperature of the devices is different, with DW2K running at $T=12.1$ mK and
DWA at $T=15.8$ mK. These scales delineate three regimes within the anneal:
when the quantum fluctuations are much stronger than the thermal ones early
on, the system closely follows the GS of $H(s)$; in the middle, where both
scales are comparable, most of the dynamics takes place, and towards the end,
once the quantum fluctuations become much smaller than the thermal ones the
dynamics are essentially frozen, and no more transitions occur. But as we can
see in Fig. 12 (bottom), while the thermal fluctuations follow a similar curve
in both devices, the same is not true for the quantum ones. $C(s)$ and $Q(s)$
cross much earlier in the anneal for DWA, meaning the region of stronger
dynamics moves earlier.

The most obvious and significant difference between the two devices is in
their number of qubits and their connectivity. DW2K has a 2048 qubit
architecture (with 2031 working qubits), arranged in a Chimera graph. See the
top panel of Fig. 13 in Appendix A for a representation. Bipartite cells with
8 qubits each are arranged in a square pattern, connected to those on either
side as well as above and below. Except for those on the edges, each qubit is
connected to 6 others. A limitation of this architecture is its inability to
natively support odd cycles. DWA has a much larger and connected architecture,
with a 5760 qubit graph (5436 working ones) arranged according to a Pegasus
graph Boothby  _et al._ (2019), as shown in the bottom panel of Fig. 13. This
graph includes Chimera as a subgraph, but adds an extra layer of connectivity,
bringing each non-edge qubit up to degree 15, and making odd cycles available.

Although decreases in noise and improvement in fabrication have an impact on
the quality of results Pokharel  _et al._ (2021), we expect that the largest
difference in performance will be due to the reduced embedding size. The
fabrication process has not changed between the two devices we are using.
(There is a low noise version of DW2K which used a different fabrication
process, but we did not have access to it DWa (2020)). The quality of
solutions depends in part on implementing the correct problem. This is
determined by Integrated Control Errors (ICE) dwa , which encompass several
sources of problem infidelity. How much the actual values of couplings
$J_{ij}$ and fields $h_{i}$ will differ from the programmed ones depends on
many factors, including the specific values being programmed, the annealing
parameter $s$ and the time at which the problem is run, so it is difficult to
determine how much these errors are affecting our results. Measurements for
most of the error sources have not yet been released for DWA adv (2020) so,
given our current knowledge, we cannot presume any significant differences in
the effect of ICE for both devices.

##  III Methods

###  III.1 Problem statements

####  III.1.1 Minimum spanning tree with bounded degree

For our comparison across devices, we solve the minimum spanning tree with
bounded degree (BD-MST) problem. The results for the older of the two devices
correspond to the results form Gonzalez Izquierdo  _et al._ (2021). Spanning
trees are important for a number of reasons, including the designing of
efficient routing algorithms and for their wide-ranging applications to
network design, cluster analysis Grygorash  _et al._ (2006) and bioinformatics
Xu  _et al._ (2002). The problem statement is as follows: given a connected,
undirected graph $G=(V,E)$ with edge weights $w_{uv}$, $(uv)\in E$, and an
integer $\Delta\geq 2$, find a minimum weight spanning tree of maximum degree
at most $\Delta$. A spanning tree of $G$ is a subgraph of $G$ that is a tree
and contains all vertices of $G$. The tree weight to be minimized is
$\sum_{(uv)\in T}w_{uv}$.

Without the bound on the degree, determining if there exists a spanning tree
of weight $W$ for a graph $G$ can be decided in polynomial time, and different
efficient algorithms exist to find a minimum weight tree Cormen  _et al._
(2001). However, once we add the constraint on the degree, the problem becomes
NP-complete Garey and Johnson (1979). We use a level-based QUBO mapping that
can be found in Appendix B.

####  III.1.2 Graph coloring

Our second problem is graph coloring (GC). Given an undirected graph $G=(V,E)$
with $n$ nodes, and $k$ colors, the objective is to find a color assignment
for the nodes such that no nodes with the same color are connected by an edge.
Our instances are randomly generated 4-regular graphs with $n=12$, 14 and 16,
and $k=5$ colors, and we use a standard one-hot encoding for the mapping to
QUBO, described in Appendix C.

####  III.1.3 Information sharing

The third problem we consider is an information sharing (INFO) problem within
the context of Collective Autonomous Mobility Idris  _et al._ that will assist
with high density Urban Air Mobility Thipphavong  _et al._ .

A set of messages are to be transmitted over a communications network. The
network will have one or more senders (e.g. an airborne vehicle), one or more
receivers (e.g a ground-based facility), as well as intermediate points (e.g.
cell towers). These elements can be represented as the nodes of a graph, with
the paths between them acting as edges. Each edge will have an associated
weight, equal to the (finite, integer) time it takes a message to traverse it.
An illustrative example of this scenario is depicted in Fig. 1.

Figure 1: Information sharing illustrative example. Three messages (black
symbols) are transmitted from an airborne vehicle (upper right) to a single
ground-based facility (lower left) via two network communication towers (blue
cell tower symbols). The vehicles, ground-based facilities, and cell towers
correspond to the nodes in the graph associated with the problem, while the
green dashed lines between them are the edges, and represent the possible
paths for message transmission. Though not relevant to our problem, the solid
gray lines depict the trajectories flown by the vehicles.

Each message has an assigned path of transmission through the network— with a
sender, a receiver, and some intermediate points—as well as a scheduled
emission time. The capacity of the network is finite (and might not be the
same at every point). This, combined with the fact that transmission takes
finite time, means that sending all the messages at their scheduled emission
times might not be possible (as capacity could be exceeded at certain
locations and times). To respect this limitation, the transmission of certain
messages can be delayed with respect to the scheduled time, with the caveat
that this delay carries some cost per unit time.

Each message has associated a cost of transmission delay per unit time. It can
be understood as a representation of the message’s priority level, and the
importance of it being delivered on time. To solve the problem, each message
transmission must be assigned a time delay in such a manner that the total
cost of delay is minimized. If every message can be transmitted at its
scheduled emission time, the total cost is zero. If delays are required, there
will be some finite cost.

Certain messages may be safety and/or time critical (which can be represented
by their cost of delay being much greater than the rest). The paths for such
messages will be hard-coded into the problem setup, by reducing the network
capacity along their path, so that no delay can be assigned to them and other
messages cannot use the network bandwidth that has been pre-assigned to the
high priority ones.

The full set of problem parameters, along with the mapping to QUBO, is
provided in Appendix D.

So far, we have considered the transmission paths of the messages to be pre-
determined. However, depending on the particular communications network,
several paths could exist between the sender and receiver of any given
message, and the total cost of delay could potentially be lowered by allowing
changing of certain paths. While we focus on solving the delay optimization
aspect with the quantum annealer in this study, both delays and paths could be
optimized by implementing a hybrid algorithm, alternating classical and
quantum routines (for an example of a similar problem, in the context of air
traffic, where paths are optimized, see Sadovsky (2014)). The full algorithm
would have the following steps:

  1. 1.

Taking as inputs the set of messages, the graph representing the
communications network, and the source and destination nodes for each message,
the shortest path for each transmission is calculated using a standard
classical algorithm.

  2. 2.

Using the calculated shortest paths, the quantum annealer finds the set of
transmission delays that result in the minimum total cost of delay.

  3. 3.

The costliest delay is identified and the path of that message changed,
avoiding the location and time where the network capacity was exceeded,
prompting a delay.

  4. 4.

With this new path, the cost optimization problem is solved again in the
quantum annealer.

  5. 5.

If the total cost is reduced, the new path is accepted, and the path-changing
routine repeated with the new costliest delay. If the total cost increases, it
can be for one of two reasons. If it is due to the new path resulting in a
later arrival time, and thus a more costly delay than the one incurred by the
original path, the path is not changed; instead, the second costliest delay is
considered and the path-changing routine repeated there. If, on the other
hand, the worse outcome is due to this new path causing a previously
nonexistent delay for a different message, the change is kept and the path-
changing routine applied to the path of the message affected by the new delay.

  6. 6.

The path-changing routine is repeated until all delays have been considered
for a path change.

###  III.2 Parameters and metrics

We use the empirical probability of success ($p_{success}$) and time to
solution ($T_{S}$) as our figures of merit for determining how likely a
problem is to be solved, defined as:

$\displaystyle p_{success}=\frac{\text{{\\#} anneals with correct
solution}}{\text{total {\\#} anneals}}$ (2) $\displaystyle
T_{S}=\frac{\log(1-0.99)}{\log(1.0-p_{success})}t_{tot}\;,$ (3)

where the total time $t_{tot}=t_{a}+t_{p}$ is the time spent on each anneal,
taking into account both the annealing time $t_{a}$ and, in the case of
schedules with a pause, the pause duration $t_{p}$.

These two measures are complementary to each other. The $T_{S}$ figure of
merit reports the expected time required to solve the problem with $99\%$
confidence. While $p_{success}$ is directly determined by and hence provides a
portal to understand the underlying physical process, $T_{S}$ gives a more
practical measure that is universal across different parameter ranges and
different solvers. A higher success probability does not necessarily mean a
lower $T_{S}$. For instance, we might get a slightly higher $p_{success}$ by
using a longer annealing time $t_{a}=100\mu$s than with a shorter one
$t_{a}=1\mu$s, yet the chance of finding the solution might be higher by
repeating the $t_{a}=1\mu$s runs 100 times than by doing the $t_{a}=100\mu$s
anneal once. It is for this reason that when optimizing the different
parameters for solving a problem, we aim to minimize $T_{S}$ rather than
maximize $p_{success}$.

We only consider that the correct solution has been found when the optimal is
returned. For GC instances, this always means that the total energy returned
by D-Wave is 0, as any non-zero contributions would come from either a node
not being assigned exactly one color or from adjacent nodes being assigned the
same color. For the BD-MST and INFO problems, each instance will have a
particular minimum energy which might be non-zero, and comes entirely from the
cost function contribution, without violating any penalty terms. All these
instances are small enough to be solved either by inspection or with a
classical routine, and we verify whether the minimum energy returned by D-Wave
is in fact that of the correct solution. We also ensure that all the
contributions from the penalty terms are 0.

We discard any potential solutions returned by D-Wave with inconsistent or
broken vertex models, i.e. where not all the qubits within a given vertex
model are aligned.

Unless otherwise noted, all our results are obtained in the following manner:
for each set of annealing parameters and each problem instance, we perform 100
gauges (or partial gauges if $|J_{F}|>1$, where the gauge is only applied to
couplings $\leq 1$), with 500 anneals each, for a total of 50,000 anneals (or
reads) per run. The number of correct solutions found in those 50,000 anneals
is then divided by the total number of anneals to obtain $p_{success}$ for
each instance. Then, a bootstrap procedure is performed over the ensemble of
instances, by drawing a number of samples equal to the number of instances
(with replacement), from which a median is obtained, and repeating this
process $10^{5}$ times to finally calculate a median of medians, which is
reported as our data point. The $35^{th}$ and $65^{th}$ percentile values are
used for the error bars.

For the standard no pause schedules, we can vary the ferromagnetic coupling
$|J_{F}|$ (sometimes referred to as chain strength in the literature) as well
as the annealing time $t_{a}$. For schedules with a pause, we consider
$|J_{F}|$ and $t_{a}$, and also the pause location within the anneal, $s_{p}$,
and its duration $t_{p}$.

After mapping to QUBO, all our instances need to be minor embedded to fit the
adjacency graph of the device. We use the standard embedding heuristic
implemented by D-Wave software, which we run ten times and keep the smallest
embedding found (i.e. with the smallest number of physical qubits). Embedding
size is one of several factors that affect performance, as we discuss in Sec.
V, and selecting a smaller embedding does not necessarily guarantee a higher
probability of success. However, it is a straightforward metric which does not
require additional running time, and we have found it can improve performance
for larger and harder problems.

##  IV Results

We present our results in this section, split into the study of a single
problem on two quantum annealers (Sec. IV.1) and that of three different
problems on a single annealer (Sec. IV.2). We limit ourselves here to laying
out our results along with some brief comments, while a more in-depth
discussion of their meaning and implications is reserved for Sec. V.

###  IV.1 Comparison across devices

We study an ensemble of 45 instances of BD-MST problems on 5 nodes that are
chosen by exhausting all connected graphs with $n=|V|=5$. Weight sets are
uniformly drawn from 1 to 7. Graphs and weight sets are combined to yield
these 45 unique instances. Once mapped to QUBO following Appendix B, they have
between 32 and 74 logical qubits. All instances require embedding on both
devices. On DW2K, once embedded, the number of physical qubits is between 83
and 485, with median vertex model size between 1.5 and 7. DWA reduces the
number of physical qubits to between 38 and 188, meaning that, for the same
number of logical qubits, the physical size decreases approximately 2.5-fold,
as expected from a 2.5-fold increase in connectivity (degree 6 vs degree 15).
The median vertex model size is between 1 and 2. Results for DW2K were
obtained from our previous study Gonzalez Izquierdo  _et al._ (2021).

We start by finding the optimal ferromagnetic coupling when using a standard
annealing schedule (without a pause). We choose an annealing time of 1 $\mu$s.
This is shown in Fig. 2. The optimal ferromagnetic couplings differs
significantly between devices. DW2K does best at $|J_{F}|=1.6$, while for DWA
it is $|J_{F}|=0.8$. Given that DW2K has much larger vertex models, keeping
them from breaking requires a stronger coupling (we discard any solutions
returned with broken vertex models). With the smaller vertex models on DWA,
the exploration of configuration space provided by a weaker coupling outweighs
the higher probability of breaking vertex models, as evidenced by a 7.5x
improvement in $p_{success}$ from DW2K to DWA, when considering them at their
respective optimal $|J_{F}|$.

Optimization of $|J_{F}|$ is crucial for getting the best possible outcome. It
is highly dependent on embedding, which makes it change not only between
devices and problems, but also between instances and even between different
embeddings of the same instance, as we will discuss in Sec. IV.2. For example,
if we use the DW2K optimal value with DWA, their $p_{success}$ would appear to
be comparable. Also, because the instances in the ensemble each have their own
optimal $|J_{F}|$, while $|J_{F}|=0.8$ is best overall on DWA, three of the 45
instances do not solve for this value; one of them solves for $|J_{F}|=0.9$
and the other two for $|J_{F}|=1.1$. We will discuss this behavior in Sec. V.
On the other hand, DW2K was unable to solve 7 of the 45 instances regardless
of $|J_{F}|$.

Figure 2: Difference in $p_{success}$ and optimal $|J_{F}|$ between the two
devices. $p_{success}$ for an ensemble of 45 BD-MST instances with $n=5$,
annealing without a pause. The best performance for DW2K is seen at
$|J_{F}|$=1.6, while DWA has its optimal at $|J_{F}|=0.8$. The maximum
$p_{success}$ attained is also much higher for the newer DWA.

We also compare the case of a schedule with a pause, which, in our previous
study, we found by carefully choosing its location and duration was able to
improve not only $p_{success}$ but also $T_{S}$. The beneficial pausing region
will be affected by the difference in the $A(s)$ and $B(s)$ functions (which
would shift the region earlier) as well as by the smaller vertex models and
weaker $|J_{F}|$ (both of which would shift it later). As shown in Fig. 3, we
find that the latter effect outweighs the former, with the optimal pause
location going from $s=0.3$-$0.32$ in DW2K to $s=0.38$-$0.4$ in DWA. The
improvement with $T_{S}$ for the ensemble is about $20\%$ for both devices.
Only two data points are shown for DW2K because data for other locations was
not obtained for the full ensemble. The pause duration for DW2K is 1 $\mu$s,
which was found to be optimal in the range $[0.25,2]\mu$s (although
$t_{p}=0.5$ and 2 $\mu$s were within margin of error), and for DWA it is 0.2
$\mu$s, with no statistically significant differences in the range
$[0.15,2]\mu$s, but worse performance observed at 5 $\mu$s (no longer
improving upon no pause results).

Figure 3: Improvement of $T_{S}$ with the introduction of a pause for the two
devices. $T_{S}$ for an ensemble of 45 BD-MST instances with $n=5$, annealing
with a pause at the location indicated on the $x$ axis, and pause duration 1
$\mu$s for DW2K and 0.2 $\mu$s for DWA. In both cases $T_{S}$ is improved by
about $20\%$ when optimizing location and duration of the pause.

Finally, we include in Appendix E additional results for an ensemble of
instances on 6 nodes (instead of 5). None of them were solvable on DW2K, so
these are only for DWA.

###  IV.2 Comparison across problems

We now explore how different problems perform on the same device, for which we
choose the newer generation DWA. In this section we limit ourselves to
presenting a summary of our results, while a more thorough discussion can be
found in Sec. V.

The ensemble of BD-MST instances is the same as described earlier in Sec.
IV.1. Their physical size is between 38 and 188, with average 81 $\pm$ 22\.
Average vertex model is $1.7\pm 0.3$, and median vertex model between 1 and 2
for all instances. The average degree is $4.6\pm 0.3$

The GC instances are random 4-regular graphs on 12, 14 and 16 nodes, with 5
colors. For each of the three sizes we generate 20 different instances, and 10
more are used for $n=16$ in some of the runs when explicitly stated, for a
total of 30. The $n=12$ instances have 60 logical qubits, and their average
number of physical qubits is $166\pm 11$, with the smallest instance at 133
qubits and the largest one at 180. The average vertex model size is $2.8\pm
0.2$, and the average degree $5.1\pm 0.2$. The $n=14$ instances have 70
logical qubits, with physical sizes between 144 and 219, being on average
$193\pm 16$, and with an average vertex model size of $2.8\pm 0.2$. Average
degree is $5.2\pm 0.3$. The $n=16$ instances have 80 logical qubits, $236\pm
20$ average physical qubits, with individual instances between 162 and 269,
and vertex models having on average $3.0\pm 0.2$ qubits. Average degree is
$5.0\pm 0.3$. For all sizes, the median vertex models are between 2 and 3
qubits.

The INFO ensemble is comprised of 9 instances, with between 2 and 4 messages
being transmitted. The graph representing the possible transmission paths has
6 nodes arranged in a 3x2 lattice, with edge weights (i.e. transmission times
between nodes) being 1 or 2 time units, and capacities of 1 or 2 messages.
Paths consist of 3 or 4 edges, cost of delays is 1 to 3 units, and number of
top priority messages 0 or 1. Once mapped to QUBO they yield between 28 and 57
logical qubits, with average 44 $\pm$ 9\. Embedded size is between 29 and 87,
with average 63 $\pm$ 19\. Average vertex model 1.4 $\pm$ 0.2, and median
vertex model between 1 and 2.

####  IV.2.1 Optimizing $|J_{F}|$

We start by optimizing $|J_{F}|$ for all three problems, as shown in Fig. 4.
The value of $|J_{F}|$ influences the location of the minimum gap for the
embedded problem and subsequently affects the general location of where a
pause will help, as will be discussed in Sec. V. The optimal $|J_{F}|$ for GC
is 0.5 (0.52 for the $n=14$ instances, although not statistically different
from the 0.5 result), and it is 0.6 for the preliminary results on INFO.

Figure 4: $p_{success}$ vs $|J_{F}|$ for the three different problems run on
DWA. Data points show median bootstrapped over instances, with error bars at
$35^{th}$ and $65^{th}$ percentiles. Left: BD-MST, 45 $n=5$ instances. Middle:
GC, 3 ensembles of 20 instances each, with legend showing their respective
$n$. Right: INFO, 9 instances.

The optimization of $|J_{F}|$ presents some differences across the three
problems. Let us turn our attention to the GC results first. Here, choosing
the right value for $|J_{F}|$ provides the largest advantage, with close to 2
orders of magnitude in improvement with respect to simply setting $|J_{F}|=1$.
Although we find $|J_{F}|=0.5$ to be optimal for the ensemble, there is a
caveat. While these weak couplings allow configurations to change easily—the
benefits of which we clearly see in the form of a much higher $p_{success}$—as
the vertex model couplings get weaker, this ease of flipping also leads to a
higher likelihood of the vertex model breaking, i.e. having spins that are not
all aligned. When this happens, we discard the returned solutions. The
fraction of solutions that contain broken vertex models increases as we
decrease $|J_{F}|$ (in principle exponentially Marshall  _et al._ (2022)). A
higher fraction of solutions with broken vertex models is not necessarily an
undesirable feature—in fact, we find that up to a point, as that fraction
increases, so does the likelihood of finding the correct minimum energy
solution among the shrinking number of configurations that do not contain any
broken vertex models. But there is a point of diminishing returns, when the
fraction of solutions with broken vertex models is so large that the higher
chances of finding the correct solution among the rest cannot compensate for
their small number—this is why $p_{success}$ decreases sharply as we reach
$|J_{F}|=0.4$. And at some point, all returned solutions contain broken vertex
models, which means that the instance does not get solved, we set its
$p_{success}$ to 0 and its $T_{S}$ to infinity. (Note that, if instead of
discarding solutions with broken vertex models we chose a different approach,
such as the commonly used majority vote, $p_{success}$ in these cases might
not be 0, and it would likely be higher in all other cases as well.) We find
that the $|J_{F}|$ at which this happens stays remarkably consistent across
instances, even those of different sizes. However, we do encounter some
exceptions, and these can lead to one or several instances not being solved
when we optimize $|J_{F}|$ for the majority.

In particular, of the 20 instances for each of the three sizes that are
included in these results, one of the $n=12$ instances did not solve for
$|J_{F}|=0.5$, nor did two of the $n=14$ and two of the $n=16$ instances. Of
these, one of the $n=16$ solved for $|J_{F}|\geq 0.6$, and the other five
solved for $|J_{F}|\geq 0.7$. Despite the unsolved instances, the median
$T_{S}$ for the ensembles is several times smaller at $|J_{F}|=0.5$ than it is
at $|J_{F}|=0.7$, when all the instances solved.

Moreover, the fact that the fraction of solutions with broken models for these
particular instances reached 1 at a $|J_{F}|$ higher than for the rest does
not appear to be related to any specific characteristics of the instances. In
fact, we find that a different embedding of the same instance might show the
same behavior as the majority of instances do. To illustrate this, we take the
$n=12$ instance that did not solve for $|J_{F}|=0.5$ with its original
embedding (henceforth called instance 1 embedding 1), and compare it to
another $n=12$ instance of similar embedded size (instance 2), as well as to
an alternative embedding to the one originally used (embedding 2). This is
shown in Fig. 5. The top figure shows that the original embedding used for
instance 1 has a worse optimal $p_{success}$ and its location is shifted to a
much larger $|J_{F}|$ compared to instance 2. When we select a new embedding
for instance 1, however, its performance is much more similar to that of
instance 2 (and of the vast majority of the instances considered). In
particular, the maximum $p_{success}$ occurs at $|J_{F}|=0.52$ for both
instances in that case. The bottom plot shows the fraction of solutions that
had broken vertex models and were thus discarded. For instance 1 with the
original embedding, we see that the fraction of solutions with broken vertex
models saturates at a much larger $|J_{F}|$ than for the other two cases. This
shift of the curve is similar to the one we observe in the $p_{success}$
results. The location of the optimal $p_{success}$ occurs at a similar
fraction in all cases, when about $60\%$ of solutions returned have broken
vertex models.

Figure 5: Optimal $|J_{F}|$ for two individual graph coloring instances (DWA).
Top: $p_{success}$ for the two individual $n=12$ instances, as we vary
$|J_{F}|$. Bottom: Fraction of solutions returned by D-Wave where at least one
vertex model was broken, and thus the solution discarded. A value of 1
indicates that all solutions found had broken vertex models.

There are no obvious differences in the structure of the embedding graphs that
would have led us to predict the disparities in behavior. Their numbers of
physical qubits and edges are similar (with embedding 2 slightly higher on
both counts, showcasing how choosing a smaller embedding does not always lead
to better performance), and so are their mean and median vertex model sizes.
It is not until we look at the embedded coefficients that we find
discrepancies that can explain the different behaviors of these cases. We
consider the ratio of maximum to minimum embedded coupling
$R_{J}=\max(|J_{ij}|)/\min(|J_{ij}|)$ and that of maximum to minimum
individual biases $R_{h}=\max(|h_{i}|)/\min(|h_{i}|)$, and find a clear
difference between their values for most of the embeddings (the ones with
optimal $T_{S}$ at or near $|J_{F}|=0.5$) and the few that require much
stronger $|J_{F}|$ to be solved. For $n=12$, all embeddings have
$R_{J}\in[15,30]$ and $R_{h}\in[3,5]$ except for the one that did not solve at
$|J_{F}|=0.5$, which has $R_{J}=7.5$ and $R_{h}=2$. The second embedding for
the same instance, which as described above showed a behavior consistent with
the rest, had $R_{J}=20$ and $R_{h}=5$. Similarly for $n=14$, all embeddings
have $R_{J}\in[15,25]$ and $R_{h}\in[3.33,5]$, except for the two that did not
solve at $|J_{F}|=0.5$, which had $R_{J}=10$ and $R_{h}=2.5$. Finally, for
$n=16$, all embeddings have $R_{J}\in[15,35]$ and $R_{h}\in[3.33,6]$, except
for the two that did not solve at $|J_{F}|=0.5$. Of those, one had $R_{J}=7.5$
and $R_{h}=2.5$, while the other one had $R_{J}=15$ and $R_{h}=3$ (note that
in the rest of the ensemble, although we find an embedding with $R_{h}=3.33$
and a few with $R_{J}=15$, these values do not occur together. Also, this
instance solved for $|J_{F}|=0.6$ while the other 4 did not).

We can see the reason that these five instances with these particular
embeddings did not solve for $|J_{F}|=0.5$ is that $|J_{F}|$ is relatively
lower for them, given that the rest of the coefficients are lower. It then
makes sense that the $p_{success}$ vs $|J_{F}|$ curve (as well as the fraction
of broken chains vs $T_{S}$ curve) are shifted to the right compared to the
rest of the instances/embeddings. Although we have left our results to include
these outlier embeddings (given that we did not explore their differences
until after performing our runs, and had initially decided that our method for
choosing embeddings would be smallest out of 10) obtaining the $R_{J}$ and
$R_{h}$ before choosing an embedding is an easy check that does not
significantly increase the time resources needed to solve the problem, and can
in fact give us a general idea, a priori, of what the optimal $|J_{F}|$ will
be.

The situation is quite different for the BD-MST instances. In that case, the
optimal $|J_{F}|$ is 0.8, but only a few instances have that value as their
individual optimal. Instead, there is a relatively even spread of optimal
$|J_{F}|$ values between 0.6 and 1.3, with the most common being 1.0 at 10
occurrences. There is one instance that does not get solved at $|J_{F}|=0.8$,
and requires $|J_{F}|\geq 1.1$ to find a valid solution, and we find that
harder instances (i.e. those with a lower $p_{success}$, considered at their
optimal $|J_{F}|$), tend to have higher optimal $|J_{F}|$. This is in contrast
with the GC case, where the instance-wise optimal $|J_{F}|$ remains unchanged
(save for the few instances discussed above) regardless of hardness.

In the INFO case, while the optimal for the ensemble is 0.6, 7 out of the 9
instances have their optimal at 0.5, one has it at 0.6, and the last one at
0.9 (this one is the hardest instance, in that its $p_{success}$ at optimal
$|J_{F}|$ is the smallest out of all instances). Given that in most cases 0.5
does better than 0.6, we choose this value for our subsequent runs.

####  IV.2.2 Optimizing $t_{a}$

Setting $|J_{F}|$ to its optimal for each ensemble (i.e. the same for all
instances), we explore a range of annealing times. Generally, the shortest
annealing time is optimal in terms of $T_{S}$, and it is likely that these
short annealing times will be particularly beneficial in combination with a
pause as discussed in Sec. V. When optimizing $|J_{F}|$, we set $t_{a}=1\mu$s,
which is the shortest allowed by the annealer. We chose this value based on
previous work that found it to be optimal in terms of $T_{S}$. We expect that
increasing $t_{a}$ leads to higher $p_{success}$, by giving the system a
longer time to thermalize. However, it is in many cases found that the
increase in $p_{success}$ is not sufficient to compensate the longer time
required to arrive at a solution. That is, given a set amount of time, it is
more beneficial to perform many short anneals than few longer ones.

Figure 6: $T_{S}$ vs $t_{a}$ for the three different problems run on DWA. Data
points show median bootstrapped over instances, with error bars at $35^{th}$
and $65^{th}$ percentiles. Left: BD-MST, 45 $n=5$ instances. Middle: GC, 3
ensembles of 20 instances each, with legend showing their respective $n$.
Right: INFO, 9 instances.

As seen in Fig. 6, the optimal $t_{a}$ is $\leq 1\mu$s in most cases (we
cannot probe the range $t_{a}<1\mu$s so it is not possible to pinpoint the
actual value). For the BD-MST ensemble there is a local minimum at
$t_{a}=4\mu$s with $T_{S}$ almost as good as that obtained at 1 $\mu$s. This
effect is consistent even if we look at the instances individually; the low
$T_{S}$ at 4 $\mu$s is not a result of some instances preferring a shorter
$t_{a}$ and others a longer one, but instead of many of them having a curve
similar to that of the ensemble.

For GC, the three ensembles separated by size tell a more complete story. The
best $T_{S}$ in the accessible range for $n=12$ is clearly at $t_{a}=1\mu$s,
with longer times doing progressively worse. This holds true for $n=14$ as
well, but a local minimum at 3 $\mu$s appears. And finally, at $n=16$, the
minimum at 3 $\mu$s is the global one, with another one at 7 $\mu$s a close
second. Both yield clearly better results than 1 $\mu$s, although it is
possible that the true global minimum could be $<1\mu$s. Because these results
were somewhat noisy, we increased the number of anneals and instances to get
more reliable data. For all $n$, the data points for $t_{a}=1$, 2 and 3 $\mu$s
are obtained with 500 gauges instead of our standard 100, for a total of
250,000 anneals. In the case of $n=16$, 500 gauges were also used for the
$t_{a}=4\mu$s point, and all the data points used 30 instances instead of the
usual 20. Interestingly, with the introduction of a short pause
($t_{p}=0.2\mu$s) at a good location ($s_{p}=0.4$, see Fig. 8 for results
through a range of $s_{p}$), $t_{a}=1\mu$s does better than either 3 or 7
$\mu$s for $n=16$ as well (although $T_{S}$ is improved by the pause for all
$t_{a}$ tested). This is shown in Fig. 7.

Figure 7: Optimal $t_{a}$ for graph coloring instances with a pause (DWA).
$T_{S}$ vs $t_{a}$ for the ensemble of $n=16$ graph coloring instances.
Without a pause, the optimal in the accessible range is $t_{a}=3\mu$s. With
the introduction of a 0.2 $\mu$s pause at $s_{p}=0.4$, $T_{S}$ improves at all
$t_{a}$, but now $t_{a}=1\mu$s does best within the range.

The optimal $t_{a}$ is also $\leq 1\mu$s for the INFO ensemble. While error
bars are very large, that is a product of there being only 9 instances of
disparate hardness, rather than lack of consistency in the results. 8 of the 9
instances have 1 $\mu$s as their optimal within the device’s range, with the
$9^{\text{th}}$ one doing better at 2 $\mu$s. For 6 out of the 9, $T_{S}$
monotonically increases with $t_{a}$.

####  IV.2.3 Improving $T_{S}$ with a pause

We test whether an appropriately located pause is able to improve $T_{S}$.
Given previous theoretical Albash and Marshall (2021b); Chen and Lidar (2020)
and experimental Marshall  _et al._ (2019); Gonzalez Izquierdo  _et al._
(2021) knowledge about the range of locations and durations where a pause is
beneficial, we explore those regions and discuss the results in Sec. V. We
have fixed $|J_{F}|$ and $t_{a}$ at the optimal found for each ensemble (in
the case of $n=16$ GC instances, we chose $t_{a}=1\mu$s rather than 3 $\mu$s
since we found the former to be better with a pause).

Figure 8: $T_{S}$ vs $s_{p}$ for the three different problems run on DWA. Data
points show median bootstrapped over instances for a pause of $t_{p}=0.2\mu$s
at the location indicated on the x axis. Horizontal lines correspond to
results without a pause. Error bars at $35^{th}$ and $65^{th}$ percentiles.
Left: BD-MST, 45 $n=5$ instances. Middle: GC, 3 ensembles of 20 instances
each, with legend showing their respective $n$. Right: INFO, 9 instances.

Fig 8 shows $T_{S}$ at a range of $s_{p}$ between 0.2 and 0.6. We have chosen
$t_{p}=0.2$, but have found no significant differences in performance in the
range $t_{p}\leq 0.5\mu$s for GC (shown in Fig. 9), and $t_{p}\leq 1\mu$s for
BD-MST. INFO instances were run with $t_{p}=0.2\mu$s, 1 $\mu$s and 2 $\mu$s,
and the shortest time was best, although, like all other results, not by a
statistically significant amount. Given the large error bars and the results
from the other two problems, we did not investigate a range of very short
$t_{p}$. We are able to improve $T_{S}$ by adding a pause for all ensembles.
The results for the INFO case are not statistically significant, due to having
too few instances with a wide range of $T_{S}$. Of the 9 instances, 8 improve
with the pause. We also repeat the runs (for the $s_{p}$ values that provide
an improvement) for $|J_{F}|=0.6$, which had the best median $p_{success}$ for
the ensemble, although it did worse than $|J_{F}|=0.5$ for 8 of the 9
instances. The medians for both $|J-F|$ values are comparable, with
$|J_{F}|=0.5$ slightly better and, its error bars considerably lower.

The BD-MST ensemble experiences the smallest improvement, of about 20$\%$,
while the GC is closer to 40 $\%$ for $n=12$, with $n=14$ around 60 $\%$ and
$n=16$ near 70 $\%$. INFO improves by almost 50 $\%$ when comparing to the
$|J_{F}|=0.6$ results, or by around 65 $\%$ compared to those with
$|J_{F}|=0.5$. The best pause location for the BD-MST ensemble was found to be
$s_{p}=0.4$, although the slightly worse data point at $s_{p}=0.38$ has much
smaller error bars. In fact, the results are too noisy to pinpoint an exact
optimal location; but we can say the region $s_{p}=0.38$ \- 0.42 does better
than earlier or later points. Similarly, all GC ensembles perform best with a
pause between 0.38 and 0.4, and the INFO between 0.4 and 0.44. The pause
becomes more beneficial as these instances increase in size and hardness.

To compare the performance across different pause durations $t_{p}$ we first
run each $t_{p}$ for a range of $s_{p}$, because the optimal location changes
for different pause durations, as we discuss in Sec. V. For clarity, we
present these results for GC in Fig. 9 by only showing the best $T_{S}$ found
for each of the pause durations $t_{p}$ (on the $x$ axis), regardless of at
what $s_{p}$ that $T_{S}$ was found. Our data indicates that no major
differences in $T_{S}$ are found in the range $t_{p}=0.1$-0.5$\mu$s, while
once we get to $t_{p}=0.75$-1$\mu$s the $T_{S}$ increases, and is clearly much
worse as we lengthen the pause to 10 $\mu$s.

Figure 9: Best $T_{S}$ found for each $t_{p}$ (DWA). Optimal $T_{S}$ found for
each $t_{p}$ for the three ensembles of 20 instances at different sizes.
$|J_{F}|=0.5$ and $t_{a}=1\mu$s always used. Note that the best $T_{S}$ at
each $t_{p}$ might be found at a different $s_{p}$.

As we increase the pause duration, the optimal location shifts later in the
anneal. This was first observed in Marshall  _et al._ (2019) (see Fig. 7
therein). We show the same effect, although for shorter pauses, in our Fig.
10. In this case, the median shown in each data point is calculated over the
full ensemble of 60 instances, taking all three sizes together, as is the no
pause result. The peak in $p_{success}$ clearly shifts to a later $s_{p}$ as
the pause duration increases by factors of 10. Note that although the pause
durations in Marshall  _et al._ (2019) were longer, we maintain the relative
difference between them, always a factor of 10. The shift in $s_{p}$ is of
similar magnitude in both cases.

We note that when performing the runs to obtain these results, certain
discrepancies were observed if the runs for different data points (e.g. with
pauses at different $s_{p}$) were performed on different days. In those cases
some results would appear more noisy. To avoid this, we chose to perform runs
for the same $t_{p}$ over a range of $s_{p}$ one right after another for
consistency. It is likely that significantly increasing the number of
instances and reads would also help smooth certain results and get rid of
“noise” due to the particulars of our ensembles, or other random factors.
However, given time limitations and the additional time resources that would
be required, we preferred to explore a wider range of parameters with fewer
instances, this way obtaining a more general picture of parameter setting.

Figure 10: Shift of optimal pause location with pause duration (DWA). Peak in
$p_{success}$ for three different pause durations $t_{p}$, for an ensemble of
60 graph coloring instances (the three ensembles of different sizes taken
together). The peak shifts later as $t_{p}$ increases.

##  V Discussion

In this section, we discuss the significance of our results and some factors
that help explain them. We first review the physical picture behind pausing,
and how our results agree with its predictions. Then we discuss how certain
characteristics of the logical and embedded problems, such as coefficient
heterogeneity and size, affect these predictions. Next, we focus more
specifically on the role of these effects when optimizing $|J_{F}|$, followed
by a discussion of optimal annealing time and some comments on the shift of
the $p_{success}$ peak location with pause duration. Finally, we provide
practical guidelines for parameter setting.

Physical picture. Let us start by reviewing the physical picture behind
pausing Marshall  _et al._ (2019), and the particular considerations when a
pause is applied to embedded problems Gonzalez Izquierdo  _et al._ (2021).

Figure 11: Diagram of the different annealing regimes. Regimes I and III (in
purple) correspond to the regions where only one Hamiltonian dominates. Regime
2, where the thermal and quantum scales are comparable, is subdivided into
three regions according to the relation between the instantaneous gap $\Delta$
and the temperature $T$.

Fig. 11 shows a cartoon diagram to understand the dynamics at different times
in the anneal. For a more detailed explanation, the reader is directed to
Sec.V of Gonzalez Izquierdo  _et al._ (2021), which we briefly summarize here:
Regime I corresponds to the earlier part of the anneal, when the driver
Hamiltonian dominates and the system stays in its GS. Regime II encompasses
the region where quantum and thermal scales are comparable. Thermal and
quantum dynamics take place, both of which can transfer population from the GS
to excited states. The minimum gap is in this region, and its relative size
with respect to temperature distinguishes three sections with different
behaviors: in section a, the gap starts out larger than the temperature and
decreases approaching it. Non-adiabaticity starts to occur, although
transitions (as well as thermal ones) are still slow compared to the anneal
evolution. In section b, the gap is smaller than the temperature and reaches
its minimum. Thermalization dominates and causes population transfer from the
GS to excited states, although quantum non-adiabatic transitions also take
place and contribute to the GS population loss. A long enough pause can allow
the system to approach thermal equilibrium. In section c the gap is again
larger than the temperature, but the quantum effects have now become weaker.
The system is unlikely to follow the thermal state unless given extra time. A
pause at the start of this section could help bring back some of the GS
population lost to thermalization in the previous one. Finally, during regime
III dynamics are essentially frozen, since the quantum scale has now become
too small.

Based on this picture, we have some information a priori about the general
location in which a pause will help. We do not know, however, exactly where
the minimum gap is. To further complicate matters, the minimum gap of the
embedded problem will not be in the same location as that of the logical one.
If instead of embedding, we increased the overall magnitude of the problem
Hamiltonian $H_{p}$, this would amount to increasing $B(s)$ with respect to
$A(s)$, hence shifting the different regions of the dynamics earlier in the
anneal, including the minimum gap. A rigorous proof of this effect can be
found in Choi (2020). What happens when we embed is not exactly the same;
additional $|J_{F}|$ couplings are added to $H_{p}$, rather than it being
simply scaled up. While no formal proof exists for this specific case that the
shift of the minimum gap still occurs, we see empirically that the optimal
pause location moves earlier, which points to the minimum gap also being
earlier. We also provide numerical evidence of the effect for a toy problem in
Sec.V.B of Gonzalez Izquierdo  _et al._ (2021). In this scenario, two factors
of the embedding would affect the shift; 1) the size of vertex models and 2)
the value of $|J_{F}|$. Larger vertex models mean more added couplings,
further adding to $H_{p}$, as would higher $|J_{F}|$ values. Moreover, these
two factors typically happen in conjunction; larger vertex models are more
likely to break since there are more moving parts, and thus usually require
stronger $|J_{F}|$ to stay intact (although this is not the case under all
circumstances, as we will discuss later). The shift is seen in our results for
the BD-MST ensemble across the two devices; the optimal pause location for
DW2K (when embeddings were large and $|J_{F}|$ strong) was found to be
$s_{p}=0.3$ \- 0.32, much earlier than was consistently found for native
problems in Marshall  _et al._ (2019) (around $s=0.5$). Vertex model sizes
were greatly reduced for DWA, and the optimal $|J_{F}|$ is much weaker. As
seen in Fig. 3, the optimal pausing region for DWA ($s_{p}=0.38$ \- 0.42) is
later than for DW2K, even though the fact that $A(s)$ decays faster for DWA
and crosses $B(s)$ earlier would shift it in the opposite direction. DWA’s
optimal pause region is still earlier than that of native problems. This all
fits well with the picture of embedding shifting the minimum gap earlier.

Problem differences: general considerations. There are significant differences
in behavior across problems, which can be related to certain characteristics
of both their QUBO mappings and embeddings. Factors that determine problem
hardness include size, density, and heterogeneity of coefficients, i.e. how
many different values of $J_{ij}$ and $h_{i}$ are present, and how close those
values are to one another. The devices we use suffer a range of integrated
control errors (ICE) dwa , due to which the coefficients implemented in
practice will differ from the user inputted ones by some $\delta J_{ij}$ and
$\delta h_{i}$. Precision is in this way limited, and having too many
different coefficients and/or some that are too close together can result in
them getting mixed up once we account for $\delta J_{ij}$ and $\delta h_{i}$.
The number of different coefficients increases after embedding and their
values change, as they get split over several qubits and couplers. Large
vertex models are detrimental since they can lead to too small coefficients,
and so is having several vertex model sizes, which add more different
coefficients, potentially closer to one another. In the standard embedding
algorithm included in the D-Wave software, the $h_{i}$ corresponding to a
single variable is evenly distributed over the physical qubits of the vertex
model (between 1 and 6 for all the problems presented in the main text,
although the largest instances of $n=6$ BD-MST in Appendix E have maximum
vertex model size going up to 10), and the $J_{ij}$ for a logical coupling
also gets evenly split over all the available physical couplers between qubits
in each vertex model (between 1 and 8 couplers for the problems in the main
text, going up to 9 for the largest $n=6$ BD-MST instances).

A more dense problem is more likely to require a larger embedding, and larger
vertex models will be needed to accommodate the extra connectivity. The more
dense embedded problem will lead to frustrations and a higher rate of broken
vertex models and discarded solutions. Larger and denser problems are in
principle more difficult, but experimentally we find this secondary to the
characteristics of the coefficients. This becomes clear when we consider the
differences in $p_{success}$ or $T_{S}$ across ensembles. For example, as we
can see in Fig. 4, the three GC ensembles have $p_{success}$ highest for the
smallest ensemble ($n=12$), subsequently decreasing for the $n=14$ and $n=16$
cases. However, all three GC ensembles have higher $p_{success}$ than the BD-
MST and INFO ones, even though they are the largest in size. The embedded size
of the smallest GC instances is similar to the larger BD-MST ones, and the
$n=14$ and $n=16$ upper range goes well beyond that of the BD-MST. This is
even more pronounced for the INFO case, which has the smallest instances, and
yet its $p_{success}$ is the lowest of the three problems. GC instances are
typically denser than BD-MST ones, with average degree $5.1\pm 0.3$ and little
variability for different $n$, while BD-MST has average degree $4.6\pm 0.3$.
INFO instances have average degree $4.6\pm 1$.

Coefficient heterogeneity, hardness and optimal ferromagnetic coupling. Based
on our results, we present the hypothesis that the main differentiating
feature across problem classes is coefficient heterogeneity. The coefficients
of GC instances are remarkably regular, and regardless of the problem
specifics they will not have very disparate values. The GC QUBO mapping
(Appendix C) does not include an objective function, but only two penalty
terms, which means that, unlike the other problems, it does not require a
penalty weight factor acting on the penalty terms to ensure that violating one
does not become advantageous. Thus, all coefficients of the GC QUBO are either
1 or 2. Once embedded, there is a wider range of values as the $h_{i}$ and
$J_{ij}$ get divided and assigned by the embedding heuristic, which leads to a
small set of coefficients that the majority of instances share. We do find a
discrepancy, discussed in Sec. IV.2.1 between the 93$\%$ of instances with
optimal $|J_{F}|=0.5$ and the rest with a different optimal $|J_{F}|$; this
difference can be traced to the fact that the minimum vertex model size in the
majority is one, while in the cases with a higher optimal $|J_{F}|$ it is two.
Yet, the physical qubits in those smallest vertex models are assigned the same
$h_{i}$ value whether there is a single physical qubit or two, resulting in
the sum of all the biases corresponding to a given vertex model, or all the
physical couplings representing a single logical one to be twice as large in
the few instances that behave differently, thus needing a roughly twice as
strong $|J_{F}|$ to obtain the best performance. This is a detail that one
should keep in mind when using this embedding routine without further
examination, as it can lead to unexpected behavior like in this case.

The level of coefficient heterogeneity can explain the greater impact that
optimizing $|J_{F}|$ has on the GC ensembles compared to the other problems.
All the individual curves (save for the few exceptions) look the same as that
of the ensemble. On the other hand, for BD-MST we have many different
individual curves, with maxima across a range of $|J_{F}|$ values, resulting
in a much less sharp peak when averaged. For INFO, although the optimal is
fairly consistent, the small number of instances and wide range of
$p_{success}$ lead to very large error bars, somewhat diluting the results.

Compare the QUBO formulation for GC with those of BD-MST or INFO. In the BD-
MST QUBO (Appendix B), we have a cost function that has each of the different
weights as coefficients. Then, the penalty terms can have coefficients 1, 2
and 4, multiplied by the penalty weight, which needs to be larger than the
maximum weight. Before embedding, we are already at quite a few different
values, depending on the number of weights. For this reason, we avoid
instances with weights that are too disparate, which will quickly run into
precision issues. Intuitively, those instances would seem easier to solve,
since it would be our first instinct to avoid the more costly parts of the
graph, while for the annealer, these are more difficult. For INFO (Appendix
D), this can become even more pronounced. The coefficients for the cost
function run through the product of each individual cost with the possible
times at which the message will be traveling, so we have two parameters coming
into play. This can result in a wide array of values, which we curb by setting
a short time horizon and keeping all other problem parameters within narrow
ranges (e.g. 2-4 messages, delay costs of 1-3, travel time for each edge 1-2).
Despite this similarity among instances, and their yielding the smallest
embedded sizes (compared to the other two problems), their median difficulty
is the highest out of the three ensembles, and their disparity even larger,
including the instance with the highest $p_{success}$ as well as some of the
most difficult ones. The penalty terms coefficients can be 1 or 2, and are
also multiplied by a penalty weight, which depends on costs and times as
described in Appendix D. Just like for BD-MST, costs that are too disparate
can take us beyond precision (and so do long times which can be a result of
long travel times in the connectivity graph, having long paths, and/or not
setting a tight time horizon bound), leading to problems that are difficult
for the annealer but would be intuitively easy to solve, by not delaying the
message(s) with the much higher cost(s) than the rest.

This picture correlates well with what we see in our instances. GC instances
have (once embedded) between 6 and 10 different $J_{ij}$ values, and 3 to 6
different $h_{i}$ ones. INFO instances have 4-11 for $J_{ij}$ and 11-30 for
$h_{i}$, and BD-MST instances have 7-14 for $J_{ij}$ and 12-27 for $h_{i}$.

Relationship between vertex model size and optimal ferromagnetic coupling. Our
results on the differences in optimal $|J_{F}|$ challenge the conventional
wisdom that smaller vertex models are less likely to break and thus will do
better with a weaker $|J_{F}|$. While this is the correct approach in some
cases (e.g. solving the BD-MST ensemble in DWA had smaller vertex models and a
much lower optimal $|J_{F}|$ compared to DW2K), things are more subtle when
other differences are present. We have seen how most of the GC instances have
a clear optimal at $|J_{F}|=0.5$, while BD-MST instances have a range of
optimal values between 0.6 and 1.3. Yet, the average vertex model size for GC
instances is $2.8\pm 0.2$ for both $n=12$ and 14, and $3.0\pm 0.2$ for $n=16$,
while for BD-MST it is significantly smaller at $1.7\pm 0.3$. We hypothesize
that the characteristics of the rest of the coefficients seem to have a
greater impact on optimal $|J_{F}|$ than simply the size of the vertex models.
Density could also play a complicated role. Similar to what happens with size,
a denser vertex model has more competing forces that can lead to breaking,
which would point towards a higher $|J_{F}|$ being needed to keep it
consistent. Yet, that density could also require more ease to change
configurations among the increased number of possible ones, thus benefiting
from weaker $|J_{F}|$. We see in Fig. 5 that the optimal $p_{success}$ is
reached when 60$\%$ of solutions returned are broken, and that number stays
consistent between the two different embeddings, despite their disparate
coefficients. When considering instances of the same BD-MST ensemble, where
there are some commonalities among coefficients, although not to the degree of
the GC case, the smaller vertex models with weaker $|J_{F}|$ correlation does
apply. Within the BD-MST ensemble, we find that larger instance size
correlates with larger vertex models, stronger optimal $|J_{F}|$, and higher
$T_{S}$.

A result that embodies these last two points is that the optimal $|J_{F}|$
remains constant through changes in hardness for GC problems, while for the
BD-MST instances it increases with hardness. This points to the fact that,
when hardness comes mostly from overall size, it does not affect optimal
$|J_{F}|$, while when it is due to the other factors a stronger $|J_{F}|$ is
required to keep up.

A more quantitative analysis to tease out the effects of each of these aspects
will be necessary to inform parameter setting across different problems and
platforms, which becomes critical as difficulty increases, as choosing the
right parameters can determine whether a valid solution is found at all.

Optimal annealing time. We see in Sec. IV.2.2 that in most cases the shortest
annealing time is optimal in terms of $T_{S}$. The exception is the GC $n=16$
ensemble, which clearly performs better at slightly longer $t_{a}$. These
instances have the largest sizes out of all our ensembles. Although D-Wave
quantum annealers can in principle be operated with shorter than 1 $\mu$s
times, this capability has not yet been offered to the public (but it has
recently been announced that shorter times of up to 0.5 $\mu$s will be made
available in the near future). It is likely that these shorter times would be
particularly beneficial in combination with a pause. When we introduce a short
pause of 0.2 $\mu$s at an appropriate location for the GC $n=16$ ensemble, the
optimal $t_{a}$ is no longer 3 $\mu$s like in the no pause case, but instead
$\leq$ 1 $\mu$s like for the other ensembles, as shown in Fig. 7. This tells
us that 0.2 $\mu$s of extra time at a well chosen location is more beneficial
than up to 9 $\mu$s of extra time spread equally over the whole anneal (since
in the no pause case, while the optimal was at 3 $\mu$s, several $t_{a}$ up to
10 $\mu$s had lower $T_{S}$ than 1 $\mu$s). This is also seen in Fig. 8, where
a pause too early or too late leads to a higher $T_{S}$. There are regions in
the anneal where going faster is unlikely to have a detrimental effect, such
as when dynamics become very slow, and that time can be better used elsewhere
to improve $T_{S}$. It is also unclear whether we have found the true optimal
$t_{a}$ at 1 $\mu$s for our ensembles, or lowering $t_{a}$ would further
decrease $T_{S}$. Access to a lower range of $t_{a}$ will help answer these
questions.

Peak shift with pause duration. We have also verified, for the ensemble of GC
instances across all three sizes, the $p_{success}$ peak shift with $t_{p}$
observed for native problems in Marshall  _et al._ (2019) (see Fig. 10). The
reason for this shift is that, to maximize our chances of obtaining the GS of
the problem Hamiltonian $H_{p}$, it would be advantageous to have the system
thermalize as late as possible in the anneal; ideally at $s=1$ when the
instantaneous Hamiltonian is $H_{p}$, so the instantaneous GS is that of
$H_{p}$. In practice, the system is not able to thermalize too late in the
anneal, after the freeze-out point Amin (2015); Marshall  _et al._ (2017), and
even before reaching this point it becomes progressively more difficult as
quantum fluctuations decrease. Introducing a pause is a way of making
thermalization more likely to happen. As we move later in the anneal, a longer
pause is required to give the system enough time to thermalize, but because
the instantaneous GS is becoming closer to that of $H_{p}$, the benefit is
larger. This is why we also see $p_{success}$ increasing as the peak shifts
later with the lengthening of $t_{p}$. However, as we have seen in Fig. 9,
that increase in $p_{success}$ is not sufficient to compensate the additional
time, and $T_{S}$ worsens for longer pauses.

Practical recommendations. We have devoted this section to examining how the
different annealing parameters relate to one another, depend on problem
characteristics and affect $T_{S}$. But a reader hoping to attack a new
problem on a quantum annealer might still be left wondering what parameters to
choose, especially if time is a constraint and exploring a range of parameters
not possible. Here, we take all our findings into consideration to provide a
set of guidelines that can be applied in practice. Given that problems vary
widely, these guidelines should be understood as a starting point and a means
to obtaining good results without the overhead of benchmarking, but cannot be
expected to provide optimal performance without additional tweaking and
exploration.

Let us start with setting $|J_{F}|$. This parameter has the largest effect on
$p_{success}$ and does not increase running time, so if additional resources
are available, it is a good idea to use them here and test more than one
value. First, consider the QUBO formulations of the set of instances to be
solved. Do they all have very similar coefficients, or is there a lot of
variety? If the former is true, they will likely all have the same (or very
similar) optimal $|J_{F}|$, which simplifies things. The best course of action
is then to explore a range of values for a single instance, identifying its
optimal value, and running the rest with that value only.

If on the other hand the coefficients vary a lot across instances, optimal
$|J_{F}|$ will too. In this case, one must decide how much extra work to put
into choosing different $|J_{F}|$ for each instance. Assigning the same to all
is the simplest and quickest way, but can lead to poor results for those at
the tail ends, in particular those with very high optimal $|J_{F}|$, for which
a value in the middle can mean that no valid solutions are returned due to all
of them having broken vertex models. In this case, because the strength of the
optimal $|J_{F}|$ correlates with problem size, vertex model size and number
of different coefficients, it can be worth splitting the ensemble into
subgroups according to one of those quantities and assigning a different
$|J_{F}|$ to each of them (logical size is a good straightforward choice,
although if there are very large differences in vertex degree across
instances, embedded size will be a better predictor). A range between 0.5 and
1.5 can be a good start. Using a majority vote approach to process solutions
with broken vertex models, rather than discarding them like we do here can
also be helpful.

It is not possible to determine whether $|J_{F}|$ is close to its optimum by
looking at the results for a single value. We can, however, get a rough idea
of whether we are in the right region by looking at the fraction of solutions
with broken vertex models, which is a simple check. A majority of the returned
solutions will contain broken vertex models when performance is best, but not
such a large fraction that obtaining valid solutions becomes almost
impossible. After testing a $|J_{F}|$ value, if the fraction is close to
100$\%$, a stronger $|J_{F}|$ should be used. If it is small, a weaker
$|J_{F}|$ is necessary.

For the annealing time $t_{a}$, we have seen that keeping it at the current
shortest possible (1 $\mu$s) combined with a pause yields the best results. We
suspect that this will hold even as shorter times become available, but there
will be a (currently unknown) limit where results worsen with shorter times,
so one should be careful if $t_{a}\ll 1\mu$s becomes available. The pause
duration $t_{p}$ should also be short. We have not found measurable
differences in the range $t_{p}\leq 0.5\mu$s, and consider $t_{p}=0.2\mu$s a
safe bet.

A pause location $s_{p}$ within the region 0.3-0.5 will generally improve
performance. The smaller the vertex models, the later the optimal pause
location and vice versa. For instance, for an embedding that only requires a
few extra qubits, and where many of the logical variables remain unembedded, a
pause close to 0.5 will likely be best. For the instances in this paper, where
vertex models are typically 1-5 qubits (when using DWA), we find the optimal
location around 0.4. The results for the older DW2K, with larger vertex
models, did best closer to 0.3. For really large vertex models this location
will probably have to be pushed even earlier (but we do not currently have the
experimental data to confirm it). If the $A(s)$ and $B(s)$ functions are
significantly different from the ones used in this work, and in particular the
region where their scales are comparable is shifted, the pause location might
need to be adjusted in the same direction as this shift.

Finally, one should keep in mind that time and other resources required to
find optimal (or even simply good) parameters must be taken into account. When
we report $T_{S}$ for a set of parameter values, the time needed to find them
is not included. But that time will be a limiting factor for practical
applications, which was considered to provide the above guidelines.

A note on our optimization strategy. In our benchmarking study, we have mainly
followed the strategy of optimizing a single parameter at a time while keeping
the rest fixed, and thus there are large regions of the multidimensional
parameter space that have been left unexplored. Although this means that we
cannot provide experimental evidence that the optimal set of parameters
resides in the region that was investigated, we believe this to be the case,
and our method to be a reasonable choice. First, it would be not just
impractical, but virtually impossible to obtain results covering the entirety
of this space, so making some assumptions is unavoidable. With this in mind,
we have based our experimental runs on the physical picture described at the
beginning of this section, whose predictions have matched our observations
well. We have then made adjustments based on the trends that the first rounds
of experimental runs revealed and, when possible, have covered a large range
of values to ensure the observed trends are not local to a small region.
Finally, in the cases where there was a physical reason for a parameter to
affect another optimal value, this was taken into account—for instance, after
finding an optimal $t_{a}>1\mu$s for certain instances with an initial no
pause schedule, several $t_{a}$ values were explored again after the
introduction of a pause.

##  VI Conclusions

We experimentally investigated the physical picture of pausing on quantum
annealers, in combination with the parameter setting problem, across multiple
devices and non-native problems. We gained insights regarding the
characteristics of the logical problem—and their transformation through
embedding—that impact the hardness of the problems, as well as the parameters
that optimize their time to solution $T_{S}$.

Pausing mid-anneal was first proposed as a strategy to improve the probability
of success $p_{success}$ for native optimization problems Marshall  _et al._
(2019), and was later confirmed to work for embedded problems Gonzalez
Izquierdo  _et al._ (2021), and to also improve $T_{S}$ with a careful choice
of pause duration Albash and Marshall (2021a). We studied an ensemble of
problem instances on two quantum annealing devices with different
architectures and annealing schedules, and also three ensembles of different
problems on a single device, leading to further confirmation that pausing can
improve $T_{S}$ in a variety of scenarios. We found the region of improvement
stays fairly robust, barring small adjustments due to large differences in
vertex model size and/or its corresponding ferromagnetic coupling $|J_{F}|$.
Smaller vertex models and smaller $|J_{F}|$ shift the location of the region
later, and vice versa.

We provide evidence that adding a short pause at a well chosen location boosts
$T_{S}$ more than a much longer annealing time overall does. Combined with the
fact that a pause outside of the beneficial region worsens $T_{S}$, this gives
experimental support to the idea that there exist rather large sections at
both ends of the anneal where annealing faster will not hinder performance (up
to a point that cannot currently be determined due to hardware limitations),
and that pauses will likely improve $T_{S}$ even when the optimal annealing
time $t_{a}$ has been reached.

Along with the pause location $s_{p}$, we also found the optimal annealing
time $t_{a}$ and pause duration $t_{p}$ to remain consistent across the
different scenarios, while the ferromagnetic coupling within vertex models,
$|J_{F}|$, requires careful examination to be chosen. Nonetheless, optimizing
$|J_{F}|$ had a much larger effect on performance than any of the other
parameters (potentially orders of magnitude greater).

We observed a correlation between coefficient heterogeneity (having many
different values for $J_{ij}$ and $h_{i}$) and problem hardness across
different problem classes, while size played a bigger role within a given
class. We also saw that when coefficient heterogeneity is low, the optimal
$|J_{F}|$ stays constant through changes in size (and hardness), while in the
opposite case (high heterogeneity), larger instance size correlated with
larger vertex models and stronger optimal $|J_{F}|$. The conventional wisdom
that large vertex models require stronger $|J_{F}|$ only held up within the
ensemble with high coefficient heterogeneity, while it proved incorrect when
applied across different problem ensembles as well as when coefficient
heterogeneity is low within an ensemble.

A quantitative analysis of the effects that these characteristics have on
hardness and parameter setting, and their interplay, should be a fruitful
field for future study and pave the way for more rigorous parameter setting.
The same can be said of the relationship between optimal pause location and
vertex model size. As new devices with more connected architectures become
available, embedding strategies such as utilizing all the available couplers
between vertex models might need to be reconsidered. These devices will also
support shorter annealing times which will allow us to experimentally verify
what we already strongly suspect—that additional time is only useful in
certain regions, so a pause lowers optimal annealing time and will be able to
improve performance even when that optimal $t_{a}$ is achievable. Theoretical
and experimental limits on how short is too short will require investigation,
as will further tailoring of the annealing schedule, such as optimizing
annealing speed at different parts of the anneal.

##  VII Acknowledgements

We are grateful for support from NASA Ames Research Center, the Autonomous
Systems discipline in the Transformational Tools and Technologies (TTT)
Project of the NASA Transformative Aeronautics Concepts Program, and from
DARPA under IAA 8839 Annexes 125 and 128. ZGI, JM, and ZW are thankful for
support from NASA Academic Mission Services, Contract No. NNA16BD14C. We thank
Alexander Sadovsky for helpful feedback on the information sharing problem.
ZGI thanks David Bernal Neira for useful discussions.

## References

  * Finnila  _et al._ (1994) A.B. Finnila, M.A. Gomez, C. Sebenik, C. Stenson, and J.D. Doll, “Quantum annealing: A new method for minimizing multidimensional functions,” Chemical Physics Letters 219, 343 – 348 (1994).
  * Farhi  _et al._ (2000) Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser, “Quantum computation by adiabatic evolution,” arXiv:quant-ph/0001106 (2000).
  * Kadowaki and Nishimori (1998) Tadashi Kadowaki and Hidetoshi Nishimori, “Quantum annealing in the transverse ising model,” Phys. Rev. E 58, 5355–5363 (1998).
  * Farhi  _et al._ (2001) Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda, “A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem,” Science 292, 472–475 (2001).
  * Aharonov  _et al._ (2007) Dorit Aharonov, Wim van Dam, Julia Kempe, Zeph Landau, Seth Lloyd, and Oded Regev, “Adiabatic quantum computation is equivalent to standard quantum computation,” SIAM Journal on Computing 37, 166–194 (2007).
  * Santoro  _et al._ (2002) Giuseppe E. Santoro, Roman Martoňák, Erio Tosatti, and Roberto Car, “Theory of quantum annealing of an ising spin glass,” Science 295, 2427–2430 (2002).
  * Harris  _et al._ (2010) R. Harris, J. Johansson, A. J. Berkley, M. W. Johnson, T. Lanting, Siyuan Han, P. Bunyk, E. Ladizinsky, T. Oh, I. Perminov, E. Tolkacheva, S. Uchaikin, E. M. Chapple, C. Enderud, C. Rich, M. Thom, J. Wang, B. Wilson, and G. Rose, “Experimental demonstration of a robust and scalable flux qubit,” Phys. Rev. B 81, 134510 (2010).
  * Harris  _et al._ (2010) R. Harris, M. W. Johnson, T. Lanting, A. J. Berkley, J. Johansson, P. Bunyk, E. Tolkacheva, E. Ladizinsky, N. Ladizinsky, T. Oh, F. Cioata, I. Perminov, P. Spear, C. Enderud, C. Rich, S. Uchaikin, M. C. Thom, E. M. Chapple, J. Wang, B. Wilson, M. H. S. Amin, N. Dickson, K. Karimi, B. MacReady, C. J. S. Truncik, and G. Rose, “Experimental investigation of an eight-qubit unit cell in a superconducting optimization processor,” Phys. Rev. B 82, 024511 (2010), arXiv:1004.1628 [cond-mat.supr-con] .
  * Johnson  _et al._ (2011) MW Johnson, MHS Amin, S Gildert, T Lanting, F Hamze, N Dickson, R Harris, AJ Berkley, J Johansson, P Bunyk, _et al._ , “Quantum annealing with manufactured spins,” Nature 473, 194–198 (2011).
  * Boixo  _et al._ (2013) Sergio Boixo, Tameem Albash, Federico M Spedalieri, Nicholas Chancellor, and Daniel A Lidar, “Experimental signature of programmable quantum annealing,” Nature communications 4, 2067 (2013).
  * Bunyk  _et al._ (2014) Paul I. Bunyk, Emile M. Hoskinson, Mark W. Johnson, Elena Tolkacheva, Fabio Altomare, Andrew J. Berkley, Richard Harris, Jeremy P. Hilton, Trevor Lanting, Anthony J. Przybysz, and Jed Whittaker, “Architectural considerations in the design of a superconducting quantum annealing processor,” IEEE Transactions on Applied Superconductivity 24, 1–10 (2014).
  * Jiang  _et al._ (2018) Shuxian Jiang, Keith A. Britt, Alexander J. McCaskey, Travis S. Humble, and Sabre Kais, “Quantum Annealing for Prime Factorization,” Scientific Reports 8, 17667 (2018), arXiv:1804.02733 [quant-ph] .
  * Hauke  _et al._ (2020) Philipp Hauke, Helmut G Katzgraber, Wolfgang Lechner, Hidetoshi Nishimori, and William D Oliver, “Perspectives of quantum annealing: methods and implementations,” Reports on Progress in Physics 83, 054401 (2020).
  * Rønnow  _et al._ (2014) Troels F. Rønnow, Zhihui Wang, Joshua Job, Sergio Boixo, Sergei V. Isakov, David Wecker, John M. Martinis, Daniel A. Lidar, and Matthias Troyer, “Defining and detecting quantum speedup,” Science 345, 420–424 (2014).
  * Boixo  _et al._ (2014) 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, 218–224 (2014).
  * Katzgraber  _et al._ (2015) Helmut G Katzgraber, Firas Hamze, Zheng Zhu, Andrew J Ochoa, and H Munoz-Bauza, “Seeking quantum speedup through spin glasses: The good, the bad, and the ugly,” Physical Review X 5, 031026 (2015).
  * Mandrà and Katzgraber (2018) Salvatore Mandrà and Helmut G. Katzgraber, “A deceptive step towards quantum speedup detection,” Quantum Science and Technology 3, 04LT01 (2018), arXiv:1711.01368 [quant-ph] .
  * Pearson  _et al._ (2019) Adam Pearson, Anurag Mishra, Itay Hen, and Daniel A. Lidar, “Analog errors in quantum annealing: doom and hope,” npj Quantum Information 5, 107 (2019), arXiv:1907.12678 [quant-ph] .
  * 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. Applied 11, 044083 (2019).
  * Venturelli and Kondratyev (2018) Davide Venturelli and Alexei Kondratyev, “Reverse Quantum Annealing Approach to Portfolio Optimization Problems,” arXiv e-prints , arXiv:1810.08584 (2018), arXiv:1810.08584 [quant-ph] .
  * Passarelli  _et al._ (2019) G. Passarelli, V. Cataudella, and P. Lucignano, “Improving quantum annealing of the ferromagnetic $p$-spin model through pausing,” Phys. Rev. B 100, 024302 (2019).
  * Winci  _et al._ (2020) Walter Winci, Lorenzo Buffoni, Hossein Sadeghi, Amir Khoshaman, Evgeny Andriyash, and Mohammad H Amin, “A path towards quantum advantage in training deep generative models with quantum annealers,” Machine Learning: Science and Technology 1, 045028 (2020).
  * Chen and Lidar (2020) H. Chen and D. A. Lidar, “Why and When Pausing is Beneficial in Quantum Annealing,” Phys. Rev. Applied 14, 014100 (2020).
  * Gonzalez Izquierdo  _et al._ (2021) Zoe Gonzalez Izquierdo, Shon Grabbe, Stuart Hadfield, Jeffrey Marshall, Zhihui Wang, and Eleanor Rieffel, “Ferromagnetically shifting the power of pausing,” Phys. Rev. Applied 15, 044013 (2021).
  * Rieffel  _et al._ (2015) Eleanor G. Rieffel, Davide Venturelli, Bryan O’Gorman, Minh B. Do, Elicia M. Prystay, and Vadim N. Smelyanskiy, “A case study in programming a quantum annealer for hard operational planning problems,” Quantum Information Processing 14, 1–36 (2015).
  * Vinci and Lidar (2016) Walter Vinci and Daniel A. Lidar, “Optimally stopped optimization,” Phys. Rev. Applied 6, 054016 (2016).
  * Albash and Lidar (2018) Tameem Albash and Daniel A. Lidar, “Demonstration of a scaling advantage for a quantum annealer over simulated annealing,” Phys. Rev. X 8, 031016 (2018).
  * Quiroz (2019) Gregory Quiroz, “Robust quantum control for adiabatic quantum computation,” Phys. Rev. A 99, 062306 (2019).
  * Grant  _et al._ (2021) Erica Grant, Travis Humble, and Benjamin Stump, “Benchmarking quantum annealing controls with portfolio optimization,” Physical Review Applied 15 (2021), 10.1103/PhysRevApplied.15.014012.
  * Weber  _et al._ (2017) Steven J. Weber, Gabriel O. Samach, David Hover, Simon Gustavsson, David K. Kim, Alexander Melville, Danna Rosenberg, Adam P. Sears, Fei Yan, Jonilyn L. Yoder, William D. Oliver, and Andrew J. Kerman, “Coherent coupled qubits for quantum annealing,” Phys. Rev. Applied 8, 014004 (2017).
  * Note (1) For definitions of the pause location and other concepts related to quantum annealing, see Sec. II.
  * Note (2) We use vertex model to refer to the set of physical qubits representing a single logical variable after minor-embedding. The word “chain” is sometimes used in the literature with the same meaning, but we prefer vertex model given that their structure is often not chain-like.
  * (33) Husni R. Idris, Quang Dao, Robert C. Rorie, and Kelley Hashemi, “A framework for assessment of autonomy challenges in air traffic management,” in  _AIAA AVIATION 2020 FORUM_.
  * (34) David P. Thipphavong, Rafael Apaza, Bryan Barmore, Vernol Battiste, Barbara Burian, Quang Dao, Michael Feary, Susie Go, Kenneth H. Goodrich, Jeffrey Homola, Husni R. Idris, Parimal H. Kopardekar, Joel B. Lachter, Natasha A. Neogi, Hok Kwan Ng, Rosa M. Oseguera-Lohr, Michael D. Patterson, and Savita A. Verma, “Urban air mobility airspace integration concepts and considerations,” in  _2018 Aviation Technology, Integration, and Operations Conference_.
  * Morita and Nishimori (2008) Satoshi Morita and Hidetoshi Nishimori, “Mathematical foundation of quantum annealing,” Journal of Mathematical Physics 49, 125210 (2008).
  * Choi (2020) V. Choi, “The effects of the problem Hamiltonian parameters on the minimum spectral gap in adiabatic quantum optimization,” Quant. Inf. Process. 19, 90 (2020).
  * Lucas (2014) Andrew Lucas, “Ising formulations of many NP problems,” Frontiers in Physics 2, 5 (2014).
  * Choi (2008) V. Choi, “Minor-embedding in adiabatic quantum computation: I. the parameter setting problem,” Quantum Inf. Process. 7, 193 (2008).
  * Fang and Warburton (2019) Yan-Long Fang and P. Warburton, “Minimizing minor embedding energy: an application in quantum annealing,” arXiv:1905.03291 (2019).
  * 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).
  * Marshall  _et al._ (2020) J. Marshall, A. Di Gioacchino, and E. G. Rieffel, “Perils of embedding for sampling problems,” Phys. Rev. Research 2, 023020 (2020).
  * Marshall  _et al._ (2022) J. Marshall, G. Mossi, and E. G. Rieffel, “Perils of embedding for quantum sampling,” Phys. Rev. A 105, 022615 (2022).
  * Albash and Marshall (2021a) Tameem Albash and Jeffrey Marshall, “Comparing relaxation mechanisms in quantum and classical transverse-field annealing,” Phys. Rev. Applied 15, 014029 (2021a).
  * DWa (2019) “QPU Properties: D-Wave 2000Q System at NASA Ames,” D-Wave User Manual 09-1151A-D (2019).
  * McGeoch and Farré (2020) Catherine McGeoch and Pau Farré,  _The D-Wave Advantage System: An Overview_, Tech. Rep. (2020).
  * Marshall  _et al._ (2017) Jeffrey Marshall, Eleanor G. Rieffel, and Itay Hen, “Thermalization, freeze-out, and noise: Deciphering experimental quantum annealers,” Phys. Rev. Applied 8, 064025 (2017).
  * Boothby  _et al._ (2019) K. Boothby, P. Bunyk, J. Raymond, and A. Roy,  _Next-Generation Topology of D-Wave Quantum Processors_, Tech. Rep. (2019).
  * Pokharel  _et al._ (2021) Bibek Pokharel, Zoe Gonzalez Izquierdo, P. Aaron Lott, Elena Strbac, Krzysztof Osiewalski, Emmanuel Papathanasiou, Alexei Kondratyev, Davide Venturelli, and Eleanor Rieffel, “Inter-generational comparison of quantum annealers in solving hard scheduling problems,” arXiv e-prints , arXiv:2112.00727 (2021), arXiv:2112.00727 [quant-ph] .
  * DWa (2020) “QPU-Specific Physical Properties: DW_2000Q_6,” D-Wave User Manual 09-1215A-B (2020).
  * (50) “Error sources for problem representation,” https://docs.dwavesys.com/docs/latest/c_qpu_ice.html.
  * adv (2020) “QPU-Specific Physical Properties: Advantage_system1.1,” D-Wave User Manual 09-1237A-A (2020).
  * Grygorash  _et al._ (2006) Oleksandr Grygorash, Yan Zhou, and Zach Jorgensen, “Minimum spanning tree based clustering algorithms,” in  _2006 18th IEEE International Conference on Tools with Artificial Intelligence (ICTAI’06)_ (2006) pp. 73–81.
  * Xu  _et al._ (2002) Ying Xu, Victor Olman, and Dong Xu, “Clustering gene expression data using a graph-theoretic approach: an application of minimum spanning trees,” Bioinformatics 18, 536–545 (2002), https://academic.oup.com/bioinformatics/articl-pdf/18/4/536/669112/180536.pdf .
  * Cormen  _et al._ (2001) Thomas H. Cormen, Charles E. Leiserson, and Ronald L. Rivestand Clifford Stein,  _”Introduction to Algorithms”_ (MIT Press, Cambridge, MA, 2001).
  * Garey and Johnson (1979) Michael R. Garey and David S. Johnson,  _Computers and Intractability: A Guide to the Theory of NP-Completeness_ (W. H. Freeman & Co., New York, NY, USA, 1979).
  * Sadovsky (2014) A. Sadovsky, “Application of the shortest-path problem to routing terminal airspace air traffic,” Journal of Aerospace Information Systems 11 (2014).
  * Albash and Marshall (2021b) Tameem Albash and Jeffrey Marshall, “Comparing relaxation mechanisms in quantum and classical transverse-field annealing,” Phys. Rev. Applied 15, 014029 (2021b).
  * Amin (2015) Mohammad H. Amin, “Searching for quantum speedup in quasistatic quantum annealers,” Phys. Rev. A 92, 052323 (2015).
  * Note (3) It can fail because of the edge capacity limitations, defined below.

##  Appendix A Device schedules and connectivity

Fig. 12 shows the annealing functions $A(s)$ and $B(s)$ for both devices used
in our experimental runs. Fig. 13 shows the graphs of their respective
architectures.

Figure 12: Top: $A(s)$ and $B(s)$ functions for the two devices. Energy scales
for the driver and problem Hamiltonian of the two devices we use in our
experimental runs, in units of GHz, and $h=1$. Bottom: $Q(s)$ and $C(s)$ for
the two devices. Dimensionless quantum ($Q(s)$) and classical ($C(s)$) scales
for the two devices we use in our experimental runs.

Figure 13: Top: Structure of the Chimera graph featured by DW2K. Shown is a
so-called $C_{3}$ graph, i.e. a square with 3 cells per side. The pattern
repeats in the same fashion to scale to larger chips. In particular, the
device we used has a $C_{1}6$. Bottom: Structure of the Pegasus graph of DWA.
Shown is a $P_{3}$ graph. As with Chimera, the pattern repeats and the actual
device has a $P_{1}6$. Chimera is a subgraph of Pegasus.

##  Appendix B Problem mapping: minimum spanning tree

Consider a graph $G=(V,E)$ with weights $w(E)$ for each edge, from which we
wish to obtain a minimal weighted spanning tree with maximum degree $\Delta$,
i.e. find its BD-MST. This will involve minimizing the sum of the weights of
the tree edges, represented by the cost function

$\displaystyle C_{0}=\sum_{p,v}w_{pv}x_{p,v},$ (4)

which we explain below. Several constraints will also be imposed to ensure
that the graph is in fact a spanning tree and its degree is bounded by
$\Delta$.

A root for the tree will be picked randomly or based on problem
structure—generally, picking a high-degree vertex as the root will result in
lower resource costs—and assigned to level $1$. Its children will be at level
$2$, their children at level $3$ and so on, leading to the ‘level-based’
designation.

The variables $x_{p,v}$ appearing in Eq. (4) represent the parent-child
relationships in the tree; $x_{p,v}=1$ if $p$ is the (adjacent) parent of $v$
(and $0$ if not). The indices $p,v$ range over $p=1,\dots,n$ and
$v=2,\dots,n$, restricted to (intersected with) pairs $(p,v)$ or $(v,p)$ that
occur in $E$. Thus there are $2$ variables for every edge not containing the
root, and one for every root edge, giving $2m-d_{r}$ total $x_{p,v}$
variables, with $m$ being the number of edges in $E$ and $d_{r}$ the degree of
the root.

Since our problem needs to be in QUBO form, the constraints will be expressed
as penalty terms. The first penalty term enforces that every node (except the
root) has exactly one parent,

$\displaystyle C_{pen1}=\sum_{v\in\\{2,\dots,n\\}}\left(\sum_{p:(pv)\in
E}x_{p,v}-1\right)^{2}\;.$ (5)

The number of terms in the sum is $2m-d_{r}$, i.e. equal to the number of
variables $x_{p,v}$.

The second penalty term enforces that each vertex exists at exactly one level
in the tree,

$\displaystyle
C_{pen2}=\sum_{v\in\\{2,\dots,n\\}}\left(\sum_{\ell=2}^{n}y_{v,\ell}-1\right)^{2}\;.$
(6)

It introduces the $y_{v,\ell}$ variables, with $y_{v,\ell}=1$ if $v$ is at
depth $\ell$ of the tree, $v=2,\dots,n$, $\ell=2,\dots,n$. There are
$(n-1)^{2}$ such variables. However, since the number of variables will
eventually determine how many logical qubits the problem requires, it is in
our interest to reduce it as much as possible. By picking the root smartly the
range of $\ell$ can be reduced. We also carry out the following pre-
processing: taking the original graph $G=(V,E)$, the distance from each node
to the one we have selected as the tree root is calculated. Given that it is
impossible for a node to be at a level smaller than its distance to the root,
we can avoid generating any $y_{v,\ell}$ for which that is the case, further
bringing down the total number of $y_{v,\ell}$ variables.

The third penalty term enforces that the tree has degree at most $\Delta$,

$\displaystyle C_{pen3}=$ $\displaystyle\sum_{p=2}^{v}\left(\sum_{v:(pv)\in
E}x_{p,v}-\sum_{j=1}^{\Delta-1}z_{p,j}\right)^{2}$ $\displaystyle+$
$\displaystyle\left(\sum_{v:(1v)\in
E}x_{1,v}-\sum_{j=1}^{\Delta}z_{1,j}\right)^{2}.$ (7)

It is separated into two terms to account for the fact that the root can have
up to $\Delta$ children, while all other nodes cannot have more than
$(\Delta-1)$, since they have a parent. To enforce the inequality
$\sum_{v:(pv)\in E}x_{p,v}\leq\Delta-1$, integer variable
$z_{p}\in[0,\Delta-1]$ is introduced as slack variable, and the inequality is
enforced as equality $\sum_{v:(pv)\in E}x_{p,v}=z_{p}$. The integer variable
is further encoded into binary variables $z_{p,j}$. In general, various
encoding methods can be applied to encode an integer into binaries, including
binary, unary, and one-hot encodings. While binary encoding is most efficient
for integers of value power of two, We use unary encoding here, which can be
applied straightforward to arbitrary value of $\Delta$.

The fourth and final penalty term enforces that the tree encoding is
consistent, i.e., that if $p$ is the parent of $v$ then its level is one less
than $v$’s,

$\displaystyle C_{pen4}=$
$\displaystyle\sum_{p,v}\sum_{\ell=3}^{n}x_{p,v}y_{v,\ell}(1-y_{p,\ell-1})$
$\displaystyle+$
$\displaystyle\sum_{v=2}^{d_{r}}x_{1,v}(1-y_{v,2})+\sum_{v=2}^{d_{r}}y_{v,2}(1-x_{1,v})\;,$
(8)

where the last two sums handle the edges connected to the root and their terms
are quadratic, while the first sum deals with the remaining edges and produces
cubic terms of the form $x_{p,v}y_{v,\ell}(1-y_{p,\ell-1})$. While the
original number of cubic terms would be

$(2m-2d_{r})*(n-2),$

thanks to the preprocessing of the $y_{v,\ell}$ variables this number is
reduced. Because cubic terms cannot be directly encoded in D-Wave, we
introduce an ancilla variable $a_{p,v,\ell}$ to encode $x_{p,v}y_{v,\ell}$,
and accordingly a penalty function $f(x,y,a)=3a+xy-2ax-2ay$ is added to raise
a penalty if $a=xy$ is violated. The term $x_{p,v}y_{v,\ell}(1-y_{p,\ell-1})$
then can be replaced by quadratic terms

$\displaystyle 4a-ay_{p,\ell-1}+x_{p,v}y_{v,\ell}-2ax_{p,v}-2ay_{v,\ell}\;.$
(9)

The total number of variables (and hence, logical qubits) without
preprocessing is at most:

$\displaystyle 2m-d_{r}+(n-1)^{2}+n(\Delta-1)+1+(2m-2d_{r})(n-2)$
$\displaystyle\simeq 2mn+n^{2}$

This would mean, for instance, that the complete graph $K_{5}$ with $\Delta=3$
would require between $86$ and $100$ logical qubits (depending on $d_{r}$).
With pre-processing, we are able to bring this number down to $74$.

Finally, we can write the overall objective function as

$\displaystyle C=C_{0}+A(C_{pen1}+C_{pen2}+C_{pen3}+C_{pen4})\;,$ (10)

and accordingly the cost Hamiltonian $H_{C}$. In Eq. (10) we have defined the
minimum penalty weight to be the maximum edge weight

$\displaystyle A=w_{max}+\varepsilon=\max_{(uv)\in E}w_{uv}+\varepsilon\;.$
(11)

##  Appendix C Problem mapping: graph coloring

We use a standard one-hot encoding for the GC problem. We have an undirected
graph $G=(V,E)$ with $n$ nodes, and a number of colors $k$. We define the
binary variables:

$x_{v,i}=\begin{cases}1&\text{if node $v$ is colored with color $i$}\\\
0&\text{otherwise},\end{cases}$ (12)

for a total of $n$x$k$ variables. Then our Hamiltonian is simply:

$H=\sum_{v=1}^{n}\left(1-\sum_{i=1}^{k}x_{v,i}\right)^{2}+\sum_{(uv)\in
E}\sum_{i=1}^{k}x_{u,i}x_{v,i}.$ (13)

The first term enforces that each node only has one color assigned, while the
second one penalizes any pair of nodes connected by an edge which have the
same color. $H=0$ if there is a $k$ coloring for graph $G$, and $H>0$
otherwise.

##  Appendix D Problem mapping: information sharing

###  D.1 Problem inputs

Given:

  * •

A discrete time variable $t$, taking the non-negative integer values.

  * •

A finite set $I=\\{0,1,\ldots,n-1\\}$, whose elements each represent a message
to be transmitted over a communications network, modeled below.

  * •

A graph $G=(V,E)$, which models a communications network. All messages $i\in
I$ are sent from the same node, called the sender node and denoted
$v_{{\footnotesize sender}}\in V$. Each sent message $i$ is intended to reach,
eventually, its recipient node $r_{i}\in V$. (Different messages generally
have different recipients.)

  * •

A transmission time function $l$ that assigns to each pair $(i,e)$, where
$i\in I$ is a message and $e\in E$ an edge in the communications graph $G$, a
time duration (a non-negative integer) $l(i,e)$, which models how long message
$i$ will take to traverse edge $e$. The latter value, when it is necessary to
indicate the constituent vertices $(v,v^{\prime})$ of $e$, can be written in
the more detailed notation $l(i,v,v^{\prime})$.

  * •

A transmission of a message along a path $P$, which is a triple

$T(i)=\big{(}i,t^{em}_{i},P=(v_{0}=v_{{\footnotesize
sender}},v_{1},\ldots,v_{N_{i}}=r_{i})\big{)},$ (14)

where $i\in I$ is a message, $t^{em}_{i}$ is the emission time of the message
(i.e., the instant when the message is sent), and $P$ is a path in $G$ from
the sender node to the recipient node of message $i$.

The set of all possible transmissions of message $i$, along all possible paths
$P$ and with all possible emission times, will be denoted by ${\cal T}^{i}$.
It follows that once a message transmission (14) has succeeded333It can fail
because of the edge capacity limitations, defined below., its total duration
is given by

$D(i,P):=\sum_{k=0}^{N_{i}-1}l(i,v_{k},v_{k+1}),$ (15)

and that message $i$ arrives at its recipient node at time

$t^{em}_{i}+D(i,P).$ (16)

The partial sums from $D(i,P)$ correspond to the times at which message $i$
enters a given edge in path $P$. In more detail, under transmission (14), edge
$(v_{j},v_{j+1})$ is entered by message $i$ at time

$t^{em}_{i}+\sum_{k=0}^{j}l(i,v_{k},v_{k+1}).$ (17)

  * •

A scheduled time of emission for each message $i\in I$,

$t^{sched}_{i}.$ (18)

In no transmission (14) of that message can the actual $t^{em}_{i}$ emission
time occur before $t^{sched}_{i}$; thus, one of the constrains is the
inequality

$t^{em}_{i}\geq t^{sched}_{i}.$ (19)

The non-negative difference

$t^{em}_{i}-t^{sched}_{i}$ (20)

is called the delay of transmission (14).

  * •

A transmission schedule; a mapping $S$ that assigns to each message $i\in I$ a
transmission:

$S:i\mapsto T(i)\in{\cal T}^{i}.$ (21)

The set of all possible transmission schedules is the Cartesian product

${\cal S}={\cal T}^{0}\times{\cal T}^{1}\times\ldots\times{\cal T}^{n-1}.$
(22)

  * •

A cost of transmission delay for each message $i\in I$, $c_{i}$, that
represents how much it costs to delay a message per unit of time, and can be
interpreted as a general way of expressing how important it is for the message
to arrive promptly.

  * •

A node capacity function $B$. Each node in $G$ is assumed capable of allowing
entry to no more than some finite number of messages. This number will be
called the capacity of the node. In some cases, a message can have such high
priority that its arrival without delay must be guaranteed (i.e. a message $i$
such that $c_{i}\gg c_{j}\forall j\neq i$). The path of such a message is
initially hardcoded into the problem through the node capacity function $B$,
with the addition of time dependence (the capacity is reduced at the nodes and
times corresponding to the message’s path). Thus, the node capacity function
$B$ (treated as a discrete control variable) assigns to each node $v$ and at
each time instant $t\in\\{0,1,\ldots\\}$ a non-negative integer

$B(t,v)\equiv B^{v}_{t}.$ (23)

The Information Sharing Problem is to find a transmission schedule

$S=(T(0),T(1),\ldots,T(n-1))\in{\cal S}$ (24)

that minimizes the total cost of delay

$\sum_{i\in I}c_{i}\left(t^{em}_{i}-t^{sched}_{i}\right),$ (25)

subject to the constraints (19) and

$\left(\begin{array}[]{c}\mbox{number of messages}\\\ \mbox{entering node
$v$}\\\ \mbox{at time $t$}\\\ \mbox{under schedule $S$}\\\
\end{array}\right)\leq B^{v}_{t}\quad\left\\{\begin{array}[]{c}\mbox{for all
node entry}\\\ \mbox{times $t$ involved in}\\\ \mbox{one or more of the}.\\\
\mbox{transmissions $S(i)$}.\\\ \end{array}\right.$ (26)

###  D.2 Mapping to QUBO

We use a quantum annealer to solve instances of the information sharing
problem where, given a path $P_{i}$ for each message $i\in I$, an emission
time $t^{em}_{i}$ is assigned to each message such that the total cost of
delays is minimized. For simplicity, the scheduled time of emission
$t^{sched}_{i}$ is taken to be 0 $\forall i\in I$.

The problem must be mapped to QUBO to be solved on the quantum annealer. The
objective function to minimize will be:

$C_{obj}=\sum_{i\in I}c_{i}a_{i},$ (27)

where $a_{i}$ is the total delay for message $i$, i.e. the actual arrival time
minus the scheduled arrival time:

$a_{i}=\sum_{t_{a}}t_{a}x_{i,{t_{a}}}^{r_{i}}-D(i,P_{i}),$ (28)

where the sum is over the possible arrival times $t_{a}$, and $D(i,P_{i})$ is
given by Eq. 15.

In the expression for the total delay, we have used the binary variables:

$x_{i,t}^{v}=\begin{cases}1&\text{if message $i$ arrives at node $v$ at time
$t$}\\\ 0&\text{otherwise},\end{cases}$ (29)

which means that for given $i$ and $v$, $x_{i,t}^{v}=1$ for a single value of
$t$ and 0 for all others.

We will have three constraints that are expressed as penalties to formulate a
QUBO. The first penalty term enforces that the capacity of the network is not
exceeded:

$C_{pen1}=\sum_{t=1}^{t_{h}}\sum_{v\in P_{i}}\left(\sum_{i\in
I}x_{i,t}^{v}-\sum_{k=1}^{B_{t}^{v}}s_{k,t}^{j}\right)^{2},$ (30)

where the $s_{k,t}^{v}\in[0,1]$ are slack variables, and $t_{h}$ is the time
horizon, a limit on the latest time we consider for all the messages to have
been transmitted, which is needed for practical purposes (the sum needs an
upper limit). If $\sum_{i\in I}x_{i,t}^{v}\leq B_{t}^{v}$ is satisfied, there
exists a variable assignment for the $s_{k,t}^{v}$ that makes $C_{pen1}=0$.

The second penalty term enforces path connectivity is respected, by ensuring
that a message does not arrive at a node faster than the travel time between
that node and the previous one in its path would allow. This penalty term
counts the number of path connectivity violations:

$C_{pen2}=\sum_{i\in I}\sum_{v\in
P_{i}}\sum_{t=1}^{t_{h}-l_{v,v+1}^{i}}\sum_{k=0}^{l_{v,v+1}^{i}-1}x_{i,t+k}^{v+1}x_{i,t}^{v}.$
(31)

The third and last penalty term makes sure that each message arrives at each
of the nodes in its path exactly once; this ensures that the message is
actually transmitted and travels through its path as it should:

$C_{pen3}=\sum_{i\in I}\sum_{v\in
P_{i}}\left(\left(\sum_{t=1}^{t_{h}}x_{i,t}^{v}\right)-1\right)^{2}.$ (32)

If we generated every variable regardless of feasibility (i.e., whether there
exists a valid scenario in which the variable is non-zero), we would have
$n$x$|E|$x$t_{h}$ $x_{i,t}^{v}$ variables. We can reduce this number with two
considerations about feasibility. First, we only need to generate variables
for the nodes in each message’s path, but not the rest. This reduces the
$n$x$|E|$ to $\prod_{i=0}^{n-1}N_{i}$. Second, we can avoid times that are too
early, by pre-calculating the earliest time that each message $i$ can arrive
at node $v_{j}\in P_{i}$:

$t_{e,v_{j}}^{i}=\sum_{j^{\prime}=0}^{j-1}l(i,v_{j^{\prime}},v_{j^{\prime}+1}).$
(33)

Then, instead of generating variables for the time interval $[1,t_{h}]$ we can
do it for $[t_{e,v_{j}}^{i},t_{h}]$ for each $i\in I$ and $v_{j}\in P_{i}$.
With this, the total number of $x_{i,t}^{v}$ variables is $\sum_{i\in
I}N_{i}t_{e,v}^{i}$.

With this, the sum over time in the second penalty term $C_{pen2}$ would start
at $t=t^{i}_{e,v^{i}_{1}}$ rather than $t=1$, where $v_{1}$ is the second node
in the path of message $i$. We can start this sum at the earliest time that
the message can be in its second path node, because the variables
corresponding to the first path node are created to respect the minimum travel
time, and thus there is no need to double check them. For the third penalty
term $C_{pen3}$, the sum over time can start at $t=t^{i}_{e,v\in P_{i}}$
instead of $t=1$.

Finally, the number of slack variables is $\prod_{v\in\bigcup\limits_{i\in
I}P_{i}}\prod_{t=1}^{t_{h}}B_{t}^{v}$, and all need to be generated. The first
penalty term $C_{pen1}$ is unaffected by the binary variable pruning described
in the previous paragraph; all times beginning at $t=1$ need to be checked for
bandwidth violations, given that these are unrelated to minimum travel times.

Taking all of the above into consideration, the full Hamiltonian for our QUBO
problem is:

$H=C_{obj}+A_{1}C_{pen1}+A_{2}C_{pen2}+A_{3}C_{pen3},$ (34)

where $A_{i}$ are the penalty weights. Penalty weights are needed because
without them it could be advantageous to violate a penalty term, by virtue of
such a penalty being smaller than the cost reduction it affords. To avoid this
scenario, we need to consider the maximum cost reduction that could be
obtained by each penalty violation, and then choose a penalty weight that will
ensure the penalty is always larger than said cost reduction.

Let us start with the first penalty term $C_{pen1}$, which enforces bandwidth.
At each time and each node in a path, surpassing bandwidth by one message
results in a $+1$ penalty before multiplying by $A_{1}$. Additional messages
would cause an even larger penalty due to the square, so it is sufficient to
consider the first one. The largest possible reduction to $C_{obj}$ without
violating any additional penalties would be the maximum, over all messages, of
a message going from arriving at the latest possible time to the earliest
possible: $\text{max}_{i}(c_{i}(t_{h}-D(i,P_{i}))$. A further reduction could
be attained by arrival at an earlier time than physically possible (i.e.
smaller than $D(i,P_{i})$), but this would result in an additional penalty by
either not arriving at all (reflected in $C_{pen3}$) or not respecting path
connectivity ($C_{pen2}$). It is then sufficient to choose
$A_{1}=\text{max}_{i}(c_{i}(t_{h}-D(i,P_{i})))+\varepsilon$.

The second penalty term $C_{pen2}$ enforces path connectivity, through pairs
of $x$ variables that cannot be 1 at the same time because they correspond to
subsequent path nodes of a given message at times closer than the minimum
travel time between said nodes. Each such pair where both variables are 1
results in a $+1$ penalty before multiplying by $A_{2}$. The maximum $C_{obj}$
reduction that can be obtained from a single pair occurs when traveling is
instantaneous, i.e. the message is early by some time $l(i,v_{j},v_{j+1})$. So
the penalty weight can be set to
$A_{2}=\text{max}_{i,j}c_{i}l(i,v_{j},v_{j+1})+\varepsilon$. Note that this is
never larger than the overall maximum reduction to $C_{obj}$ described in the
previous paragraph ($\text{max}_{i}(c_{i}(t_{h}-D(i,P_{i})))$), given that
$t_{h}\geq D(i,P_{i})$ and $D(i,P_{i})\geq l(i,v_{j},v_{j+1})$ (see Eq. 15).

Figure 14: Optimal $|J_{F}|$ for ensemble of 50 $n=6$ BD-MST instances (DWA).
Figure 15: Improvement of $T_{S}$ with pause for ensemble of 50 $n=6$ BD-MST
instances (DWA).

Finally, the third penalty term $C_{pen3}$ involves arrivals. A penalty of +1
(before multiplying by $A_{3}$) is incurred every time a message does not
arrive to one of its path nodes. A $C_{pen3}$ violation has the potential for
the greatest reduction in $C_{obj}$, in the case where a message should arrive
at its destination at time $t_{h}$, and it simply does not arrive, reducing
$C_{obj}$ by $c_{i}t_{h}$. Then, we need to set
$A_{3}=\text{max}_{i}c_{i}t_{h}+\varepsilon$.

Because this last penalty weight is the highest of the three, if we want to
set a single penalty weight for all penalties that will suffice in any
situation, $A_{1}=A_{2}=A_{3}=\text{max}_{i}c_{i}t_{h}+\varepsilon$ is the
appropriate choice.

For our particular set of instances, we empirically found a lower penalty
weight of $A_{1}=A_{2}=A_{3}=1/2\sum_{i\in
I}c_{i}\sum_{t=t_{a_{min}}}^{t_{h}}t+\varepsilon$ to be sufficient, while
helping keep coefficients at a similar scale to avoid running into precision
issues with the hardware.

##  Appendix E Results for $n=6$ BD-MST instances

We also extend our results to an ensemble of 50 $n=6$ instances. These are
much harder than the $n=5$ and, when a few were attempted on DW2K, none of
them could be solved. Instead of 50,000 reads, we use $10^{6}$ reads to obtain
good statistics. As shown in Fig. 14, $|J_{F}|=0.8$ still appears to yield the
best $T_{S}$, although it is not as clear as for the $n=5$ ensemble. Due to
the costly nature of solving these instances, we have not exhaustively
explored different schedules as we did for the previous problems, but were
able to verify that a pause of $t_{p}=0.2\mu$s at $s_{p}=0.4$ improves upon
the no pause results, which can be seen in Fig. 15. Without a pause, 11 of the
50 instances did not solve after $10^{6}$ reads. With the pause, 10 out of 50
did not solve. It is possible that some of them would solve at a different
$|J_{F}|$, given that a limited range was tested and, as we have seen for the
$n=5$ ensemble, particularly difficult instances tend to require a stronger
$|J_{F}|$.

We were able to confirm that for these difficult instances, our embedding
method of picking the smallest (in number of physical qubits) out of 10 leads
to a significant improvement in the results. For $|J_{F}|=1.0$, if we simply
use the first embedding we generate, we find $p_{success}=$
$7\text{\times}{10}^{-6}\text{\,}\mathrm{,}$ with $35^{th}$ and $65^{th}$
percentile values at $6\text{\times}{10}^{-6}\text{\,}\mathrm{a}$nd
$8\text{\times}{10}^{-6}\text{\,}\mathrm{r}$espectively. By generating 10
embeddings and picking the smallest, we increase to $p_{success}=$
$1.4\text{\times}{10}^{-5}\text{\,}\mathrm{,}$ with $35^{th}$ and $65^{th}$
percentiles at $1.2\text{\times}{10}^{-5}\text{\,}\mathrm{a}$nd
$1.7\text{\times}{10}^{-5}\text{\,}\mathrm{,}$ a $2\times$ improvement.