# Quantum error mitigation in quantum annealing

Mohammad H. Amin<EMAIL_ADDRESS>D-Wave Systems, Burnaby, British Columbia,
Canada  Department of Physics, Simon Fraser University, Burnaby, British
Columbia, Canada  Andrew D. King  D-Wave Systems, Burnaby, British Columbia,
Canada  Jack Raymond  D-Wave Systems, Burnaby, British Columbia, Canada
Richard Harris  D-Wave Systems, Burnaby, British Columbia, Canada  William
Bernoudy  D-Wave Systems, Burnaby, British Columbia, Canada  Andrew J. Berkley
D-Wave Systems, Burnaby, British Columbia, Canada  Kelly Boothby  D-Wave
Systems, Burnaby, British Columbia, Canada  Anatoly Smirnov  D-Wave Systems,
Burnaby, British Columbia, Canada  Fabio Altomare  D-Wave Systems, Burnaby,
British Columbia, Canada  Michael Babcock  D-Wave Systems, Burnaby, British
Columbia, Canada  Catia Baron  D-Wave Systems, Burnaby, British Columbia,
Canada  Jake Connor  D-Wave Systems, Burnaby, British Columbia, Canada  Martin
Dehn  D-Wave Systems, Burnaby, British Columbia, Canada  Colin Enderud  D-Wave
Systems, Burnaby, British Columbia, Canada  Emile Hoskinson  D-Wave Systems,
Burnaby, British Columbia, Canada  Shuiyuan Huang  D-Wave Systems, Burnaby,
British Columbia, Canada  Mark W. Johnson  D-Wave Systems, Burnaby, British
Columbia, Canada  Eric Ladizinsky  D-Wave Systems, Burnaby, British Columbia,
Canada  Trevor Lanting  D-Wave Systems, Burnaby, British Columbia, Canada
Allison J. R. MacDonald  D-Wave Systems, Burnaby, British Columbia, Canada
Gaelen Marsden  D-Wave Systems, Burnaby, British Columbia, Canada  Reza Molavi
D-Wave Systems, Burnaby, British Columbia, Canada  Travis Oh  D-Wave Systems,
Burnaby, British Columbia, Canada  Gabriel Poulin-Lamarre  D-Wave Systems,
Burnaby, British Columbia, Canada  Hugh Ramp  D-Wave Systems, Burnaby, British
Columbia, Canada  Chris Rich  D-Wave Systems, Burnaby, British Columbia,
Canada  Berta Trullas Clavera  D-Wave Systems, Burnaby, British Columbia,
Canada  Nicholas Tsai  D-Wave Systems, Burnaby, British Columbia, Canada  Mark
Volkmann  D-Wave Systems, Burnaby, British Columbia, Canada  Jed D. Whittaker
D-Wave Systems, Burnaby, British Columbia, Canada  Jason Yao  D-Wave Systems,
Burnaby, British Columbia, Canada  Niclas Heinsdorf  Department of Physics and
Astronomy and Quantum Matter Institute, The University of British Columbia,
Vancouver, British Columbia, Canada  Nitin Kaushal  Department of Physics and
Astronomy and Quantum Matter Institute, The University of British Columbia,
Vancouver, British Columbia, Canada  Alberto Nocera  Department of Physics and
Astronomy and Quantum Matter Institute, The University of British Columbia,
Vancouver, British Columbia, Canada  Marcel Franz  Department of Physics and
Astronomy and Quantum Matter Institute, The University of British Columbia,
Vancouver, British Columbia, Canada

###### Abstract

Quantum Error Mitigation (QEM) presents a promising near-term approach to
reduce error when estimating expectation values in quantum computing. Here, we
introduce QEM techniques tailored for quantum annealing, using Zero-Noise
Extrapolation (ZNE). We implement ZNE through zero-temperature extrapolation
as well as energy-time rescaling. We conduct experimental investigations into
the quantum critical dynamics of a transverse-field Ising spin chain,
demonstrating the successful mitigation of thermal noise through both of these
techniques. Moreover, we show that energy-time rescaling effectively mitigates
control errors in the coherent regime where the effect of thermal noise is
minimal. Our ZNE results agree with exact calculations of the coherent
evolution over a range of annealing times that exceeds the coherent annealing
range by almost an order of magnitude.

##  I Introduction

Quantum computation has emerged as a computational paradigm with the potential
to solve complex problems that remain beyond the capabilities of classical
computers. Nevertheless, its performance is substantially hampered by
environmental noise and hardware imperfections. Quantum Error Correction (QEC)
[1, 2] is regarded as the ultimate solution to eliminate the impact of these
errors. However, the significant overhead of QEC limits its practicality only
to very large scales, far beyond the current state of technology [3].
Recently, Quantum Error Mitigation (QEM) has been proposed [4, 5, 6, 7] as a
near-term solution that can be used to estimate error-free expectation values
when the impact of noise is small. Among the various QEM techniques, Zero-
Noise Extrapolation (ZNE) [4, 5] stands out as one of the most practical
methods. In ZNE, one systematically varies the noise amplitude experienced by
the system. By observing the system’s response to this controlled change, it
becomes possible to make predictions about how the system would behave under
noise-free conditions.

While QEM was initially developed and tested for circuit model quantum
computing [8, 9, 10, 11, 12], the same principles can be applied to other
quantum computing protocols, including quantum annealing (QA) [13]. Quantum
Anneling Correction (QAC), founded on repetition codes, has been proposed and
examined as a means to enhance the performance of QA at the expense of
requiring a larger number of qubits [14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24]. Over the past few years, quantum simulation has evolved into an important
application of QA for the exploration of exotic states of condensed matter
systems [25, 26, 27, 28] and their quantum phase transitions [29, 30]. QEM can
enhance the accuracy of expectation values in quantum simulation experiments
without incurring any additional overhead in qubit count. Indeed, the first
experimental attempt to extract information about noise-free evolution via
extrapolation was reported in Ref. [31], in which QA performance was examined
as a function of temperature for a problem with a very small spectral gap. By
extrapolating the final ground state probability to zero temperature, the
noise-free Landau-Zener behavior was replicated (see Fig. 3 and 4 of Ref.
[31]).

In this paper, we provide a theoretical description of ZNE in quantum
annealing, with a particular emphasis on extracting information about quantum
phase transitions. We then employ ZNE in experimental investigations into the
critical dynamics of a 1D quantum Ising spin chain. The extrapolated results
compare well with exact solutions as well as time dependent density matrix
renormalization group (DMRG) simulations for the closed system Schrödinger
evolution across a wide range of control parameters.

##  II Theory

The Hamiltonian of an annealing-based quantum processing unit (QPU) is written
as

