Quantum Annealing: An

Overview
[Link] Atanu Rajak1 , Sei Suzuki2 , Amit Dutta3
and Bikas K. Chakrabarti4,5

Research 1 Presidency University, Kolkata -700073, India

2 Saitama Medical University, Moroyama, Saitama

350-0495, Japan
arXiv:2207.01827v4 [[Link]-mech] 17 Jan 2023

Article submitted to journal 3 Indian Institute of Technology Kanpur, 208016, India

4 Saha Institute of Nuclear Physics, 1/AF Bidhannagar,

Subject Areas: Kolkata-700064, India

quantum annealing and optimisation, 5 Indian Statistical Institute, 203 B. T. Road,
quantum phase transitions, Kolkata-700108, India
disordered systems, open quantum
systems
In this review, after providing the basic physical
concept behind quantum annealing (or adiabatic
Keywords: quantum computation), we present an overview
quantum tunnelling, Kibble-Zurek of some recent theoretical as well as experimental
Scaling, transverse Ising models, developments pointing to the issues which are still
quantum spin glass, Griffiths debated. With a brief discussion on the fundamental
ideas of continuous and discontinuous quantum
singularities, p-spin models,
phase transitions, we discuss the Kibble-Zurek scaling
decoherence, antiferromagnetic
of defect generation following a ramping of a
interactions, NP-complete and quantum many body system across a quantum
NP-hard problems critical point. In the process, we discuss associated
models, both pure and disordered, and shed light
on implementations and some recent applications
Author for correspondence:
of the quantum annealing protocols. Furthermore,
Atanu Rajak
we discuss the effect of environmental coupling on
e-mail: [Link]@[Link] quantum annealing. Some possible ways to speed
up the annealing protocol in closed systems are
elaborated upon: We especially focus on the recipes
to avoid discontinuous quantum phase transitions
occurring in some models where energy gaps vanish
exponentially with the system size.

© The Authors. Published by the Royal Society under the terms of the
Creative Commons Attribution License [Link]
by/4.0/, which permits unrestricted use, provided the original author and
source are credited.
1. Introduction 2
Following the recent technological advance in manipulation of a quantum state, the notion of

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
quantum computation and simulation which initially stemmed from pure theoretical concepts
has now spread to flourish an industry with an immense possibility of technological applications.
In particular, studies of quantum annealing (QA) have gained a tremendous momentum since
programmable QA machines, dubbed as quantum annealers, with more than thousands of qubits
have been realized and commercialized. In this review, having provided an overview of QA
protocol, we discuss some recent theoretical and experimental developments of the QA exploiting
the advantage of utilizing quantum tunneling in finding the minimum of a classical energy
function.
QA is usually aimed at seeking the ground state of a generic Ising model, which may contain
random biases and/or random many-body interactions [1–3]. Many optimization problems
including the traveling salesman problem, job scheduling problem, knapsack problem, and so
on, are shown to reduce to this problem. Therefore the application of QA extends from physics
to our daily life. This broadness of application is another reason why QA is attracted much
attention in industry. Now, let us consider an Ising model denoted by the Hamiltonian HP , where
the subscript P stands for the problem Hamiltonian. We assume that HP is a classical many-
body Ising Hamiltonian described in terms of the z components of the Pauli operator {σjz }. We
further introduce a driver Hamiltonian HD which is not commutative with HP and has the trivial
ground state. A simple choice for HD is the transverse field: HD = − j σjx , so that HD does not
P

commute with HP . The total Hamiltonian of QA is given as

H(t) = A(t)HD + B(t)HP , (1.1)

where A(t) and B(t) are the scheduling function satisfying A(ti )  B(ti ) at the initial time ti and
A(tf )  B(tf ) at the final time tf so that H(t) interpolates between HD at t = ti and HP at t = tf .
The initial state at t = ti is set at the ground state of HD ≈ H(ti )/A(ti ). If the change in H(t) with
t is “sufficiently” small, the spin state evolves adiabatically (i.e., stays in the ground state of the
instantaneous Hamiltonian) and arrives at the ground state of HP at t = tf which we seek. This
constitutes the basic notion of the QA, also known as the adiabatic quantum computation [4–10].
Throughout this paper, we shall employ QA scheme using the transverse Ising Hamiltonian (if not
otherwise mentioned). To illustrate, we consider the following Hamiltonian with ferromagnetic
nearest neighbour interactions in one dimension:
X z z X x
H = −J σj σj+1 − Γ σj . (1.2)
j j

where J denotes the strength of the interaction and Γ is the strength of the non-commuting
transverse field. Here HD = − j σjx and HP = − j σjz σj+1 z
P P
. The transverse field Γ is annealed
to reach the ground state of HP from the ground state of HD .
The success of QA is determined by how slowly the Hamiltonian changes with time.
According to the adiabatic theorem of quantum mechanics, the criterion of the adiabatic time
evolution is given by [11] h i
dH(t)
max |h1(t)| dt |g(t)i|
 1, (1.3)
min[∆(t)]2
where |g(t)i and |1(t)i are the instantaneous ground and first-excited states at time t, respectively,
and ∆(t) denotes the instantaneous energy gap above |g(t)i. The min and max functions are taken
with respect to the variable t. Thus, roughly speaking, QA works better for larger ∆(t) [12].
As a classical counterpart to QA, simulated annealing (SA) is a known method of computation
for optimization problems [13]. In this method, we prepare the Gibbs-Boltzmann distribution of
HP at sufficiently high temperature by means of the Monte-Carlo method and literally anneal
the system down to zero temperature. If annealing is sufficiently slow, then we are expected
to arrive at the ground state of HP with high probability. SA utilizes the thermal fluctuation
 3
SA

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
Energy
h

QA

w
Spin configuration
Figure 1. Schematic picture of the thermal fluctuation and quantum tunneling in a system with local energy minima
separated by an energy barrier with the height h and the width w.

for optimization, which induces the thermal (Arhenius) jump from a local energy minimum to
another separated by an energy barrier. The escape rate from a local minimum over the energy
barrier with height h is given by e−h/kB T , where kB and T denote the Boltzmann constant and
the temperature. Assuming that h is proportional to the system size N , this suggests that an
exponentially long time in N is necessary to reach the global energy minimum by SA.
In contrast to SA, quantum tunneling induces an escape from a local minimum through an
energy√barrier as shown schematically in Fig. 1. The tunneling probability is approximately given
by e− hw/g [14,15], where g denotes the strength of quantum fluctuation which corresponds to
the transverse field Γ in transverse Ising models. Therefore, assuming the height h ∼ O(N ) and
the width w < O(N 1/2 ), the time necessary to escape from a local minimum due to quantum
tunneling is subexponential in N . For such a system, quantum tunneling helps the system to
equilibrate even though the system is glassy, i.e., non-ergodic in the absence of the quantum
fluctuation, leading to a potential advantage of QA over SA in glassy systems. This role of
quantum tunneling was first discussed by Ray et al., in 1989 [16] (see discussions in [17]) in this
regard) in the context of the restoration of the replica symmetry or ergodicity due to quantum
fluctuation in the quantum version of the Sherrington-Kirkpatrick model [18], which is detailed
in the next section. Although the existence of an energy landscape with thin and high barriers
in specific models is still an issue of debate, it must be a foundation for the speedup of QA over
SA [17,19]. In addition, several numerical and experimental studies have provided evidences for
such an advantage of QA over SA in some specific models as shown in Secs. 1(b)i and 2. We show
a brief time-line for the development of QA in Fig. 2.
The review is organized in the following fashion: Having discussed the basic idea behind the
QA scheme and the results for various models in the context of annealing and defect generation
especially for annealing across a quantum critical point in Sec. 1, we move to discuss various
implementations of annealing protocols in Sec. 2. In Sec. 3. , we probe how does coupling to an
external environment influence the QA process. In Sec. 4, we again refer to the close systems and
discuss possible ways to speed QA processes especially in the context of avoiding discontinuous
phase transitions. Some recent applications are discussed in Sec. 5.

(a) Quantum phase transition and quantum annealing

The minimum gap min[∆(t)] appearing in Eq. (1.3) often decreases with increasing the number
of spins. In general the energy gap vanishes at a quantum phase transition (QPT) because a QPT
separates disordered and ordered phases and the ground state is degenerate at a transition. Let
us consider the transverse Ising Hamiltonian introduced in Eq. (1.2). The initial state, for Γ  J,
 Chakrabarti [47] Wu et al. [84] Farhi et al. [3] proposed the Santoro et al.
introduced the reported supporting methodology of quantum [79] reported on
transverse field Ising evidences for rapid adiabatic computation in the the study of
spin glass model and decrease in context of NP-hard problems. quantum
reported the first characteristic annealing in 4
study of the quantum relaxation times to Ising spin
phase transition reach the ground Brooke et al. [85] glasses and showed
behavior of the state in LiHoYF reported their first that the residual energy

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
Mattis (unfrustrated) dipolar (Ising) glass experimental decreases faster with a
and the Edwards- due to the effective demonstration of larger power (than in
Anderson (frustrated) quantum tunneling in quantum annealing the classical case) of
Ising Spin glass the system. and its advantages the inverse annealing
models in transverse in LiHoYF Ising time.
field. glass magnets.

1981 1989 1991 1994 1998 1999 2000 2002 2011

Ray et al. [16] proposed Finnila et al. [1] Kadowaki & Johnson et al. [88]
that quantum fluctuations reported a successful Nishimori [2] first reported on the
could help explore the computational search formulated and remarkable
rugged (free) energy of the minima (ground numerically development and
landscapes of the state) of a demonstrated clearly functioning of a
Sherrington-Kirkpatrick multidimensional the computational Josephson-Junction-
spin glasses and search energy landscape (in advantages in Ising Coupled circuit
the ground state(s) by the context of large glass-like mixed quantum Ising spin
escaping from local molecules) using magnetic models. glass annealing
minima (having tall but quantum annealing. machine, built and
thin barriers) using later marketed by D-
tunneling. Wave Systems.

Figure 2. A brief Time-line for the development of Quantum Annealing.

is given by the ground state of the transverse field, in which all spins are aligned along the x
axis of spin. This is a disordered state where the ground-state averaged magnetization in the z
direction of spin is zero, i.e., hσiz i = 0. The targeted ground state of the Ising Hamiltonian for
Γ = 0, however, is an ordered state in the sense that it has a fixed magnetization +1 or −1 for
each σjz . This implies that the system encounters a QPT during QA. Indeed, the model in Eq. (1.2)
has QPTs at Γ/J = ±1. The finite size scaling of the energy gap at QPT depends on the character
of the associated QPT, and the latter is determined by the property of the Ising Hamiltonian.
The character of a conventional continuous QPT is specified by critical exponents [20–23]. The
size scaling of the energy gap at a quantum critical point is given as ∆c ∼ L−z where L denotes the
linear size of the system and the exponent z, known as the dynamical exponent, characterizes the
associated quantum critical point (QCP). Therefore the time for QA to work scales polynomially
with the system size. However, apart from this simple situation, the polynomial scaling of the
energy gap at QPT is not always true. In fact, a discontinuous QPT usually gives rise to an
exponential scaling with the system size. This can be understood phenomenologically as follows.
Consider a quantum many-body system and focus on the two lowest energy levels. We assume
that higher energy levels are highly separated from them. The effective Hamiltonian is then
written as " #
εA ∆
H= ,
∆ εB

where εA and εB corresponds to the energies of the two local minima, and ∆ denotes the
tunneling energy between these two states. Figure 3 shows the energy levels of this Hamiltonian
schematically. The discontinuous QPT corresponds to the change of the lowest energy level
between A and B. The transition takes place where the bare energies εA and εB of two levels
are degenerate. The energy gap at the transition is given by the twice of the tunneling energy ∆,
and ∆ is given by an exponential of the Hamming distance between the states A and B. Note
 (1) (3) 5

(2)

Energy

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
(1) (2) (3)

Figure 3. Schematic picture for an interchange of two energy levels. εA and εB corresponds to the energies of the two
local minima. (1), (2), and (3) shows the situations with εA  εB , εA = εB , and εA  εB , respectively. In the case (2)
with εA = εB , the energy gap is given by the twice of the tunneling energy ∆ between the states A and B .

that the Hamming distance is the number of sites at which the spin orientation along z axis is
different. Usually this distance increases linearly with the system size. Therefore the energy gap
at a transition decays exponentially with the system size. Since the discontinuous QPT hinders
QA, several ways to avoid the discontinuous QPT have been proposed. We will mention some of
them in Sec. 4.
QA across a QPT is closely related to the Kibble-Zurek mechanism of defect generation
following an annealing across a QCP. [24–30]. The system starting from the initial disordered
ground state evolves adiabatically as far as the characteristic time of the instantaneous ground
state (i.e., the inverse of the gap) is shorter than the annealing speed. However, on approaching a
QPT, the characteristic time grows and hence the dynamics becomes non-adiabatic in the vicinity
of the QCP. The state of the system after the passage through the QCP is no longer the ground
state, rather a state with topological defects. The residual energy density, εres , i.e., the excess
energy over the expected final ground state at the end of the QA is a monotonically decreasing
function of the annealing duration τ . In case of a linear annealing through a conventional
continuous QPT in a d-dimensional many-body system with the critical exponent ν for the
correlation length and the dynamical exponent z, Kibble-Zurek scaling of the residual energy is
given by εres ∼ τ −dν/(zν+1) as far as the system after annealing is in a gapped phase. The scaling
of the residual energy density is modified from the Kibble-Zurek scaling for other unconventional
continuous QPTs or discontinuous QPTs and when the annealing protocol involves a non-linear
variation of the tuning parameter [30–32]. The scaling of the residual energy together with the
scaling of the energy gap at QPT is an important measure that characterizes the property of
QA [33,34].

(b) Transverse Ising models