$H_{S}(s)=-\Gamma(s)\sum_{i}\sigma_{i}^{x}+\mathcal{J}(s)H_{P},$ (1)

where

$H_{P}=\sum_{i}h_{i}\sigma_{i}^{z}+\sum_{\langle
i,j\rangle}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ (2)

is the problem Hamiltonian and $s=t/t_{a}$, with $t$ being time and $t_{a}$
being the annealing time. The energy scales $\Gamma(s)$ and $\mathcal{J}(s)$
evolve with time in such a way that $\Gamma(0)\gg\mathcal{J}(0)$ and
$\Gamma(1)\ll\mathcal{J}(1)$. They are predetermined at the time of
calibration (see Fig. 4 in Appendix C). The dimensionless biases ($h_{i}$) and
coupling coefficients ($J_{ij}$) are programmable.

To incorporate all sources of error, whether thermal, $1/f$ noise, or
parameter misspecification, we consider the qubits to be in contact with their
environment. The total Hamiltonian is written as

$H=H_{S}+H_{B}+H_{\rm int},$ (3)

where $H_{B}$ is the environment Hamiltonian, and

$H_{\rm int}=-\sqrt{\gamma}\sum_{\alpha}Q^{\alpha}{\cal O}^{\alpha}$ (4)

describes the interaction between system and environment. The operators ${\cal
O}^{\alpha}$ are local operators acting on the qubits. They are typically
Pauli matrices (e.g., $\sigma^{z}_{i}$ for flux noise and $\sigma^{y}_{i}$ for
charge noise), but such details are irrelevant for the purpose of this work.
The operators $Q^{\alpha}$ act on the corresponding environments, with
$\alpha$ spanning all sources of noise. We include control error within the
same framework by taking $Q^{\alpha}$ to be a constant random number. For
example, to incorporate static errors $\delta h_{i}$ ($\delta J_{ij}$) in the
Hamiltonian parameters, we use ${\cal O}^{\alpha}=\sigma^{z}_{i}\
(\sigma^{z}_{i}\sigma^{z}_{j})$ and $\sqrt{\gamma}Q^{\alpha}=\delta h_{i}\
(\delta J_{ij})$. The corresponding noise spectral density in that case is a
$\delta$-function at zero frequency. The prefactor $\sqrt{\gamma}$ is
introduced to keep track of the perturbation order; $\gamma$ can be thought of
as an overall decay rate, in units of inverse time, multiplying all relaxation
and dephasing rates. Its precise definition is unimportant for our discussion.

The state of an open quantum system is commonly described by the reduced
density matrix $\rho(t)$, which is the solution to a master equation, e.g.,
Bloch-Redfield [32] or Lindblad [33, 34]. As we show in Appendix A, the
density matrix at the end of the evolution can be expanded in powers of
$\gamma$ as

$\rho(t_{a})=\sum_{\mu=0}^{\infty}(\gamma t_{a})^{\mu}\rho_{\mu}(t_{a})$ (5)

where $\rho_{0}(t_{a})$ is the noise-free contribution to the density matrix
and $\rho_{\mu>0}(t_{a})$ take care of perturbative corrections due to the
environment. As such, any observable $A$ at the end of annealing can be
expanded as

$\langle A(t_{a})\rangle=\sum_{\mu=0}^{\infty}\langle A(t_{a})\rangle_{\mu},$
(6)

where

$\langle A(t_{a})\rangle_{\mu}=\text{Tr}[(\gamma
t_{a})^{\mu}\rho_{\mu}(t_{a})A].$ (7)

The zeroth order term, $\langle A(t_{a})\rangle_{0}$, is the error-free
expectation value for coherent evolution. Our goal is to extract $\langle
A(t_{a})\rangle_{0}$ by extrapolation from a set of noisy evolutions.

ZNE is implemented by amplifying the effect of noise by a controllable factor
$\lambda$ in such a way that $\rho_{0}(t_{a})$ and therefore $\langle
A(t_{a})\rangle_{0}$ remain unchanged. This can be easily achieved by
increasing the coupling coefficient:

$\gamma\to\lambda\gamma.$ (8)

Substituting (8) into (7) and (6) and expanding up to order $M$, we obtain

$\langle A(t_{a},\lambda)\rangle\approx\sum_{\mu=0}^{M}C_{\mu}\lambda^{\mu},$
(9)

where $C_{\mu}=\langle A(t_{a})\rangle_{\mu}$. We can now estimate the error-
free expectation value, $\langle A(t_{a})\rangle_{0}=C_{0}$, by measuring
$\langle A(t_{a},\lambda)\rangle$ for a range of $\lambda$-values, fitting the
result into (9), and extrapolating back to $\lambda\to 0$. When the effect of
the environment is small, only a few terms in the expansion (usually $M=1$ or
2) suffice to yield a good estimate. In practice, it is not feasible to
amplify noise by increasing $\gamma$ because it requires manipulating the
qubits’ environment at the microscopic level. However, it is possible to mimic
such an effect by other means, as we shall discuss next.

##  III Noise amplification by temperature

The simplest way to increase the influence of a thermal environment is by
increasing its temperature $T$. Clearly, the coherent contribution to the
density matrix $\rho_{0}(t_{a})$ remains not affected by changing $T$. As we
show in Appendix A, in the regime of weak coupling to the environment, all
decay processes, including relaxation and dephasing, are functions of the
thermal noise spectral density, which we denote by $S(\omega)$ [35]. For
example, the relaxation rate between states $\ket{n}$ and $\ket{m}$ is
proportional to $S(\omega_{nm})$, where $\omega_{nm}=E_{n}-E_{m}$ is the
energy difference between the two states, and the dephasing rates are
functions of $S(0)$ or spectral density at low frequencies. For an ohmic
environment in thermal equilibrium at temperature $T$, the noise spectral
density is given by ($\hbar=k_{B}=1$)

$S(\omega)={\eta\omega\over 1-e^{-\omega/T}},$ (10)

where $\eta$ is a constant. When the evolution of the system is slow (near
adiabatic) only low energy eigenstates of the instantaneous Hamiltonian get
occupied. Thus, only transitions between these states, for which
$\omega_{nm}<T$, affect the density matrix. We can therefore expand
$e^{-\omega/T}$ in the denominator to first order to obtain
$S(\omega)\approx\eta T$. This means that

$T\to\lambda T\quad\Longrightarrow\quad S(\omega)\to\lambda S(\omega).$ (11)

Since the coefficient $\lambda$ multiplies all decay rates, it can be factored
out and absorbed into the coupling coefficient $\gamma$, leading to (8).
Factoring out $\lambda$ would keep $\rho_{\mu}(t_{a})$ unchanged in (5),
resulting in expansion (9), thus allowing ZNE.