Assuming the transverse field Hamiltonian as HD , the total Hamiltonian of QA forms the
transverse Ising models (TIMs). We briefly review properties of some representative TIMs in this
subsection [24–27].
(i) Pure and disordered transverse Ising chain 6
As the simplest case we first consider the pure ferromagnetic one-dimensional TIM (1dTIM) given

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
by Eq. (1.2). There are two phases of the ground state in this model separated by a quantum phase
transition at Γ/J = 1. The ground state is disordered for Γ/J > 1, while it is ferromagnetically
ordered for Γ/J < 1. At the critical point Γ/J = 1, the energy gap above the ground level vanishes
as the system size L → ∞. The scaling of the energy gap with the linear size L at the critical point
is given by ∆ ∼ L−1 . Critical exponents of the correlation length and the dynamical exponent are
ν = 1 and z = 1, respectively.
The scaling of the defect density following a QA of the pure 1dTIM was solved exactly by
Dziarmaga [35]. Let us assume Γ = −tJ/τ with the parameter τ which denotes the inverse of the
annealing speed. Using the periodic boundary condition in Eq. (1.2), and applying the Jordan-
Wigner transformation followed by the Fourier transformation, the quantum time evolution
of the spin state is reduced to decoupled Landau-Zener models of two-level systems for each
momentum mode [22,23]. When the time t is varied from −∞ to 0, with the initial state chosen to
be the ground state of the initial Hamiltonian, the residual energy per spin at the final time t = 0
is found to to be of the form
1 1
εres = p (1.4)
π 2Jτ /~

in the thermodynamic limit L → ∞. We recall that the residual energy is defined as the excess
energy that is the difference of the energy expectation value of H(t = 0) with respect to the
evolved state at t = 0 from the ground energy of H(t = 0). According to the Kibble-Zurek scaling
with ν = z = 1, one has εres ∼ τ −1/2 consistent with Eq. (1.4).
The disordered version of 1dTIM is given by

Jj σjz σj+1
z
hj σjx ,
X X
H =− −Γ (1.5)
j j

where Jj is a random ferromagnetic coupling and hi is a random transverse field obeying

distributions πJ (J) and πh (h), respectively. The phase transition of this model happens
when [log J]av = log Γ + [log h]av , where [· · · ]av denotes the random average. The ground
state is ferromagnetically ordered for Γ < exp([log J]av − [log h]av ), and disordered for Γ >
exp([log J]av − [log h]av ) [36–39]. The phase transition is characterized by the infinite randomness
fixed point, where the dimensionless parameter appearing in √ the distribution of the energy gap
√ −C L
∆ is (− log ∆)/ L, implying that the energy gap scales as e , where C is a positive constant.
Therefore, even though the ground state of this model with Γ = 0 is trivial, the dynamics of QA to
this state across the quantum phase transition is highly nontrivial. The size scaling of the typical
gap suggests that the time to arrive at the target state by QA scales as subexponential with L. In
connection to the Kibble-Zurek scaling, Dziarmaga and Caneva et al., reported that the density of
kinks produced after QA scales as

[ρkink ]av ∼ 1/ log2 ατ (1.6)

with an O(1) constant α [40,41]. Thus, the defect density decays in a logarithmically slow
fashion with the annealing time τ which we reiterate makes QA difficult. However, it has been
reported that SA for the one-dimensional disordered Ising model (i.e., Eq. (1.5) with Γ = 0) yields
[ρkink ]av ∼ 1/ log α0 τ , where α0 is a constant, which decays slower than Eq. (1.6) [42]. Therefore,
this model reveals an evident advantage of QA over SA.

(ii) Pure transverse Ising model in higher dimensions

The two-dimensional TIM (2dTIM) may be the simplest model next to 1dTIM, though unlike
the 1d case, the 2d model is not integrable. The equilibrium properties of the 2dTIM has been
studied numerically and some of thermodynamic properties including the character of quantum
and thermal phase transitions are available. Recently Schmitt et al., carried out numerical study of
QA in 2dTIM in the context of the Kibble-Zurek scaling using state-of-the-art numerical methods
[43]. Their results are consistent with the Kibble-Zurek prediction. Study of out-of-equilibrium 7
dynamics of a two-dimensional quantum system will be a direction of study in the near future.
The situation becomes simpler in the infinite dimension. The pure TIM in the infinite

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
dimension is written as
J X z z X x
H =− σj σk − Γ σj , (1.7)
N
j<k j

where N denotes the number of spins. Note that each spin interacts with all the other spins with
~ = (1/2) PN ~σj , this Hamiltonian can be
an equal strength. Defining the total spin operator as S j=1
arranged into
J 2 J
H = −2 Sz − 2Γ Sx + . (1.8)
N 2

In the thermodynamic limit, this model undergoes a continuous QPT at Γ = J. The energy gap
above the ground state behaves as ∆ [(Γ − J)Γ ]1/2 for Γ ≥ J and the size scaling of the energy
gap at Γ = 1 is ∆c ∼ N −1/3 [44]. Introducing the effective dimension deff so that the system size
N is tied to the linear size L by Ldeff = N , one has relations between critical exponents as zν = 1/2
and z/deff = 1/3. Then, assuming z = 1 as in pure TIMs with finite dimension, one obtains ν = 1/2
and deff = 3. Caneva et al., studied QA of the present model and obtained εres ∼ τ −1/3 . This
scaling is inconsistent with the Kibble-Zurek scaling, since the latter predicts τ −1 . Acevedo et al.,
revealed that there is an anomaly in the transition amplitude between the ground and excited
states in the present model [45]. Therefore the naive phenomenological argument to derive the
Kibble-Zurek scaling does not apply to the system in infinite dimension. We shall also discuss an
extension of the Hamiltonian (1.7) to the p-body interacting model in Sec. 4 and argue that the QA
does not work in this model with odd p.

(iii) Transverse Ising spin glass

The Hamiltonian of the Edwards-Anderson [46] version of the transverse Ising spin glass,
introduced by Chakrabarti in 1981 [47], is written as

Jjk σjz σkz − Γ

X X x
H =− σj , (1.9)
hjki j

where hjki stands for nearest neighbour pairs and Jjk are independent random variables. The
order parameter of a spin glass is defined in terms of the spin overlap between different replicas.
Supposing that σjα,a denotes the spin operator for a replicated system labeled by a, the overlap
z,1 z,2
operator between replicas a = 1 and a = 2 is defined by R1,2 = (1/N ) N
P
i=1 σj σj . The order
parameter is then given by q = [hR1,2 i]av . The spin glass order is characterized by q > 0 with
zero magnetization m = 0, where the magnetization is defined by m = [h(1/N ) N z
P
i=1 σj ]av . This
means that the spin configuration is spatially random but frozen. Rieger et al., and Guo et al.,
investigated the character of QPTs of this model with the Gaussian distribution of Jjk with zero
mean and unit variance in square and cubic lattices, respectively, by means of the quantum Monte
Carlo simulation [48,49]. Singh and Young studied ±J model where Jkj takes +1 or −1 with
equal probability for dimensions up to d = 8 using the linked cluster expansion to determine the
location of the QCP [50]. Subsequently, QPTs of these models were reconsidered by Miyazaki
and Nishimori [51] and by Matoz-Fernandez and Romá [52] using the real-space renormalization
group and the quantum Monte-Carlo with parallel-tempering, respectively. They concluded that
the QPTs in transverse Ising spin glasses in two and three dimensions were compatible with
the infinite randomness fixed point with the critical exponents ν and ψ, where ψ specifies the
activation type of size scaling of the energy gap as [log ∆]av ∼ N ψ/d [53–57]. The estimated
exponents for the Gaussian model were ν ≈ 1.2 and ψ ≈ 0.44 in two dimension [51,52] and
ν ≈ 0.94 in three dimension [51].
 The Hamiltonian of the transverse Ising spin glass in infinite dimension, i.e., the quantum 8
Sherrington-Kirkpatrick (SK) model, is written as [16]

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
N
1
Jjk σjz σkz − Γ σjx ,
X X
H = −√ (1.10)
N 1≤j<k≤N j=1

The classical SK model in the absence of the transverse field unveiled the existence of so-called
replica symmetry breaking (RSB) in the spin glass phase [58,59], where the overlap R1,2 has a
dispersed continuous distribution in the thermodynamic limit. Ray et al., conjectured on the basis
of the quantum Monte-Carlo simulation the collapse of a continuous distribution for the classical
SK model into a delta function in the presence of any amount of the transverse field [16], which
paved the way for using quantum tunneling in finding the global minimum or ground state of
SK spin glass model. In the classical model, due to random interactions between spins at different
lattice sites, such systems have many local minima in free energy which are separated by large
energy barriers of order O(N ), where N is the system size [59]. This induces non-ergodicity in
the system and eventually breaks the replica symmetry of the system. As a result, finding the
ground state or global minimum of such systems is a very hard problem; for SK spin glass model,
it turns out to be NP (non-deterministic polynomial-time) hard. The system indeed gets trapped
into one of the local minima inside the spin glass phase, due to the highly rugged nature of free-
energy landscape. This leads to a broad order parameter distribution in the spin glass phase [58].
In addition to a peak value of the order parameter distribution, it is extended up to the zero value
of the order parameter even in the thermodynamic limit.
It seems that the scenario may change drastically, when a transverse field is applied on the
SK spin glass [16]. The presence of quantum fluctuations induces ergodicity in the system, since
quantum tunneling becomes possible between the local minima separated by tall and narrow
free-energy barriers. This indicates the restoration of replica symmetry breaking for quantum SK
spin glass model. As a result, the order parameter distribution should be sharply peaked at a
point for quantum SK model in the thermodynamic limit. This ergodic behavior of quantum SK
model is responsible for advantage in quantum annealing in comparison to simulated annealing.
This conjecture was criticized by Young [60] by solving numerically the effective one-
dimensional model to which the quantum SK model can be mapped in the N → ∞ limit; this
work predicted that the replica symmetric solution is unstable down to zero temperature. On
the contrary, Mukherjee et al., [61] explored the behavior of the order parameter distribution of
the quantum SK model in the spin glass phase using Monte Carlo technique for the effective
Suzuki-Trotter Hamiltonian at finite temperatures (see Eq. (2.1) discussed later) and the exact
diagonalization method at zero temperature. It has been found that there exists a low temperature
regime in the spin glass phase, where the order parameter distribution becomes peaked around
its most probable value in thermodynamic limit, thus suggesting the ergodic behavior. On the
other hand, the order parameter distribution remains Parisi type in high temperature regime,
which indicates the non-ergodic behavior of the system in this part of the spin glass phase.
These two regions of the spin glass phase are separated by a boundary, connecting the zero
temperature-zero transverse field point and the quantum-classical crossover point on the phase
boundary [61,62]. In addition, quantum annealing has also been investigated for quantum SK
model using Suzuki-Trotter Hamiltonian dynamics in both the ergodic and non-ergodic regimes.
The average annealing time was estimated, when both the temperature and the transverse-field
were annealed down to some fixed low values, starting from the paramagnetic phase. It was
found that the average annealing time is independent of the system size, when the annealing
is performed through the ergodic (quantum fluctuation dominated) region, whereas it grows
strongly with the system size, when the annealing is carried out through the non-ergodic (classical
fluctuation dominated) region. This suggests that the quantum annealing has potential to detect
whether a phase is ergodic or non-ergodic. Also, the average annealing time to approach a
same ground state is small for annealing through ergodic regime compared to that through the
non-ergodic regime. The QA for SK spin glass is also studied by tuning both transverse and
longitudinal fields, and it has been shown that this protocol exhibits some effectiveness compared 9
to the QA by varying the transverse field only [63].
Recently, Leschke et al., proved rigorously nonzero variance of the overlap in the

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
thermodynamic limit of the quantum SK model at sufficiently low temperature with small but
finite transverse field [64]. Their study reveals a dispersed distribution of the overlap, therefore
the existence of RSB. However the controversy about the continuous distribution of the overlap
in the quantum SK model is still an open problem [65].

(c) Transverse Ising model for satisfiability

The satisfiability problem is known as one of the basic combinatorial optimization problem
in computer science. Given the number of bits and constraints among bits, the problem is to
determine whether the bit configuration satisfying all the constraints exists or not. In the case of
3-satisfiability problems, or 3-SATs in short, each constraint involves three bits drawn randomly.
Using the spin language, a constraint for three spins Si , Sj and Sk taking the values ±1 can be
represented as (Si + Sj + Sk − 1)2 , for instance. This vanishes and is called satisfied when two
of the three spins are +1 and the other spin is −1, otherwise it gives a nonzero and positive value.
The Exact Cover, a variant of the 3-SAT, consists of M such constraints for N spins, and thus the
corresponding quantum model that needs to be annealed is given by
M N
(σizm + σjzm + σkzm − 1)2 − Γ σjx .
X X
H= (1.11)
m=1 j=1

If one arrives at the exact ground state for Γ = 0, by annealing the field Γ , then one can solve the
Exact Cover. However, the Exact Cover is an NP-complete problem which no known algorithm
can solve in a time polynomial in N . Young et al., studied Eq. (1.11) by means of the quantum
Monte-Carlo method [66]. They found that some instances of the model show a discontinuous
first-order QPT with an exponentially small energy gap and the fraction of such instances grows
toward unity with increasing N . Jörg et at. also reported occurrence of a first-order QPT in the
random 3-XORSAT problem, which is another variant of the 3-SAT [67].

2. Implementation of quantum annealing

Implementing QA is a challenging task, since one needs to evolve a many-spin state under a
quantum many-body Hamiltonian. In this section, we review results from numerical simulation
using real-time dynamics as well as Monte-Carlo dynamics, and from quantum simulations using
hardwares.

(a) Results from numerical simulation for real time dynamics

The real-time evolution of a quantum system governed by the Schrödinger equation can be
computed in general by solving a linear differential equation. However, since the number of
unknown functions of time increases as 2N with the number of spins, the system size acceptable
to a conventional computer is limited to N ∼ O(10). Kadowaki and Nishimori reported in the
seminal paper that QA yields better solutions than the classical simulated annealing on the basis
of their simulation for 8 spin systems of a frustrated model and the SK model [2]. Farhi et al.,
reported the numerical result for the Exact Cover problem with the number of spins up to 20 [3].
Exact Cover is one of the NP-complete problems, which no known classical algorithm can solve
in a time polynomial in the number of spins. The numerical result suggested a quadratic scaling
of the runtime in QA with respect to the number of spins. However, this scaling should turn into
an exponential one for larger size as shown by the quantum Monte-Carlo study for the energy
gap [66].
Restricted to systems in one dimension, there are efficient methods of numerical simulation
for real-time evolution. In case of 1dTIMs without the longitudinal field, irrespective of disorder,
 S1,M
10

Temporal direction

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
S1,3
S1,2
Spatial direction

S1,1 S2,1 S3,1 SN,1

1dTIM
Spatial direction
2d Ising model

Figure 4. Schematic picture of the Suzuki-Trotter mapping. A 1dTIM is mapped to a two-dimensional classical Ising model
on the square lattice. The additional dimension corresponds to the time. Sj,m denotes an Ising-spin variable at spatial
site j and temporal site m.

the Schrödinger equation of the spin state reduces to the Bogoliubov-de Gennes equation of
2N unknown functions of time through the Jordan-Wigner’s fermionization and the Bogoliubov
transformation [68]. Generic 1dTIMs with longitudinal fields cannot be mapped to free fermion
models. However, time evolution of generic 1dTIMs can be simulated using the time-dependent
density matrix renormalization group (tDMRG) proposed by White and feiguin [69] or the time
evolving block decimation (TEBD) by Vidal [70]. In addition, the infinite system of 1dTIMs with
no disorder can be simulated using an infinite method of TEBD (iTEBD) [71]. These methods serve
the study of QA in 1dTIM with a uniform or disordered longitudinal field [72,73].

(b) Results from numerical simulation for Monte-Carlo dynamics

A d-dimensional TIM in finite temperature (denoted by β −1 ) with spin-spin coupling Jjk and
transverse field Γ can be mapped to a (d + 1)-dimensional classical Ising model by the Suzuki-
Trotter mapping [74,75]. The resulting Hamiltonian Heff is given by
M X N M
X Jjk M βΓ X X
Heff = − Sj,m Sk,m − log coth Sj,m Sj,m+1 , (2.1)
M 2β M
m=1 j<k j=1 m=1

where β is the inverse temperature, M is the Trotter number, Sj,m denotes the spin variable
with the spatial site j and temporal site m taking values ±1, and we defined the sign of Jjk
according to Eq. (1.9). In Figure 4, we schematically illustrate the mapping of 1dTIM into a
(1+1)-dimensional classical Ising model. One can simulate in principle any TIM in and out of
equilibrium using this effective Hamiltonian and the Monte-Carlo method. This method is called
the quantum Monte-Carlo method (QMC). Although the number β/M controls the accuracy of
QMC, the cluster-update method invented by Swendsen and Wang along the temporal direction
enables to have β/M → 0 [76,77]. QMC is known to give rise to the sign problem and fail when the
model involves the frustration. However, QMC for TIM is free from the sign problem. Therefore
QMC is a powerful method of classical computation in simulating TIM.
QA can be implemented in QMC by regarding the Monte-Carlo step as time. The dynamics
realized by QMC is not the quantum dynamics governed by the Schrödinger equation but the
stochastic one. However, QA with QMC serves the purpose of solving an optimization problem
using a classical computer [78]. Several works have shown so far that QA with QMC works in
variety of optimization problems, such as two-dimensional Ising spin glass [79,80], travelling
salesman problem [81], and 3-SATs [82]. Figure 5 shows comparison between QA by QMC and
 0.1 11

Residual energy per spin

[Link] Phil. Trans. R. Soc. A 0000000

SA

..................................................................
(ln τ)-2

0.01
QA

(ln τ)-3

0.001
100 1000 10000 100000
Monte-Carlo step

Figure 5. Comparison of the residual energy between QA by QMC and SA for the two-dimensional spin glass model
with 99×99 spins with random coupling Jjk from the uniform distribution between -2 and 2. The cluster-flip algorithm in
the imaginary-time direction was used in QMC. SA was started from the initial temperature T = 5.0, while QA was from
Γ = 5.0 with T = 0.01. The average was taken over 100 runs for a single instance in SA, and for 16 instances in QA.
The decay of the residual energy in SA is well fitted by (log τ )−2 . As for QA, it is approximated by (log τ )−3 , except for
long annealing time where the decay rate is smaller. (Taken from ref. [83].)

SA in the two-dimensional spin glass model [83]. This result implies outperformance of QA over
SA. Although an opposite result has been reported for harder 3-SAT problems [82], numerical
studies using QMC suggest that there are problems for which QA is potentially advantageous
over SA due to the restoration of ergodicity by quantum fluctuation [16].

(c) Results from quantum simulation using hardware

The most efficient way to perform QA is use a quantum magnetic material which realizes a
TIM with a temporally controllable transverse field. LiHox Y1−x F4 is a material that models a
disordered Ising model and it also realizes a disordered TIM by application of the magnetic field
perpendicular to the easy axis of magnetization [84]. Brooke et al., investigated two protocols, QA
and thermal annealing (TA), using this material with x = 0.44. In QA, the transverse field Γ was
strengthened at high temperature, the system was cooled, and then Γ is weakened to Γf at low
temperature. On the other hand, in TA, the temperature was lowered with keeping Γ = 0, and
then Γ was raised to Γf . In both protocols, the initial and final sets of Γ and temperature were the
same and the durations were the same as well. The ground state at Γf as the target state is a glassy
ferromagnetic state. Brooke et al., reported that the state after QA much is closer to an equilibrium
state than the state after TA [85]. This result implies that QA brings the state to the target faster
than TA. Quite recently Säubert et al., studied QA and TA of the same material and detailed the
dynamical behavior of the energy landscape during QA. They showed that the transverse field
applied in QA induced random longitudinal fields, implying that the energy landscape of the
problem Hamiltonian HP evolved as QA proceeded [86]. This evolving landscape may be an
issue of future work related to QA.
Progress in Rydberg atom experiments enables to use Rydberg atoms as a quantum simulator.
Keesling et al., performed a quantum simulation of a sort of QA using an array of 51
Rydberg aroms. In this simulation, the system is described by the many-body Hamiltonian H =
P P P
(Ω/2) i (|gi ihri | + h.c.) − ∆ i ni + j<k Vjk nj nk , where |gj i and |rj i denote the ground
and the excited Rydberg states, respectively, of atoms, nj = |rj ihrj |, and Vjk is the van der Waars
 12
85th

QMC and SA single-core annealing time ( )

75th

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
50th

D-Wave annealing time ( )

QMC

SA

D-Wave

Problem size (bits)

Figure 6. Comparison of the time to reach the ground state with 99% success probability as a function of the problem
size in D-Wave 2X, simulated annealing (SA), and quantum Monte-Calro (QMC). The runtime in SA and QMC is defined
by nsweeps N Tupdate , where nsweeps is the number of sweeps (one sweep attempts to update all spins). Tupdate
is the single-spin update time for SA and the update time of a spin-cluster along the temporal direction. It is set as
Tupdate = 1/5 ns for SA and 10 × 870 ns for QMC. Data for 50th, 75th, and 85th percentile taken from a set of 100
instances are shown. The error bars represent 95% confidence interval from bootstrapping. Taken from ref. [89].

interaction with the strength which decays as 1/|j − k|6 . The model exhibits a QPT, achieved by
tuning the parameter ∆, belonging to the same universality class as that of 1dTIM . Keesling et
al., observed that the Kibble-Zurek scaling for the kink density arising due to the sweeping of the
parameter ∆, turns out to be the very same as that in 1dTIM [87].
In order to apply QA as a computation to an optimization problem in practice, spin-
spin interactions and longitudinal fields in addition to the transverse field need to be locally
controllable. A Canadian venture company, D-Wave Systems, has developed a quantum
annealing machine named as a quantum annealer, which consists of programmable coupled
superconducting flux qubits and performs QA to various Ising models [88]. The number of qubits
in the latest machine is beyond five thousands. This is 100 times larger than the number of qubits
in the current gate-based quantum computer. Denchev et al., benchmarked D-Wave 2X using
100 instances of the weak-strong cluster model with up to 945 spins [89]. Qubits in D-Wave 2X
form the so-called chimera graph with unit cells consisting of 8 qubits. In the weak-strong cluster
model, there are all-all ferromagnetic couplings inside the cell, and half of the spins in a cell
ferromagnatically couple with those in neighboring cells. In addition, weak longitudinal fields
are applied to spins in randomly chosen cell, while strong fields anti-parallel to the weak ones
are applied to spins in the other cells. Figure 6 shows the time to reach the ground state with
99% success probability. For D-Wave, this time is defined by 20 µs [log(1 − 0.99)/ log(1 − p)] for
an instance, where the annealing time is fixed at 20 µs and p denotes the success probability to
obtain the ground state estimated from many runs. As for SA and QMC, it is the runtime on a
single processor. Regarding the median from 100 instances, D-Wave 2X is 108 and 107 times faster
than SA and QMC, respectively.
Boixo et al., tested D-Wave’s quantum annealer to a spin glass model HP = − hjki Jjk σjz σkz ,
P

where Jjk is chosen randomly from J = ±1, with N = 108 spins and reported that the results of
quantum annealer correlated well with those obtained by QA with QMC [90]. Figure 7 shows
the comparison of the histogram of the success probability between D-Wave’s quantum annealer
(DW) and QA with QMC (named as Simulated QA). The bimodal distribution which is common
in D-Wave and QMC could be an evidence that the system embedded in D-Wave’s quantum
 500 500 13
Number of instances (a) DW (b) SQA

Number of instances
400 400

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
300 300

200 200

100 100

0 0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Success probability Success probability

Figure 7. Comparison of the histogram of the probability of finding the ground state between D-Wave One and QA with
QMC (simulated QA) over 1000 instances of the spin glass model with N = 108 spins. Taken from ref. [90].

annealer was a quantum system. However, Shin et al., reported that the classical spin vector
model along with the Monte-Carlo dynamics, named as the spin vector Monte-Carlo (SVMC)
model, provided as strong correlation with D-Wave’s data as QMC [91]. The classical spin vector
model is represented by the Hamiltonian
N
X X
H(t) = −A(t) cos θj − B(t) Jjk sin θj sin θk , (2.2)
j=1 hjki

where θj denotes the angle of the unit vector at site j in the xz plane. These works have raised
a problem to identify the model of D-Wave’s quantum annealer [92]. Bando et al., studied the
Kibble-Zurek scaling in 1dTIM using D-Wave 2000Q and found that the kink density defined
by n = (1/2N ) N z z −α
P
j=1 (1 − hσj σj+1 i) = εres /2 scaled with the annealing time ta as n ∼ ta with
α ≈ 0.20 by the device at NASA and α ≈ 0.34 by the one at D-Wave Systems [93]. As mentioned in
−1/2
Sec. i, the scaling of the kink density is predicted as n ∼ ta for an isolated system belonging to
one-dimensional Ising universality class. The authors in ref. [93] compared numerical simulations
for 1dTIM with coupling to an environment and for SVMC with the experiment, and concluded
that the quantum model agreed better with the experiment. Recently, King et al., studied QA of
1dTIM using D-Wave 2000Q [73], focusing on shorter annealing times than those in the previous
works. For short annealing times, the system in the device is less affected by environment as
we shall discuss in the next section. Comparing analytic and numerical computation for the
Schrödinger dynamics of the isolated 1dTIM, the QMC simulation, and the SVMC simulation
with the experiment by D-Wave 2000Q, King et al., reported that only the Schrödinger dynamics
of the isolated 1dTIM with a small amount of disorder can explain all the experimental results
[73]. Also, in ref. [94], fully connected Sherrington-Kirkpatrick model with random couplings
was programmed using D-Wave TwoTM annealer, where optimal parameter setting allowed
better performance of the quantum annealer when compared to those obtained using optimized
simulated annealing algorithms.

3. Effects of environment on QA
Although QA is ideally performed in an isolated system, any real system is always coupled to
an environment and hence susceptible to decoherence. In fact, the system in D-Wave’s device is
believed to be affected considerably by an environment when the annealing time is longer than a
few µs. Therefore it is very important to study an effect of environment in QA.
There is a variety of models representing an environment. Caldeira and Leggett, in their
seminal work analyzed the dynamics of flux state in a SQUID and constructed a simple model of
a two-level system coupled to a boson bath, where bosons are attributed to the electro-magnetic 14
field coming from the fluctuating current [95]. Leggett et al., elaborated on the single-spin model
coupled a boson bath [96]. Thus, considering QA performed with superconducting flux qubits, the

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
model which includes the effect of an environment should be an extension of the Caldeira-Leggett
model to many spins. Hamiltonian is written as H(t) = HS (t) + HB + Hint , where the system
Hamiltonian HS (t) is given by Eq. (1.1). The bath is represented by the collection of harmonic
oscillators. Hence the bath Hamiltonian HB is given, using the boson operators bj,a and b†j,a for
site j and mode a, by

~ωa b†j,a bj,a ,

X
HB = (3.1)
j,a

where ωa is the frequency of the harmonic oscillator of mode a. Hint represents the interaction
between the system and the bath, and written as
N
σjx Qx z z
X

Hint = j + σj Qj (3.2)
j=1

where Qγj (γ = x, z) is the bath operators given by Qγj = a λγa (b†j,a + bj,a ). The spectral density
P

of the boson bath is assumed as Jγ (ω) = a (λγa )2 δ(ω − ωa ) = ηγ ω s e−ω/ωc , where ηγ denotes the
P