While the above explanation relied on low energy expansion of an Ohmic
spectral density, the conclusion goes beyond such assumptions. As long as the
change in $T$ can be treated as a perturbation, an expansion similar to (9) is
expected, but with different coefficients $C_{\mu}$ and with $C_{0}$
representing the zero-temperature expectation value, which is not necessarily
equivalent the error-free one (see below).

We experimentally tested zero-$T$ extrapolation using a prototype D-Wave
Advantage2TM processor, with 232 flux qubits coupled in a one-dimensional
periodic chain. We applied zero bias $h_{i}=0$ to all qubits and used uniform
ferromagnetic coupling $J_{ij}=-1.8$. We annealed the QPU using fast anneal to
simulate critical dynamics of a transverse-field Ising spin chain. Details of
the experiment are provided in Appendix C and Ref. [29]. We measured kink
density defined as

$n={1\over L}\sum_{i=1}^{L}\langle K_{i}\rangle$ (12)

where $L$ is the length of the chain and

$K_{i}=(1-\sigma^{z}_{i}\sigma^{z}_{i+1})/2$ (13)

is the kink operator. Figure 1 shows $n$ as a function of annealing time
$t_{a}$. The triangle symbols represent QPU results obtained at different
temperatures and the solid lines represent the theoretical prediction,
$n\propto t_{a}^{-1/2}$, according to Jordan-Wigner transformation and DMRG
simulations using the time dependent variational principle (TDVP). (details
about the DMRG simulations are provided in Appendix D). The experimental data
follow the theoretical predictions up to $t_{a}\lesssim 40\,$ns, above which
thermal excitations generate additional kinks [29, 36]. The inset shows data
corresponding to the three vertical cross sections in the main panel,
distinguished by three colors in both plots. Linear extrapolations to $T=0$
for these three cross sections are indicated by black circles in both panels.
The extrapolated points follow the power law $n\propto t_{a}^{-1/2}$, as
predicted by coherent theory, up to $t_{a}\lesssim 200\,$ns and then deviate
sharply from the theoretical prediction. Beyond this point thermal excitations
cannot be treated as perturbation.

Figure 1:  Error mitigation by temperature extrapolation. Main panel: Kink
density as a function of annealing time for four different temperatures. The
gray straight line represents the theoretical prediction for error-free
coherent dynamics $n\propto t_{a}^{-1/2}$. The $T$-dependent deviation from
the prediction is due to thermal excitations. The solid black circles
represent ZNE results at each $t_{a}$. The extrapolated results follow the
theoretical prediction for an extended range of $t_{a}$ compared to
unmitigated data. Inset: extrapolations in $T$ for three annealing times
$t_{a}$ indicated by the three color-coded vertical lines in the main plot.

While zero temperature extrapolation seems to reproduce correct results for
the kink density, it is not expected to correct miscalibration or other
nonthermal errors that do not depend on $T$. These errors can affect other
physical observables such as higher order correlations, as we shall discuss
later. Moreover, our QPU architecture currently does not support rapid changes
in temperature and therefore ZNE in $T$ cannot be regarded as an efficient way
of improving simulation accuracy.

##  IV Noise amplification by energy-time rescaling

Figure 2:  Error mitigation by energy extrapolation. (a) Kink density is
measured for various annealing times $t_{a}$ and coupling energies
$J=-1.8/\kappa$, with $\kappa$ listed in the legend of panel b. An empirical
rescaling of time is performed by collapsing fits of $\rho\propto
t_{a}^{-1/2}$. (b) Collapse of panel a data plotted with rescaled x-axis.
Circles represent the extrapolated points. Inset: ZNE is performed in
$\lambda$, using a quadratic fit, Eq. (16). (c) $\lambda(\kappa)$ obtained
from collapse (circles) compared to those obtained from Eq. (15) (triangles).

A more effective way of amplifying noise is via energy-time rescaling. Since
$\gamma$ is multiplied by $t_{a}$ in (5), it is possible to mimic the effect
of increasing $\gamma$ by increasing $t_{a}$. However, in order to keep the
noiseless unitary evolution $\rho_{0}(t_{a})$ from being affected, one needs
to rescale both anneal time and energy. This can be achieved by substituting
$t_{a}\to\lambda t_{a}$ and $H_{S}\to H_{S}/\lambda$, which keeps
$\int_{0}^{t_{a}}H_{S}(t)dt$ and therefore the unitary time evolution operator
unchanged. In gate model quantum computation, this approach is referred to as
analog-ZNE, in contrast to digital-ZNE, in which noise is amplified by
inserting additional (identity) gates [12]. Intuitively, extending the
annealing time would increase the effect of relaxation and dephasing on the
qubits. Reducing the energy scale, on the other hand, would increase the
relative importance of control errors as long as they remain unchanged by the
rescaling. These all lead to a $\lambda$-dependent reduction of accuracy. It
should be noted that $\rho_{\mu>0}(t_{a})$ may also have some residual
dependence on $\lambda$, but that can be easily incorporated into the Taylor
expansion (9) by redefining the expansion coefficients $C_{\mu}$. This means
that only the invariance of $\rho_{0}(t_{a})$ is required for ZNE to work.

In QA processors, rescaling the system Hamiltonian $H_{S}$ in (1) is
challenging because both $\Gamma(s)$ and ${\cal J}(s)$ are fixed functions.
However, one can easily scale down $H_{P}$ by reprogramming $h_{i}$ and
$J_{ij}$ in (2). Since only one part of $H_{S}$ is now rescaled, making sure
the coherent evolution remains unaffected is nontrivial. Let us introduce an
energy-time rescaling as

$t_{a}\to\lambda t_{a},\qquad H_{P}\to H_{P}/\kappa.$ (14)

We need to choose $\lambda$ and $\kappa$ in such a way that $\rho_{0}(t_{a})$
remains unchanged. In general, this may not be exactly possible, but can be
achieved approximately. When the evolution is dominated by a quantum phase
transition, one can obtain an approximate $\lambda(\kappa)$ in such a way that
the critical dynamics remains unaffected by the rescaling (see Appendix B for
details). This is possible because critical dynamics occur within a narrow
region close to the critical point $s_{c}$ and the statistics only depend on
the speed of passing this point. Rescaling $H_{P}$ would move the critical
point to a new point $s_{c}^{\kappa}\equiv s_{c}(\kappa)$ with a new effective
speed. This change has to be compensated by rescaling the annealing time. For
the case of the transverse-field Ising spin chain studied here, one can show
that (see Appendix B)

$\lambda(\kappa)={\mathcal{J}(s_{c})^{2}[\mathcal{J}^{\prime}(s_{c}^{\kappa})-\Gamma^{\prime}(s_{c}^{\kappa})]\over\mathcal{J}(s_{c}^{\kappa})^{2}[\mathcal{J}^{\prime}(s_{c})-\Gamma^{\prime}(s_{c})]}.$
(15)

where primes indicate first derivatives with respect to $s$. It is also
possible to determine $\lambda(\kappa)$ experimentally as we discuss below.
Once $\lambda$ is determined, we can use (9) to fit the data and extrapolate
to $\lambda=0$. It should be emphasized that $\lambda\to 0$ corresponds to
very short annealing times [37].

Figure 3:  Mitigation of kink-kink correlation by energy extrapolation. We
compute the observable ${\cal C}(r)$ by (a) fermionizing the system via
Jordan-Wigner transformation, and (b) using numerics with time dependent DMRG
using the TDVP algorithm (see Appendix D for details). Points fall
approximately on a non-universal curve for a range of $t_{a}$. Red line
indicates a degree-10 polynomial that is replicated on the other panels to the
right to facilitate comparison between theory and experiment. (c) Unmitigated
QA data show correct qualitative behavior, but with a significantly flattened
peak speculated to arise primarily due to unintended Hamiltonian disorder. (d)
Results obtained via linear ZNE to $\lambda=0$ show close agreement with the
fermionized solution. (e) We show example extrapolations for varying lattice
distances $r$, with $t_{a}=7.34$ ns.

For the case of thermal noise, the environment does not affect the system when
the annealing time is shorter than the thermal relaxation time, hence no
$\lambda$-dependence. As such, we should expect a flat curve with zero slope
at small $\lambda$, i.e., $C_{1}=0$. This means that for thermal noise, the
lowest order polynomial is quadratic:

$\langle A(t_{a},\lambda)\rangle\approx C_{0}+C_{2}\lambda^{2}.$ (16)

Figure 2a shows measurements of kink density for different values of
$J=-1.8/\kappa$. As expected, reducing the energy scale increases the number
of kinks. Moreover, it increases the relative influence of the thermal
environment as evidenced by the increasingly obvious deviation from $n\propto
t_{a}^{-1/2}$ at high $t_{a}$. Nonetheless, it is possible to rescale the
horizontal axis $t_{a}\to t_{a}/\lambda$, using $\lambda$ as a free parameter,
to collapse all of the data within the coherent annealing regime. The
resulting collapse is shown in Fig. 2b and the values of $\lambda(\kappa)$
obtained by this collapse are plotted in Fig. 2c together with theoretical
predictions using (15). Note that the collapse method cannot mitigate static
errors, such as Hamiltonian misspecification because they affect both coherent
and thermally activated evolution. The agreement between the empirical
collapse and the annealing schedule inferred theoretical values of
$\lambda(\kappa)$ in the inset suggests that kink-density is insensitive to
such errors.

The next step involves using extrapolation to obtain a noise-free kink density
in the regime where the thermal environment is effective. Since the dominant
noise is thermal, one should extrapolate using quadratic fits per Eq. (16).
The inset in Fig. 2b illustrates this extrapolation for three different values
of $t_{a}/\lambda$. The black circles in the main panel show the result of
this extrapolation. As in Fig. 1, the extrapolated points agree with the exact
$n\propto t_{a}^{-1/2}$ behavior for up to $t_{a}/\lambda\lesssim 200\,$ns.
Thus one can successfully infer noise-free expectation values over a much
broader range of annealing time than the ostensible coherent annealing regime
that is nominally bounded by $t_{a}\lesssim 40\,$ns.

While kink density is unaffected by control error, as the inset of Fig. 2a
suggests, other quantities may not exhibit this resilience. One such quantity
is kink-kink correlation [38, 39]

${\cal C}(r)={1\over L}\sum_{i=1}^{L}{\langle K_{i}K_{i+r}\rangle-n^{2}\over
n^{2}},$ (17)

with $K_{i}$ defined in Eq. (13). Figure 3 shows ${\cal C}(r)$ as a function
of $r$ normalized to the correlation length $\xi=1/n$. Figures 3a and 3b
illustrate results obtained from the exact solution using Jordan-Wigner
transformation and DMRG, respectively. The measured ${\cal C}(r)$ depicted in
Fig. 3c exhibits a systematic deviation from the exact results in Fig. 3a,
independent of $t_{a}$. Such a deviation was first observed in Ref. [29] and
was attributed to the disorder due to control error or $1/f$ noise. Disorder
in the potential energy can attract kinks toward local minima, thus affecting
the kink-kink correlation while keeping the density of kinks unchanged. As
discussed before, the effect of disorder and control error can be mitigated by
ZNE through energy-time rescaling. Figure 3d illustrates the ZNE result,
showing a significant improvement in alignment with the exact theoretical
predictions (depicted by the solid lines), compared to unmitigated results. In
Fig. 3e, a few examples of extrapolation at different values of r are
displayed and, unlike the case of thermal noise (Fig. 2), linear extrapolation
works well.

##  V Conclusion

We have successfully demonstrated the implementation of Zero-Noise
Extrapolation (ZNE) for quantum annealing. Through the analysis of critical
dynamics in a 1D quantum spin chain, we mitigated the impact of thermal noise
on kink density, extending the effective decoherence time by almost an order
of magnitude. Additionally, we mitigated the influence of static errors
possibly due to mis-calibration or $1/f$ noise, on coherent evolution, as
evidenced by having successfully estimated kink-kink correlations. All
extrapolated results are in good agreement with the exact solution to the
Schrödinger equation obtained from the Jordan-Wigner transformation as well as
DMRG.

It is important to emphasize that ZNE primarily enhances the accuracy of
expectation values but does not directly impact sample quality. Nevertheless,
one can employ extrapolation of probabilities to gain insights into error-free
probability distribution. This allows for the reconstruction of the
distribution by adjusting sample weights or employing Monte Carlo techniques.
For some problems this may require a significant amount of sampling.
Extrapolating probabilities is particularly useful when evaluating measures
such as time-to-solution [44] for increasingly lower noise and more coherent
systems, hence constructing a performance model for the QPU and providing a
means to specify requirements on coherence and noise levels. When the
objective involves optimizing for a lower energy solution, QAC techniques [14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24] may represent more appropriate
approach to enhance overall performance. The technique described in this work
can be applied to tackle more intricate challenges, such as quantum
simulations of exotic magnetic materials. Leveraging ZNE makes it feasible to
estimate expectation values in situations where classical computation
approaches become intractable.

### Acknowledgment

We would like to thank Gonzalo Alvarez, Jacek Dziarmaga, Daniel Lidar,
Hidetoshi Nishimori, and Marek Rams for fruitful discussions and comments on
the manuscript. Work at UBC was supported by NSERC Alliance Quantum Program
(Grant ALLRP-578555), CIFAR and the Canada First Research Excellence Fund,
Quantum Materials and Future Technologies Program.