coupling strength of the system-bath interaction and ωc is the energy cutoff of the bath spectrum.
The Ohmic bath refers to s = 1, while the super-Ohmic and sub-Ohmic baths refer to s > 1 and
s < 1, respectively.
Let us now move to study the time evolution assuming that the state of the composite system
is described by the density operator ρ(t) at the instant t. The initial state is assumed to be a direct
product state of the form ρ(0) = |ψ(0)ihψ(0)| ⊗ e−HB /T /ZB , where |ψ(0)i denotes a state vector
of the system, T is the temperature, and ZB is the partition function of the bath. Since we are
interested in the behavior of the system, we consider the reduced density operator describing the
system, ρS (t) = TrB ρ(t), where TrB stands for the trace with respect to the bosonic degrees of
freedom.
Amin explored the success probability of QA for a range of annealing time ta obtained by
solving numerically the quantum Redfield master equation for an instance of 16 spins of random
Ising model in a random longitudinal field with nonzero ηz and ηx = 0 [97]. The obtained success
probability is a nonmonotonic function of ta . For short ta , the spin system is not influenced by
the bath and hence the success probability increases with increasing ta . In a middle range of
ta , the thermal environment disturbs the system’s adiabatic evolution more for longer ta , hence
leading to decreasing success probability. For very long ta , finally, the system evolves keeping
the thermal equilibrium with the bath until it is frozen near the end of QA. The freezing happens
because HS (ta ) and Hint are commutable when ηx = 0 and hence the relaxation time diverges as
t → ta . Thus the success probability in this regime goes to the probability of the ground state at
the thermal distribution as ta → ∞.
QA of 1dTIM in the presence of an environment has been attracted a lot of attention in
the context of the Kibble-Zurek scaling. Assuming Qzj = 0, namely, the boson bath coupled
to σ x , 1dTIM with coupling to the boson bath is mapped to a noninteracting fermion model
with a fermion-boson coupling through the Jordan-Wigner transformation. Then the problem is
significantly tractable, compared to the situation with ηz 6= 0. Patané et al., studied the density
of excitation following QA using the Keldish technique. Based on the ansatz that the density of
excitation E is given by the sum of the coherent part Ecoh and the incoherent part Einc due to
the environment, Patané et al., obtained Einc ∼ ηx T 4 τ for the Ohmic bath with temperature T
when QA ends near the quantum critical point [98,99]. The incoherent part increases with τ in
contrast to the coherent part, hence its scaling is called the anti-Kibble-Zurek scalng. Nalbach et
al., studied the model with the spatially correlated bath where all spins are coupled to a single
bath, i.e., Hint = Qx j σjx and HB = a ~ωz b†a ba with Q = a λa (b†a + ba ). In this situation,
P P P

the correlation length is the largest of the Kibble-Zurek length scale ξKZ ∼ τ 1/2 and the thermal
length scale ξT ∼ T −1 . When 1  ξT < ξKZ , the thermal effect comes into play in the density 15
√
of excitation. Thus, it is suggested that E ∼ ηT (τ T 2 ) ∼ ηx T 3 τ for τ  T  1 in this model.
Nalbach et al., proposed this scaling relation and confirmed using the dissipative Landau-Zener

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
theory for the two-level system [100]. Dutta et al., studied 1dTIM in the presence of a spatially
homogeneous Gaussian white noise on the transverse field, instead of considering the coupling
to a boson bath. This stochastic perturbation yields an effective dynamics of the noise-averaged
density operator for an open quantum system. Then the noise-averaged density of excitation has
an incoherent part which scales as Einc ∼ η̃x τ , where η̃x is the strength of the noise [101]. Weinberg
et al., studied a similar model but with a spatially uncorrelated and temporally correlated noise.
The numerical result showed the same scaling as ref. [101]. Weinberg et al., also performed
quantum simulation for 2dTIM using D-Wave 2000Q and obtained scaling for the residual energy
εres ∼ aτ −α + bτ β with α ≈ 0.74 and β ≈ 0.456 [102]. The first term is consistent with the Kibble-
Zurek scaling for 2dTIM with the exponent dν/(zν + 1) ≈ 0.77, where d = 2, z = 1 and ν ≈ 0.63
[103]. Note that the scaling of the second term corresponding to the incoherent part is different
from that in 1dTIM, implying that Einc ∼ τ is specific to 1dTIM. The anti-Kibble Zurek mechanism
in the presence of non-thermal bath has been also discussed in the framework of the Lindblad
formalism in refs. [104–108]
The Kibble-Zurek and anti-Kibble-Zurek scalings implies the existence of a global or local
minimum of the density of excitation. Based on the numerical study using the Redfield master
equation in the momentum space for 1dTIM with the boson bath, Eqs. (3.1) and (3.2) with Qzj = 0,
Arceci et al., identified the region in the T − ηx -plane where there exists the local or global
minimum of the density of defects after QA [109]. Interestingly, a global minimum of the density
of defects appears in the digitized QA, in which the time-evolution operator inducing QA is split
into slices with a finite time width and further split into those involving only HD and those
involving only HP [110], suggesting that the decomposition of the time-evolution operator has
an influence as a decoherence on the dynamics of a closed system.
Many studies of 1dTIM with a boson bath focusing on the Kibble-Zurek physics have assumed
Qzj = 0, namely, coupling between the system operators σjx and bosons. However, in experimental
systems such as those made of coupled superconducting flux qubits, coupling between σjz and
QzJ is rather important [96]. Recently, Suzuki et al., developed new matrix-product-state-based
methods for 1dTIM with a boson bath with Hint = j Qzj σjz , which enables the simulation of
P

finite pure and disordered systems with O(102 ) spins [111,112] or an infinite translationally
invariant system [113]. Using the infinite-system method, Oshiyama found modified Kibble-
Zurek scaling in 1dTIM coupled to the bath at zero temperature [112]. Oshiyama et al., also
studied QA of 1dTIM with the bath at finite temperatures. In the thermal environment at finite
temperature T , the infinitely slow QA (ta → ∞) can be regarded as the quasistatic and isothermal
process, hence the final energy should be identical to the thermal average of B(ta )HP at T . When
ta is sufficiently long but finite, the energy of the system has an excess from the thermal average.
−1/3
Oshiyama et al., found and numerically confirmed that this excess energy scales with ta as ta ,
for the linear annealing protocol [113].

4. Avoiding first-order phase transitions: Closed systems

As discussed in Sec. 1(a), it has been generally observed (with a few exceptions) that the minimum
energy gap decreases exponentially with the system size for a first-order phase transition, whereas
it shows a polynomial decrement in system size for a second-order phase transition. Therefore,
the order of phase transitions is an important factor in determining the efficiency of the quantum
annealing algorithm.
(a) Quantum ferromagnetic model 16
We consider here a ferromagnetic p-spin model in the transverse field. The Hamiltonian for such

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
a system is given by
1 XN p N
σiz − Γ σix ,
X
H = −N (4.1)
N
i=1 i=1

where σiz and σix are usual Pauli matrices at the lattice site i, Γ is the magnetic field in transverse
direction and N is the number of spins and p is an integer. These type of models were initially
introduced in the context of spin glasses. The ground state of the classical model at zero
temperature with Γ = 0, corresponds to all spins aligned in the same direction. For even p, all
the spins in up or down states are valid ground states, whereas odd p has a unique ground state
when all the spins are in up states. Therefore, for simplicity, we will concentrate here on the odd
p cases. For p = 2, the Hamiltonian in Eq. (4.1) reduces to an infinite-range Ising model which can
be mapped to the usual mean field Curie-Weiss model exhibiting continuous phase transitions.
On the other hand, for p > 2, both classical and quantum phase transitions of the system are
discontinuous.
Using Suzuki-Trotter formalism and “static” approximation, the phase diagram of the p-spin
ferromagnetic model can be found in Γ − T plane for different values of p (see Fig. 8) [114]. In
the limit of p → ∞, using perturbation theory, the minimum energy gap of the system can be
calculated as ∆min = 2N 2−N/2 [114]. This indicates that the energy gap between the ground and
excited states closes exponentially fast with the system size at the transition point. For a general
p, an explicit form of the energy gap is not available so that one can comment about its scaling
with the system size, however, the same scaling can be inferred from numerical calculations.
The energy gap of the system can be calculated numerically using two complementary
methods as a function of the transverse-field Γ [114]. Using these numerical methods, we can
find the transition point Γc where the energy gap shows a minima that scales with the system
size. In the present case, the energy spectrum of the system has been studied for 3 ≤ p ≤ 31. The
Hamiltonian in Eq. (4.1) is represented by a sparse matrix of dimension 2N . For such systems,
Lanczos method provides nearly exact extreme eigenvalues of the Hamiltonian for the system
size N ≤ 21. From the results of the Lanczos method for N ≤ 21, it has been found that the
transition happens between two states with the maximum possible angular momentum l = N/2.
The efficiency of the numerical simulation can be improved by exploiting the fact that the
total angular momentum L2 commutes with the Hamiltonian H in Eq. (4.1), (where L is the
total angular momentum of N spins). Therefore the transition occurs mainly in the subspace of
dimension 2l + 1 = N + 1. In this subspace, the Hamiltonian assumes a tri-diagonal form and the
resulting tri-diagonal matrix can be diagonalized efficiently for a system with size N ∼ 100 in
just a few seconds. The energy gap has been shown in the left panel of Fig. 9 as a function of Γ
for p = 3 and different N . The gap becomes minimum at the critical value of Γ that agrees with
analytically predicted value. One can observe that the region where the gap closes gets narrower
as the value of N is increased. The minimum energy gap ∆min is further plotted as a function of
N for different values of p to find its dependence on N (see right panel of Fig. 9). It has been found
that the minimum energy gap decays exponentially as ∆min ∝ N 2−N α for p ≥ 3. The minimum
energy gap closes exponentially fast as expected for the first order phase transition. The value
of exponent α can be computed numerically from the right panel of Fig. 9. These exponents are
also calculated analytically using instantonic approach. A comparison of values of α for different
values of p are given in Table 1 of Ref. [ [114]].
Due to an exponential decay of energy gap with the system size, the running time increases
exponentially for the case of a first-order phase transition, and thus reducing the efficiency of
QA process. Therefore, it is an important issue to investigate whether one can avoid first-order
phase transitions in the annealing path to solve the optimization problem efficiently using QA
algorithm. Below we discuss various methods to speed up a quantum annealing process.
 1.8
p=Infinity
p=3
17
1.6 p=5
Paramagnet p=9
1.4 p=15

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
p=25
1.2

T
0.8

0.6
Ferromagnet
0.4

0.2

0
0.2 0.4 0.6 0.8 1 1.2 1.4
Γ

Figure 8. Phase diagram of the ferromagnetic p-spin model in the T − Γ plane for different values of p. The ferromagnetic
and quantum paramagnetic phases are separated by first-order phase transitions. Taken from [114].

6 1

5 p=3 0.01

0.0001
4
1e-06
∆min/N
∆min

3
1e-08
p=3
2 N=10 p=5
N=20 1e-10 p=7
N=30 p=11
1 N=60 p=15
N=90 1e-12 p=21
N=120 p=31
N=150 2-N/2
0 1e-14
0 0.5 1 1.5 2 0 20 40 60 80 100 120 140
Γ N

Figure 9. Left Panel: Variation of the energy gap as a function of Γ for p = 3 computed using an exact diagonalization
method as described in the text. The gap vanishes exponentially fast with the system size N near the critical point Γc
(the black vertical line). It also shrinks near the critical regime as N increases. Right Panel: Minimum energy gap as a
function of N for a few values of p on semi-log scale. It shows an exponential fall of the minimum gap with N for all values
of p according the scaling relation ∆min ∝ N 2−N α . Taken from [114].

(b) Application of antiferromagnetic fluctuations

In the context of speed-up of QA, Seki and Nishimori in 2012 proposed a method [115] to
overcome issues related to first order phase transitions, by studying quantum annealing in the
presence of antiferromagnetic fluctuations in addition to the transverse-field term. They applied
the method to the infinite-range ferromagnetic p-spin model (see Eq. (4.1)) and showed that there
exists a quantum path that avoids first-order transitions for some intermediate values of p. The
Hamiltonian for p-spin is given by

N
1 X p
H0 = −N σiz . (4.2)
N
i=1

This is indeed the classical counterpart of the Hamiltonian as in Eq. (4.1) with zero transverse field.
Here H0 is the target Hamiltonian HP , whose ground state is the optimal solution of the problem.
The QA for this model is studied before with the transverse-field as a driver Hamiltonian HD ,
which takes an exponentially long time to reach the ground state of the target Hamiltonian
due to the presence of a quantum first-order phase transition during the time evolution. The
ferromagnetic p-spin model reduces to the Grover problem when p → ∞ and there is no known
algorithm that can solve the problem efficiently in a polynomial time.
 0 1 18
p=3
p=5
p=11

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
p=21

mx
-0.2 0.5
f

p=3
p=5
p=11
p=21
-0.4 0
0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
s s

Figure 10. Left panel: Free energy of the system with the Hamiltonian in Eq. (4.3), as a function of s for some values
of p with λ = 0.3. The free energy of the QP phase, Eq. (A.9), is represented by the dash-dotted line in light green. The
thin solid line with blue color represents the free energy of the ferromagnetic phase (F), Eq. (A.10), and the thick solid line
with red color is for the QP2 phase, Eq. (A.11). The lower limit of the QP2 domain (s = 1/(3 − 2λ)) is denoted by the
vertical dashed line. Although it is hard to see in this present scale, all the data for finite p studied here, have lower values
than that of fQP2 in the QP2 regime. Right panel: Magnetization mx vs. s for λ = 0.3 and the same values of p as in the
case of free energies. The verical dashed line is the same as shown in the left panel. The solid line in red exhibits the
x component of magnetization of the QP2 phase. The magnetization decreases to zero as s is increased with a jump at
the boundary of the QP2 domain for p ≥ 5. Taken from [115].

We discuss here how the inclusion of an antiferromagnetic fluctuation term can improve the
performance of QA of the model when both the transverse-field term and the antiferromagnetic
term are tuned. The total Hamiltonian of the problem is then given by

H(s, λ) = s{λH0 + (1 − λ)V̂AF } + (1 − s)VTF , (4.3)

  2