##  Appendix A Open quantum model

In this appendix we use open quantum modeling to explain ZNE. We write the
total Hamiltonian as

$H=H_{S}+H_{B}+H_{\rm int},$ (18)

where $H_{S}$ and $H_{B}$ are system and bath Hamiltonians, respectively. We
write the interaction Hamiltonian in a general form as

$H_{\rm int}=-\sqrt{\gamma}\sum_{\alpha}Q^{\alpha}{\cal O}^{\alpha}$ (19)

where ${\cal O}^{\alpha}$ is an operator acting on the system and $Q^{\alpha}$
is the noise operator acting on the corresponding environment. The operators
${\cal O}^{\alpha}$ are typically Pauli matrices, e.g., $\sigma^{z}_{\alpha}$
for flux noise and $\sigma^{y}_{\alpha}$ for charge noise. For static noise,
such as control error, $Q^{\alpha}$ can be considered as a constant random
number.

In open quantum systems, the state of the system is described by the reduced
density matrix $\rho(t)$. We choose the preferred basis to be the eigenstates
$\ket{n}$ of the system Hamiltonian $H_{S}$ with eigenvalues $E_{n}$. The
Bloch-Redfield equation for the reduced density matrix $\rho(t)$ is written as
[32]

$\displaystyle\dot{\rho}_{nm}(t)=-i\omega_{nm}\
\rho_{nm}(t)-\gamma\sum_{k,l}R_{nmkl}\ \rho_{kl}(t)$ (20)

where $\omega_{nm}=E_{n}-E_{m}$ and the Redfield tensor is defined as

$\displaystyle R_{nmkl}$ $\displaystyle=$
$\displaystyle\delta_{lm}\sum_{r}\Gamma^{(+)}_{nrrk}+\delta_{nk}\sum_{r}\Gamma^{(-)}_{lrrm}$
(21) $\displaystyle-\Gamma^{(+)}_{lmnk}-\Gamma^{(-)}_{lmnk}$

with

$\displaystyle\Gamma^{(+)}_{lmnk}={1\over 2}\sum_{\alpha,\beta}{\cal
O}^{\alpha}_{lm}{\cal O}^{\beta}_{nk}S_{\alpha\beta}(-\omega_{nk}),$ (22)
$\displaystyle\Gamma^{(-)}_{lmnk}={1\over 2}\sum_{\alpha,\beta}{\cal
O}^{\alpha}_{lm}{\cal O}^{\beta}_{nk}S_{\alpha\beta}(\omega_{lm}),$ (23)

Here, ${\cal O}^{\alpha}_{nm}\equiv\bra{n}{\cal O}^{\alpha}\ket{m}$ and the
noise spectral density is defined as [41]

$\displaystyle S_{\alpha\beta}(\omega)=\int_{-\infty}^{\infty}dt\ e^{i\omega
t}\langle Q^{\alpha}(t)Q^{\beta}(0)\rangle.$ (24)

For the case of uncorrelated heat baths we have

$S_{\alpha\beta}(\omega)=\delta_{\alpha\beta}S_{\alpha}(\omega).$ (25)

For a time-dependent system Hamiltonian, the generalized Bloch-Redfield
equation becomes [42]

$\displaystyle\dot{\rho}_{nm}$ $\displaystyle=$
$\displaystyle-i\omega_{nm}\rho_{nm}-\sum_{kl}\left(\gamma
R_{nmkl}+M_{nmkl}\right)\rho_{kl}.$ (26)

where

$\displaystyle M_{nmkl}=-\delta_{nk}\langle
l|\dot{m}\rangle-\delta_{ml}\langle\dot{n}|k\rangle.$ (27)

In the energy basis, we write the density matrix as a linear vector with
$2^{2N}$ elements:

$\hat{\rho}=\left[\begin{array}[]{c}\rho_{11}\\\ .\\\ .\\\ \rho_{nm}\\\ .\\\
.\end{array}\right]$ (28)

The master equation becomes a matrix equation:

${d\over
dt}\hat{\rho}(t)=[-i\hat{\omega}(t)-\hat{M}(t)-\gamma\hat{R}(t)]\hat{\rho}(t)$
(29)

where

$\hat{\omega}(t)=\left[\begin{array}[]{ccccccc}\omega_{11}(t)&&&&&&\\\
&.&&&&&\\\ &&.&&&&\\\ &&&\omega_{nm}(t)&&&\\\ &&&&.&&\\\
&&&&&.&\end{array}\right],$ (30)

is a diagonal matrix made of energy differences and $\hat{M}$ and $\hat{R}$
are the basis-rotation and relaxation matrices corresponding to $M_{nmkl}$ and
$R_{nmkl}$, respectively. A formal solution to this equation is

$\hat{\rho}(t)=\hat{C}_{1}+{\cal
T}e^{-\int_{0}^{t}[i\hat{\omega}(\tau)+\hat{M}(\tau)+\gamma\hat{R}(\tau)]d\tau}\hat{C}_{2}.$
(31)

The time-ordering operator ${\cal T}$ is introduced to take care of non-
commuting $\hat{\gamma}(t)$, $\hat{M}(t)$, and $\hat{R}(t)$ at different
times. Changing the integration variable to $s=\tau/t_{a}$, we get

$\hat{\rho}(t_{a})=\hat{C}_{1}+{\cal
T}e^{-t_{a}\int_{0}^{1}[i\hat{\omega}(s)+\hat{M}(s)+\gamma\hat{R}(s)]ds}\hat{C}_{2}.$
(32)

The term $\gamma\hat{R}(s)$ in the exponent take care of all relaxation and
dephasing processes during the evolution. In the limit of small $\gamma$ when
those decay processes are not very strong it is possible to assume this term
is small and expand the exponent. Using Taylor expansion, we can write

$\hat{\rho}(t_{a})=\hat{\rho}_{0}(t_{a})+\sum_{\mu=1}^{\infty}(\gamma
t_{a})^{\mu}\,\hat{\rho}_{\mu}(t_{a}),$ (33)

where

$\hat{\rho}_{0}(t_{a})=\hat{C}_{1}+{\cal
T}e^{-t_{a}\int_{0}^{1}[i\hat{\omega}(s)+\hat{M}(s)]ds}\hat{C}_{2}$ (34)

is the noise-free contribution to the density matrix and

$\hat{\rho}_{\mu}(t_{a})=\Bigg{[}{d^{\mu}\over\mu!(t_{a}d\gamma)^{\mu}}{\cal
T}e^{-t_{a}\int_{0}^{1}[i\hat{\omega}(s)+\hat{M}(s)+\gamma\hat{R}(s)]ds}\bigg{]}_{\gamma=0}\hat{C}_{2},$
(35)