1 PN x
where VAF = +N N i=1 σi is an antiferromagnetic interaction term, whereas VTF is the
conventional transverse-field term. The parameters s and λ are functions of time, assumed to
lie between 0 and 1, which are chosen appropriately for a QA process. The initial Hamiltonian is
defined by s = 0 and an arbitrary λ, and the final one is given by s = λ = 1.

(i) Numerical results

We now focus on analyzing the phase diagram of the model on the s − λ plane for finite values of
p, using the saddle point method and the static approximation, as elaborated in Appendix A. The
method we adopt to construct the phase diagram as follows. The self-consistent equations (A.6)
and (A.7) are initially solved numerically for a particular value of p and a set of values of s and λ
to find out corresponding free energy. By comparing these free energies and all possible solutions,
the stable phases of the system are identified with smallest value of free energy. The variation of
free energy with s for some values of p and λ = 0.3 is shown in Fig. 10. It can be seen that the free
energies for different values of p lie below fQP2 in the QP2 regime and, therefore, the QP2 phase
is not a stable phase. As we vary λ, the system undergoes a quantum phase transition from the
QP phase for small s to the Ferromagnetic (F) phase for large s.
To determine the type of phase transitions, i.e., first or second order, the magnetization mx
is numerically calculated as a function of s for the same parameter values as in the case of free
energy. We observe a change in mx around s = 0.4167 and mx decreases continuously to zero
from its unit value for p ≥ 5. Equivalently, it indicates that mz increases continuously from zero
to a finite value as we increase s for p ≥ 5. This identifies that for p ≥ 5 there exists a second
order phase transition at the boundary of QP and QP2 phases. An interesting scenario arises for
some parameter values (e.g., λ = 0.3, p = 11) where the magnetization shows a jump within the
ferromagnetic phase. This suggests the existence of first-order transition within the F phase and
the energy gap at the transition point decreases exponentially with the system size. Nevertheless,
this peculiar behavior does not appear for smaller values of λ for any non-zero p, except p = 3.
 1 1 1 19

0.8 0.8 0.8

[Link] Phil. Trans. R. Soc. A 0000000

F F F

..................................................................
0.6 0.6 0.6
s

0.4 0.4 0.4

0.2 QP 0.2 QP 0.2 QP

0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
λ

Figure 11. Phase diagrams on the s-λ plane for the system with Hamiltonian in Eq. (4.3) for p = 3 (left), p = 5 (middle),
and p = 11 (right). The boundary of the QP2 domain (s = 1/(3 − 2λ)), where a transition occurs between the QP and
F phases, is represented by the dash-dotted line. The red lines are for first-order transitions and the light green lines
represent second-order transitions. For p = 5 and 11, the magnetization shows sudden jumps on the dashed blue line
(F-F boundary) within the F phase. Taken from [115].

Therefore, for smaller λ, there exists only a second order transition when one increase s from zero
to a value near unity.
Using these results, phase diagrams of the system for p = 3, 5, and 11 are drawn on the s − λ
plane (see Fig. 11). We can see that a boundary of second-order transitions between F and QP
phases exists for small λ and p ≥ 5. As a consequence, there are possibilities to find a path to
reach the F phase from the QP phase by avoiding a first-order transition provided the first-order
F-F boundary does not reach the λ = 0 axis, that occurs probably in the limit of p → ∞ [115].
Let us now focus on analyzing the behavior of the energy gap across the phase transition points
of the system. The energy gap of the system is calculated numerically using perturbation theory
as described in Ref. [114]. The variation of energy gap with s for λ = 0.3 and p = 11 is shown
in Fig. 12. If the range of s where the energy gap has minimum value is zoomed, it can be seen
that the gap shows wiggly behavior throughout the range. This behavior starts at s ' 0.4184 for
λ = 0.3, which indeed corresponds to the second-order transition point between the QP and F
phases. The wiggly behavior ends at s ' 0.4676 for λ = 0.3, which corresponds to the first-order
transition point at the F-F boundary. The dashed vertical lines in Fig. 12 indicate two transition
points that are evaluated analytically using Eqs. A.6 and A.7. The analytical results show nearly
a good agreement with the numerical results in the interval where the gap is very small. It has
been found that the rightmost local minimum of the energy gap in Fig. 12 corresponding the F-F
boundary shows different scaling relation with the system size N compared to the other local
minima. The rightmost minimum energy gap decays exponentially with the system size, which
is expected from discontinuous behavior of the magnetization in Fig. 10 at the F-F boundary
implying the first-order transition. Although, for the present case, the above mentioned energy
gap is not a global minimum, this will affect the efficiency of QA for larger systems where the
rightmost gap can become a global minimum since the other local minima decay ploynomially
with the system sizes (see Figs. 6 and 7 of [115] for details).
These analytical and numerical results suggest that it is possible to increase the efficiency of
QA by choosing a path around λ = 0.1, which avoids first-order transition to reach the F phase
from the QP phase. For this process, s is the tuning parameter and the value p needs to be chosen
within the range 5 ≤ p ≤ 21 achieving maximum efficiency.
 1
N=20 1.0×100 20
N=80
N=140
Energy gap

Energy gap

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
0.5 1.0×10
-2

0 1.0×10-4
0.36 0.39 0.42 0.45 0.48 40 80 120 160
s N

Figure 12. Left panel: Energy gap as a function of s for p = 11 and λ = 0.3. The positions of minima of the energy gap
are shown by the vertical dashed lines at the QP2 domain with s ' 0.4167 and the F-F boundary at s ' 0.4701. Right
panel: The rightmost local minimum of the energy gap with the system size N for p = 11 and λ = 0.3 on a semi-log
scale. The energy gap closes exponentially fast with N . Taken from [115].

(c) Inhomogenous transverse field driving

In Sec. 4(b), we have discussed a method to speed up a QA process in the presence of first order
phase transitions by adding an antiferromagnetic fluctuation term. Here we consider a relatively
simpler approach of inhomogenous driving of the transverse field to overcome the issue of first
order phase transitions. In this case, the strength of transverse field is turned off sequentially from
one site to the next according to the annealing schedule. Using both analytical and numerical
calculations, it has been shown that inhomogenous driving can completely remove QPTs from
the annealing path and thus, it speeds up the annealing process exponentially [116,117].
The total Hamiltonian for the inhomogenous driving is given by

N (1−τ )
σix ,
X
H(s, τ ) = sH0 − (4.4)
i=1

where H0 is the Hamiltonian for p-spin model in Eq. (4.2). The parameters s and τ both are time-
dependent, where s = τ = 0 at t = 0 and s = τ = 1 at t = t0 . This shows that the initial Hamiltonian
has only transverse field and the final one has only p-spin interacting term with the Hamiltonian
H0 . Both the initial and final Hamiltonians are in agreement with the traditional QA protocol.
We note that the transverse field in Eq. (4.4) is applied only to N (1 − τ ) spins, where τ increases
from 0 to 1 as time proceeds from 0 to t0 . This indicates that the transverse field is turned off at
neighbouring sites one by one as time increases, starting from site i = N to ending with site i = 1
at τ = 1. This is the process of how the transverse field is driven inhomogeneously. It can be noted
that the parameter τ can have only discrete values for a finite N , since the upper limit N (1 − τ )
in Eq. (4.4) should be an integer.

(i) Results
Using Trotter decomposition and the static approximation in Hamiltonian (4.4), the free energy of
the system can be calculated analytically for both finite and zero temperatures (see appendix B).
By minimizing the zero-temperature free energy with respect to magnetization m produces a
ground state phase diagram as depicted in Fig. 13.
For a fixed value of p, a line of first-order phase transitions originated from a point on the
s-axis, terminates before approaching to any point of the τ -axis. Remarkably, all these lines for
different values of p end before they reach one of the axes, τ = 1 or s = 0. Therefore, there exists
a path starting from s = τ = 0 to s = τ = 1, that does not encounter any kind of phase transitions.
This leads to an exponential speedup of QA, since the energy gap always remains finite even
 1 21
p=3
p=4

[Link] Phil. Trans. R. Soc. A 0000000

0.8

..................................................................
p=5
p=6
p=7
0.6 p=8
τ

p=9
0.4

0.2

0
0 0.2 0.4 0.6 0.8 1
s

Figure 13. Phase diagram of the p-spin ferromagnetic model under inhomogenous field. All the lines represent first-order
phase transitions for different values of p, which are extended up to the middle of the phase diagram from the axis τ = 0
corresponding the homogenous model. Taken from [116].

6 ∆a 6 ∆a 6 N=5
5 ∆b 1 5 ∆b 1 5 N = 10
N = 100
Energy gap

Energy gap

Energy gap

4 4 4
3 3 3
2 2 2
1 1 1
0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
s s s

Figure 14. Plots of two types of energy gap ∆a1 and ∆b for p = 3 as functions of s with (left) τ = s (i.e., away from
the transition line) and (middle) τ = s2.366 (i.e., just touches the critical point). The final energy gap is defined by the
smaller of these two gaps. Right: The energy gap for different system sizes with τ = s computed by direct numerical
diagonalization. The location of the minimum energy gap is shown by an arrow for each N . Taken from [117].

for a system with large size. The positions of the critical points on the τ − s plane where the first-
order transitions terminate for different p values, are calculated analytically by using the standard
Landau theory of phase transitions (see Eq. (B.3)). In this calculation, it has been considered that
the coefficients of the expansion of the free energy (B.2) around its minimum at m = mc vanish to
third order [118].
To strengthen the above conclusion, the energy gap of the system has been calculated both
analytically and numerically. Since our system is mean-field-type, the semi-classical treatment
can be applied to evaluate the energy gap [119,120]. In this context, the parameterization of a
path τ = sr is considered to connect s = τ = 0 and s = τ = 1 with a parameter r that determines
the shape of the path. Figures i and i exhibit two energy gap candidates, ∆a1 and ∆b , for the
system with p = 3 along the paths τ = s, that does not encounter phase transitions, and τ = s2.366 ,
that just touches the critical point where the first-order line ends. The smaller one between these
two candidates is the actual energy gap of the system.
 As shown in Fig. i, ∆b is found to be smaller one and, it monotonically increases with s. On the 22
other hand, as expected, the energy gap ∆a1 vanishes at the critical point sc ≈ 0.52. To investigate
the effect of finite-size systems, the energy gap is calculated by a direct numerical diagonalization

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
method along the τ = s path. The result is shown in Fig. i, which shows a very good agreement
with the asymptotic behavior in the N → ∞ limit as observed in in ∆b of Fig. i. As seen in Fig. (i),
the energy gap becomes minimum when the transverse field is turned off at the first site as shown
by the arrows, thus implying the location of the minimum gap at s = 0 in the N → ∞ limit. It is
important to note that the the energy gap becomes minimum near the origin τ = s = 0, when
the annealing path (τ = s) does not encounter any transitions (see Figs. i and i ), whereas the
minimum of the gap occurs at the critical point when such a transition exists along the path ( see
Fig. i). In addition, a series of paths is considered to examine the inhomogenous driving protocol,
and it is found that the minimum energy gap shows an exponential decrement with the system
sizes when the paths cross the first order transitions (for details, refer to the discussion around
Eq. (12) of Ref. [117]).
The problem that we discussed so far, for inhomogeneous driving, considers ideal situations,
i.e., zero temperature and a complete turning off of the transverse field at each site. This
problem also has been studied at finite temperature and zero temperature with different types
of inhomogenity. It has already been studied that the first-order transitions that exists under
homogeneous driving at zero temperature can be circumvented by inhomogeneous driving
with complete turning off the transverse field at each site. For nonideal situations, with a finite
temperature or a non-zero value of the final transverse field, one can not avoid new first-order
transition lines like the ideal case. Nevertheless, it has been observed that the new first-order
transitions are weaker than the original one, since the free energy barrier between two arbitrary
local minima is smaller than the original homogeneous case. This leads to an increase in the
quantum tunneling rate. Therefore one can infer that the inhomogeneous driving of the transverse
field has the potential to provide a better performance in quantum annealing. In addition,
Matsuura et al. [121] studied analytically the p-body ferromagnetic infinite-range Ising model
in transverse-field using a mean-field analysis and demonstrated that for p ≥ 3, where the phase
transition is of first order, Quantum Annealing Correction softens the closing of the gap for small
energy penalty values and prevents its closure for sufficiently large energy penalty values, thereby
providing from excitations that occur near the quantum critical point. It also has been shown
analytically that nested quantum annealing correction can suppress errors effectively in Ising
models with infinite-range interactions and their analysis revealed that the nesting structure can
significantly weaken or remove the first-order phase transitions, where the energy gap closes
exponentially [122].

(d) Suppression of Griffiths singularities

In a recent work by Knysh et al. [123], it has been shown that a QA process can be accelerated
using an embedded spin chain system with random interactions. A randomly interacting spin
chain exhibits Griffiths-McCoy singularities [124,125], since different parts of the system can not
reach criticality simultaneously for random fluctuations. This leads to the diverging dynamical
exponent z and a stretched exponential scaling of the energy gap [36]. Therefore the presence of
Griffiths singularities increases the annealing time for such systems.
On the other hand, quantum annealing has been studied for a embedded spin chain problem,
where logical qubits were replaced by ferromagnetic Ising spin chains [123]. For this study, an
ansatz is considered to find a balanced choice of coupling parameters based on renormalization
group intuition for the better performance of QA. This results an exponential improvement of
annealing time, which is also confirmed numerically. It indicates that this protocol prevents to
occur randomly oriented domains in the system by ensuring a simultaneous criticality of spatially
separated regions.
5. Application 23
In this section, we discuss some recent results on QA of a initeracting bosonic system coupled

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
with cavity field modes for both classical and quantum limits of the system. In addition, we
will mention about the recent development of QA in context of parallel computation, and its
effectiveness compared to the existing methods.

(a) Quantum annealing vs. semi-classical annealing

Starchl and Ritsch [17], have used the idea of quantum tunneling to show the success of QA
over semi-classical annealing for an interacting bosonic model in presence of cavity modes. The
model that is considered for this study is described by a tight binding Bose-Hubbard lattice model
with four sites, which are filled by two interacting bosons. The tunable non-local interactions are
introduced in the model via collective light scattering to two independent cavity modes. The
Hamiltonian for such system is given by