captures the effect of noise to order $\mu$ in perturbation. Turning back the
density matrices from vector form to matrix form we obtain (5).

In some cases, only a small number of energy states get occupied during the
annealing process. Therefore, we only need to incorporate relaxation and
dephasing processes between those states. The experiment reported in Ref. [31]
is one such case. The evolution is limited to the lowest two energy states
between which the relaxation process is extremely slow. This allows expansion
in powers of those relaxation rates only, although other decay processes could
be much faster. As a result zero-T extrapolation worked for time scales
millions of times longer than the expected single qubit decoherence time of
the qubits.

##  Appendix B Energy-time rescaling

The goal of this section is to introduce a relation between energy and time
rescaling so that the critical dynamics remains unchanged. We first write the
time-dependent Hamiltonian in a dimensionless form with one dimensionless
parameter $\tau_{Q}$, which measures the speed of passing the critical point.
The same quantity is used in some theory papers (see e.g. [38, 39]), hence our
formalism would allow direct comparison with these analytical results.

We start by writing the Hamiltonian (1) as

$H(s)=\mathcal{J}(s)\bigg{(}\\!\\!-g(s)\sum_{i}\sigma_{i}^{x}+H_{P}\bigg{)},$
(36)

where

$g(s)={\Gamma(s)\over\mathcal{J}(s)}$ (37)

is the dimensionless transverse field. The corresponding Schrodinger equation
reads

$i{d\over
dt}\ket{\psi(t)}=\mathcal{J}(s)\bigg{(}\\!\\!-g(s)\sum_{i}\sigma_{i}^{x}+H_{P}\bigg{)}\ket{\psi(t)}.$
(38)

Introducing a dimensionless time

$\tilde{t}=\int_{0}^{t}\mathcal{J}(t^{\prime})dt^{\prime}=t_{a}\int_{0}^{s}\mathcal{J}(s^{\prime})ds^{\prime},$
(39)

we can write (38) in dimensionless form

$i{d\over
d\tilde{t}}\ket{\psi(\tilde{t})}=\tilde{H}(\tilde{t})\ket{\psi(\tilde{t})}$
(40)

where

$\tilde{H}(\tilde{t})=-g(\tilde{t})\sum_{i}\sigma_{i}^{x}+H_{P}$ (41)

is the dimensionless Hamiltonian. Note that $\tilde{t}\in[0,\tilde{t}_{a}]$
with $\tilde{t}_{a}=t_{a}\int_{0}^{1}\mathcal{J}(s)ds$ being the dimensionless
annealing time.

For fast annealing, only the speed of passing the critical point determines
the critical dynamics. We define the critical point $s=s_{c}$ by

$\Gamma(s_{c})/\mathcal{J}(s_{c})=g_{c}$ (42)

where $g_{c}$ is the critical value of $g$. Using linear expansion near the
critical point we write

$g(\tilde{t})\approx g_{0}-{\tilde{t}\over\tau_{Q}}$ (43)

where $\tau_{Q}$ is the dimensionless quench time scale. It is possible to
express Kibble-Zurek dynamics only in terms of $\tau_{Q}$. In other words, two
systems with the same $\tau_{Q}$ are expected to yield the same statistics as
long as only critical dynamics affect their evolution. For a slower evolution
where the dynamics outside the critical region affect the results, details of
the schedule become important.

We expand the schedule near the critical point as

$\displaystyle\Gamma(s)$ $\displaystyle=$
$\displaystyle\Gamma(s_{c})+\Gamma^{\prime}(s_{c})(s-s_{c})$ (44)
$\displaystyle\mathcal{J}(s)$ $\displaystyle=$
$\displaystyle\mathcal{J}(s_{c})+\mathcal{J}^{\prime}(s_{c})(s-s_{c}).$ (45)

Expanding $g(s)$ near $s_{c}$, we obtain

$\displaystyle g(s)$ $\displaystyle\approx$ $\displaystyle
g_{c}+{\Gamma^{\prime}(s_{c})\mathcal{J}(s_{c})-\Gamma(s_{c})\mathcal{J}^{\prime}(s_{c})\over\mathcal{J}^{2}(s_{c})}(s-s_{c})$
(46) $\displaystyle\approx$ $\displaystyle
g_{c}+g_{c}\bigg{[}{\Gamma^{\prime}(s_{c})\over\Gamma(s_{c})}-{\mathcal{J}^{\prime}(s_{c})\over\mathcal{J}(s_{c})}\bigg{]}(s-s_{c})$
$\displaystyle\approx$ $\displaystyle
g_{c}\bigg{[}{\Gamma^{\prime}(s_{c})\over\Gamma(s_{c})}-{\mathcal{J}^{\prime}(s_{c})\over\mathcal{J}(s_{c})}\bigg{]}{t\over
t_{a}}+const.$

On the other hand

$\tilde{t}\approx\mathcal{J}(s_{c})t+const.$ (47)

Therefore

$g(\tilde{t})\approx
g_{c}\bigg{[}{\Gamma^{\prime}(s_{c})\over\Gamma(s_{c})}-{\mathcal{J}^{\prime}(s_{c})\over\mathcal{J}(s_{c})}\bigg{]}{\tilde{t}\over\mathcal{J}(s_{c})t_{a}}+const.$
(48)

Comparing with (43), we obtain

$\tau_{Q}={t_{a}\over t_{Q}},$ (49)

where

$t_{Q}={g_{c}\over\mathcal{J}(s_{c})}\bigg{[}{\mathcal{J}^{\prime}(s_{c})\over\mathcal{J}(s_{c})}-{\Gamma^{\prime}(s_{c})\over\Gamma(s_{c})}\bigg{]}$
(50)

is a time scale that only depends on the annealing schedule. For 1D spin
chain, we have $g_{c}=1$ and the critical point is the point where
$J(s_{c})=\Gamma(s_{c})$. On the other hand, for 2D and 3D problems we have
$g_{c}>1$ and therefore we expect the critical point to occur earlier in the
schedule.

In order for the two systems to have the same critical dynamics, they need to
have the same $\tau_{Q}$. Equation (49) requires that $t_{a}$ should be scaled
the same way as $t_{Q}$ in (50). If both $\Gamma(s)$ and ${\cal J}(s)$ are
scaled the same way such that $H_{S}\to H_{S}/\kappa$, we obtain
$t_{Q}\to\kappa t_{Q}$, therefore we need $\lambda=\kappa$, as expected for
simple rescaling.

For the more practical case of only rescaling $H_{P}$, it is equivalent ${\cal
J}(s)\to{\cal J}(s)/\kappa$ without changing $\Gamma(s)$. Obviously, the
critical point will be shifted to a new point $s_{c}^{\kappa}$ which is the
solution to $\Gamma(s_{c}^{\kappa})=g_{c}\kappa{\cal J}(s_{c}^{\kappa})$.
Thus, one needs to calculate the new $t_{Q}^{\kappa}$, obtained from (50) at
$s_{c}^{\kappa}$. The correct time rescaling coefficient will therefore be:

$\lambda(\kappa)=t_{Q}^{\kappa}/t_{Q}.$ (51)

For case of the 1D transverse field Ising problem, we can write

$\lambda(\kappa)={\mathcal{J}(s_{c})^{2}[\mathcal{J}^{\prime}(s_{c}^{\kappa})-\Gamma^{\prime}(s_{c}^{\kappa})]\over\mathcal{J}(s_{c}^{\kappa})^{2}[\mathcal{J}^{\prime}(s_{c})-\Gamma^{\prime}(s_{c})]}.$
(52)

##  Appendix C Experimental details

Figure 4: QPU annealing schedule. $\Gamma(s)$ and $\mathcal{J}(s)$ are
estimated using single-qubit measurements, then both are rescaled by an
overall fit parameter $1/1.15$ to account for any deviation between single-
qubit measurements and many-body flux qubit behavior. The y-axis units is GHz
$=10^{-9}$ Joules/$h$, where $h$ is the Planck constant. For equations with
$\hbar=1$ convention, these numbers should be multiplied by $2\pi$.

Quantum annealing experiments were performed using a prototype D-Wave
Advantage2 processor, with 232 flux qubits coupled in a one-dimensional
periodic chain. At each temperature, coupling strength, and annealing time,
calibration refinement was performed to homogenize coupler frustration,
balance each qubit at zero magnetization, and synchronize annealing schedules.
These refinements are performed by tuning, respectively, individual couplings,
flux-bias offsets, and anneal offsets as previously detailed in the
supplementary materials of [29] and as described in [40] (the latter does not
employ or describe anneal offset tuning, as it is only impactful for anneals
faster than 100 ns).

Average kink densities are measured over 20 programmings; 100 samples are
taken in each programming. Kink-kink correlations require more robust
statistics, and are measured over 1000 programmings; 1000 samples are taken in
each programming. All error bars represent 95% confidence intervals generated
from bootstrap resampling of the 20 or 1000 programmings.

Qubit temperatures are measured by the population slope for uncoupled qubits
subjected to a varying longitudinal field.

Zero-noise extrapolation was performed using linear extrapolation in $T$ (Fig.
1), quadratic extrapolation in $\lambda$ (Fig. 2), and linear extrapolation in
$\lambda$ (Fig. 3) to zero. Relative values of $\lambda$ were determined from
an estimated annealing schedule shown in Fig. 4, which was derived from
single-qubit measurements; rather than refining the single-qubit schedule with
detailed many-body modeling, as in previous works [29, 30], we simply
introduced a single parameter and multiplied both $\Gamma(s)$ and
$\mathcal{J}(s)$ by $1/1.15$ to align measured and simulated kink densities in
the 5-10 ns region.

##  Appendix D Details about DMRG simulations

For kink density as well as kink-kink correlations, we compared QA output
against results from exact solution via Jordan-Wigner transformation and
against DMRG, with dynamics simulated by the time-dependent variational
principle (TDVP) with a two site update [43, 45] using ITensor library [46].
We also compared the results against the $W_{\text{II}}$ method introduced in
Ref. [47] using the TenPy library [48] and the local Krylov method [49] (for
details about this method in MPS language see Ref. [50]) implemented by the
DMRG++ library [51]. Overall, we found that TDVP emerges as the most efficient
method, providing converged results using the largest time step ($dt$=0.01
ns). For the lattice geometry, we adopted the same mapping from periodic to
open chain used in Ref. [29] as shown in Fig. 5a. We found a much better
performance of the numerical simulations when using a local Hilbert space of
two qubits (see Fig. 5b) for the MPS wave-functions.

Figure 5:  (a) Site ordering for MPS simulations of a periodic chain with one-
qubit local Hilbert space. (b) Site ordering with a two-qubit local Hilbert
space.

For the results shown in the main part of the manuscript, TDVP simulations
were performed with a bond dimension $D=32$ and a time step of 0.01 ns, with a
maximum truncation error of $10^{-10}$.