P † U P
H =J kP BC (bk bk+1 + h.c.) + 2 k nk (nk − 1)

∆(a†1 a1 + a†2 a2 ) + J˜ M̂1 (a1 + a†1 ) + M̂2 (a2 + a†2 ) ,

− (5.1)

where bk and b†k are bosonic annihilation and creation operators, respectively. The photonic
annihilation operators a1 and a2 are associated with two independent cavity modes. The
interactions between the bosons and cavity field modes are represented by the fourth term
of Eq. (5.1), where M̂1 and M̂2 are called effective scattering operators. Using mean-field
approximation of the field operators, the Hamiltonian in Eq. (5.1) can be written in semiclassical
form
∆J˜2
H sc = 2 ˆ M̂1 i2 + 2M̂2 hM̂2 i − Ih
ˆ M̂2 i2 ,


2M̂1 hM̂1 i − Ih (5.2)
κ + ∆2

where κ determines the strength of cavity loss, where I denotes the identity operator. Both
the Hamiltonians in Eqs. (5.1) and (5.2) with periodic boundary conditions are translationally
invariant and thus provide approximately degenerate ground states. In order to create an unique
target ground state for the annealing process, a certain amount of impurity of strength V is added
in the Hamiltonian.

Figure 15. Color density plot of the fidelity calculated as the overlap between the final state after an adiabatic evolution
using (a) the Hamiltonian with semi-classical mean-field approximation and as well as (b) for the full quantum Hamiltonian,
and the desired target state on the U − V plane. Here the parameter values are: tf = 1000. For (a): ∆ = −1, J˜ = 1
√
and for (b): ∆ = −5, J˜ = 5. Taken from [17].
 The dynamics of the system is started with the ground state at zero pump J˜ = 0, and the pump
√ 24
strength is increased linearly towards J˜ = 5, following an adiabatic schedule: J˜ ≈ t/tf , where
tf is the final time. The results of the study of annealing for this system are summarized in Fig. 15.

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
It shows the density plot of the fidelity on U − V plane after an adiabatic sweep using both the
semi-classical (5.1) and the full quantum (5.2) Hamiltonian. In this case, the fidelity is calculated
as the overlap between the final state after an adiabatic sweep and the desired target state. By
comparing two density plots for both semi-classical and quantum cases, i.e., when the dynamics
is driven by the semi-classical Hamiltonian and the full quantum Hamiltonian, respectively, one
can identify a clear quantum improvement in the success rate. It is observed that the semi-classical
approximation provides a reliably correct solution for small onsite interaction strength. For this
scenario, one can see a sudden fall of the success rate, when the interaction strength is increased.
If the repulsive interaction is further increased, the gap between the final ground state energy and
the first excitation becomes very small and the classical model effectively never succeeds. On the
other hand, the adiabatic evolution with full quantum Hamiltonian (5.1) provides almost correct
solution with 99% probability, even for higher U values as shown in Fig. 15(b). Therefore, for this
system, a large parameter region is found where quantum annealing is highly successful, whereas
the semi-classical approach largely fails. In addition, for quantum scenario, a direct connection is
found between atom-field entanglement in the dynamics and a high probability to find the correct
solution at end of annealing process.

(b) Parallel quantum annealing

There are a few recent studies on QA in the context of parallel quantum computation [126–128].
Recently, Pelofske et al. [128] propose a method named as parallel quantum annealing that
has potential to solve many independent problems on a quantum annealer during the same
annealing process. The authors applied their proposed method of parallel quantum annealing
on both D-Wave 2000Q at Los Alamos National Laboratory (refereed to as D-W 2000Q) and the
newer D-Wave Advantage_System 1.1 (referred to as D-W Advantage). The results of parallel
quantum annealing have been compared with those found from sequential quantum annealing,
i.e., when the same problems are solved sequentially on D-Wave machine. It has been observed
that, although, there is a slight decrement in the accuracy of the solution for simultaneously solved
problems, parallel quantum annealing can provide a considerable speedup of up to two orders of
magnitude [128].

6. Summary and outlook

We have provided an overview of recent developments of QA which is based on the possible
advantage of utilizing quantum tunneling. When the energy landscape of an Ising Hamiltonian,
where the corresponding ground state is the target state of an optimisation problem, consists
of high but thin barriers surrounding local minima, quantum tunneling has an advantage over
thermal fluctuation in overcoming barriers and thus getting the system equilibrated. This nature
of quantum tunneling provides a foundation which asserts that QA can outperform SA in a glassy
system with a rugged energy landscape. Indeed, we have focused on analytical and numerical
evidences that QA yields a better solution than SA in several glassy systems. The question of the
restoration of ergodicity due to quantum tunneling in the quantum SK model is still unresolved.
Nevertheless, recent studies do indicate the possibility of the existence of an ergodic phase at
least in low temperature region [62,64]: this is expected to lead to a remarkable possibility of the
success of the annealing scheme in those systems [61].
We have also discussed quantum phase transitions (QPTs) in connection with QA. Generally
speaking, a QPT hinders the adiabatic time evolution underlying QA, since the energy gap above
the ground state closes at a QPT; this leads to inevitable generation of defects and excess energy in
the final evolved state. In this context, it is worth noting that the short-cut to adiabaticity [129] or
counter-diabatic driving protocols [130] have been proposed as a method to realize the adiabatic
time evolution with a finite time, providing a possible route to avoid a continuous QPT in 25
the process of annealing. On the contrary, a discontinuous QPT involves an exponentially fast
closure of the energy gap with the system size. Therefore, it is desirable and at the same time far

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
more challenging to circumvent a discontinuous QPT for the success of QA. We have reviewed
methods using additional antiferromagnetic interactions [123] or inhomogeneous transverse
fields [116,117] proposed for avoiding discontinuous QPTs which would lead to acceleration
of QA. Random interactions in spin chain systems induce Griffiths singularities that eventually
increase annealing time for such systems. These singularities can be suppressed for the embedded
spin chain problem, where ferromagnetic Ising spin chains are used as logical qubits. It has
been shown that the performance of QA for such embedded systems improves exponentially
in the context of annealing time [123]. To be precise, the performance advantage of QA is still
model specific and a generic prescription has not been known for discontinuous phase transitions
involved by practical optimization problems. In the cases of spin glass models, however, major
advantages of the standard QA have now even been established using QMC simulation [79] and
hardware implementations [85,131,132].
The idea of quantum tunneling has been further used by Starchl and Ritsch [17] to establish
the superiority of QA over semi-classical annealing for a realistic system of interacting bosons
in presence of cavity modes. Using numerical results, it has been shown that there exists a large
parameter regime where the QA provides a better performance than semi-classical annealing.
In addition, we have noted the recent development of parallel quantum computation using an
annealing algorithm [128]. The idea of parallel QA is to solve many independent problems on a
quantum annealer during the same annealing schedule. The authors have checked their method
of parallel QA on D-wave quantum machines and indicated the effectiveness of the same [128].
Finally, in recent years the progress in developing hardware that performs QA using physical
qubits has gained a tremendous momentum. So far, devices with more than 5000 qubits have been
made available, and employed to study a wide gamut of fields that include condensed matter
systems in and out of equilibrium [93,131,133–135], high-energy physics [136–138], quantum
chemistry [139], and biology [140,141]. Application to various optimization problems has been
developing as well. We reviewed some of the experimental studies using a QA hardware. From
the viewpoint of application as well as gaining theoretical rigor, decoherence inherent in a device
coupled to an environment is a fundamental issue of interest. We have briefly reviewed effects
of thermal and non-thermal environments on QA. In order to perform QA ideally, coupling to
an environment leading to decoherence should be reduced. However, environment assisted QAs
have been proposed for specific situations [142–144]. Utilising specially engineered environments
to accelerate QA would be an important direction of future study.

A. Static approximation and low-temperature limit

Using Suzuki-Trotter formula, the partition function for the Hamiltonian in Eq. (4.3) can be
written as

Z = lim ZM
M →∞
β β
M
= lim Tr e− M sλH0 e− M {s(1−λ)VAFF +(1−s)VTF }
M →∞
 N
h βsλN  1 X p i
h{σiz }| exp σiz
X
= lim
M →∞ M N
{σiz } i=1

N
h βs(1 − λ)N  1 X N
2 β(1 − s) X iM
× exp − σix + σix |{σiz }i, (A.1)
M N M
i=1 i=1
 configurations in the z basis, and |{σiz }i ≡
P
where {σiz } denotes the summation over all spin 26
NN z z z
i=1 i i. The state |σi i is the eigenstate of σi ,
|σ having the eigenvalue σiz (= ±1). Similar
notations will be used for the x basis.

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
Following saddle point method in the limit N → ∞ and static approximation (i.e., neglecting
the imaginary-time dependence of the partition function in Eq. (A.1) [114,115,145]), the partition
function of the system can be written as
ZZ
Z = dmz dmx exp[−N βf (β, s, λ; mz , mx )], (A.2)

where f (β, s, λ; mz , mx ) is the pseudo free energy defined as follows:

f (β, s, λ; mz , mx ) = (p − 1)sλ(mz )p − s(1 − λ)(mx )2

q
1  2  2
− ln 2 cosh β psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx . (A.3)
β
The saddle point equations are thus

psλ(mz )p−1
mz = q
2  2
psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx
q
 2  2
× tanh β psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx , (A.4)
x
1 − s − 2s(1 − λ)m
mx = q
2  2
psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx
q
 2  2
× tanh β psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx . (A.5)

To examine quantum phase transitions of the model, we consider low-temperature limits of

the above self-consistent equations. For a finite value of the square root in Eq. (A.4) and Eq. (A.5),
the hyperbolic tangent tends to unity in the limit of β → ∞. Then the equations are given as

psλ(mz )p−1
mz = q , (A.6)
2  2
psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx
1 − s − 2s(1 − λ)mx
mx = q . (A.7)
2  2
psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx

In that case, the pseudo free energy (A.3) becomes

f (s, λ; mz , mx ) = (p − 1)sλ(mz )p − s(1 − λ)(mx )2

q
 2  2
− psλ(mz )p−1 + 1 − s − 2s(1 − λ)mx . (A.8)

Equations (A.6) and (A.7) provide a ferromagnetic (F) solution with mz > 0 and a quantum
paramagnetic (QP) solution for mz = 0 and mx 6= 0. Using these properties of a quantum
paramagnetic phase, the regions of QP phases can be found on the s − λ plane. It appears that
there exists two types of QP phases in this problem and we call them QP and QP2 phases to
distinguish from each other.
The regions of the different phases in terms of system parameters can be calculated using the
above conditions of those phases in Eqs. (A.6) and (A.7). It has been found that the QP phase
exists in the region 0 ≤ s < 1/(3 − 2λ), and its free energy is given by

fQP (s, λ) = −sλ + 2s − 1, (A.9)

which is independent of p. The free energy of the F phase can not be calculated analytically for a
general p from Eqs. (A.6) and (A.7). However, in the limit of p → ∞, the free energy of the F phase
is given as 27

fF (s, λ)|p→∞ = − sλ. (A.10)

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
The free energy of the QP2 phase is given by

(1 − s)2
fQP2 (s, λ) = − , (A.11)
4s(1 − λ)

with the domain of applicability is restricted by 1/(3 − 2λ) ≤ s < 1.

B. Inhomogenous driving of the transverse field

(a) Free energy and first-order transitions
By applying Trotter decomposition and the static approximation on the Hamiltonian in Eq. (4.4),
the resulting free energy at finite temperature is given by [116]

f (m; s, τ )
 q 
p p−1 2
= (1 − τ ) (p − 1)sm − T log 2 cosh β (spm ) +1
n o
+ τ (p − 1)smp − T log 2 cosh(βspmp−1 ) , (B.1)

where m is the magnetization of the system along the z axis. In the limit of zero temperature, the
free energy takes form
 q 
f0 (m; s, τ ) = (1 − τ ) (p − 1)smp − (spmp−1 )2 + 1
n o
+ τ (p − 1)smp − spmp−1 . (B.2)

For the calculation of zero-temperature free energy it has been assumed that m ≥ 0.
Using the standard Landau theory of phase transitions, the locations of the critical points sc ,
τc (see Fig. 13) where the first-order transition lines terminate for different p, can be found as
1 1
τc = , sc = , (B.3)
pmp−1
s p
27(p − 1) c 1 − m1c 2 /m1c
1+
4(p − 2)3
p
where m1c = (p − 2)/(3(p − 1)) and mc = τc + (1 − τc )m1c .

(b) Semiclassical theory of energy gap

One can rewrite the Hamiltonian (4.4) in terms of two macroscopic spin operators (for details see
Refs. [116,117]),
N (1−sr ) N
1 1
S1z,x σiz,x , S2z,x = σiz,x
X X
= (B.4)
2 2
i=1 i=N (1−sr )+1

as

p
 
2 z
H(s, τ ) = −sN S1 + S2z − 2S1x . (B.5)
N

These giant operators can be considered as classical vectors for sufficiently large N , and the
quantum fluctuations are subsequently applied around the classically stable directions through
an expansion of the Holstein-Primakoff transformation to the quadratic order in terms of boson
operators, as done in Refs. [119,120]. The result is given as 28
δ p
H(s, τ ) =N e + γ + ( 1 − 2 − 1)

[Link] Phil. Trans. R. Soc. A 0000000

2

..................................................................
+ ∆a1 ã†1 ã1 + ∆a2 ã†2 ã2 + ∆b b† b, (B.6)

where ã1 and ã2 are bosonic annihilation operators, and e is the energy per spin of the classical
ground-state. The parameters ∆a1 , ∆a2 and ∆b represent quantum fluctuations, where
p
∆a1 = δ 1 − 2 , ∆a2 = δ. (B.7)

Because ∆a2 ≥ ∆a1 , the minimum energy gap of the system is the smaller of ∆a1 and ∆b :

∆ = min(∆a1 , ∆b )
p
∆a1 = δ 1 − 2 , ∆b = 2sp{τ + (1 − τ ) cos θ0 }p−1 , (B.8)

where

θ0 = arg min −s[τ + (1 − τ ) cos θ]p − (1 − τ ) sin θ


θ
2γ
=− ,
δ
1
γ = − sp(p − 1)(1 − τ ) sin2 θ0 {τ + (1 − τ ) cos θ0 }p−2 ,
2
δ = ∆b cos θ0 + 2 sin θ0 + 2γ. (B.9)

Ethics. NA.
Data Accessibility. The article has no additional data over those given in different figures taken from
published papers.

Authors’ Contributions. AD and BKC conceptualized the review. AR and SS contributed to materials and
organised it. All authors contributed to editing and finalizing the review.

Competing Interests. We declare, we do not have any competing interests.

Funding. A.R. acknowledges UGC, India for Start-up Research Grant No. F. 30-509/2020(BSR). The work
of S.S. was supported by JSPS KAKENHI Grant No. 22K03455. A.D. acknowledges support from SPARC
program, MHRD, India and SERB, DST, New Delhi, India. BKC is grateful to the Indian National Science
Academy for their Senior Scientist Research Grant.

Acknowledgements. We acknowledge our collaborations with Gabriel Aeppli, Arunava Chakrabarti,

Arnab Das, Uma Divakaran, Jun-ichi Inoue, Sudip Mukherjee, Hidetoshi Nishimori, Masato Okada,
Purusattam Ray, Thomas F. Rosenbaum, Diptiman Sen, Parongama Sen, Robin B. Stinchcombe, Ryo Tamura
and Shu Tanaka on this development. AD acknowledges Souvik Bandyopadhyay and Sourav Bhattacharjee
for comments. We are grateful to the anonymous referees for their useful comments and important
suggestions.

Disclaimer. This review is limited by our personal knowledge and also by the size limit (which we have
already crossed). We do not claim any completeness of discussions on even some important contributions in
this incredibly active field of research.

References
1. A. B. Finnila, M. A. Gomez, C. Sebenik, et al.
Quantum annealing: A new method for minimizing multidimensional functions.
Chem. Phys. Lett., 219, 343 (1994).
2. T. Kadowaki and H. Nishimori.
Quantum annealing in the transverse Ising model.
 Phys. Rev. E, 58, 5355 (1998). 29
3. E. Farhi, J. Goldstone, S. Gutmann, et al.
A Quantum Adiabatic Evolution Algorithm Applied to Random Instances of an NP-

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
Complete Problem.
Science, 292, 472 (2001).
4. G. E. Santoro and E. Tosatti.
Optimization using quantum mechanics: quantum annealing through adiabatic evolution.
J. Phys. A: Math. Gen., 39, R393 (2006).
5. A. Das and B. K. Chakrabarti.
Colloquium: Quantum annealing and analog quantum computation.
Rev. Mod. Phys., 80, 1061 (2008).
6. S. Morita and H. Nishimori.
Mathematical foundation of quantum annealing.
J. Math. Phys., 49, 125210 (2008).
7. S. Tanaka, R. Tamura, and B. K. Chakrabarti.
Quantum spin glasses, annealing and computation.
Cambridge University Press (2017).
8. T. Albash and D. A. Lidar.
Adiabatic quantum computation.
Rev. Mod. Phys., 90, 015002 (2018).
9. P. Hauke, H. G. Katzgraber, W. Lechner, et al.
Perspectives of quantum annealing: Methods and implementations.
Rep. Prog. Phys., 83, 054401 (2020).
10. A. Das, B. K. Chakrabarti, and R. B. Stinchcombe.
Quantum annealing in a kinetically constrained system.
Phys. Rev. E, 72, 026701 (2005).
11. A. Messiah.
Quantum Mechanics.
North Holland, Amsterdam (1961).
12. S. Suzuki and M. Okada.
Residual energies after slow quantum annealing.
J. Phys. Soc. Jpa., 74, 1649 (2005).
13. S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi.
Optimization by simulated annealing.
science, 220, 671 (1983).
14. R. Shankar.
Principles of quantum mechanics.
Springer (2012).
15. S. Mukherjee and B. K. Chakrabarti.
Multivariable optimization: Quantum annealing and computation.
Eur. Phys. J. ST, 224, 17 (2015).
16. P. Ray, B. K. Chakrabarti, and A. Chakrabarti.
Sherrington-Kirkpatrick model in a transverse field: Absence of replica symmetry breaking
due to quantum fluctuations.
Phys. Rev. B, 39, 11828 (1989).
17. E. A. Starchl and H. Ritsch.
Unraveling the origin of higher success probabilities in quantum annealing versus semi-
classical annealing.
J. Phys. B: Atomic, Molecular and Optical Physics (2022).
18. D. Sherrington and S. Kirkpatrick.
Solvable model of a spin-glass.
Phys. Rev. Lett., 35, 1792 (1975).
19. R. Yaacoby, N. Schaar, L. Kellerhals, et al.
Comparison between a quantum annealer and a classical approximation algorithm for
computing the ground state of an Ising spin glass.
Phys. Rev. E, 105, 035305 (2022).
20. S. L. Sondhi, S. Girvin, J. Carini, et al.
Continuous quantum phase transitions.
Rev. Mod. Phys., 69, 315 (1997).
21. S. Sachdev. 30
Quantum Phase Transitions.
Cambridge University Press, Cambridge (1999).

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
22. S. Suzuki, J.-i. Inoue, and B. K. Chakrabarti.
Quantum Ising phases and transitions in transverse Ising models, volume 862.
Springer (2012).
23. A. Dutta, G. Aeppli, B. K. Chakrabarti, et al.
Quantum phase transitions in transverse field spin models: from statistical physics to quantum
information.
Cambridge University Press (2015).
24. T. W. B. Kibble.
Topology of cosmic domains and strings.
J. Phys. A: Math. Gen., 9, 1387 (1976).
25. W. H. Zurek.
Cosmological experiments in superfluid helium?
Nature, 317, 505 (1985).
26. W. H. Zurek, U. Dorner, and P. Zoller.
Dynamics of a Quantum Phase Transition.
Phys. Rev. Lett., 95, 105701 (2005).
27. A. Polkovnikov.
Universal adiabatic dynamics in the vicinity of a quantum critical point.
Phys. Rev. B, 72, 161201 (2005).
28. B. Damski.
The simplest quantum model supporting the Kibble-Zurek mechanism of topological defect
production: Landau-Zener transitions from a new perspective.
Phys. Rev. Lett., 95, 035701 (2005).
29. V. Mukherjee, U. Divakaran, A. Dutta, et al.
Quenching dynamics of a quantum X Y spin-1 2 chain in a transverse field.
Phys. Rev. B, 76, 174303 (2007).
30. D. Sen, K. Sengupta, and S. Mondal.
Defect production in nonlinear quench across a quantum critical point.
Phys. Rev. Lett., 101, 016806 (2008).
31. R. Barankov and A. Polkovnikov.
Optimal nonlinear passage through a quantum critical point.
Phys. Rev. Lett., 101, 076801 (2008).
32. V. Mukherjee and A. Dutta.
Adiabatic multicritical quantum quenches: Continuously varying exponents depending on
the direction of quenching.
EPL (Europhysics Lett.), 92, 37004 (2010).
33. J. Dziarmaga.
Dynamics of a quantum phase transition and relaxation to a steady state.
Adv. Phys., 59, 1063 (2010).
34. A. Polkovnikov, K. Sengupta, A. Silva, et al.
Colloquium: Nonequilibrium dynamics of closed interacting quantum systems.
Rev. Mod. Phys., 83, 863 (2011).
35. J. Dziarmaga.
Dynamics of a Quantum Phase Transition: Exact Solution of the Quantum Ising Model.
Phys. Rev. Lett., 95, 245701 (2005).
36. D. S. Fisher.
Random transverse field Ising spin chains.
Phys. Rev. Lett., 69, 534 (1992).
37. D. S. Fisher.
Critical behavior of random transverse-field Ising spin chains.
Phys. Rev. B, 51, 6411 (1995).
38. A. P. Young and H. Rieger.
Numerical study of the random transverse-field Ising spin chain.
Phys. Rev. B, 53, 8486 (1996).
39. D. S. Fisher and A. P. Young.
 Distributions of gaps and end-to-end correlations in random transverse-field Ising spin 31
chains.
Phys. Rev. B, 58, 9131 (1998).

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
40. J. Dziarmaga.
Dynamics of a quantum phase transition in the random Ising model: Logarithmic
dependence of the defect density on the transition rate.
Phys. Rev. B, 74, 064416 (2006).
41. T. Caneva, R. Fazio, and G. E. Santoro.
Adiabatic quantum dynamics of a random Ising chain across its quantum critical point.
Phys. Rev. B, 76, 144427 (2007).
42. S. Suzuki.
Cooling dynamics of pure and random Ising chains.
J. Stat. Mech. Theory Exp., 2009, P03032 (2009).
43. M. Schmitt, M. M. Rams, J. Dziarmaga, et al.
Quantum phase transition dynamics in the two-dimensional transverse-field Ising model
(2021).
44. R. Botet and R. Jullien.
Large-size critical behavior of infinitely coordinated systems.
Phys. Rev. B, 28, 3955 (1983).
45. O. L. Acevedo, L. Quiroga, F. J. Rodríguez, et al.
New Dynamical Scaling Universality for Quantum Networks Across Adiabatic Quantum
Phase Transitions.
Phys. Rev. Lett., 112, 030403 (2014).
46. S. F. Edwards and P. W. Anderson.
Theory of spin glasses.
J. Phys. F: Metal Phys., 5, 965 (1975).
47. B. K. Chakrabarti.
Critical behavior of the Ising spin-glass models in a transverse field.
Phys. Rev. B, 24, 4062 (1981).
48. H. Rieger and A. P. Young.
Zero-temperature quantum phase transition of a two-dimensional Ising spin glass.
Phys. Rev. Lett., 72, 4141 (1994).
49. M. Guo, R. N. Bhatt, and D. A. Huse.
Quantum critical behavior of a three-dimensional Ising spin glass in a transverse magnetic
field.
Phys. Rev. Lett., 72, 4137 (1994).
50. R. R. P. Singh and A. P. Young.
Critical and Griffiths-McCoy singularities in quantum Ising spin glasses on d-dimensional
hypercubic lattices: A series expansion study.
Phys. Rev. E, 96, 022139 (2017).
51. R. Miyazaki and H. Nishimori.
Real-space renormalization-group approach to the random transverse-field Ising model in
finite dimensions.
Phys. Rev. E, 87, 032154 (2013).
52. D. A. Matoz-Fernandez and F. Romá.
Unconventional critical activated scaling of two-dimensional quantum spin glasses.
Phys. Rev. B, 94, 024201 (2016).
53. O. Motrunich, S.-C. Mau, D. A. Huse, et al.
Infinite-randomness quantum Ising critical fixed points.
Phys. Rev. B, 61, 1160 (2000).
54. D. Karevski, Y.-C. Lin, H. Rieger, et al.
Random quantum magnets with broad disorder distribution.
Eur. Phys. J. B, 20, 267 (2001).
55. Y.-C. Lin, F. Iglói, and H. Rieger.
Entanglement Entropy at Infinite-Randomness Fixed Points in Higher Dimensions.
Phys. Rev. Lett., 99, 147202 (2007).
56. F. Iglói and C. Monthus.
Strong disorder RG approach of random systems.
 Phys. Rep., 412, 277 (2005). 32