## References

  * [1] P. W. Shor, Scheme for reducing decoherence in quantum computer memory, Phys. Rev. A 52, 2493 (1995).
  * [2] A. M. Steane, Error correcting codes in quantum theory, Phys. Rev. Lett. 77, 793 (1996).
  * [3] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * [4] Temme, K., Bravyi, S., and Gambetta, J. M. Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * [5] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * [6] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Phys. Rev. X 8, 031027 (2018).
  * [7] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum Error Mitigation, eprint: arXiv:2210.00921.
  * [8] A. Kandala, K.Temme, E. Novais, S. Woerner, and J. M. Gambetta, Error mitigation extends the computational reach of a noisy quantum processor, Nature 567 491 (2019).
  * [9] E. van den Berg et al., Probabilistic error cancellation with sparse Pauli–Lindblad models on noisy quantum processors, Nat. Phys. 19, 1116 (2023).
  * [10] Y. Kim et al., Scalable error mitigation for noisy quantum circuits produces competitive expectation values, Nat. Phys. 19, 752 (2023).
  * [11] Y. Kim et al., Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * [12] Ritajit Majumdar, Pedro Rivero, Friederike Metz, Areeq Hasan, and Derek S. Wang, Best practices for quantum error mitigation with digital zero-noise extrapolation, eprint: arXiv:2307.05203.
  * [13] M. W. Johnson et al., M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, _Quantum annealing with manufactured spins_ , Nature 473, 194 (2011).
  * [14] K.C. Young, R. Blume-Kohout, D.A. Lidar, Adiabatic Quantum Optimization with the Wrong Hamiltonian, Phys. Rev. A 88, 062314 (2013).
  * [15] K.P. Pudenz, T. Albash, and D.A. Lidar, Error Corrected Quantum Annealing with Hundreds of Qubits, Nature Communications 5, 3243 (2014).
  * [16] K. Pudenz, T. Albash, and D. Lidar, Quantum Annealing Correction for Random Ising Problems, Phys. Rev. A 91, 042302 (2015).
  * [17] W. Vinci, T. Albash, G. Paz-Silva, I. Hen, and D. A. Lidar, Quantum Annealing Correction with Minor Embedding, Phys. Rev. A 92, 042310 (2015).
  * [18] A. Mishra, T. Albash and D.A. Lidar, Performance of two different quantum annealing correction codes, Quant. Info. Proc. 15, 2, pp. 609 (2016).
  * [19] S. Matsuura, H. Nishimori, T. Albash, D.A. Lidar, Mean Field Analysis of Quantum Annealing Correction, Phys. Rev. Lett. 116, 220501 (2016).
  * [20] W. Vinci, T. Albash, and D. A. Lidar, Nested Quantum Annealing Correction, Nature Quant. Info. 2, 16017 (2016).
  * [21] S. Matsuura, H. Nishimori, W. Vinci, T. Albash, and D. A. Lidar, Quantum annealing correction at finite temperature: ferromagnetic p-spin models, Phys. Rev. A 95, 022308 (2017).
  * [22] W. Vinci and D. A. Lidar, Scalable effective temperature reduction for quantum annealers via nested quantum annealing correction, Phys. Rev. A 97, 022308 (2018).
  * [23] S. Matsuura, H. Nishimori, W. Vinci, D. A. Lidar, Nested Quantum Annealing Correction at Finite Temperature: p-spin models, Phys. Rev. A 99, 062307 (2019).
  * [24] D. A. Lidar, Arbitrary-Time Error Suppression for Markovian Adiabatic Quantum Computing Using Stabilizer Subspace Codes, Phys. Rev. A 100, 022326 (2019).
  * [25] R. Harris et al., Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi,R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, J. Yao, _Phase transitions in a programmable quantum spin glass simulator,_ Science 361, 162 (2018).
  * [26] A. D. King et al., J. Carrasquilla, J. Raymond, I. Ozfidan, E. Andriyash, A. Berkley, M. Reis, T. M. Lanting, R. Harris, F. Altomare, K. Boothby, P. I. Bunyk, C. Enderud, A. Frechette, E. Hoskinson, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich, Y. Sato, A. Yu. Smirnov, L. J. Swenson, M. H. Volkmann, J. Whittaker, J. Yao, E. Ladizinsky, M. W. Johnson, J. Hilton, and M. H. Amin, _Observation of topological phenomena in a programmable lattice of 1,800 qubits_ , Nature 560, 456 (2018).
  * [27] A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, A. Lopez-Bezanilla, Qubit spin ice Science 373, 576 (2021).
  * [28] A. Lopez-Bezanilla, J. Raymond, K. Boothby, J. Carrasquilla, C. Nisoli, and A. D. King, Kagome qubit ice, arXiv:2301.01853.
  * [29] Andrew D. King, Sei Suzuki, Jack Raymond, Alex Zucca, Trevor Lanting, Fabio Altomare, Andrew J. Berkley, Sara Ejtemaee, Emile Hoskinson, Shuiyuan Huang, Eric Ladizinsky, Allison MacDonald, Gaelen Marsden, Travis Oh, Gabriel Poulin-Lamarre, Mauricio Reis, Chris Rich, Yuki Sato, Jed D. Whittaker, Jason Yao, Richard Harris, Daniel A. Lidar, Hidetoshi Nishimori, and Mohammad H. Amin, Coherent quantum annealing in a programmable 2,000 qubit Ising chain, Nat. Phys. 18, 1324 (2022).
  * [30] Andrew D. King, Jack Raymond, Trevor Lanting, Richard Harris, Alex Zucca, Fabio Altomare, Andrew J. Berkley, Kelly Boothby, Sara Ejtemaee, Colin Enderud, Emile Hoskinson, Shuiyuan Huang, Eric Ladizinsky, Allison J. R. MacDonald, Gaelen Marsden, Reza Molavi, Travis Oh, Gabriel Poulin-Lamarre, Mauricio Reis, Chris Rich, Yuki Sato, Nicholas Tsai, Mark Volkmann, Jed D. Whittaker, Jason Yao, Anders W. Sandvik, and Mohammad H. Amin, Quantum critical dynamics in a 5000-qubit programmable spin glass, Nature 617, 61 (2023).
  * [31] N. G. Dickson et al., Thermally assisted quantum annealing of a 16-qubit problem, Nat. Commun. 4 1903 doi: 10.1038/ncomms2920 (2013).
  * [32] K. Blum, Density Matrix Theory and Applications, (Plenum, New York, 1981).
  * [33] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976).
  * [34] Evgeny Mozgunov1 and Daniel Lidar, Completely positive master equation for arbitrary driving and small level spacing, Quantum 4, 227 (2020).
  * [35] In general, different qubits may be coupled to different heat baths with different spectral densities, but these details don’t affect our argument. See Appendix A for more details.
  * [36] Bando, Y. et al. Probing the universality of topological defect formation in a quantum annealer: Kibble-Zurek mechanism and beyond, Phys. Rev. Res. 2, 033369 (2020).
  * [37] We stress that experimentally attainable annealing times are away from the regime of fast quenches, where deviations from Kibble-Zurek mechanism have been predicted in: H.-B. Zeng, C.-Y. Xia, and A. del Campo, Phys. Rev. Lett. 130 060402 (2023).
  * [38] R.J. Nowak and J. Dziarmaga, Quantum Kibble-Zurek mechanism: Kink correlations after a quench in the quantum Ising chain, Phys. Rev. B 104, 075448 (2021).
  * [39] Jacek Dziarmaga and Marek M. Rams, Kink correlations, domain size distribution, and the emptiness formation probability after the Kibble-Zurek quench in the quantum Ising chain, Phys. Rev. B 106, 014309 (2022)
  * [40] Kevin Chern, Kelly Boothby, Jack Raymond, Pau Farré, Andrew D. King, Tutorial: Calibration refinement in quantum annealing, Frontiers in Computer Science 5 (2023).
  * [41] In the master equation, environmental normalization of energy differences (Lamb Shift) is ignored. One can incorporate the effect of high frequency noise in the renormalization of the Hamiltonian parameters at the single qubit level. Starting from the renormalized Hamiltonian, based on measured parameters, no additional Lamb shift due to high frequency noise should be introduced. See also Ref. [34].
  * [42] M. H. S. Amin, C. J. S. Truncik, and D. V. Averin, Role of single-qubit decoherence time in adiabatic quantum computation, Phys. Rev. A 80, 022303 (2009).
  * [43] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107 070601 (2011).
  * [44] 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 (2014).
  * [45] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94 165116 (2016).
  * [46] M. Fishman and S. R. White and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases, (2022).
  * [47] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015).
  * [48] J. Hauschild and F. Pollmann, Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), SciPost Phys. Lect. Notes, 5 (2018).
  * [49] G. Alvarez, L. G. G. V. Dias da Silva, E. Ponce, and E. Dagotto, Time evolution with the density-matrix renormalization-group algorithm: A generic implementation for strongly correlated electronic systems, Phys. Rev. E 84, 056706 (2011).
  * [50] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019).
  * [51] G. Alvarez, The density matrix renormalization group for strongly correlated electron systems: A generic implementation, Computer Physics Communications 180, 1572 (2009).