57. F. Iglói and C. Monthus.
Strong disorder RG approach – a short review of recent developments.

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
Eur. Phys. J. B, 91, 290 (2018).
58. G. Parisi.
The order parameter for spin glasses: a function on the interval 0-1.
J. Phys. A: Math. Gen., 13, 1101 (1980).
59. K. Binder and A. P. Young.
Spin glasses: Experimental facts, theoretical concepts, and open questions.
Rev. Mod. Phys., 58, 801 (1986).
60. A. P. Young.
Stability of the quantum Sherrington-Kirkpatrick spin glass model.
Phys. Rev. E, 96, 032112 (2017).
61. S. Mukherjee, A. Rajak, and B. K. Chakrabarti.
Possible ergodic-nonergodic regions in the quantum Sherrington-Kirkpatrick spin glass
model and quantum annealing.
Phys. Rev. E, 97, 022146 (2018).
62. S. Mukherjee, A. Rajak, and B. K. Chakrabarti.
Classical-to-quantum crossover in the critical behavior of the transverse-field Sherrington-
Kirkpatrick spin glass model.
Phys. Rev. E, 92, 042107 (2015).
63. A. Rajak and B. Chakrabarti.
Quantum annealing search of Ising spin glass ground state (s) with tunable transverse and
longitudinal fields.
Indian J. Phys., 88, 951 (2014).
64. H. Leschke, C. Manai, R. Ruder, et al.
Existence of Replica-Symmetry Breaking in Quantum Glasses.
Phys. Rev. Lett., 127, 207204 (2021).
65. P. M. Schindler, T. Guaita, T. Shi, et al.
A Variational Ansatz for the Ground State of the Quantum Sherrington-Kirkpatrick Model.
arXiv preprint arXiv:2204.02923 (2022).
66. A. P. Young, S. Knysh, and V. N. Smelyanskiy.
First-Order Phase Transition in the Quantum Adiabatic Algorithm.
Phys. Rev. Lett., 104, 020502 (2010).
67. T. Jörg, F. Krzakala, G. Semerjian, et al.
First-order transitions and the performance of quantum algorithms in random optimization
problems.
Phys. Rev. Lett., 104, 207206 (2010).
68. E. Barouch, B. M. McCoy, and M. Dresden.
Statistical mechanics of the XY model. I.
Phys. Rev. A, 2, 1075 (1970).
69. S. R. White and A. E. Feiguin.
Real-Time Evolution Using the Density Matrix Renormalization Group.
Phys. Rev. Lett., 93, 076401 (2004).
70. G. Vidal.
Efficient simulation of one-dimensional quantum many-body systems.
Phys. Rev. Lett., 93, 040502 (2004).
71. R. Orús and G. Vidal.
Infinite time-evolving block decimation algorithm beyond unitary evolution.
Phys. Rev. B, 78, 155117 (2008).
72. F. Pollmann, S. Mukerjee, A. G. Green, et al.
Dynamics after a sweep through a quantum critical point.
Phys. Rev. E, 81, 020101 (2010).
73. A. D. King, S. Suzuki, J. Raymond, et al.
Coherent quantum annealing in a programmable 2000-qubit Ising chain.
arXiv preprint arXiv:2202.05847 (2022).
74. H. F. Trotter.
On the product of semi-groups of operators.
Proc. Am. Math. Soc., 10, 545 (1959).
75. M. Suzuki. 33
Quantum statistical Monte Carlo methods and applications to spin systems.
J. Stat. Phys., 43, 883 (1986).

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
76. R. H. Swendsen and J.-S. Wang.
Nonuniversal critical dynamics in Monte Carlo simulations.
Phys. Rev. Lett., 58, 86 (1987).
77. T. Nakamura and Y. Ito.
A quantum Monte Carlo algorithm realizing an intrinsic relaxation.
J. Phys. Soc. Jpn., 72, 2405 (2003).
78. T. Kadowaki.
Study of Optimization Problems by Quantum Annealing (2002).
79. G. E. Santoro, R. Martonák, E. Tosatti, et al.
Theory of quantum annealing of an Ising spin glass.
Science, 295, 2427 (2002).
80. R. Martoňák, G. E. Santoro, and E. Tosatti.
Quantum annealing by the path-integral Monte Carlo method: The two-dimensional random
Ising model.
Phys. Rev. B, 66, 094203 (2002).
81. R. Martoňák, G. E. Santoro, and E. Tosatti.
Quantum annealing of the traveling-salesman problem.
Phys. Rev. E, 70, 057701 (2004).
82. D. A. Battaglia, G. E. Santoro, and E. Tosatti.
Optimization by quantum annealing: Lessons from hard satisfiability problems.
Phys. Rev. E, 71, 066707 (2005).
83. S. Suzuki.
Kibble-Zurek mechanism in simulated annealing and quantum annealing.
In J. Phys.: Conf. Ser., volume 302, page 012046. IOP Publishing (2011).
84. W. Wu, B. Ellman, T. F. Rosenbaum, et al.
From classical to quantum glass.
Phys. Rev. Lett., 67, 2076 (1991).
85. J. Brooke, D. Bitko, Rosenbaum, et al.
Quantum annealing of a disordered magnet.
Science, 284, 779 (1999).
86. S. Säubert, C. L. Sarkis, F. Ye, et al.
Microscopics of Quantum Annealing in the Disordered Dipolar Ising Ferromagnet
LiHox Y1−x F4 .
arXiv preprint arXiv:2105.03408 (2021).
87. A. Keesling, A. Omran, H. Levine, et al.
Quantum Kibble–Zurek mechanism and critical dynamics on a programmable Rydberg
simulator.
Nature, 568, 207 (2019).
88. M. W. Johnson, M. H. Amin, S. Gildert, et al.
Quantum annealing with manufactured spins.
Nature, 473, 194 (2011).
89. V. S. Denchev, S. Boixo, S. V. Isakov, et al.
What is the Computational Value of Finite-Range Tunneling?
Phys. Rev. X, 6, 031015 (2016).
90. S. Boixo, T. F. Rønnow, S. V. Isakov, et al.
Evidence for quantum annealing with more than one hundred qubits.
Nat. Phys., 10, 218 (2014).
91. S. W. Shin, G. Smith, J. A. Smolin, et al.
How" quantum" is the D-Wave machine?
arXiv preprint arXiv:1401.7087 (2014).
92. S. Suzuki and A. Das.
Quantum annealing: The fastest route to quantum computation? (2015).
93. Y. Bando, Y. Susa, H. Oshiyama, et al.
Probing the universality of topological defect formation in a quantum annealer: Kibble-
Zurek mechanism and beyond.
Phys. Rev. Research, 2, 033369 (2020).
 94. D. Venturelli, S. Mandrà, S. Knysh, et al. 34
Quantum optimization of fully connected spin glasses.
Physical Review X, 5, 031040 (2015).

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
95. A. O. Caldeira and A. J. Leggett.
Quantum tunnelling in a dissipative system.
Ann. Phys. (N Y), 149, 374 (1983).
96. A. J. Leggett, S. Chakravarty, A. T. Dorsey, et al.
Dynamics of the dissipative two-state system.
Rev. Mod. Phys., 59, 1 (1987).
97. M. H. Amin.
Searching for quantum speedup in quasistatic quantum annealers.
Phys. Rev. A, 92, 052323 (2015).
98. D. Patanè, A. Silva, L. Amico, et al.
Adiabatic Dynamics in Open Quantum Critical Many-Body Systems.
Phys. Rev. Lett., 101, 175701 (2008).
99. D. Patanè, L. Amico, A. Silva, et al.
Adiabatic dynamics of a quantum critical system coupled to an environment: Scaling and
kinetic equation approaches.
Phys. Rev. B, 80, 024302 (2009).
100. P. Nalbach, S. Vishveshwara, and A. A. Clerk.
Quantum Kibble-Zurek physics in the presence of spatially correlated dissipation.
Phys. Rev. B, 92, 014306 (2015).
101. A. Dutta, A. Rahmani, and A. del Campo.
Anti-Kibble-Zurek Behavior in Crossing the Quantum Critical Point of a Thermally Isolated
System Driven by a Noisy Control Field.
Phys. Rev. Lett., 117, 080402 (2016).
102. P. Weinberg, M. Tylutki, J. M. Rönkkö, et al.
Scaling and Diabatic Effects in Quantum Annealing with a D-Wave Device.
Phys. Rev. Lett., 124, 090502 (2020).
103. M. Hasenbusch, K. Pinn, and S. Vinti.
Critical exponents of the three-dimensional Ising universality class from finite-size scaling
with standard and improved actions.
Phys. Rev. B, 59, 11471 (1999).
104. M. Karl and T. Gasenzer.
Strongly anomalous non-thermal fixed point in a quenched two-dimensional Bose gas.
New J. Phys., 19, 093014 (2017).
105. M. Keck, S. Montangero, G. E. Santoro, et al.
Dissipation in adiabatic quantum computers: lessons from an exactly solvable model.
New J. Phys., 19, 113029 (2017).
106. S. Bandyopadhyay, S. Bhattacharjee, and A. Dutta.
Dynamical generation of Majorana edge correlations in a ramped Kitaev chain coupled to
nonthermal dissipative channels.
Phys. Rev. B, 101, 104307 (2020).
107. D. Rossini and E. Vicari.
Dynamic Kibble-Zurek scaling framework for open dissipative many-body systems crossing
quantum transitions.
Phys. Rev. Research, 2, 023211 (2020).
108. R. Puebla, A. Smirne, S. F. Huelga, et al.
Universal Anti-Kibble-Zurek Scaling in Fully Connected Systems.
Phys. Rev. Lett., 124, 230602 (2020).
109. L. Arceci, S. Barbarino, D. Rossini, et al.
Optimal working point in dissipative quantum annealing.
Phys. Rev. B, 98, 064307 (2018).
110. G. B. Mbeng, L. Arceci, and G. E. Santoro.
Optimal working point in digitized quantum annealing.
Phys. Rev. B, 100, 224201 (2019).
111. S. Suzuki, H. Oshiyama, and N. Shibata.
Quantum Annealing of Pure and Random Ising Chains Coupled to a Bosonic Environment.
 J. Phys. Soc. Jpn., 88, 061003 (2019). 35
112. H. Oshiyama, N. Shibata, and S. Suzuki.
Kibble–Zurek Mechanism in a Dissipative Transverse Ising Chain.

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
J. Phys. Soc. Jpn., 89, 104002 (2020).
113. H. Oshiyama, S. Suzuki, and N. Shibata.
Classical simulation and theory of quantum annealing in a thermal environment.
Phys. Rev. Lett., 128, 170502 (2022).
114. T. Jörg, F. Krzakala, J. Kurchan, et al.
Energy gaps in quantum first-order mean-field–like transitions: The problems that quantum
annealing cannot solve.
EPL (Europhysics Lett.), 89, 40004 (2010).
115. Y. Seki and H. Nishimori.
Quantum annealing with antiferromagnetic fluctuations.
Phys. Rev. E, 85, 051112 (2012).
116. Y. Susa, Y. Yamashiro, M. Yamamoto, et al.
Exponential speedup of quantum annealing by inhomogeneous driving of the transverse
field.
J. Phys. Soc. Jpn., 87, 023002 (2018).
117. Y. Susa, Y. Yamashiro, M. Yamamoto, et al.
Quantum annealing of the p-spin model under inhomogeneous transverse field driving.
Phys. Rev. A, 98, 042326 (2018).
118. H. Nishimori and G. Ortiz.
Elements of phase transitions and critical phenomena.
Oup Oxford (2010).
119. B. Seoane and H. Nishimori.
Many-body transverse interactions in the quantum annealing of the p-spin ferromagnet.
J. Phys. A: Math. Theor., 45, 435301 (2012).
120. M. Filippone, S. Dusuel, and J. Vidal.
Quantum phase transitions in fully connected spin models: An entanglement perspective.
Phys. Rev. A, 83, 022327 (2011).
121. S. Matsuura, H. Nishimori, T. Albash, et al.
Mean field analysis of quantum annealing correction.
Physical review letters, 116, 220501 (2016).
122. S. Matsuura, H. Nishimori, W. Vinci, et al.
Nested quantum annealing correction at finite temperature: p-spin models.
Physical Review A, 99, 062307 (2019).
123. S. Knysh, E. Plamadeala, and D. Venturelli.
Quantum annealing speedup of embedded problems via suppression of Griffiths
singularities.
Phys. Rev. B, 102, 220407 (2020).
124. R. B. Griffiths.
Nonanalytic behavior above the critical point in a random Ising ferromagnet.
Phys. Rev. Lett., 23, 17 (1969).
125. B. M. McCoy.
Incompleteness of the critical exponent description for ferromagnetic systems containing
random impurities.
Phys. Rev. Lett., 23, 383 (1969).
126. K. Jałowiecki, A. Wi˛eckowski, P. Gawron, et al.
Parallel in time dynamics with quantum annealers.
Sci. Rep., 10, 1 (2020).
127. N. A. Aadit, A. Grimaldi, M. Carpentieri, et al.
Massively parallel probabilistic computing with sparse Ising machines.
Nat. Electron., pages 1–9 (2022).
128. E. Pelofske, G. Hahn, and H. N. Djidjev.
Parallel quantum annealing.
Sci. Rep., 12, 1 (2022).
129. D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, et al.
Shortcuts to adiabaticity: Concepts, methods, and applications.
Rev. Mod. Phys., 91, 045001 (2019).
130. D. Sels and A. Polokovnikov. 36
Minimizing irreversible losses in quantum systems by local counterdiabatic driving.
Proc. Natl. Acad. Sci., 114, E3909 (2017).

[Link] Phil. Trans. R. Soc. A 0000000

..................................................................
131. A. D. King, J. Raymond, T. Lanting, et al.
Quantum critical dynamics in a 5000-qubit programmable spin glass (2022).
132. N. Mohseni, P. L. McMahon, and T. Byrnes.
Ising machines as hardware solvers of combinatorial optimization problems.
Nature Reviews Physics, 4, 363 (2022).
133. R. Harris, Y. Sato, A. J. Berkley, et al.
Phase transitions in a programmable quantum spin glass simulator.
Science, 361, 162 (2018).
134. P. Kairys, A. D. King, I. Ozfidan, et al.
Simulating the shastry-sutherland ising model using quantum annealing.
Phys. Rev. X Quantum, 1, 020320 (2020).
135. A. D. King, C. D. Batista, J. Raymond, et al.
Quantum annealing simulation of out-of-equilibrium magnetization in a spin-chain
compound.
Phys. Rev. X Quantum, 2, 030317 (2021).
136. A. Mott, J. Job, J.-R. Vlimant, et al.
Solving a Higgs optimization problem with quantum annealing for machine learning.
Nature, 550, 375 (2017).
137. S. Das, A. J. Wildridge, S. B. Vaidya, et al.
Track clustering with a quantum annealer for primary vertex reconstruction at hadron
colliders.
arXiv preprint arXiv:1903.08879 (2019).
138. S. Abel and M. Spannowsky.
Quantum-Field-Theoretic Simulation Platform for Observing the Fate of the False Vacuum.
Phys. Rev. X Quantum, 2, 010389 (2021).
139. A. Teplukhin, B. K. Kendrick, and D. Babikov.
Calculation of molecular vibrational spectra on a quantum annealer.
J. Chem. Theory Comput., 15, 4555 (2019).
140. A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, et al.
Finding low-energy conformations of lattice protein models by quantum annealing.
Sci. Rep., 2, 1 (2012).
141. R. Y. Li, R. Di Felice, R. Rohs, et al.
Quantum annealing versus classical machine learning applied to a simplified computational
biology problem.
NPJ quantum information, 4, 1 (2018).
142. M. H. S. Amin, P. J. Love, and C. J. S. Truncik.
Thermally Assisted Adiabatic Quantum Computation.
Phys. Rev. Lett., 100, 060503 (2008).
143. N. G. Dickson, M. W. Johnson, M. H. Amin, et al.
Thermally assisted quantum annealing of a 16-qubit problem.
Nat. Commun., 4, 1903 (2013).
144. A. Mishra, T. Albash, and D. A. Lidar.
Finite temperature quantum annealing solving exponentially small gap problem with non-
monotonic success probability.
Nat. Commun., 9, 2917 (2018).
145. V. Bapst and G. Semerjian.
On quantum mean-field models and their quantum annealing.
J. Stat. Mech.: Theory Exp., 2012, P06007 (2012).