Optimization using quantum mechanics: quantum annealing through adiabatic evolution

Santoro et al., doi:10.1088/0305-4470/39/36/r01 的机器可读全文页。

文献信息

Santoro et al., doi:10.1088/0305-4470/39/36/r01，“Optimization using quantum mechanics: quantum annealing through adiabatic evolution”. DOI:10.1088/0305-4470/39/36/r01 原文为 arXiv 预印本版本的机器可读转换，公式与图注以原文为准；本页仅作站内索引与全文查阅，引用请以正式出版物为准。

全文

TOPICAL REVIEW

Optimization using quantum mechanics: quantum annealing through adiabatic evolution

To cite this article: Giuseppe E Santoro and Erio Tosatti 2006 J. Phys. A: Math. Gen. 39 R393

View the article online for updates and enhancements.

You may also like

Optimization using quantum mechanics:- quantum annealing through adiabatic evolution Giuseppe E Santoro and Erio Tosatti

Diabatic quantum annealing for thefrustrated ring model Jeremy Côté, Frédéric Sauvage, Martín Larocca et al.

Mathematical aspects of quantumannealing H Nishimori and S Morita

TOPICAL REVIEW

Optimization using quantum mechanics: quantum annealing through adiabatic evolution

Giuseppe E Santoro^1,2 and Erio Tosatti^1,2

¹ International School for Advanced Studies (SISSA) and CNR-INFM Democritos National Simulation Center, Via Beirut 2-4, I-34014 Trieste, Italy

² International Center for Theoretical Physics (ICTP), Trieste, Italy

E-mail: santoro@sissa.it and tosatti@sissa.it

Received 22 May 2006, in final form 21 July 2006

Published 18 August 2006

Online at stacks.iop.org/JPhysA/39/R393

Abstract

We review here some recent work in the field of quantum annealing, alias adiabatic quantum computation. The idea of quantum annealing is to perform optimization by a quantum adiabatic evolution which tracks the ground state of a suitable time-dependent Hamiltonian, where $\cdot {\hslash }^{\prime }$ is slowly switched off. We illustrate several applications of quantum annealing strategies, starting from textbook toy-models—double-well potentials and other one-dimensiona examples, with and without disorder. These examples display in a clear way the crucial differences between classical and quantum annealing. We then discuss applications of quantum annealing to challenging hard optimization problems, such as the random Ising model, the travelling salesman problem and Boolean satisfiability problems. The techniques used to implement quantum annealing are either deterministic Schrodinger’s evolutions, for the toy models, or path-¨ integral Monte Carlo and Green’s function Monte Carlo approaches, for the hard optimization problems. The crucial role played by disorder and the associated non-trivial Landau–Zener tunnelling phenomena is discussed and emphasized.

PACS numbers: 03.67.Lx, 75.10.Nr, 03.65.Xp, 02.70.Ss, 05.10.Ln, 07.05.Tp

(Some figures in this article are in colour only in the electronic version)

1. Introduction

Optimization has to do with finding the minimum of some given cost function; it is a very common task in science and a self-standing field of research. Quite often, the cost (or energy) function one is seeking to optimize depends on a large number N of variables (either discrete or continuous), and because of multiple constraints—which may lead to frustration—it will possess many local minima, a feature that makes the search for the optimal solution generically ‘hard’. As non-trivial examples of hard energy functions, one can mention (i) finding the ground-state energy configuration of an assembly of N Ising spins with frustrating interactions (a spin glass) [1]; (ii) finding the minimal length itinerary of a travelling salesman problem (TSP) visiting N cities [2]; (iii) finding the satisfying bit assignment for an N-bit Boolean formula with many logical clauses each involving k bits (k-SAT) [3]. These examples, taken from the fields of physics (i) and computer science (ii, iii), are just a small repertoire from the vast world of disordered/frustrated physical systems, and of combinatorial optimization problems [3]; we will use them, in the present work, as our favourite illustration and benchmark problems. Common to all of them is the fact that the dimension of the configuration space of the N input variables increases exponentially with N—it is $2^N$ in the Ising and k-SAT cases, (N − 1)! in the TSP case—and that no obvious rule is available to find the optimal solution in such a space. On the other hand, the space is simply too large for an exhaustive enumeration to be feasible beyond a relatively small value of N $(N\sim 30)$

Theoretical computer science has classified the complexity of a given problem according to the time (and space) resources needed to solve it, in its worst-case instance, on a classica computer [4]. Without entering the details of such a classification (that will not be relevant in our work), we just mention that a problem belongs to the P or to the NP complexity classes if a polynomial (in N) algorithm is or is not known, respectively, which is able to solve worst-case instances. The examples quoted above all belong to the class of the NP-complete problems (the hardest of all the NP problems [3]), except for the Ising case on some two-dimensional lattices [5] or the 2-SAT problem, which are both polynomial.

Returning to physics, classical statistical mechanics has contributed a crucial input of concepts and methods to the field of optimization, starting from the pioneering idea of thermal simulated annealing by Kirkpatrick et al [6], down to the recent advances in the applications of replica methods and cavity techniques to general combinatorial optimization problems [7]. In simulated annealing, the optimization problem is tackled by the introduction of a (real or fictitious) temperature variable, which is slowly lowered (annealed) in the course of a Monte Carlo or molecular dynamics simulation [6]. This device allows exploration of the large configuration space of the problem at hand, effectively avoiding trapping at unfavourable local minima through thermal hopping above energy barriers. It makes for a very simple, robust and effective minimization tool [8, chapter 10], often much more effective than standard, gradient-based, minimization methods.

The idea of quantum annealing (QA) is an elegant and fascinating alternative to classical thermal simulated annealing (CA); it consists in helping the system escape the local minima using quantum mechanics—by tunnelling through the barriers rather than thermally overcoming them—with an artificial and appropriate source of quantum fluctuations (the counterpart of the temperature) initially present and slowly (adiabatically) switched off. To the best of our knowledge, this idea was first explicitly formulated, and tested in numerical simulations, in the early 1990s [9–11]. In the Ising spin glass context—more precisely, for the infinite range Sherrington–Kirkpatrick model [12]—the idea that the addition of a transverse field quantum term $-\Gamma {\sum }_i{\sigma }_i^x({\sigma }_i^x$ being Pauli matrices at site i) might help the
system in tunnelling through the infinitely high classical barriers separating the infinitely many metastable states was indeed put forward even earlier, in [13, 14]. More recently, experimental evidence in the disordered Ising ferromagnet ${\mathrm{LiHo}}_{0.44}{\mathrm{Y}}_{0.56}{\mathrm{F}}_4$ in a transverse magnetic field showed that the QA strategy is not only feasible but presumably winning in certain cases [15, 16].

A parallel advance, again during the 1990s and with quantum mechanics as a protagonist, has been the blooming idea of quantum computation (QC) [17]. A remarkable result in this field, due to Schor [18], was the proof that the problem of factorizing a large integer of N bits—believed to be classically NP (although not NP-complete)—can be solved with a number of quantum operations which are polynomial in N $(O(N^2))$ ), in the framework of the standard approach to QC (the so-called circuit theory) [17].

The remainder of this Introduction will be devoted to a brief list of basic concepts and issues.

1.1. The basic idea ofquantum annealing

Let us consider $\mathrm{QA},$ and concentrate discussion on classical optimization problems, where the energy function to be optimized is represented by a classical Hamiltonian, denoted by $H_{\mathrm{cl}}$ In QA, one supplements the classical Hamiltonian $H_{\mathrm{cl}}$ with a suitable, in principle arbitrary, time-dependent quantum kinetic term, $H_{\mathrm{kin}}(t)$ , which is initially very large, for $t⩽0$ , and then gradually reduced to zero in a certain annealing time τ. The total Hamiltonian one aims at studying is therefore a time-dependent one given by $H(t)=H_{\mathrm{cl}}+H_{\mathrm{kin}}(t)$ , where $H_{\mathrm{kin}}(0)$ dominates initially—that is, H(0) ≈ $H_{\mathrm{kin}}(0)\cdot$ —while, at the end of the annealing, $H_{\mathrm{kin}}(\tau )=0$ so that we simply recover $H(\tau )=H_{\mathrm{cl}}$ . The kinetic energy $H_{\mathrm{kin}}(t)$ can be chosen quite freely, with the only requirement of being represented by a local Hamiltonian. Using the Ising glass case as an illustrative example, we could take

$$H(t)=H_{\mathrm{cl}}+H_{\mathrm{kin}}(t)=-\sum _{⟨ij⟩}{J}_{ij}{\sigma }_i^z{\sigma }_j^z-\Gamma (t)\sum _i{\sigma }_i^x,\text{\hspace{1em}(1)}$$

where $\begin{array}{r}{H}_{\mathrm{cl}}=-{\sum }_{⟨ij⟩}{J}_{ij}{\sigma }_i^z{\sigma }_j^z\end{array}$ is an Edward–Anderson disordered Ising model—with ${\sigma }_i^z$ $\mathrm{spin}-1/2$
Ising variables (z-component Pauli matrices) at each site i of a lattice, and $J_{ij}$ random couplings—and represents the classical problem for which we are trying to find the ground state. Due to the fact that the nearest-neighbour spin–spin couplings $J_{ij}$ are frustrating $(i.e.$ both positive and negative) and random, finding the ground state of $H_{\mathrm{cl}}$ is generally hard. The choice of the kinetic term as a simple time-dependent transverse field $\begin{array}{r}{H}_{\mathrm{kin}}(t)=-\Gamma (t){\sum }_i{\sigma }_i^x\end{array}$ (with ${\mathbf{\sigma }}_i^x$
the x-component Pauli matrix at site i) is in this case very natural, as is also directly realized in the experiment [15]. Indeed this form of $H_{\mathrm{kin}}$ is quite natural for all those problems that can be classically formulated in terms of spin- $\cdot 1/2$ Ising variables (or bits, as in the case of k-SAT). The condition that $H_{\mathrm{kin}}(0)$ dominates initially is realized by $\Gamma (0)\gg \mathrm{Max}\{|J_{ij}|\};$ for $\Gamma (0)\infty$ , the ground-state $|{\Psi }_0⟩$  of $H(0)=H_{\mathrm{cl}}+H_{\mathrm{kin}}(0)$ approaches the simple product state $\begin{array}{r}{2}^{-N/2}{\prod }_i(|\uparrow ⟩+|\downarrow ⟩)\end{array}$ , with all spins aligned in the x-direction; at the end of the annealing $\Gamma (\tau )=0$ , and the lowest energy of $H(\tau )$ coincides with the classical ground state of $H_{\mathrm{cl}}$ we are looking for. At zero temperature, the quantum state of the system $\left|\Psi (t)\right|$ , initially prepared in the fully quantum ground state $|{\Psi }_0$ of $H(t=0),|\Psi (t=0)⟩=|{\Psi }_0⟩$ , evolves according to the time-dependent Scrodinger equation ¨

$$\mathrm{i}\hslash \frac{\mathrm{d}}{\mathrm{d}t}|\Psi (t)⟩=H(t)|\Psi (t)⟩,\text{\hspace{1em}(2)}$$

to reach a final state $|\Psi (t=\tau )⟩$ The key question is then how the residual energy $E_{\mathrm{res}}(\tau )=E_{\mathrm{fin}}(\tau )-E_{\mathrm{opt}}$ _t decreases for increasing annealing time τ , where $E_{\mathrm{opt}}$ is the absolute minimum of $H_{\mathrm{cl}}.$ , and $E_{\mathrm{fin}}(\tau )$ is the average energy attained by the quantum system after evolving for a time τ, i.e.,

$$E_{\mathrm{res}}(\tau )=E_{\mathrm{fin}}(\tau )-E_{\mathrm{opt}}=\frac{⟨\Psi (\tau )|(H_{\mathrm{cl}}-E_{\mathrm{opt}})|\Psi (\tau )⟩}{⟨\Psi (\tau )|\Psi (\tau )⟩}.\text{\hspace{1em}(3)}$$

Generally speaking, the question posed is obviously related to the adiabaticity of the quantum evolution, i.e., whether the system is able, for sufficiently slow annealing (i.e., sufficiently long τ), to follow the instantaneous ground state of $H(t)$ . Quite generally, also, the resul for $E_{\mathrm{res}}(\tau )$ will depend on the specific instance of the problem taken (given N, the specific choice of the couplings $\{J_{ij}\}$ , in the Ising case), so that we have a distribution $P[E_{\mathrm{res}}(\tau )]$ of residual energies, and we can meaningfully talk about average (av) and typical (typ) residual energies $[E_{\mathrm{res}}(\tau ){]}_{av/typ}$ (which might differ if the distribution is wide and non-Gaussian [19]). A crude, but useful, classification of the ‘annealing complexity’ of a problem might then be the following: A problem is on average/typically easy if, asymptotically for large τ ,

$${\left[\frac{E_{\mathrm{res}}(\tau )}{N}\right]}_{av/typ}\propto {\tau }^{-\alpha },\text{\hspace{1em}(4)}$$

where the power-law exponent α does not vanish for $N\to \infty ;$ if, in contrast, $[E_{\mathrm{res}}(\tau ){]}_{av/typ}$ decreases slower than a power law, for instance as

$${\left[\frac{E_{\text{ res }}(\tau )}{N}\right]}_{\text{ av / typ }}\propto (\log \tau {)}^{-\zeta },\text{\hspace{1em}(5)}$$

then we say that the problem is on average/typically hard. In essence, if a problem as a powerlaw decrease of $[E_{\mathrm{res}}(\tau ){]}_{av/typ}$ then we can increase the accuracy of our annealing calculation, reducing the residual energy by a given factor, say 2, by increasing the annealing time from τ to $2^{1/\alpha }\tau$ , while, for instance, a logarithmically slow decrease would require, for the same improvement, an annealing time ${\tau }^{2^{1/\zeta }}$

1.2. Quantum annealing and adiabatic quantum computation

As previously mentioned, the idea of annealing, in a quantum mechanics framework, is inherently related to the idea of adiabaticity. Computing by adiabatic evolution of a quantum system has become quite popular in the quantum computing (QC) community, where it is commonly known as adiabatic quantum computation (AQC) [20]. The standard reference for the QC community is [20], where the idea is formulated in the following alternative way: take $H_{\mathrm{kin}}$ to be time independent, and write $H(t)=(t/\tau )H_{\mathrm{cl}}+(1-t/\tau )H_{\mathrm{kin}}.$ so that $H(0)=H_{\mathrm{kin}}$ and, once again, $H(\tau )=H_{\mathrm{cl}}$ . We stress, therefore, that quantum annealing and adiabatic quantum computation are two names given by two different scientific communities to the very same idea; we will therefore use QA and AQC interchangeably, or even together, QA–AQC. We mention, finally, that a generalized form of AQC has been recently shown to be polynomially equivalent [21] to the standard paradigm of QC, i.e., circuit theory, where computations are performed by the application of a sequence of universal gates (unitary operators) to a system of qubits [17]. In essence, if a problem can be efficiently solved by standard QC, then one can construct $H_{\mathrm{init}}$ _t and $H_{\mathrm{final}}$ (both local, the latter in general non-classical) such that a QA–AQC with $H(t)=(1-t/\tau )H_{\mathrm{init}}+(t/\tau )H_{\mathrm{final}}$ leads to an equivalent result, at the price of an extra computational effort that scales polynomially with N. Therefore, the computational power of standard QC and that of such a generalized QA–AQC is the same [21]. The important point, however, is that $H_{\mathrm{final}}$ is not necessarily a classical Hamiltonian. Therefore, it is still not obvious whether the computational power of QA–AQC for classical optimization problems, where one insists in taking $H_{\mathrm{final}}=H_{\mathrm{cl}}$ , is the same as that of standard QC. Indeed, we wil suggest that even polynomial classical optimization problems can be hard for QA–AQC.

1.3. Adiabaticity and Landau–Zener tunnelling

Returning to the problem of adiabaticity, it is quite common that for a given instance of a problem of input size $N,$ the Hamiltonian $H(t)$ will have an instantaneous spectrum of eigenvalues $E_n(t)$ with a strictly positive minimum gap $\begin{array}{r}{\Delta }_{\min }=\mathbf{M}{\mathrm{in}}_{0⩽t<\tau }[E_1(t)-E_0(t)]>0\end{array}$ above the ground state. For instance, in the Ising case—more generally, in all cases in which the matrix elements of $H_{\mathrm{kin}}(t)$ are non-positive, and any two configurations can be connected by the application of a suitable number of $H_{\mathrm{kin}}.$ —due to the Perron–Frobenius theorem [22, p 54ff], one can immediately conclude that the ground state of (1) for any non-vanishing transverse field $\Gamma (t)>0$ has an everywhere-positive wavefunction and is unique, hence ${\Delta }_{\min }>0.$ . Ideally, therefore, the adiabatic theorem of quantum mechanics [23] guarantees that the system will adiabatically follow the ground state if the annealing time τ is larger than a certain characteristic time, ${\tau }_{\mathrm{LZ}}\propto {\Delta }_{\min }^{-2}$ , beyond which $E_{\mathrm{res}}(\tau )$ decreases very fast (see, for instance, [24], where a behaviour $E_{\mathrm{res}}(\tau )\propto {\tau }^{-2}$ is argued). (The earliest realization of a such a phenomenon is just the textbook $2\times 2$ Landau–Zener problem [25, 26], see section 3.) In order for this argument to be of any real use, however, one should make sure that ${\Delta }_{\min }$ does not get very small over the ensemble of instances one is considering, and, more generally, when $N\to \infty$ . In practice, unfortunately, due to disorder and the associated Landau–Zener tunnelling between broad barriers, such a minimum gap ${\Delta }_{\min }$ can be so tiny (even in systems with a relatively small value of N) that the annealing time needed to adiabatically follow the ground state is astronomically large. Later on, we will present a simple case, based on a disordered Anderson model, where this happens (see section 3). As a general remark, therefore, beware of QA–AQC results for cases where N is so small that, typically, the system will have a relatively large ${\Delta }_{\min }.$ : they are not really representative of the actual computational power of QA–AQC on that problem.

1.4. Real-time versus imaginary-time Schrodinger evolution ¨

In practical implementations of QA on an ordinary classical computer, the task of following directly the time-dependent Schrodinger evolution in (¨ 2) is feasible only for toy models with a sufficiently manageable Hilbert space [11, 27, 28]; for spin-1/2 problems, this restricts N to the range $N<20-30$ . Actual optimization problems of practical interest usually involve astronomically large Hilbert spaces, a fact that calls for alternative quantum Monte Carlo (QMC) approaches. These QMC techniques, in turn, are usually based, in some way or another, on an imaginary-time (IT) framework. For instance, it is much easier to implement, stochastically, an imaginary-time Schrodinger evolution,¨

$$-\hslash \frac{\mathrm{d}}{\mathrm{d}t}|\Psi (t)⟩=H(t)|\Psi (t)⟩,\text{\hspace{1em}(6)}$$

rather than the standard real-time (RT) evolution in (2). Obviously, such an evolution is non-unitary; hence the norm of |(t) is not conserved. It is still perfectly legitimate to ask if, at the end of the annealing, the residual energy $E_{\mathrm{res}}^{\mathrm{IT}}(\tau )$ , which is formally given by the same expression (3), is better or worse than the corresponding one obtained in a standard real-time evolution. We have recently addressed such a question in the context of simplified toy problems [28], finding that, as far as annealing is concerned, IT is essentially equivalent to RT, and, as a matter of fact, it can be quantitatively better [28]; we conjecture that, quite generally, $E_{\mathrm{res}}^{\mathrm{IT}}(\tau )<E_{\mathrm{res}}(\tau )$ for a given instance of the problem. This aspect will be directly addressed in section 2.3.

1.5. Quantum Monte Carlo studies

A number of recent studies have applied strategies related to path-integral Monte Carlo (PIMC) [29] or to diffusion quantum Monte Carlo (DMC) [30], to perform QA. PIMC is based on an imaginary-time framework, but now an equilibrium one, rather than a time-dependent one. Instead of solving (6), one simulates the equilibrium properties of the quantum partition function $Z={\mathrm{Tre}}^{-\beta \overline{H(}\Gamma )}$ for a given fixed value of the quantum coupling $\Gamma (t)$ appearing in

$H_{\mathrm{kin}}$ and sufficiently small temperature $T=1/\beta$ , reducing the value of $\Gamma (t)$ as a function of the Monte Carlo ‘time’, in much the same way as temperature is reduced in standard Metropolis Monte Carlo-simulated annealing. A certain success has been obtained in a number of optimization problems, such as the folding of off-lattice protein models [31, 32], the ground-state structure of classical Lennard–Jones clusters [9, 10, 33, 34], the random Ising model [35, 36] and the random field Ising model ground-state problems [37] (see section 4.2.1), the travelling salesman problem [38] (see section 4.2.2). On the other hand, for the interesting case of Boolean satisfiability (SAT)—more precisely, a prototypical NPcomplete problem such as 3-SAT—a recent study of our group shows that PIMC annealing performs worse than simple CA [39] (see section 4.2.3). The same conclusion might apply to the one-dimensional (1D) Ising ANNNI model case [40].

In view of these results, it is fair to stress that it is a priori not obvious or guaranteed that a QA approach should do better than, for instance, CA, on a given problem. Evidently, the comparative performance of QA and CA depends in detail on the energy landscape of the problem at hand, in particular on the nature and type of barriers separating the different local minima, a problem about which very little is known in many practical interesting cases [41]. That in turn depends crucially on the type and effectiveness of the kinetic energy chosen. Unfortunately, there is still no reliable theory predicting the performance of a QA algorithm, in particular correlating it with the energy landscape of the given optimization problem.

In order to gain understanding on these problems, we have moved, more recently, one step back and concentrated attention on the simplest textbook problems where the energy landscape is well under control: essentially, one-dimensional potentials, starting from a double-well potential, the simplest form of barrier. On these well-controlled landscapes we carried out a detailed and exhaustive comparison between quantum adiabatic Schrodinger¨ evolution, both in real and in imaginary time, and its classical deterministic counterpart, i.e., Fokker–Planck evolution [28]. This work will be illustrated in section 2.2. On the same double well potential, we also studied [42] the performance of different stochastic approaches, both classical Monte Carlo and path-integral Monte Carlo. Some of this work, which turns out to be quite instructive, is briefly presented in section 4.3.

The scope of this review is to give an overview of some of the recent work in the field of QA, from the authors’ perspective. No systematic attempt has been made in being exhaustive in quoting references in the field, and we apologize for omissions: the reader can find a more extensive source of references in the recent book on the subject by Das and Chakrabarti [43]. A review of the experimental work on the disordered Ising ferromagnet ${\mathrm{LiHo}}_{0.44}{\mathrm{Y}}_{0.56}{\mathrm{F}}_4$ is given in [16].

The rest of this review is organized as follows: section 2 illustrates the deterministic annealing approaches applied to toy problems, essentially the minimization of a function of a continuous coordinate. Section 3 discusses the crucial role played by disorder and the issue of Landau–Zener tunnelling in QA. Section 4 introduces the path-integral Monte Carlo techniques, and illustrates some of the recent applications, notably on the random Ising model, on the travelling salesman problem, and on Boolean satisfiability problems. Section 5 introduces the Green’s function Monte Carlo technique, and illustrates its recent application to the same random Ising model instance previously studied by PIMC–QA as well as CA. Section 6, finally, contains a brief summary of the main points, and some concluding remarks.

2. Deterministic approaches on toy problems

The deterministic annealing approaches (both classical and quantum) described in the present section are directly related to the work of Kadowaki and Nishimori on Ising models on small clusters $(N⩽8)$ [11]. Further references along the same lines include [27], where a Schrodinger evolution¨ $\mathrm{QA}-\mathrm{AQC}$ is applied to small $(N⩽20)$ random instances of a Boolean satisfiability problem (exact cover), and recent work on the disordered Anderson model (with $N⩽20$ sites) and on the $\pm J$ Ising model on a $3\times 3$ lattice [24].

Conceptually, one of the simplest optimization problems to illustrate is that of finding the global minimum of an ordinary function of several continuum variables with many minima. Suppose the classical Hamiltonian $H_{\mathrm{cl}}$ mentioned in the introduction is just a potential energy $V(\mathbf{x})$ , (with x being a Cartesian vector of arbitrary dimension), of which we need to determine the absolute minimum $({\mathbf{x}}_{\mathrm{opt}},E_{\mathrm{opt}}=V({\mathbf{x}}_{\mathrm{opt}}))$ . Assume, generally, a situation in which a steepest descent approach, i.e., the strategy of following the gradient of $V,$ would lead to trapping into one of the many local minima of $V,$ , and would thus not work. Classically, as an obvious generalization of a steepest descent approach, one could imagine performing a stochastic (Markov) dynamics in x-space according to a Langevin equation

$$\dot{\mathbf{x}}=-\frac{1}{\eta (T)}\vec{\nabla }V(\mathbf{x})+\xi (t),\text{\hspace{1em}(7)}$$

where the strength of the noise term ξ is controlled by the squared noise correlations $\overline{{\xi }_i(t){\xi }_j(t^{\prime })}=2D(T){\delta }_{ij}\delta (t-t^{\prime })$ , with $\bar{\xi }=0$ . Both $D(T)$ and $\eta (T)$ —with dimensions of a diffusion constant and of a friction coefficient and related, respectively, to fluctuations and dissipation in the system—are temperature-dependent quantities which can be chosen, for the present optimization purpose, with a certain freedom. The only obvious constraint is in fact that the correct thermodynamical averages will be recovered from the Langevin dynamics only if $\eta (T)D(T)=k_{B}T$ , an equality known as Einstein’s relation [44]. Physically, $D(T)$ should be an increasing function of $T,$ so as to lead to increasing random forces as $T$ increases, with $D(T=0)=0,$ , since noise is turned off at $T=0,$ . Classical annealing can in principle be performed through this Langevin dynamics, by slowly decreasing the temperature $T(t)$ as a function of time, from some initially large value $T_0$ down to zero. Instead of working with Langevin equation—a stochastic differential equation—one might equivalently address the problem by studying the probability density $P(\mathbf{x},t)$ of finding a particle at position x at time t. The probability density is well known to obey a deterministic time-evolution equation given by the Fokker–Planck (FP) equation [44]

$$\frac{\partial }{\partial t}P(\mathbf{x},t)=\frac{1}{\eta (T)}\mathrm{div}(P\vec{\nabla }V)+D(T){\nabla }^{2}P.\text{\hspace{1em}(8)}$$

Here, the second term on the right-hand side represents the well-known diffusion term, proportional to the diffusion coefficient $D(T)$ , whereas the first term represents the effect of the $drift$ force $-\vec{\nabla }V$ , inversely proportional to the friction coefficient $\eta (T)=k_{B}T/D(T)$ [44]. Annealing can now be performed by keeping the system for a long enough equilibration time at a large temperature $T_0,$ and then gradually decreasing $T$ to zero as a function of time, $T(t)$ , in a given annealing time τ . We can model this by assuming $T(t)=T_{0}f(t/\tau )$ , where $f(y)$ is some assigned monotonically decreasing function for $y\in [0,1]$ , with $f(y⩽0)=1$ and $f(1)=0$ . In this manner the diffusion constant $D$ in equation (8) becomes a timedependent quantity, $D_t=D(T(t))$ . The FP equation should then be solved with an initial condition given by the equilibrium Boltzmann distribution at temperature $T(t=0)=T_0$ , i.e., $P(\mathbf{x},t=\bar{0})={\mathrm{e}}^{-\dot{V}(\mathbf{x})/k_B{\bar{T}}_0}$ . The final average potential energy after annealing, in excess of the true minimum value, will then be simply given by

$${ϵ}_{\text{ res }}(\tau )=\int \mathrm{d}\mathbf{x}V(\mathbf{x})P(\mathbf{x},t=\tau )-E_{\text{ opt }}⩾0,\text{\hspace{1em}(9)}$$

where $E_{\mathrm{opt}}$ is the actual absolute minimum of the potential $V.$

In a completely analogous manner, we can conceive using Schrodinger’s equation to ¨ perform a deterministic quantum annealing $(\mathbf{QA})$ evolution of the system, by introducing quantum fluctuations through a standard non-relativistic kinetic term $H_{\mathrm{kin}}(t)=-({\hslash }^2/2m_t){\nabla }^2$ with a fictitious time-dependent mass $m_t$ . We are therefore led to studying the time-dependent Schrodinger problem¨

$$\xi \frac{\partial }{\partial t}\psi (\mathbf{x},t)=[-\Gamma (t){\nabla }^2+V(\mathbf{x})]\psi (\mathbf{x},t),\text{\hspace{1em}(10)}$$

where $\xi =\mathrm{i}\hslash$ for a real-time (RT) evolution, while $\xi =-\hslash$ for an imaginary-time (IT) evolution. Here $\Gamma (t)={\hslash }^2/2m_t$ will be our annealing parameter, playing the role that the temperature $T(t)$ had in classical annealing. Once again we may take $\Gamma (t)$ varying from some large value ${\Gamma }_0$ at $t⩽0$ —corresponding to a small mass of the particle, hence to large quantum fluctuations—down to $\Gamma (t=\tau )=0.$ , corresponding to a particle of infinite mass, hence without quantum fluctuations. Again, we can model this with $\Gamma (t)={\Gamma }_{0}f(t/\tau )$ , where $f$ is a preassigned monotonically decreasing function. The initial condition here will be $\psi (\mathbf{x},t=0)={\psi }_0(\mathbf{x})$ , where ${\psi }_0(\mathbf{x})$ is the ground state of the system at $t⩽0.$ , corresponding to the large value $\Gamma (t)={\Gamma }_0$ and hence to large quantum fluctuations. The residual energy after annealing will be similarly given by equation (9), where now, however, the probability $P(\mathbf{x},t=\tau )$ should be interpreted, quantum mechanically, as

$$P(\mathbf{x},t)=\frac{|\psi (\mathbf{x},t){|}^2}{\int \mathrm{d}{\mathbf{x}}^{\prime }|\psi ({\mathbf{x}}^{\prime },t){|}^2}.$$

In general, the residual energy will be different for RT or IT Schrodinger evolutions. We will ¨ comment further on RT versus IT Schrodinger evolution later on.¨

In the remaining part of this section, we will present some of the results obtained along the previous lines on simple one-dimensional potentials [28].

2.1. The harmonic potential: a warm-up exercise

Preliminary to any further treatment of a potential with barriers, and as a warm-up exercise which will be useful later on, we start here with the simple case of a parabolic potential in one dimension, $V(x)=k{x}^2/2$ , which has a trivial minimum in $x=0,$ with $E_{\mathrm{opt}}=0$ , and no barriers whatsoever.

Let us consider classical FP annealing first. As detailed in [28], it is a matter of simple algebra to show that, for the harmonic potential, one can write a simple closed linear differential equation [45] for the average potential energy ${ϵ}_{\mathrm{pot}}(t)$ , which has the form

$$\frac{\mathrm{d}}{\mathrm{d}t}{ϵ}_{\text{ pot }}(t)=k{D}_t\left[1-\frac{2}{k_{B}T(t)}{ϵ}_{\text{ pot }}(t)\right],\text{\hspace{1em}(11)}$$

the initial condition being simply given by the equipartition value ${ϵ}_{\mathrm{pot}}(t=0)=k_B{T}_0/2$ As with every one-dimensional linear differential equation, equation (11) can be solved by quadrature for any choice of $T(t)$ and $D_t=D(T(t))$ ). Assuming the annealing schedule to be parameterized by an exponent ${\alpha }_T>0,T(t)=T_0(1-t/\tau {)}^{{\alpha }_T}$ , τ being the annealing time, and the diffusion constant $D(T)$ to behave as a power law of temperature, $D(T)=D_0(T/T_0{)}^{{\alpha }_D}$ with ${\alpha }_D⩾0$ , we can easily extract from the analytical solution for ${ϵ}_{\mathrm{pot}}(t)$ the large-τ asymptotic behaviour of the final residual energy ${ϵ}_{\mathrm{res}}(\tau )={ϵ}_{\mathrm{pot}}(t=\tau )$ . That turns out to be

$${ϵ}_{\text{ res }}(\tau )\approx {\tau }^{-{\Omega }_{\text{ CA }}}\text{ with }{\Omega }_{\text{ CA }}=\frac{{\alpha }_T}{{\alpha }_T({\alpha }_D-1)+1}.\text{\hspace{1em}(12)}$$

Trivial as it is, annealing proceeds here extremely fast, with a power-law exponent ${\Omega }_{\mathrm{CA}}$ that can increase without bounds (for instance if ${\alpha }_D=1)$ upon increasing the exponent ${\alpha }_T$ of the annealing schedule.

Consider now the Schrodinger evolution problem for this potential, ¨

$$\xi \frac{\partial }{\partial t}\psi (x,t)=\left[-\Gamma (t){\nabla }^2+\frac{k}{2}{x}^2\right]\psi (x,t)\psi (x,t=0)={\psi }_0(x),\text{\hspace{1em}(13)}$$

where ${\psi }_0(x)\propto \exp (-B_0{x}^2/2)$ is the ground-state Gaussian wavefunction corresponding to the initial value of the Laplacian coefficient $\begin{array}{r}\vec{\nabla }(t=0)={\Gamma }_0,\end{array}$ and $\xi :=:\mathrm{i}\hslash$ or $\xi =-\hslash$ for RT or IT evolution, respectively. As shown in [28], a Gaussian ansatz for $\psi (x,t)$ , of the form $\psi (x,t)\propto \exp (-B_t{x}^2/2)$ with Real $\left(B_t\right)>0$ , satisfies the time-dependent Schrodinger¨ equation as long as the inverse variance $B_t$ of the Gaussian satisfies the following ordinary nonlinear first-order differential equation:

$$-\xi {\dot{B}}_t=k-2\Gamma (t)B_t^2{B}_{t=0}=B_0=\sqrt{\frac{k}{2{\Gamma }_0}}.\text{\hspace{1em}(14)}$$

Contrary to the classical case, there is no simple way of recasting the annealing problem in terms of a closed linear differential equation for the average potential energy ${ϵ}_{\mathrm{pot}}(t)$ . The final residual energy ${ϵ}_{\mathrm{res}}(\tau )={ϵ}_{\mathrm{pot}}(t=\tau )$ ) is still expressed in terms of $B_{\tau }$ (or better, of its real part 3 $\mathrm{i}(B_{\tau }))$ ,

$${ϵ}_{\mathrm{res}}(\tau )=\frac{\int \mathrm{d}xV(x)|\psi (x,t=\tau ){|}^2}{\int \mathrm{d}x|\psi (x,t=\tau ){|}^2}=\frac{k}{4\Re (B_{\tau })},\text{\hspace{1em}(15)}$$

but the behaviour of $B_t$ must be extracted from the study of the nonlinear equation (14). The properties of the solutions of equation (14) are studied in detail in [28], where it is shown that

(i) ${ϵ}_{\mathrm{res}}(\tau )$ cannot decrease faster than $1/\tau$ , for large τ, i.e., a power-law exponent ${ϵ}_{\mathrm{res}}(\tau )$ ≈ ${\tau }^{-{\Omega }_{\mathrm{QA}}}$ is bounded by ${\Omega }_{\mathrm{QA}}⩽1$

(ii) Adopting a power-law annealing schedule $\Gamma (t)={\Gamma }_0(1-t/\tau {)}^{{\alpha }_{\Gamma }}$ , the exponent ${\Omega }_{\mathrm{QA}}$ for the IT case is

$${\Omega }_{\mathrm{QA}}=\frac{{\alpha }_{\Gamma }}{{\alpha }_{\Gamma }+2},\text{\hspace{1em}(16)}$$

increasing towards the upper bound 1 as ${\alpha }_{\Gamma }$ is increased towards ∞. In particular, a linear annealing schedule, ${\alpha }_{\Gamma }=1$ , leads to ${\Omega }_{\mathrm{QA}}=1/3$

(iii) RT quantum annealing proceeds with exactly the same exponent ${\Omega }_{\mathrm{QA}}$ as IT quantum annealing—although ${ϵ}_{\mathrm{res}}^{\mathrm{RT}}(\tau )\gtrsim {ϵ}_{\mathrm{res}}^{\mathrm{IT}}(\tau )$ in general—except that the limit ${\alpha }_{\Gamma }=\infty$ (abrupt switch-off of the Laplacian coefficient) is singular in the RT case.

Note that, in the present continuum context, the residual energy cannot decrease faster $1/\tau$ contrary to the $1/\overline{{\tau }^2}$ decrease observed in small-sized discrete problems (disordered Anderson and random Ising models) [24].

Summarizing, we have learned that, for a single parabolic valley in configuration space, both CA and QA proceed with power laws, but CA can be much more efficient than QA, with an arbitrarily larger power-law exponent. We underline however that this is merely an academic matter at this point, steepest descent being much more efficient than both CA and QA in such a simple case. The power of annealing, in particular of QA, shows up only when potentials with barriers are considered.

2.2. The simplest barrier: a double-well potential

Consider, as a potential $V(\mathbf{x})$ to be optimized, a slightly generalized double-well potential in one dimension

$$V_{\text{asym}}(x)=\{\begin{array}{ll}{V}_0\frac{(x^2-a_{+}^2{)}^2}{a_{+}^4}+\delta x& \text{for}x⩾0\\ V_0\frac{(x^2-a_{-}^2{)}^2}{a_{-}^4}+\delta x& \text{for}x<0,\end{array}\text{\hspace{1em}(17)}$$

with, in general, $a_{+}\ne a_{-}$ , both positive, $V_0.$ , and $\delta$ real constants. (The discontinuity in the second derivative at the origin is of no consequence in our discussion.) In the absence of the linear term $(\delta =0)$ , the potential has two degenerate minima located at $x_{-}=-a_{-}$ and $x_{+}=a_{+}$ , separated by a barrier of height $V_0$ . When a small linear term δ > 0 is introduced, with $\delta a_{\pm }\ll V_0$ , the two degenerate minima are split by a quantity ${\Delta }_V\approx \delta (a_{+}+a_{-})$ , the minimum at $x\approx -a_{-}$ becoming slightly favoured. For reasons that will be clear in a moment, it is useful to consider the situation, which we will refer to as ‘asymmetric double well’, in which the two wells possess definitely distinct curvatures at the minimum $(i.e.$ , their widths differ too), realized by taking $a_{+}\ne a_{-}$ . (To lowest order in $\delta ,$ we have $V^{\prime \prime }(x=x_{\pm })=8V_0/a_{+}^2.)$ In particular, we shall examine the case in which the metastable ‘valley’ at $x_{+}$ is ‘wider’ than the absolute minimum at $x_{-}$ , a situation realized by choosing $a_{+}>a_{-}$ . This curvature asymmetry will have a rather important effect on the quantum evolution, since, as we shall see, for small and intermediate values of the mass of the particle, the wavefunction of the system will be predominantly located on the metastable minimum. Obviously, if we set $a_{+}=a_{-}=a$ , and $\delta =0$ we recover the standard double-well potential. QA work on the latter problem, via an approximate solution of the imaginary-time Schrodinger problem with an Gaussian ansatz for ¨ the wavefunction, is discussed in [9].

We now present the results obtained by the annealing schemes introduced in section (2). The Fokker–Planck and the Schrodinger equation (both in RT and in IT) can be integrated¨ numerically using a fourth-order adaptive Runge–Kutta method [8], after discretizing the x variable in a sufficiently fine real-space grid [28]. For the FP classical annealing, the results shown are obtained with a linear temperature schedule, $T(t)=T_0(1-t/\tau )$ , and a diffusion coefficient simply proportional to $T(t),D_t=D_0(1-t/\tau )$ . (Consequently, the friction coefficient is kept constant in $t,{\eta }_t=k_{B}T(t)/D_t=k_B{T}_0/D_0.)$ Similarly, for the Schrodinger¨ quantum annealing we show results obtained with a coefficient of the Laplacian $\Gamma (t)$ vanishing linearly in a time $\tau ,\Gamma (t)={\Gamma }_0(1-t/\tau )$

Figure 1 shows the results obtained for the final annealed probability distribution $P(x,t=\tau )$ at different values of τ, for both the Fokker–Planck (FP-CA, panel (a)) and the Schrodinger imaginary-time case (IT-QA, panel (b)), for an ‘asymmetric’ double-well¨ potential $V_{\mathrm{asym}}(x)$ , with $V_0=1$ (our unit of energy), $a_{+}=1.25,a_{-}=0.75,\delta =0.1$ Figure 1(c) summarizes the results obtained for the residual energy ${ϵ}_{\mathrm{res}}(\tau )$ in equation (9).

We note immediately that QA wins over FP-CA, giving a better annealing for large enough value of τ . The RT-QA, which behaves as its IT counterpart for a symmetric double-well $(a_{+}=a_{-},\mathsf{see}$ [28]), shows a slightly different behaviour from IT-QA in the asymmetric case (see below for comments). We discuss first the FP-CA data (panel (a) and (c) of figure 1). Starting from an initially broad Boltzmann distribution at a high $T=T_0=V_0,P(x,t=0)$ (solid line), the system quickly sharpens the distribution $P(x,t)$ into two well-defined and quite narrow peaks located around the two minima $x_{\pm }$ of the potential. This agrees very well with what CA does for a harmonic potential [28]. If we denote by $p_{\pm }$ the integral of each of the two narrow peaks, with $p_{-}+p_{+}=1$ , it is clear that the problem has effectively been reduced

10

10

Figure 1. (a, b): the annealed final probability distribution $P(x,t=\tau )$ at different values of the annealing time τ , for both the Fokker–Planck classical annealing $(\mathrm{FP}-\mathrm{CA},$ panel (a)), and the imaginary-time Schrodinger quantum annealing (IT-QA, panel ( ¨ b)). (c): final residual energy ${ϵ}_{\mathrm{res}}(\tau )$ versus annealing time τ for quantum annealing in real-time (RT-QA) and imaginary-time $(\mathrm{IT}-\mathrm{QA})$ , compared to the Fokker–Planck classical annealing (FP-CA). The solid line in (c) is a fit of the FP-CA data. Here the double-well potential (dashed line in (a, b), inset of (c)) is given by equation (17) with $a_{+}=1.25,a_{-}=0.75,V_0=1,\delta =0.1$

to a discrete two-level problem. The time evolution of $p_{\pm }$ _±, therefore, obeys a discrete master equation which involves the thermal promotion of particles over the barrier, $V_0,$ of the form presented and discussed by Huse and Fisher in [46]. They show that, apart from logarithmic corrections, the leading behaviour of the residual energy is of the form ${ϵ}_{\mathrm{res}}\sim {\tau }^{-{\Delta }_V/B}$ , with the power-law exponent controlled by the ratio ${\Delta }_V/B$ between the energy splitting of the two minima ${\Delta }_V$ and the barrier $B=V_0-V(x_{+})$ . As shown in figure 1(c) (solid line through solid circle), the power-law asymptotic behaviour anticipated by Huse and Fisher fits nicely our $\mathrm{FP}-\mathrm{CA}$ residual energy data (solid circles), as long as the logarithmic corrections are accounted for in the fitting procedure [28]. Obviously, we can make the exponent as small as we wish by reducing the linear term coefficient δ, and hence the ratio ${\Delta }_V/B$ , leading to an exceedingly slow classical annealing.

The behaviour of the QA evolution is remarkably different. Observe, as a first point, that the final annealed wavefunctions narrow only slowly around the minimum of the potential, although the residual energy asymptotics of QA are clearly winning over FP-CA. The asymptotic behaviour of the QA residual energy is ${ϵ}_{\mathrm{res}}(\tau )\propto {\tau }^{-1/3}$ , indicated by the dashed line in figure 1(c). This rather strange exponent was anticipated by equation (16) for the Schrodinger annealing with a linear schedule¨
(t) within a harmonic potential (the lower minimum valley). Going back to figure 1(b), the initial wavefunction squared $|\psi (x,t=0){|}^2$ corresponds to a quite small mass (a large ${\Gamma }_0=0.5)$ , and is broad and delocalized over both minima (solid line). As we start annealing, and if the annealing time τ is relatively short— that is, if $\tau <{\tau }_c$ , with a characteristic time ${\tau }_c$ which depends on which kind of annealing,

Figure 2. (a, b): instantaneous eigenvalues (a) and ground-state wavefunctions (b) of the Schrodinger problem¨ $H\psi =E\psi$ for different values of
, for the potential in equation (17) with $a_{+}=1.25,a_{-}=0.75,V+0=1,\delta =0.1$ 1. Note the clear Landau–Zener avoided crossing in (a), indicated by the arrow and magnified in the inset. (c)–(e): instantaneous Fokker–Planck eigenvalues (c) as a function of temperature T, and the corresponding eigenstates for $T=1\left(d\right)$ and $T=0.1(e)$ , for the same potential used in (a, b).

RT or IT, we perform—the final wavefunction becomes mostly concentrated on the wrong minimum, roughly corresponding to the ground state with a still relatively large ${\Gamma }_1<{\Gamma }_0$ (see also figure 2 and accompanying discussion). The larger width of the wrong valley is crucial in causing that, giving a smaller quantum kinetic energy contribution, so that tunnelling to the other (deeper) minimum does not immediately occur, and the system is trapped in the wrong but broader valley. By increasing τ, there is a crossover: the system finally recognizes the presence of the other minimum, and effectively tunnels into it, with a residual energy that, as previously mentioned, decays asymptotically as ${ϵ}_{\mathrm{res}}(\tau )\propto {\tau }^{-1/3}$ (dashed line in figure 1(c)). There is a characteristic annealing time τ_c—different in the two Scrodinger cases, RT and¨ IT—above which tunnelling occurs, and this shows up as the clear crossover in the residual energy behaviour of both IT and RT, shown in figure 1(c).

These findings can be rationalized by looking at the instantaneous (adiabatic) eigenvalues and eigenstates of the associated time-independent Schrodinger problem, which we show in ¨ figures 2(a) and (b). Looking at the instantaneous eigenvalues shown in figure 2(a) we note a clear avoided crossing occurring at $\Gamma ={\Gamma }_{\mathrm{LZ}}\approx 0.038$ , corresponding to a resonance condition between the states in the two different valleys of the potential. For $\Gamma >{\Gamma }_{\mathrm{L}Z}$ the ground-state wavefunction is predominantly concentrated in the wider but metastable valley, while for $\Gamma <{\Gamma }_{\mathrm{LZ}}$ it is mostly concentrated on the deeper and narrower global minimum valley. In the full time-dependent RT evolution, transfer to the lower valley is a Landau–Zener problem [25, 26]: the characteristic time ${\tau }_c$ for the tunnelling event is given by ${\tau }_{\mathrm{LZ}}=\hslash \alpha {\Gamma }_0/2\pi {\Delta }^2$ where α is the relative slope of the two crossing branches as a function of $\Gamma ,2\Delta$ is the gap at the avoided crossing point and ${\Gamma }_0$ is the initial value of the annealing parameter. (For the case shown in figure 2, we have $2\Delta =0.0062,\alpha =2.3$ , hence ${\tau }_{\mathrm{LZ}}\approx 18980.$ , see rightmost arrow in figure 1(c).) The Landau–Zener probability of jumping, during the evolution, from the ground state onto the ‘wrong’ (excited) state upon fast approaching of the avoided level crossing is $P_{\mathrm{ex}}={\mathrm{e}}^{-\tau /{\tau }_{\mathrm{LZ}}}$ , so that adiabaticity applies only if the annealing is slow enough, $\tau >{\tau }_{\mathrm{L}Z}$ Note that the gap 2, and hence the probability of following adiabatically the ground state, can be made arbitrarily small by increasing the asymmetry of the two well, $i.e.$ by making $a_{+}\gg a_{-}$ The IT characteristic time is smaller, in the present case, than the RT one. This point is discussed in some detail in [28]. In a nutshell, the reason for this is the following. After the system has jumped into the excited state, which occurs with a probability $P_{\mathrm{ex}}={\mathrm{e}}^{-\tau /{\tau }_{\mathrm{LZ}}}$ , the residual IT evolution will filter out the excited state; this relaxation towards the ground state is controlled by the annealing rate as well as by the average gap seen during the residual evolution. Numerically, the characteristic time ${\tau }_c$ seen during the IT evolution is of the order of $\hslash /(2\Delta )$ , see leftmost arrow in figure 1(c), rather than being proportional to $1/{\Delta }^2$ as ${\tau }_{\mathrm{L}Z}$ would imply.

Obviously, instantaneous eigenvalues/eigenvectors can be studied for the Fokker–Planck equation as well; their properties, however, are remarkably different from the Landau–Zener scenario just described for the Schrodinger case. Figure ¨ 2(c) shows the first four low-lying eigenvalues of the FP equation as a function of $T,$ while figures $2(d)$ and (e) show the corresponding eigenstates for two values of the temperature, $T/V_0=1$ and $T/V_0=0.1$ for the same asymmetric potential used in figures (1) and (2). The lowest eigenvalue of the FP operator is identically $0,$ and the corresponding eigenvector [44] is the Boltzmann distribution ${\mathbf{e}}^{-V(x)/k_{B}T}$ , with roughly symmetric maxima on the two valleys. The first excited state corresponds to a function peaked on the two valleys but with a node at the origin, and is separated from the ground state by an exponentially small Arrhenius-like ${\mathrm{gape}}^{-B/k_B\overline{T}}$ . Higher excited states are separated by a very large gap, so that, effectively, only the two lowest lying states dominate the dynamics at small temperature. The reduction of a continuum double-wel FP classical dynamics onto a discrete effective two-level system, previously noted, is quite evident from this form of the spectrum. In contrast, the true quantum case does not allow a discrete two-level system description to hold for small enough
. Indeed, when $\Gamma <{\Gamma }_{\mathrm{LZ}}$ the tower of oscillator states within the valley at x₋ is always very close and collapsing in energy onto the actual ground state, and the quantum annealing evolution reduces effectively to a particle in a single harmonic well. This explains the rather large width of the final distributions $P(x,\tau )$ observed in the quantum case.

Summarizing, we have found that QA and CA proceed in a remarkably different way. CA is sensitive to the height of the barrier, more precisely to the ratio ${\Delta }_V/B$ between the energy offset ${\Delta }_V$ of the two minima, and the barrier height B. In contrast, QA crucially depends on the tunnelling probability between the two valleys, which is reflected in a Landau–Zener (avoided crossing) gap: a wide tunnelling barrier is obviously bad for QA. Finally, we noted that RT and IT proceed with different characteristic times: we discuss this issue a bit more in the following section.

2.3. Real- versus imaginary-time Schrodinger evolution¨

A Schrodinger dynamics in imaginary time (IT) is clearly much more convenient than that¨ in real time (RT) for QA simulations on current classical computers: but does it make a difference in the final results? The answer to this question is, we believe: no, it does not make a difference, in essence, although IT does give quantitatively better results.

To qualify this statement, let us denote by $|{\Psi }^{(\xi )}(t)⟩$  the solution of the Scrodinger equation¨

$$\xi \frac{\mathrm{d}}{\mathrm{d}t}|{\Psi }^{(\xi )}(t)⟩=[H_{\mathrm{cl}}+H_{\mathrm{kin}}(t)]|{\Psi }^{(\xi )}(t)⟩|{\Psi }^{(\xi )}(t_0)⟩=|{\Psi }_0⟩,\text{\hspace{1em}(18)}$$

where we assume that $|{\Psi }_0$ is the ground state of the initial Hamiltonian at time $t_0,H_{\mathrm{cl}}+$ $H_{\mathrm{kin}}(t_0)$ , while $\xi =\mathrm{i}\hslash$ for RT or $\xi =-\hslash$ for IT. By definition, the final residual energy after annealing up to time $t_f=t_0+\tau$ , where the kinetic energy is finally turned off, is given by

$${ϵ}_{\mathrm{res}}^{(\xi )}(\tau )=\frac{⟨{\Psi }^{(\xi )}(t_0+\tau )|H_{\mathrm{cl}}|{\Psi }^{(\xi )}(t_0+\tau )⟩}{⟨{\Psi }^{(\xi )}(t_0+\tau )|{\Psi }^{(\xi )}(t_0+\tau )⟩}-E_{\mathrm{opt}}.\text{\hspace{1em}(19)}$$

We conjecture that the residual energies for the two alternative ways of doing a Schrodinger¨ evolution verify the following: (i) the IT residual energy is not larger than the RT one, that is

$${ϵ}_{\text{ res }}^{\text{(IT)}}(\tau )⩽{ϵ}_{\text{ res }}^{\text{(RT)}}(\tau ),\text{\hspace{1em}(20)}$$

and (ii) in many problems, the leading asymptotic behaviour, for $\tau \to \infty$ , might be identica for ${ϵ}_{\mathrm{res}}^{(\mathrm{IT})}(\tau )$ and ${{\dot{ϵ}}_{\mathrm{res}}}^{(\mathrm{RT})}(\tau )$

Expectation (i) seems very reasonable, and is inspired by the time-independent case, where it is well known that the IT Schrodinger dynamics tends to ‘filter the ground state’¨ out of the initial trial wavefunction, as long as the gap between the ground state and the first excited state is non-zero. However, we have here a time-dependent situation, and the result is a priori not guaranteed. We do not have a proof of this statement, but we have verified it in all the cases where an explicit integration of the Schrodinger equation was possible (see, for¨ instance, the results of the previous section). Needless to say, we have no proof of (ii) either, but, again, it never failed in all our tests.

The simplest time-dependent problem where one can test our conjectures is the discrete two-level system (TLS) problem. Here, in terms of Pauli matrices, $H_{\mathrm{cl}}=\Delta {\sigma }^z$ , while $H_{\mathrm{kin}}(t)=-\Gamma (t){\sigma }^x$ , with $\Gamma (t)=-vt$ . The full $H(t)$ is therefore

$$H(t)=\Delta {\sigma }^z-\Gamma (t){\sigma }^x.\text{\hspace{1em}(21)}$$

The annealing interpretation is very simple: the classical optimal state is |↓, with energy $E_{\mathrm{opt}}=-\Delta$ , separated from the excited state |↑ by a gap 2. The kinetic term induces transitions between the two classical states. Starting from the ground state of $H(t_0)$ at time $t_0=-\tau$ we let the system evolve up to time $t_f=t_0+\tau =0$ , at which point the Hamiltonian is entirely classical, $H(t_f=0)=H_{\mathrm{cl}}=\Delta {\sigma }^z$ The probability of missing the instantaneous final ground-state |↓, ending up instead with the excited-state |↑, is $P_{\mathrm{ex}}(0)=|⟨\uparrow |{\Psi }^{(\xi )}(0)⟩{|}^2/⟨{\Psi }^{(\xi )}(0)|{\Psi }^{(\xi )}(0)⟩$ . In principle, $P_{\mathrm{ex}}$ depends, for given , both on the initial $\Gamma (t_0)={\Gamma }_0=\upsilon \tau$ and on the annealing time τ . The really important parameter, however, turns out to be the ratio v between these two quantities, which determines the ‘velocity of annealing’: taking $\tau \infty (i.e.,t_0-\infty )$ , and ${\Gamma }_0\infty$ with $\Gamma (t)=-vt$ for every t, the problem can be solved analytically (in terms of parabolic cylinder functions, see for instance [11] for the RT case) for both RT and IT. The probability $P_{\mathrm{ex}}(0)$ of ending into the excited state can be expressed in terms of the variable $\gamma ={\Delta }^2/4v$ . The explicit expressions, in terms of Gamma functions, are

$$P_{\mathrm{ex}}(0)=\frac{|R+1{|}^2}{2(1+|R{|}^2)}\text{ with }R={\mathrm{e}}^{\mathrm{i}{\varphi }_0}\frac{1}{\sqrt{\gamma }}\frac{\Gamma (1+z)}{\Gamma (1/2+z)},\text{\hspace{1em}(22)}$$

where ${\varphi }_0=3\pi /4$ and $z=\mathrm{i}\gamma$ for RT, while ${\varphi }_0=\pi$ and $z=\gamma$ for IT. A plot of $P_{\mathrm{ex}}$ for both RT and IT is shown in figure $3(a)$ as a function of $\gamma ={\Delta }^2/4v$ . Note that (i) the IT result for $P_{\mathrm{ex}}$ (dashed line) is always below the RT result, (ii) the difference between the two curves is only quantitative: one can verify analytically that the leading behaviour for large $\gamma$ is the same in both cases, i.e., $P_{\mathrm{ex}}\approx 1/(256{\gamma }^2)$ . Similar results are obtained by direct numerical integration of the Schrodinger equation for finite¨ ${\Gamma }_0$ and $\tau ,$ and with other forms of
(t).

Figure 3. (a) The probability $P_{\mathrm{ex}}(0)$ of ending up into the excited state, given by equation $(22)$ for the discrete two-level system problem in equation $(21).$ , for both imaginary-time (IT, dashed line) and real-time (RT, solid line) Schrodinger annealing. The large¨ γ behaviour of $P_{\mathrm{ex}}(0)$ is, in both cases, given by $P_{\mathrm{ex}}(0)\approx 1/(256{\gamma }^2)$ . (b) Comparison between the RT (solid lines) and the IT (dashed lines) evolution of a Landau–Zener problem, equation (23), for several values of the tunnelling gap $2\Delta$ (the values of  shown are $\dot{\Delta }=0.4,0.\dot{2},0.1,10^{-2},10^{-3},10^{-4},10^{-5},10^{-6}$ while $v=1)$ . The inset shows the two instantaneous eigenvalues of the problem, $E_{\pm }(t)$ , as a function of t.

With the same toy model, we can illustrate another point raised in the previous section: what happens to the IT evolution after a Landau–Zener avoided crossing gap is encountered? The Hamiltonian we consider is essentially that of equation (21), simply rotated in spin space,

$$H(t)=-vt{\sigma }^z-\Delta {\sigma }^x.\text{\hspace{1em}(23)}$$

In the absence of the tunnelling amplitude $\Delta$ , the two energy levels would cross at $t=0$ , while for $\Delta >0$ the two instantaneous eigenvalues are simply $E_{\pm }(t)=\pm \sqrt{(vt{)}^2+{\Delta }^2}$ (see inset in figure 3(b)). Starting with the system in the ground state at $t=-\infty$ , we can monitor the probability ofjumping onto the excited state at any time $t,$ which we plot in figure 3(b) for both the RT and the IT evolution and for several values of $\Delta$ (taking $v=1)$ . The RT data provide an illustration of the well-known Landau–Zener result: after a (relatively short) tunnelling time, and possibly a few oscillations, the probability of getting onto the excited state saturates to a value given by $P_{\mathrm{ex}}(t=\infty )={\mathrm{e}}^{-\pi {\Delta }^2/\hslash v}$ . As for the IT data, the initial (tunnelling) part and the subsequent plateau of the curves are similar to the RT case: the plateau value attained, call it $P_{\mathrm{ex}}^{\ast }$ , is indeed very close to the RT saturation value (in fact, asymptotically the same for $\Delta 0);$ ; after that, the IT evolutions start to filter out the ground-state component— initially present in the state with a small amplitude $1-P_{\mathrm{ex}}^{\ast }$ —through the usual mechanism of suppression of excited states, leading to $P_{\mathrm{ex}}(t)$ which is nicely fit by the curve

$$P_{\mathrm{ex}}(t)=\frac{P_{\mathrm{ex}}^{\ast }\exp \left(-2{\int }_0^t\mathrm{d}{t}^{\prime }[E_{+}(t^{\prime })-E_{-}(t^{\prime })]\right)}{(1-P_{\mathrm{ex}}^{\ast })+P_{\mathrm{ex}}^{\ast }\exp \left(-2{\int }_0^t\mathrm{d}{t}^{\prime }[E_{+}(t^{\prime })-E_{-}(t^{\prime })]\right)},\text{\hspace{1em}(24)}$$

which asymptotically goes to zero as $t\infty$ . This rather trivial effect of filtering, if on one hand explains the discrepancy between the IT and the RT evolution observed in the asymmetric double-well case of the previous section, is, on the other hand, of no harm at all: in contrast, it provides a quantitative improvement of IT over RT.

Figure 4. Parabolic washboard potential resulting in a logarithmically slow classical annealing. The minima are regularly located at positions $x_i=ia.$ and the dashed line shows the parabolic envelope potential.

In summary, the essential equivalence of IT and RT Schrodinger annealing (with,¨ moreover, a quantitative improvement of IT over RT) justifies practical implementations of quantum annealing based on imaginary-time quantum Monte Carlo schemes.

2.4. Other simple one-dimensional potentials with many minima

Moving on to multi-minima problems, we would like to mention one interesting onedimensional potential (see [28] for details) which shows a remarkably different behaviour of CA and QA. The problem was proposed and solved, for CA, by Shinomoto and Kabashima in [45, 47], and consists in a parabolically shaped washboard potential, sketched in figure 4. This example will display a logarithmically slow classical annealing, showing that CA may run into trouble even in simple models with no complexity whatsoever, whereas quantum mechanics can do much better in this case. The problem consists in a wiggly one-dimensional potential with barriers of individual height ≈B separating different local minima, regularly located a distance a apart one from each other, i.e., at positions $x_i=ai,i=0,\pm 1,\pm 2$ The ith local minimum is at energy ${ϵ}_i=k{a}^2{i}^2/2$ , so that the resulting envelope is parabolic. By writing the appropriate master equation governing the probability $P_i(t)$ that the particle is found in the ith valley at time t, and taking the continuum limit $a0$ , Shinomoto and Kabashima [45, 47] showed that the equation governing the evolution of $P(x,t)$ turns out to be a Fokker–Planck (FP) equation, equation (8), with an effective diffusion constant of the form

$$D_{\mathrm{eff}}(T)=\gamma a^2{\mathrm{e}}^{-B/k_{B}T},\text{\hspace{1em}(25)}$$

η $T)=k_{B}T/D_{\mathrm{eff}}(T)$ , and an effective drift potential $V(x)=k{x}^2/2$ given by the macroscopic parabolic envelope potential. This exponentially activated $D_{\mathrm{eff}}(T)$ makes the annealing behaviour of the $P(x,t)$ exceedingly slow. In fact, the surprising result of this exercise [45, 47] is that the optimal annealing schedule $T(t)$ is logarithmic and the residual energy converges to 0 at best as ${ϵ}_{\mathrm{res}}(t)\sim \log (t{)}^{-1}$ . The physical reason behind such a slow CA annealing is that the relaxation time $t_{\mathrm{relax}}=k_{B}T/(2\dot{\gamma }k{\dot{a}}^2){\mathrm{e}}^{B/k_{B}T}$ for the system to thermalize at any temperature T diverges exponentially at low T. As a result, the system will never be able to follow the decreasing T till the end of the annealing, by maintaining roughly the equilibrium value ${ϵ}_{\mathrm{pot}}=k_{B}T/2$ . Indeed, if we assume for instance $T(t)=T_0(1-t/\tau )$ , the relaxation of the systems will cease to be effective—i.e., the system will fall out of equilibrium—at a time $t^{\ast }$ , and temperature $T^{\ast }=T(t^{\ast })$ , at which $t_{\mathrm{relax}}\approx \tau ,i.e.$ ., when $k_B{T}^{\ast }$ ≈ B/ log γτ . The residual energy at this point cannot be smaller than the equipartition value $k_B{T}^{\ast }/2$ , hence ${ϵ}_{\mathrm{res}}\approx B/\log \gamma \tau$ as well. This freezing and falling out of equilibrium for classical systems with barriers seem to provide an ubiquitous source of logarithms in classical annealing [46].

The quantum mechanical approach to the same problem has been illustrated in [28]. In essence, starting from a tight-binding description in which the on-site energies ${ϵ}_i$ are supplemented by a time-dependent nearest-neighbour hopping term which contains the inverse mass $\Gamma ={\hslash }^2/2m$ in the typical semi-classical (WKB) form $\sim {\mathrm{e}}^{-\sqrt{V_h/\Gamma }}(V_h$ being an energy related to the details of the barrier), one can take, once again, the continuum limit $a0$ The dynamics for the $\psi (x,t)$ reduces, in strict analogy with the classical case, to an effective Schrodinger equation for a particle moving in the parabolic envelop potential¨ $V(x)=k{x}^2/2$ with an effective Laplacian coefficient ${\Gamma }_{\mathrm{eff}}(t)\propto {\mathsf{e}}^{-\sqrt{V_h/\Gamma (t)}}$ , which plays here the role that the effective diffusion constant in equation (25) played in the FP case. Contrary to the classical case, however, where an exponentially activated behaviour of the diffusion constant $D_{\mathrm{eff}}$ was strongly detrimental to the annealing (turning a power law into a logarithm), here the exponential WKB-like behaviour of ${\Gamma }_{\mathrm{eff}}$ does no harm at all: surprisingly, it improves the annealing. Indeed, as shown in [28], the power-law exponent ${\Omega }_{\mathrm{QA}}$ determining the decrease of the residual energy for a particle in a harmonic well, ${ϵ}_{\mathrm{res}}(\tau )\propto {\tau }^{-{\Omega }_{\mathrm{QA}}}$ , increases as one switches off the Laplacian coefficient more and more rapidly, tending to the value ${\Omega }_{\mathrm{QA}}1$ for an infinitely fast switching off.

We believe that one of the important points that makes QA so different from CA in this case is that the spectrum of the instantaneous eigenvalues of the quantum problem does not show any dangerous Landau–Zener avoided crossing, and, correspondingly, the ground-state wavefunction is always more peaked in the central valley (the minimum at $x_i=0)$ than elsewhere. As in the two-level case, a disorder in the width of the different valleys would drastically change this result, and likely ruin the power-law behaviour of $\mathrm{QA}$

3. Role of disorder and Landau–Zener tunnelling

Despite their disarming simplicity, the cases illustrated above turn out to be extremely informative in qualifying the profound difference of $\mathrm{QA}$ from CA, and their surprising consequences. Of course the cases studied, although instructive, do not possess the real ingredient which makes annealing, especially $\mathrm{QA}$ , difficult, i.e., some form of disorder in the distribution of the minima.

As an explicit demonstration of the effect of disorder in a very simple framework, consider the quantum walk of a particle in a given lattice, say cubic in D dimensions, of N sites. If an on-site energy ${ϵ}_i$ is associated with each site of the lattice, the Hamiltonian describing the quantum walk is just a disordered Anderson model

$$H(t)=\sum _i{ϵ}_i{\hat{n}}_i-\Gamma (t)\sum _{⟨i,j⟩}\left[c_j^{†}{c}_i+H.c.\right],\text{\hspace{1em}(26)}$$

where $\begin{array}{r}{H}_{\mathrm{cl}}={\sum }_i{ϵ}_i{\hat{n}}_i\end{array}$ represents the disordered on-site energy term, while the hopping term,
proportional to $-\Gamma (t)$ , provides the quantum kinetic term. The problem is to search for the lowest ${ϵ}_i$ in the lattice—with the on-site energy ${ϵ}_i$ randomly distributed in a given energy interval, say [0, 1]—by adiabatic evolution of a single-electron wavefunction from a very large ${\Gamma }_0=\Gamma (t⩽0)$ (when the wavefunction is extended over all sites) down to $\Gamma (\tau )=0$ Recent work on this problem has been reported in [24], where small instances $(N=20$ sites)

Figure 5. Instantaneous low-lying eigenvalues of a three-dimensional disordered Anderson model, on a lattice of size $10\times 10\times 10,$ as a function of the hopping integral
. Note the Landau–Zener avoided crossings, particularly the one occurring at smaller
, which can have a tremendously small gap 2. The inset shows a schematic of the Landau–Zener process.

of the problem are shown to be annealed rather fast by QA–AQC, with $E_{\mathrm{res}}(\tau )\sim 1/{\tau }^2$ . This disordered Anderson model is clearly related to the problem of searching for an item in an unsorted database of N entries, Grover’s problem [17, 48], one of the cases where QC is known to provide a speed-up, from the classical exhaustive search, $O(N)$ , to $O(\sqrt{N})$ . Indeed, the spatial search Grover’s problem associated with a structured database can be formulated in terms of equation (26) with ${ϵ}_{\bar{i}}=0$ while ${ϵ}_{j\ne \bar{i}}=1$ , where ^¯i is the desired entry of the lattice (database), and $⟨ij⟩$ denote neighbouring lattice sites (database entries) connected by the hopping in the given lattice structure; in this case, for instance, a quadratic speed-up has been demonstrated for $D⩾4[49]$ . However, the disorder associated with the on-site energies ${ϵ}_i.$ in a more general circumstance, makes the spectrum of the Hamiltonian in (26) much more complicated than in the spatial search Grover’s case. Obviously, the properties of the instantaneous spectrum of the Hamiltonian (26) depend on the dimensionality D of the lattice, and on the nearest-neighbour kinetic energy terms included (the so-called adjacency matrix of the quantum walker, in other words, on the topography of the landscape of the problem, specifying who is neighbour of whom). On quite general grounds, Anderson’s localization [50, 51] would predict that instantaneous wavefunctions of $H(t)$ are localized, for a genuinely disordered potential and for large enough mass (i.e., small enough $\Gamma (t)$ and hence kinetic energy bandwidth) in any $D>2$ (this localization occurs for all values of the mass in $D=1,2)$ . Therefore, quantum annealing should always, via a cascade of Landau–Zener events [35], end up into some localized state which has, $apriori$ , nothing to do with the actual potential minimum, i.e., the smallest ${ϵ}_i$ . Figure 5 shows the instantaneous eigenvalues of a disordered three-dimensional Anderson model, on a lattice of size $10\times 10\times 10$ , as a function of the nearest-neighbour hopping integral
. In the process of reducing
, the ground state encounters several ‘identity crisis’ associated with tunnelling from one region of the lattice to another. Some of these tunnelling amplitudes can be tremendously small, so that one has to wait for an astronomically large Landau–Zener characteristic time ${\tau }_{\mathrm{L}Z}$ (see section (2.2))

Figure 6. (Right) Instantaneous low-lying eigenvalues of a one-dimensional disordered Anderson model, on a lattice of size $N=100$ , as a function of the hopping integral
. Note the Landau– Zener avoided crossings, particularly the one occurring at smaller
, which has a very small gap $2\Delta .$ . (Left) The residual energy obtained by a Schrodinger imaginary-time evolution versus the ¨ annealing time τ. The inset shows the particular realization of the disordered on-site energies ${ϵ}_i,$ with the two arrows signalling the lowest minimum, at $i=60$ and the competing one a $i=72$

in order for the probability of ‘adiabatically following the ground state’ to be non-negligible. An even more concrete illustration of this problem is obtained for the $D=1$ case shown in figure 6. Here, for the realization of the disorder shown in the left inset, there are two competing minima: the absolute minimum at $i=60$ , and a higher one, very close in energy, at $i=72$ . For reasons that have to do with the quantum kinetic energy—the valley around $i=72$ being broader—the energy of the state localized around $i=72$ is the lowest one, down to very small values of $\Gamma >{\Gamma }_{\mathrm{LZ}}\approx 0.02$ . At such a value of
, a Landau–Zener crossing with an exceedingly small gap occurs, see right inset of figure 6. The system is, however, unable to follow the instantaneous ground state, for all practical purposes, and the final annealed energy of an IT–QA Schrodinger evolution remains trapped at ¨ ${ϵ}_{i=72}:$ no $E_{\mathrm{res}}\sim 1/{\tau }^2$ behaviour [24] is seen until astronomically large τ values are reached.

Moving to slightly more complicated models, let us consider the crucial role that disorder plays in the one-dimensional $(D=1)$ disordered Ising ferromagnet:

$$H=-\sum _i{J}_i{\sigma }_i^z{\sigma }_j^z-\Gamma \sum _i{\sigma }_i^x,\text{\hspace{1em}(27)}$$

where $J_i⩾0$ are non-negative random variables in the interval [0, 1], and $\Gamma$ is the transverse field inducing quantum fluctuations. Obviously, the ground state of the classical model at $\Gamma =0$ is the ferromagnetic state with all spins aligned up (or down). So, once again we are considering a trivial minimization problem, without any frustration at all (frustration requires competition between different couplings, and is therefore absent in $D=1$ nearest-neighbour models). However, even at $\Gamma =0$ there are important low-energy defects—domain walls between up and down ferromagnetic regions—which are typically pinned at the bonds where the couplings $J_i\mathbf{s}$ are the smallest. (For a finite system with periodic boundary conditions, domain walls appear in pairs, and separate sections of the system with alternating ↑ and $\downarrow$ ferromagnetic ground states.) Given two domain walls, pinned at weak $J_i<ϵ\ll 1$ bonds, typically at a distance $L\sim {ϵ}^{-1}\gg 1$ apart, healing the system via single spin flip moves requires flipping L spins, which can be a formidable barrier to tunnel through (or to go over, thermally). The situation is even worse at the critical point [19], ln $\begin{array}{r}{\Gamma }_c=\overline{\ln J}={\int }_0^1\mathrm{d}J\end{array}$ ln J , where it is known that the typical gap  between the ground state and the first relevant excited state is exponentially small, in the lattice size N [52],

$$[\Delta {]}_{\mathrm{typ}}\sim {\mathrm{e}}^{-\mathrm{const}\sqrt{N}}.\text{\hspace{1em}(28)}$$

The critical point ${\Gamma }_c-\mathrm{separating}$ the gapped large
quantum paramagnetic phase with all spins aligned along x, from the gapless low
classically disordered phase—is therefore a crucial source of small Landau–Zener gaps. We expect that, due to such exponentially small gaps, the typical characteristic Landau–Zener time ${\tau }_{\mathrm{LZ}}\propto {\Delta }_{\min }^{-2}$ beyond which the system will be able to adiabatically follow the ground state will be exponentially large in $N,$ resulting in a residual energy per spin ${ϵ}_{\mathrm{res}}(\tau )=E_{\mathrm{res}}(\tau )/N$ which, on average, will decrease only logarithmically as a function of the annealing time τ,

$$[{ϵ}_{\mathrm{res}}(\tau ){]}_{\mathrm{av}}\propto (\log \tau {)}^{-\zeta },\text{\hspace{1em}(29)}$$

for $N\to \infty$ . Work on this problem is currently in progress [53]. In conclusion, the system will have a very slow annealing while showing, at the same time, no complexity whatsoever: simple disorder is enough to yield slow annealing logarithmically.

4. Path-integral Monte Carlo quantum annealing

In order to move from toy problems with a manageable Hilbert space to real optimization problems, stochastic approaches are mandatory. As discussed in the introduction and in section 2.3, imaginary-time stochastic approaches are perfectly suitable to the goal: there is no gain in doing, on a classical computer, a Schrodinger evolution in real time [¨ 28].

A very simple quantum Monte Carlo approach, suitable to the proposed goal, is the PIMC approach. PIMC work in the QA context has been pioneered by Berne and collaborators, who have combined clever thermal and quantum annealing schemes within a path-integral framework, and applied these strategies to finding the minimal energy configurations of protein models on the continuum (with up to $N=46$ monomers) [31, 32], and the minimal energy configurations of Lennard-Jones clusters of up to $N=100$ atoms [33]. More recent PIMC-based QA work on the Lennard-Jones problem is reported in [34].

We start by briefly sketching the idea behind this approach with two introductory examples: the Ising case, as representative of discrete optimization problems, and the particle in a potential, as representative of continuum problems.

4.1. Path-integral Monte Carlo: introduction

The first, crucial, observation is that PIMC is intended to simulate the equilibrium behaviour of a system atfinite temperature T. Both these features are potential limitations of the method. To clarify this point, consider, for instance, the Edward–Anderson Ising glass in a transverse field: $\begin{array}{r}H=H_{\mathrm{cl}}+H_{\mathrm{kin}}=-{\sum }_{⟨ij⟩}{J}_{ij}{\sigma }_i^z{\sigma }_j^z-\Gamma (t){\sum }_i{\sigma }_i^x\end{array}$ . Strictly speaking, in the quantum annealing

context, this is a time-dependent Hamiltonian, of which we would like to follow the low-lying states (ideally, the ground state) as a function of time, by turning off the transverse field
(t). PIMC allows you to simulate the thermodynamics (at fixed strictly positive temperature T ) for a fixed value of
(t), by an approximate sampling of the quantum partition function

$$Z(T,\Gamma )=\mathrm{Tr}{\mathrm{e}}^{-\beta \left(H_{\mathrm{cl}}+H_{\mathrm{kin}}\right)}=\sum _{s^1}⟨s^1|{\mathrm{e}}^{-\beta \left(H_{\mathrm{cl}}+H_{\mathrm{kin}}\right)}|s^1⟩,\text{\hspace{1em}(30)}$$

where $s^1$ denotes a generic configuration of all the N spins. The idea behind the path integral is to reduce equation (30) to a classical partition function which is then sampled in the usual way using, for instance, a Metropolis Monte Carlo. In order to do that, one needs to split the exponential of the Hamiltonian, appearing in equation (30), into products of exponentials. This is allowed by the Trotter theorem, stating that

$${\mathrm{e}}^{-\beta (H_{\mathrm{cl}}+H_{\mathrm{kin}})}=\underset{P\to \infty }{\lim }{\left({\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{cl}}}{\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{kin}}}\right)}^P.\text{\hspace{1em}(31)}$$

Using this relationship, and inserting identities between the various exponentials, we get

$$Z(T,\Gamma )=\underset{P\to \infty }{\lim }\sum _{s^1\dots s^P}\exp \left(-\frac{\beta }{P}\sum _{k=1}^P{H}_{\mathrm{cl}}\left(s^k\right)\right)⟨s^1|{\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{kin}}}|s^2⟩\dots ⟨s^P|{\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{kin}}}|s^1⟩.\text{\hspace{1em}(32)}$$

The various configurations $s^k(k=1\cdots P)$ are often referred to as Trotter replicas of the original configuration $s^1$ . The next thing one needs to do is to calculate explicitly the relevant exponential of the kinetic term, $⟨s^k|{\mathbf{e}}^{-\stackrel{⌣}{H_{\mathrm{kin}}}/PT}|s^{k+1}⟩$ , between two generic configurations. This is sometimes very easy to do—like in the Ising transverse field case, or in the Laplacian case (see below)—but can also be very difficult, for other choices of $H_{\mathrm{kin}}$ (like in the travelling salesman case, see section (4.2.2)). In the Ising case, the problem factorizes into N independent sites, each of which involves a simple Pauli matrix expectation value, yielding

$$⟨s^k|{\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{kin}}}|s^{k+1}⟩=C^N\exp \left(\frac{\beta }{P}{J}^{\perp }\sum _i{s}_i^k{s}_i^{k+1}\right),\text{\hspace{1em}(33)}$$

where the transverse coupling $J^{\perp }$ is given by

$$J^{\perp }=-\frac{P}{2\beta }\ln (\tanh \beta \Gamma /P)>0\text{\hspace{1em}(34)}$$

(the constant $C$ is not relevant for our discussion). This kinetic term has a very transparent interpretation: $J^{\perp }$ gives a ferromagnetic Ising-like coupling between nearest-neighbour (k and $k+1)$ Trotter replicas of the same spin. In order to implement the problems numerically, a finite number of Trotter replicas $P$ is mandatory. This leads to an approximation, the error of which is proportional to the square of the Trotter break-up time, $O((\beta /P{)}^2)$ [54]. (Better Trotter break-ups, for given finite values of $P,$ can lead to smaller errors, see section (4.3), but we will concentrate here on the basic form proposed above.) For the full partition function we thus finally get, in the Ising case

$$Z(T,\Gamma )\approx C^{\mathrm{NP}}\sum _{s^1}\cdots \sum _{s^P}{\mathrm{e}}^{-\frac{\beta }{P}{S}_{D+1}}\text{\hspace{1em}(35)}$$ $$S_{D+1}=-\sum _{k=1}^P\left(\sum _{⟨ij⟩}{J}_{ij}{s}_i^k{s}_j^k+J^{\perp }\sum _i{s}_i^k{s}_i^{k+1}\right),\text{\hspace{1em}(36)}$$

which represents the partition function of a classical $(D+1)$ -dimensional anisotropic Ising system at temperature $P/\beta =PT$ The system has couplings $J_{ij}$ along the original D-dimensional lattice bonds (same for all Trotter replicas), and $J^{\perp }$ (same for all sites $i)$ along the extra Trotter dimension where the system has a finite length $P.$

Similar expressions hold, for instance, for the problem of a particle in a potential $V(\mathbf{x})$ where $H_{\mathrm{cl}}=V(\mathbf{x}),H_{\mathrm{kin}}=-\Gamma {\nabla }^2$ , and sums over configurations $\sum {}_{s^k}$ transform into integrals over the variables ${\mathbf{x}}^k$ . Similarly, the kinetic term contribution

$$⟨{\mathbf{x}}^k|{\mathrm{e}}^{-H_{\text{kin}}/PT}|{\mathbf{x}}^{k+1}⟩={\left(\frac{K^{\perp }}{2\pi PT}\right)}^{D/2}{\mathrm{e}}^{-\frac{1}{PT}\frac{K^{\perp }}{2}({\mathbf{x}}^k-{\mathbf{x}}^{k+1}{)}^2},\text{\hspace{1em}(37)}$$

where $D$ is the dimension of the x-space and $K^{\perp }=(PT{)}^2/(2\Gamma )$ admits a perfectly transparent interpretation: the transverse coupling $K^{\perp }$ between different Trotter replicas has the form of a spring coupling neighbouring configurations ${\mathbf{x}}^k$ and ${\mathbf{X}}^{k+1}$

In all cases, QA, in the present context, consists in externally controlling, during the PIMC dynamics, the value of the transverse field
—leaving T untouched—in much the same way as one externally controls T in classical simulated annealing [6]. This approach, which we will refer to as PIMC-QA, does not lead, therefore, to the simulation of a true quantum mechanical annealing dynamics, of the type implied by equation (2), but only to a MC annealing dynamics [42]. A clever mixing of thermal (reduction of T ) and quantum annealing (reduction of
)—possibly combined with decimation techniques on the number of Trotter replicas P as
is reduced—has been successfully applied in [31–33] to various continuum optimization problems (protein folding models and Lennard-Jones clusters). In the following, however, we will stick to the simple PIMC-QA scheme outlined above, whereby T is untouched during a QA simulation and the number of Trotter replicas P is also kept fixed, in order to more clearly discriminate possibly genuine quantum effects from thermal ones. We now move on to describe some of the results obtained so far with this technique on discrete optimization problems.

4.2. PIMC-QA applied to combinatorial optimization problems

In a quite general way, one could define a combinatorial optimization problem as the algorithmic task of minimizing any given cost function which depends on the configuration of variables assuming discrete values [4]. Further classifications are of course possible, but they somehow hide the fundamental fact that it is straightforward to map such problems over the search for the ground state of some Hamiltonian depending on Potts (or Ising) spin degrees of freedom [7, 55]. This is the case, for instance, of the travelling salesman problem [56–59], Boolean satisfiability [60, 61], vertex covering [62], graph colouring [63] and many others.

Random problem instances are of particular interest, because they can be investigated resorting to powerful techniques developed in the context of disordered statistical mechanics systems [1]. The physical approach to combinatorial optimization has often allowed the derivation of phase diagrams, telling us in which range of some control parameters hard instances are expected to be found [7, 55]. These analyses provide insight about the typicalcase complexity of problem solving, in contrast to the more rigorous but less informative worstcase complexity theory, which constitutes one of the corner stones of theoretical computer science [2]. The basic distinction between the P and NP complexity classes (that is, between problems for which a polynomial algorithm able to solve worst-case instances is or is not known) can sometimes be misleading. Easy instances of NP-complete problems (the hardest of all the NP problems, [3]) can easily be found (see e.g. [64]), while sometimes the optimization of instances of problems in P can take an exponential time using local search techniques (see e.g. [65]).

In the following, we shall briefly illustrate the results obtained on three classes of problems which have been recently addressed using PIMC-QA.

4.2.1. Random ising spin models. Determining the ground state of an Ising spin glass model can be an extraordinarily difficult task. To get the picture, it is enough to think that the number of possible configurations of a very small 32 × 32 square lattice Ising model is of order $10^{308}$ while the number of electrons in the universe is ${\dot{\mathrm{just}}}^{\prime }$ of the order $1{\bar{0}}^{80}{!}^3$ It can be rigorously shown that, in the 3D lattice case, the ground-state determination belongs to the NP-complete complexity class [5], but here we shall report results on the simpler 2D lattice case, where $E_{\mathrm{GS}}$ can be calculated up to sufficiently large lattice sizes [66]. The Hamiltonian of an Ising spin glass has been already discussed in section 4.1 and is given by equation (1).

Figure 7. Residual energy per site for an $80\times 80$ disordered 2D Ising model after CA and path integral Monte Carlo (PIMC) QA. We show the PIMC-QA data for the optimal value of $PT=1$ with $T=0.05$ and $P=20$ Trotter replicas. The actual inverse annealing rate τ used in QA has been rescaled (multiplied by P) for fair comparison with CA. Still, QA is faster than $\mathbf{CA}$

For a given 2D lattice size $L\times L(L$ up to 80), and for various quenched realizations of the random couplings $J_{ij}$ , drawn from a flat distribution in the interval $(-2,2)$ , we carried out several repeated classical and quantum annealings (for more details, see [35, 36]). At the end of both QA and CA, the system remains generally trapped at energy $E_{\mathrm{final}}=E_{\mathrm{GS}}+{ϵ}_{\mathrm{res}}$ and the efficiency of each protocol is monitored by considering the average residual energy ${ϵ}_{\mathrm{res}}(\tau )$ as a function of τ.

The annealing parameters T (CA) or
(QA) were decreased linearly from the initial value of $T_0=3$ or ${\Gamma }_0=2.5$ down to zero, with a total of τ MC steps per spin. In QA we used fixed values of $T_q=PT=1,1.5,2$ at several P values and prepared the initial state (same for all replicas) by a classical pre-annealing stage. The computational cost scales linearly with P, but increasing P beyond a certain characteristic length (see inset in figure 7) does not produce any further improvement. The choice of $P=20$ was found to be optimal. The moves proposed in both CA and QA are single-spin flip moves, but QA also attempts slightly more ‘global moves’ by proposing a spin flip for the spins $s_i^k,k=1,\dots ,P$ of all the Trotter replicas of a given site i. Figure 7 shows that QA is definitely superior to CA in the case of the Ising spin glass, in agreement with the experimental observation of significantly faster frequency-dependent relaxations during QA of the disordered magnet [15]. Theoretical arguments, however, can be given in favour of the fact that the residual energy will decrease, in both cases, logarithmically slow [35]

$${ϵ}_{\mathrm{res}}(\tau )\propto (\log \tau {)}^{-\zeta },$$

albeit with a possibly different annealing exponent ζ , which might favour QA over SA, ${\zeta }_{\mathrm{QA}}>{\zeta }_{\mathrm{CA}}$

Very recently, PIMC-QA and CA have been extensively applied [37] to the problem of finding the ground state of a random field Ising model (once again, a P-problem)

$$H=-J\sum _{⟨ij⟩}{\sigma }_i^z{\sigma }_j^z-\sum _i{h}_i{\sigma }_i^z,$$

for 1D, 2D and 3D cubic lattices, via the usual strategy, for QA, of adding and annealing a transverse field term. The results of [37] seem to confirm that the log behaviour is present in both cases, but that QA is indeed quantitatively better, with a larger $\zeta .$

Finally, let us mention that there is an Ising model case where PIMC-QA might be possibly worse than CA, according to [40]. This is the 1D Ising ANNNI model

$$H=-J_1\sum _i{\sigma }_i^z{\sigma }_{i+1}^z+J_2\sum _i{\sigma }_i^z{\sigma }_{i+2}^z,$$

with the point $J_2=J_1/2$ representing a highly frustrated case. While CA slows down only close to the frustrated case $J_2=J_1/2$ , PIMC-QA seems to have a uniformly poor performance even for small values of $J_2.$ deep in the ferromagnetic phase of the model [40].

4.2.2. Travelling salesman problem. Given N cities and their tabulated distances $d_{ij}$ , TSP consists in finding the shortest path connecting them, visiting each city only once and returning to the starting point. An account of the vast literature about algorithms for TSP can be found e.g. in [67], while three classical papers analysing physics approaches to the problem are [57–59].

As a first step to a QA optimization we have to choose a representation for the classical potential energy $H_{\mathrm{pot}}$ of a given configuration (in our case, the length of a tour), and, most crucially, a suitable source of quantum fluctuations $H_{\mathrm{kin}}$ . TSP can be mapped to a highly constrained Ising-like model—in a way similar to [56, 68]—in which each configuration of the system (a valid tour) is associated with an $N\times N0/$ 1-matrix $\hat{T}.$ . For every ordered sequence of cities, ${\hat{T}}_{i,j}=1$ if the tour visits city i immediately after city $j.$ , and ${\hat{T}}_{i,j}=0$ otherwise. For the symmetric TSP problem we wish to consider (a TSP with symmetric distance matrix $d_{ij}=d_{ji})$ the direct tour represented by a $\hat{T}.$ , and the reversed tour, represented by the transposed matrix ${\hat{T}}^t$ , have exactly the same length. It is then convenient to introduce the symmetric matrix $\hat{U}=\hat{T}+{\hat{T}}^t$ as representative of undirected tours. The length of a tour can now be written as

$$H_{\mathrm{pot}}(\hat{U})=\frac{1}{2}\sum _{(ij)}{d}_{ij}{\hat{U}}_{i,j}=\sum _{⟨ij⟩}{d}_{ij}{\hat{U}}_{i,j},\text{\hspace{1em}(38)}$$

where ij signifies counting each link only once. $H_{\mathrm{kin}}$ should be chosen in order to induce fluctuations generating the important elementary ‘moves’ of the problem, changing a tour into another one. Deciding which configurations are to become direct neighbours of a given configuration is indeed a crucial step, because it determines the problem’s effective landscape [69]. A very important move that is often used in heuristic TSP algorithms is the so-called two-opt move, which consists in eliminating two links in the current tour, $(c_1{c}_2)$ and $(c_{1^{\prime }}{c}_{2^{\prime }})$ , and rebuilding a new tour in which the connections are exchanged, $(c_1{c}_{1^{\prime }})$ and $(c_2{c}_{2^{\prime }})$ (see figure 8). Associating a spin variable +1 (−1) with each entry 1 (0), the whole two-opt move, when working with $\hat{U}$ matrices, can be represented by just four spin-flip operators

$$S_{⟨c_{1^{\prime }},c_1⟩}^{+}{S}_{⟨c_{2^{\prime }},c_2⟩}^{+}{S}_{⟨c_2,c_1⟩}^{-}{S}_{⟨c_{2^{\prime }},c_{1^{\prime }}⟩}^{-},$$

where, by definition, each $S_{⟨i,j⟩}^{\pm }$ flips an Ising spin variable (defined as $S_{⟨i,j⟩}^z=(2{\hat{U}}_{i,j}-1)=$ $\pm 1)$ at position $(i,j)$ and at the symmetric position $(j,i),i.e.,S_{⟨i,j⟩}^{\pm }=S_{i,j}^{\pm }{\bar{S}}_{j,i}^{\pm }$ . However, this kinetic Hamiltonian does not allow for an obvious Trotter discretization of the path integral (see discussion in section 4), and the PIMC scheme cannot deal with it (for this purpose, Green’s function MC methods, that do not use a Trotter break-up, would be in principle more effective, see section 5). We then introduce a drastic simplification to our kinetic energy term, replacing it altogether with a standard transverse Ising form, arriving finally at the Hamiltonian

Figure 8. (a) (Left) Representation of an eight-city tour, with the corresponding matrices ${\hat{T}}_{\mathrm{in}}$ and ${\hat{U}}_{\mathrm{in}}={\hat{T}}_{\mathrm{in}}+{\hat{T}}_{\mathrm{in}}^t$ . (Right) The final tour obtained when a two-opt move is performed, with a whole section reversed (dotted line). The matrices ${\hat{T}}_{\mathrm{fin}}$ and ${\hat{U}}_{\mathrm{fin}}$ are shown, the circles indicating the entries that have been switched $(01)$ by the two-opt move. The dotted circles in ${\hat{T}}_{\mathrm{fin}}$ are entries related to the trivial reversal of a section of the tour. (b) Average residual excess length found afte CA and PIMC-QA for a total time τ (in MC steps), for the $\bar{N}=1002$ instance pr1002 of the TSPLIB. PIMC-QA is once again faster than CA. (c), (d) the pr1002 instance considered, with its optimal tour and best tours obtained in the QA (c) and CA (d) cases.

$${\tilde{H}}_{\mathrm{TSP}}=\sum _{⟨ij⟩}{d}_{ij}\frac{\left(S_{⟨i,j⟩}^z+1\right)}{2}-\Gamma (t)\sum _{⟨ij⟩}\left[S_{⟨j,i⟩}^{+}+H.c.\right].\text{\hspace{1em}(39)}$$

This simplified form of kinetic energy no longer fulfils the constraint to take a valid tour to another valid tour, but this problem is avoided by proposing exclusively two-opt moves in the MC algorithm [38].

We tested our QA algorithm against CA [38] on a standard benchmark TSP problem, namely the printed circuit board instance pr1002 of the TSPLIB [70]. It is a structured TSP problem with $N=1002$ cities whose optimal tour length $L_{\mathrm{opt}}$ is known exactly. For CA, we chose an optimal initial temperature $T_0$ by first performing several CA with various short cooling times τ and starting from sufficiently high temperatures. The point where the cooling curves for different τs start to differ identifies an approximate ‘dynamical temperature’ $T_{\mathrm{dyn}}$ For pr1002, we obtained $T_{\mathrm{dyn}}\sim 100$ . As expected [67], the optimal $T_0$ for CA approximately coincides with $T_{\mathrm{dyn}}$ . Not surprisingly, for QA the same choice $PT\sim T_{\mathrm{dyn}}$ yields the optimal results, together with the choice ${\Gamma }_0=300$ . Figure 8 shows the results obtained [38] for the average percentage best-tour excess length ${ϵ}_{\mathrm{exc}}(\tau )=({\bar{L}}_{\mathrm{best}}(\tau )-L_{\mathrm{opt}})/L_{\mathrm{opt}}.$ , both with CA (filled squares) and with QA (open circles). As a reference, the best out of 1000 runs of the Lin–Kernighan algorithm (one of the standard local-search algorithms for TSP [67]) is also plotted (dashed line in figure 8). The results show that, once again, QA anneals more efficiently, even accounting for the extra factor $P$ in the total CPU time (rightmost open circles), reducing the error at a much steeper rate than CA.

4.2.3. Random Boolean satisfiability. In order to state the problem, consider a set of $N$ Boolean variables $z_1,\dots ,z_N$ , where z = 1 or 0 (‘true’ or ‘false’). Denoting by ${\zeta }_i$ the variable $z_i$ or its negation z¯_i, one then considers the disjunction (logical OR) of three variables $C=({\zeta }_i\vee {\zeta }_j\vee {\zeta }_k)$ , which is called a three-clause. The random 3-SAT problem consists in deciding if the conjunction (logical AND) of M different clauses $C_1\wedge C_2\cdot \cdot \cdot \wedge C_M-$ each clause being formed by three variables extracted at random among the N available, and appearing negated or direct with uniform probability—can be simultaneously satisfied by a truth value assignment $\{z_i\}$ . (In the exact cover variant studied in [27], a clause C is satisfied if one of the three variables is 1, the other two being 0.) If we associate an Ising spin variable $S_i=(-1{)}^{z_i}$ with each Boolean variable $z_i$ , we can assign to any clause $C_a$ involving three variables $z_i,z_j,z_k$ an energy $E_a$ given by

$$E_a=\frac{(1+J_{a,i}{S}_i)(1+J_{a,j}{S}_j)(1+J_{a,k}{S}_k)}{8},\text{\hspace{1em}(40)}$$

where the coupling $J_{a,i}$ assumes the value −1 if the variable $z_i$ appears negated in clause ${}^{a,}$ +1 otherwise. Evidently, $E_a=0$ if the corresponding clause is satisfied, $E_a=1$ otherwise.

As in the case of TSP, archives of hard structured instances exist [71]. In addition, statistical mechanics techniques can be used to determine the phase diagram of the random 3-SAT problem [7, 60, 61]. The main parameter determining the hardness of a formula is the ratio $\alpha =M/N$ between the number, M, of clauses and the number, $N,$ of variables. For $\alpha <{\alpha }_c\simeq 4.26$ it is typically possible to find satisfying assignments, but instances particularly hard to solve are expected to be found if $\alpha >{\alpha }_G\simeq 4.15[72]$ . It is expected that, due to the proliferation of an exponential number of metastable states acting as dynamical traps, local search gets trapped at an energy close to some finite threshold level, lower bounded by the so-called Gardner energy [73]. The trapping effect induced by the threshold states cannot be neglected when the instancesize is large $(N⩾10000)$ and large statistical fluctuations become sufficiently rare [72]. Smaller random formulae are, on the other hand, often easily solvable by classical simulated annealing and cannot be used as significant benchmarks.

We performed a first set of annealings over a single hard 3-SAT random instance with $N=10^4\mathrm{and}\alpha =4.24[39]$ . The kinetic term was given by a simple transverse field inducing single spin-flip fluctuations, like in the Ising case, since no clever sets of moves are known for 3-SAT, unlike the TSP case [74]. Using an efficient ad hoc algorithm (described in [72]), we verified that the chosen formula was actually satisfiable, as expected from theory for $\alpha <{\alpha }_c$ As in the case of the TSP optimization, we set both $T_0$ for CA and PT for QA equal to $T_{\mathrm{dyn}}=0.3$ . The optimal field-ramp range was found to be between ${\Gamma }_0=0.7$ and ${\Gamma }_f\simeq 10^{-3}$

A comparison between the performance of the optimal CA and the optimal QA at $P=50$ both with and without global (i.e., all $s_i^k,k=1,\dots ,P$ , are flipped) moves [39], is shown in figure 9. For each point, an average has been taken over 50 different realizations of the same experiment; in the case of QA, a second average was performed among the energies of the P replicas, which are in general different. It can be seen that the linear-schedule CA always performs better than the linear-schedule QA. No further improvement can be obtained for $P⩾100$ , see inset of figure 9—a much larger value than in the case of the Ising spin glass and the TSP instance—but we chose $P=50$ in order to extend the simulation time as much as possible. The asymptotic slope of the linear-schedule QA curves seems indeed to be definitely less steep than that of CA, independently of the number of replicas involved in the simulation and of the use of global moves.

Figure 9. Comparison between optimal linear-schedule classical (CA) and path-integral Monte Carlo quantum annealing (PIMC-QA) for a 3-SAT problem with $N=10^4$ and $\alpha =M/N=4.24$ CA always performs better than PIMC-QA simulated with $P=50$ Trotter replicas. The average performance of linear PIMC-QA is worse than that of CA, even if an improvement in the result can be obtained by introducing global moves (G) and by increasing P (in the inset the final average energy found by PIMC-QA after 2000 iterations for increasing P is plotted and compared with the average result of a CA of the same length, dashed line). The solid triangles are the data obtained by the field cycling PIMC-QA hybrid strategy described in [39].

The sobering message converged by this failure is that superiority of QA over CA is not universal, and is only achieved when we can use some understanding of the problem, especially when building the kinetic energy operator.

4.3. PIMC-QA ofa double-well: lessonsfrom a simple case

We would like to finish our discussion about path-integral Monte Carlo -based QA by mentioning recent results on a very simple case from which one can learn much about the limitations of the method [42]. Suppose we want to perform a QA optimization of the simple double-well potential which was investigated in section 2.2 using PIMC. One is then lead to simulate the behaviour of a closed polymer made up of P Trotter replicas $\{x^k\}(k=1,\dots ,P)$ of the original particle, held at temperature $\beta /P$ and moving in the potential $V_{\mathrm{asym}}$ with a nearest-neighbour spring coupling, as shown in equation (37). One can actually be more sophisticated than that, and perform a higher order Trotter break-up, correct to $O(\beta /P{)}^4$ instead of $O(\beta /P{)}^2$ , using, for instance, the Takahashi–Imada approximation [75]. Moreover, instead of performing single-bead moves, i.e., moves involving a single $x^k$ at a time, one can reconstruct, during the move, entire sections of the polymer, using the bisection method [29]. We have applied this rather sophisticated PIMC to our textbook double-well problem, working with a temperature $T=0.03V_0$ , a number of Trotter slices up to $P=160$ , and a bisection of level up to 4, i.e., involving up to $2^4+1$ replicas $x^k$ at each move. The initial temperature value of $\Gamma ={\hslash }^2/2m$ was taken to be ${\Gamma }_0=0.5$ , and its value was reduced linearly to 0 in a certain total number τ of Monte Carlo moves. The results, shown in figure 10 by solid circles, are altogether disappointing: the data barely start to cross below the level ${ϵ}_{\mathrm{res}}\approx 2\delta =0.2$ , which corresponds to the metastable minimum of the potential, for the largest values of τ simulated. This means that the system had no occasion, up to these large values of τ, of realizing that there was another minimum available through tunnelling. Moreover, the overall slope of the data is definitely less steep than that of direct Schrodinger annealing, shown for comparison¨ by solid triangles in figure 10. (The absolute values of τ are not comparable between the two sets of data, because they refer to different quantities: a real Schrodinger dynamics ¨ versus a Monte Carlo dynamics.) The situation improves substantially (see solid squares in figure 10) if we introduce, as candidate Monte Carlo moves, also the so-called instanton moves, i.e., basically classical trajectories that move from one minimum to the other (plus fluctuations) [76]. This is, however, not a fair game: we have substantially exploited a crucial information on the landscape which is generally not available for a complicated optimization problem!

Figure 10. Comparison between the Schrodinger annealing data, solid triangles and different types ¨ of PIMC-QA, on the double-well potential of section 2.2.

One other lesson we can learn, in the present context, is the role of the kinetic energy operator $H_{\mathrm{kin}}$ on which the quantum fluctuations are based. Up to now, we were using as $H_{\mathrm{kin}}$ the usual non-relativistic kinetic energy $H_{\mathrm{kin}}=p^2/2m$ , and annealed the system by increasing the mass m of the particle. (The propagator of this kinetic term is just the Gaussian, as shown in equation (37).) Imagine now we pretend that the particle has a relativistic ‘photon-like’ dispersion

$$H_{\mathrm{kin}}=\Gamma |p|\to ⟨{\mathbf{x}}^k|{\mathrm{e}}^{-\frac{\beta }{P}{H}_{\mathrm{kin}}}|{\mathbf{x}}^{k+1}⟩=\frac{1}{\pi }\frac{\hslash \Gamma \beta /P}{(x^k-x^{k+1}{)}^2+(\hslash \Gamma \beta /P{)}^2},\text{\hspace{1em}(41)}$$

and that we anneal the system by reducing to 0 the velocity
of the dispersion. The bisection method can be generalized for this kind of kinetic energy, and the results obtained are shown by solid diamonds in figure 10. As one notes, the residual energy versus τ is now considerably lower than in the non-relativistic (Gaussian) case, and even levels-off, for large τ , to the thermal limit $k_{B}T/2=0.015$ set by our finite temperature T. This example shows the important role played by the choice of the kinetic energy $H_{\mathrm{kin}}$ (as well as the limitations imposed by the unavoidable finite temperature T ).

Figure 11. Sketch of the annealing of the transverse field
, linearly in a time $\tau ={\tau }_1+{\tau }_2+\cdot \cdot \cdot {\tau }_n,$ or in a stepwise fashion.

Summarizing, PIMC-QA suffers evidently of a number of problems: (i) it is only a fake Monte Carlo annealing dynamics, in principle not fully representative of the true imaginarytime Schrodinger dynamics; (ii) the sampling of relevant ‘action’ might be highly inefficient¨ (recall the instanton problem above) and the cure for that might not be obvious at all; (iii) the finite temperature T imposes a lower thermal limit for the residual energy ${ϵ}_{\mathrm{res}}(\tau )$ below which we can never possibly go; (iv) the calculation of the propagator of ${\mathrm{e}}^{-{\overline{H}}_{\mathrm{kin}}/PT}$ might be very difficult for a kinetic term which we would like to implement (see section 2.2). These various problems clearly call for an alternative to PIMC in order to implement a stochastic QA approach. In principle, a useful alternative is given by a Green’s function Monte Carlo (GFMC) method, which we now turn to discuss.

5. Green’s function Monte Carlo quantum annealing

Green’s function Monte Carlo is a very effective quantum Monte Carlo approach for correlated quantum Hamiltonians on a lattice, the method of choice, for instance, for studying zero temperature properties of Heisenberg spin models [77]. It is naturally a method that works at zero temperature, by projecting out the ground-state wavefunction from an arbitrary initial state. We now illustrate the general ideas behind a GFMC-based quantum annealing approach.

5.1. Green’sfunction Monte Carlo quantum annealing: ideas

As outlined in the introduction, the ideal scope of a QA approach is to take some initial-state |(0) and let it evolve according to the Schrodinger dynamics, see equation ( ¨ 2), associated with a time-dependent Hamiltonian H(t) interpolating between an extreme quantum regime and the classical problem one is interested in. For the random Ising problem, for instance, a very natural choice is given by $\begin{array}{r}H(t)=-{\sum }_{⟨i,j⟩}{J}_{i,j}{\sigma }_i^z{\sigma }_j^z-\Gamma (t){\sum }_i{\sigma }_i^x\end{array}$ , where $\Gamma (t)$ is initially

very large, and slowly decreased towards zero in a given time τ , as sketched in figure 11. We will use the Ising case, throughout this section, to exemplify the ideas of the method, which can be adapted essentially to any other model in a straightforward way.

We already argued, see section 2.3, that an imaginary-time Schrodinger evolution is, for ¨ optimization purposes on a classical computer, equally good, and most likely even superior, to the standard real-time evolution. With this in mind, we can aim at solving the imaginary-time evolution in equation (6), instead of the real-time one in (2). If we also imagine the gradual decrease of
to be made stepwise, as sketched in figure 11, then the solution for $|\Psi (t)⟩$ is obtained by repeated applications of an imaginary-time evolution operator

$$|\Psi (\tau )⟩={\mathrm{e}}^{-H({\Gamma }_n){\tau }_n}\dots {\mathrm{e}}^{-H({\Gamma }_1){\tau }_1}|\Psi (0)⟩,\text{\hspace{1em}(42)}$$

where $\tau ={\tau }_1+\cdot \cdot \cdot {\tau }_n$ is the total annealing time, ${\Gamma }_1>{\Gamma }_2>\cdots >{\Gamma }_n\sim 0$ is a decreasing sequence of transverse fields and $H({\Gamma }_i)$ is a shorthand for $H(t)$ with a value ${\Gamma }_i$ of the transverse field. Each application of the imaginary-time evolution operator ${\mathbf{e}}^{-H({\Gamma }_i){\tau }_i}$ effectively tends to $\#lter$ out the corresponding ground state of $H({\Gamma }_i)$ from the state to which it is applied. In turn, ${\mathbf{e}}^{-H({\Gamma }_i){\tau }_i}$ can be obtained by a repeated application of many infinitesimal propagators of the type $[1-\Delta tH({\Gamma }_i)]$ , i.e.,

$$|\Psi (t+{\tau }_i)⟩={\mathrm{e}}^{-H({\Gamma }_i){\tau }_i}|\Psi (t)⟩=\underset{m\to \infty }{\lim }[1-\Delta tH({\Gamma }_i){]}^m|\Psi (t)⟩,\text{\hspace{1em}(43)}$$

where $\Delta t={\tau }_i/m$ . The Green’s function Monte Carlo (GFMC) is just a stochastic technique which implements equation (43). More precisely, if we define, recursively,

$${\psi }_{n+1}(x^{\prime })=\sum _x{G}_{x^{\prime },x}^{(\Gamma )}{\psi }_n(x),\text{\hspace{1em}(44)}$$

where ${\psi }_n(x)=⟨x|{\psi }_n⟩,|x⟩$ is a shorthand for a generic spinconfiguration describing the Hilbert space of the problem, and the Green’s function $G_{x^{\prime },x}^{(\Gamma )}$ is given by

$$\begin{array}{c}{G}_{x^{\prime },x}^{(\Gamma )}=⟨x^{\prime }|G^{(\Gamma )}|x⟩=⟨x^{\prime }|1-\Delta t[H(\Gamma )-E_T]|x⟩\\ =(1+\Delta t{E}_T){\delta }_{x^{\prime },x}-\Delta t⟨x^{\prime }|H(\Gamma )|x⟩,\end{array}\text{\hspace{1em}(45)}$$

one can show that—for large n—the iterated state ${\psi }_n$ converges (apart from a normalization constant) to the ground state ${\psi }_{\mathrm{GS}}^{(\Gamma )}(x)$ $H(\Gamma )$ , if $\Delta t$ is chosen to be suitably small [77]. $(E_T$ is an estimate of the ground-state energy which allows us to reduce the statistical fluctuations, see also [78].)

The problem in equation (44) looks superficially similar to an ordinary Markovian master equation [44], with a few very crucial differences: (i) ${\psi }_n(x)$ are not probabilities, but amplitudes; (ii) the Green’s function $G_{x^{\prime },x}^{(\Gamma )}$ in equation (45), unlike the transition probability of a master equation, is not necessarily made of non-negative elements, and is, in general, not columnnormalized, $\begin{array}{r}{\sum }_{x^{\prime }}{G}_{x^{\prime },x}^{(\Gamma )}\ne 1\end{array}$ , unlike a Markov transition probability [44]. In summary,
the process underlying the iterated-power method is not a properly defined Markov chain, and, therefore, it cannot be immediately simulated, as it stands, with a Monte Carlo approach.

Problem (ii) above can be quite serious: if some of the matrix elements of $G_{x^{\prime },x}^{(\Gamma )}$ are negative, no possible interpretation of it as a ‘transition probability’ is possible. This is at the heart of the so-called sign problem [78] in quantum Monte Carlo. In the following, we will assume that a choice of basis is possible in which no sign problem exists, $i.e.$ , all matrix elements of $G^{(\Gamma )}$ are non-negative, ${\dot{G}}_{x^{\prime },x}^{(\Gamma )}⩾0$ . This is certainly true for the Ising glass in a transverse field. More generally, since the choice of the kinetic energy to be used in $\mathrm{QA}$ is at our disposal, it is wise to choose the signs of $H_{\mathrm{kin}}$ in equation (1) such that no sign problem occurs. Still, we miss the correct column-normalization:

$$\sum _{x^{\prime }}{G}_{x^{\prime },x}^{(\Gamma )}=1+\Delta t{E}_T-\Delta t\sum _{x^{\prime }}{H}_{x^{\prime },x}\stackrel{\text{ def }}{=}{b}_x\ne 1.\text{\hspace{1em}(46)}$$

A way out of this difficulty is to factorize $G^{(\Gamma )}$ in terms of a stochastic matrix $p_{x^{\prime },x}-\mathsf{by}$ definition [44], a matrix with all positive elements $p_{x^{\prime },x}⩾0$ , and with the normalization condition ${\sum }_{x^{\prime }}{p}_{x^{\prime },x}=1$ for all columns—times the scale factor $b_x$ defined above. Indeed,
with the previous definition (46) of $b_x$ , the matrix

$$p_{x^{\prime },x}=G_{x^{\prime },x}^{(\Gamma )}/b_x\text{\hspace{1em}(47)}$$

is trivially positive and column normalized and, therefore, it is a suitable transition matrix for a Markov chain in x-space.

The crucial idea is then to extend the configuration space where the Markov process is defined, adding to x a non-negative weight factor (hereafter, the weight) w. This extended configuration space will be labelled by $(x,w)$ . The pair (x, w) is often called a walker, because it will be the basic entity in the Markov chain ‘random walk’. The weight part will take care of $b_x$ , while x will be taken care of by $p_{x^{\prime },x}$ . More precisely, if $(x_n,w_n)$ indicates a walker at iteration time $n,$ in this extended configuration space, we set up the following Markov process:

(a) Generate $x_{n+1}=x^{\prime }$ with probability $p_{x^{\prime },x_n}$

(b) Update the weight with $w_{n+1}=w_n{b}_x$

(48)

In words, the walker performs a random walk in the Hilbert space (x) of the system and in the weight space (w); such a random walk is composed of a standard Markov chain in x-space, associated with $p_{x^{\prime },x}$ , and a multiplicative process for the weight $w_n{w}_{n+1}=w_n{b}_x$ By moving in this way, the walkers visit every point in the (x, w)-space with a probability $P_n(x_n,w_n)$ whose first w-moment can be shown to be proportional to the wavefunction $\begin{array}{r}{\psi }_n,{\psi }_n(x)\propto \int \mathrm{d}{w}_n{w}_n{P}_n(x,w_n)\end{array}$ [79].

The algorithm in (48) is the basic version of a GFMC [77]. In this form, however, the algorithm simply does not work in practice. The reason for this failure is not difficult to grasp. While $x_n\to x_{n+1}$ is a plain Markov process, the weight update $w_n{w}_{n+1}=w_n{b}_{x_n}$ is a multiplicative process with random factors $b_{x_n},$ , which is prone to very large fluctuations $[77,79]:w_n$ might grow large, or become negligibly small, in just a few iterations of the algorithm, and the whole method would go wild, because error bars in the calculations of the averages grow in an uncontrolled way. The cure to this disease goes through the introduction of many walkers, simulated in parallel, and through performing occasional ‘reconfigurations’ of their weights, via the so-called branching [78]. In practice, one propagates simultaneously a set of M walkers defined by weights $w_i$ and configurations $x_i$ , for $i=1,\dots ,M$ . Before the variance of the weights $w_i$ becomes too large, one appropriately redefines the set of walkers— by reproducing some of them and deleting some others—in such a way as to drop those with excessively small weight, and to generate copies of the more important ones⁴.

The analogy of such a many-walker GFMC with a genetic-like algorithm [80] is worth noting. Each walker (x, w) plays the role of an individual that propagates (mutates) increasing or decreasing itsfitness—represented by the accumulated weight $w,$ related to the wavefunction amplitude $\psi (x)$ . A mutation is here simply a step of the algorithm, which attempts a single spin flip of a certain site in the configuration, this is what the off-diagonal matrix elements $⟨x^{\prime }|H|x⟩$ do for an Ising model in the transverse field. At certain times, branching occurs, which modifies the population of individuals by favouring the survival of those with highest fitness (largest $w)$ . The only genetic feature that is missing in the quantum mechanical case is the possibility of cross-breeding (mixing of genetic codes of two configurations, to give rise to new configurations); this would correspond to non-local moves which change the configurations in a global way.

The final, important, ingredient that makes the algorithm work is the so-called importance sampling [30]. It can be seen, in the genetic analogy proposed before, as a way of proposing mutations (single spin flips) that instead of being equally probable, with a common matrix element
, are biased by a function which guides the system towards the most representative configurations. More precisely, suppose we have a reasonable guess ofthe ground-state ${\psi }_{\mathrm{GS}}(x)$ in the form of some nodeless wavefunction ${\psi }_T(x)$ , known as trial (or guiding) wavefunction. It is then enough to substitute the Green’s function $G^{(\Gamma )}$ with the so-called importance-sampling Green’s function,

$${\bar{G}}_{x^{\prime },x}={\psi }_T(x^{\prime })G_{x^{\prime },x}^{(\Gamma )}{\psi }_T^{-1}(x),\text{\hspace{1em}(49)}$$

which just rescales by an extra factor ${\psi }_T(x^{\prime })/{\psi }_T(x)$ a transition from $x\mathrm{to}{x}^{\prime }.$ thus favouring those transitions $x{x}^{\prime }$ where ${\psi }_T(x^{\prime })/{\psi }_T(x)$ is largest. In general, ${\hat{G}}_{x^{\prime },x}$ is not symmetric, but one can still apply to it the same decomposition in equation (47), defining the corresponding Markov chain (48) with

$$\begin{array}{l}{p}_{x^{\prime },x}={\bar{G}}_{x^{\prime },x}/{\bar{b}}_x,\\ {\bar{b}}_x={\sum }_{x^{\prime }}{\bar{G}}_{x^{\prime },x}=1+\Delta t{E}_T-\Delta t\frac{{\sum }_{x^{\prime }}{\psi }_T(x^{\prime })H_{x^{\prime },x}}{{\psi }_T(x)}=1-\Delta t(E_L(x)-E_T),\end{array}\text{\hspace{1em}(50)}$$

where the local energy $E_L(x)$ is defined as

$$E_L(x)=\frac{⟨{\psi }_T|H|x⟩}{⟨{\psi }_T|x⟩}.\text{\hspace{1em}(51)}$$

Quite amusingly, if the guessed trial wavefunction ${\psi }_T$ coincides with the ground-state wavefunction ${\psi }_T(x)={\psi }_{\mathrm{GS}}(x)$ , then $E_L(x)=E_{\mathrm{GS}}$ and ${\bar{b}}_x=1-\Delta t(E_{\mathrm{GS}}-E_T)$ are constant, and one can show that statistical fluctuations in the calculation vanish exactly. (This is the so-called zero variance property [78].) Therefore, by variationally improving the quality of the guiding wavefunction ${\psi }_T(x)$ ) one can substantially improve the quality of the calculation (e.g., reducing the error bars).

We will soon show that the use of a good importance-sampling guiding wavefunction is indeed a crucial ingredient of a GFMC–QA application; without importance sampling, GFMC–QA is simply worthless. Therefore, we start illustrating our application of GFMC on the random Ising model with a discussion of possible trial wavefunctions we have tested, and the difficulties encountered.

5.2. GFMC–QA on random ising model

Finding a good trial wavefunction for a random Ising model in a transverse field is a highly non-trivial task. The first idea that comes to mind is a kind of ‘mean-field’ wavefunction, made up of a product of single-site factors as

$$|{\psi }_T^{(\mathrm{MF})}⟩=\prod _{i=1}^N\left(\frac{{\mathrm{e}}^{+\frac{h_i}{2}}|\uparrow {⟩}_i+{\mathrm{e}}^{-\frac{h_i}{2}}|\downarrow {⟩}_i}{\sqrt{2\cosh (h_i)}}\right),\text{\hspace{1em}(52)}$$

where $\{h_i\}$ , the local fields on each site $i,$ are variational parameters to be optimized for each given value of the transverse field
. The optimization of $\{h_i\}$ amounts to finding the minimum of the variational energy,

$$\begin{array}{rl}& E_T^{(\mathrm{MF})}(\{m_i\};\Gamma )=⟨{\psi }_T^{(\mathrm{MF})}\left|H\right|{\psi }_T^{(\mathrm{MF})}⟩\\ & =-{\sum }_{⟨i,j⟩}{J}_{i,j}{m}_i{m}_j-\Gamma {\sum }_i\sqrt{1-m_i^2},\end{array}\text{\hspace{1em}(53)}$$

where $m_i=\tanh (h_i)$ are the local magnetizations. The stationary conditions required by the minimization read

$$0=\frac{\partial E_T^{(\mathrm{MF})}}{\partial h_i}=-\left(1-m_i^2\right)\left(\sum _{j\in \mathcal{N}(i)}{J}_{i,j}{m}_j\right)+\Gamma m_i\sqrt{1-m_i^2}\forall i,\text{\hspace{1em}(54)}$$

where $\mathcal{N}(i)$ indicates the set of nearestneighbours of site i. As it turns out, finding solutions of equation (54) with optimal variational energies is simple only for large enough
, where the quantum paramagnetic solution $h_i=m_i=0$ is found. Such a solution, representing the $\Gamma =+\infty$ ground state with all spins aligned along the +xˆ direction, survives down to some value ${\Gamma }_{\mathrm{cr}}$ of the transverse field, below which non-trivial solutions of equation (54)— with non-vanishing local magnetizations $m_i\ne 0⟶\mathrm{start}$ to appear. However, the number of those solutions (local minima)—found for instance by a straightforward conjugate gradients algorithm [8]—is large [81]. In the low-
region, our minimization problem is just the quantum counterpart of the well-known Weiss mean-field approach for the classical random Ising model [82], which is known to run into difficulties in the classical glassy phase. In a sense, minimizing equation (53) in the low-
glassy phase is not much simpler than finding the classical ground state of the problem for $\Gamma =0$ We transformed a minimization task in a discrete space of variables, $S_i=\pm 1$ , into one where the variables are continuous, $m_i\in (-1,1)$ , but the task itself is of comparable difficulty. We refer the reader to [81] for more details on this aspect.

A second, quite natural, choice of trial wavefunction that comes to mind is a Boltzmannlike wavefunction of the form

$${\psi }_T^{(\beta )}(\{S_i\})=\mathcal{N}(\beta ){\mathrm{e}}^{-(\beta /2)H_{\mathrm{cl}}(\{S_i\})},\text{\hspace{1em}(55)}$$

where $1/\beta$ plays the role of an effective temperature, with $\beta$ a variational parameter to be optimized, and $H_{\mathrm{cl}}(\{S_i\})$ (see equation (1)) is the classical energy of a given configuration $\{S_i\}.N(\beta )$ is an appropriate normalization factor, which we will not need to calculate. Once again, for large $\Gamma$ we expect to find $\beta =0$ (the exact $\Gamma =+\infty$ solution), while, by decreasing
, larger and larger values of $\beta$ will favour regions where the ‘potential energy $H_{\mathrm{cl}}(\{S_i\})$ has a local minimum, until we get, for $\Gamma =0$ , to the asymptotic limit $\beta \to \infty (\mathrm{ideally})$ , required by a wavefunction which is perfectly localized in the global minimum (see below for a discussion of this point).

To calculate the expectation value of energy with the Boltzmann-like trial wavefunction in equation (55), as a function of the single parameter $\beta ,$ , we used a standard variational Monte Carlo (VMC) algorithm [78]. Figure 12 shows (left panel) the optimal value ${\beta }_{\mathrm{opt}}$ of $\beta$ which minimizes the variational energy $E_{\mathrm{tot}}^{(\mathrm{Boltz})}=⟨{\psi }_T^{(\beta )}\left|H\right|{\psi }_T^{(\beta )}⟩$ , for several values of the transverse field $\Gamma .$ Note that ${\beta }_{\mathrm{opt}}$ saturates for small $\Gamma$  to about ${\dot{\beta }}_{\mathrm{opt}}\approx 2$ , somewhat surprisingly at first sight, since, for $\Gamma 0$ , one would expect ${\beta }_{\mathrm{opt}}+\infty$ , in such a way that the classical ground state dominates (i.e., equation (55) becomes a delta-like function localized in the exact classical ground state). This is clearly an effect of a severe ergodicity loss of the VMC algorithm, which is not difficult to understand. For a given $\beta$ , indeed, the VMC samples a probability distribution ${\left|{\psi }_T^{(\beta )}(x)\right|}^2={\mathrm{e}}^{-\beta H_{\mathrm{cl}}(x)}$ with single spin-flip moves; its efficiency  in exploring the phase space, therefore, is exactly identical to that of a classical Metropolis Monte Carlo at the temperature $T=1/\beta$ . Finding an optimal $\beta$ for a given
is therefore totally equivalent to asking what is the effective temperature of a classical Ising spin glass which provides the best approximation to the wavefunction of a quantum Ising glass at zero temperature and non-zero
. Now, from classical spin-glass physics [1, 83] we know that a threshold energy $E_{\mathrm{th}}$ exists below which the system has a finite complexity, i.e., it displays an exponentially (∼ exp N) large number of metastable minima. Close to this threshold energy, the relaxation of any local algorithm towards equilibrium becomes exceedingly slow—the algorithm gets stuck for a long time in every minimum visited—and the average quantities measured are not representative of their true thermodynamical values. Evidently, for $\Gamma 0.$ the VMC algorithm does not visit the regions near the true minima of the classical energy, but wanders in a high-energy band of metastable states, separated by moderate energy barriers. In such a case, a small and finite value of $\beta$ allows us to (still) overcome such barriers, so as to find slightly more favourable local minima, while perfect localization $(\beta \to +\infty )$ in a wrong excited state would lead to an average bigger residual energy.

Figure 12. (Left) Plot of the optimal $\beta ,{\beta }_{\mathrm{opt}},$ , for the ‘Boltzmann’ trial wavefunction $|{\psi }_T^{(\beta )}⟩$  defined in equation (55), for several values of
. The dashed line is a fit to the data. (Right) The variation residual diagonal energy ${ϵ}_{\mathrm{res}}^{(\mathrm{Boltz})}=⟨{\psi }_T^{({\beta }_{\mathrm{opt}})}|H_{\mathrm{cl}}|{\psi }_T^{({\beta }_{\mathrm{opt}})}⟩/N$ − _GS corresponding to the ${\beta }_{\mathrm{opt}}$ is shown in the left panel, for several
. The dashed horizontal line represents the best residual energy ever achieved, for $\Gamma >0.01$ , by employing the mean-field trial wavefunction in [52].

The right panel in figure 12 shows the optimal variational residual energy ${ϵ}_{\mathrm{res}}^{(\mathrm{Boltz})}=$ $⟨{\psi }_T^{({\beta }_{\mathrm{opt}})}\left|H_{\mathrm{cl}}\right|{\psi }_T^{({\beta }_{\mathrm{opt}})}⟩/N-{ϵ}_{\mathrm{GS}}$ corresponding to the optimal $\beta$ shown in the left panel, for   several values of the transverse field
. For large
values, the variational residual energy per spin is of order 1, as a random spin configuration would do. By decreasing
, we note that the variational residual energy saturates, for small
, to finite non-zero values, of order 0.03 per spin, in agreement with the previously noted saturation in the optimal ${\beta }_{\mathrm{opt}}$ , due to ergodicity breaking. (A closer inspection shows that the variational residual energy is actually non-monotonic for $\Gamma <0.25$ , again an artefact of sampling difficulties.) Note, however, that this saturation value is definitely below the best results provided by the previously discussed ${\psi }_T^{(\mathrm{MF})}$ , which is of order 0.035 per spin (shown for comparison by a dashed horizontal line). Therefore, with all its pitfalls, the Boltzmann-like trial wavefunction in equation (55) provides, at low
, a marginally better approximation of the true GS, than that obtained by the mean-field Ansatz in equation (52). Moreover, ${\psi }_T^{(\beta )}$ is also much better behaved, and simpler to optimize, as far as the minimization problem is concerned. For these reasons, we decided to work out our GFMC results using the Boltzmann-like wavefunction only.

We finally present the results of a GFMC-based QA approach, where the transverse field $\Gamma$ is decreased stepwise during the simulation, while, at the same time, the importance sampling Boltzmann-like wavefunction is changed according to the corresponding value of the variational parameter ${\beta }_{\mathrm{opt}}(\Gamma )$ . As a benchmark, we will compare GFMC-QA outcomes with the path-integral Monte Carlo quantum annealing (PIMC-QA) and classical simulated annealing (CA) results described in [35, 36]. We reduce the coupling
in equation (1) at each Monte Carlo step (MCS) in a linear way. We start from an initial large enough value of the transverse field, ${\Gamma }_0=2.5$ , and then we set ${\Gamma }_n={\Gamma }_0(1-n/\tau )$ during the nth MCS $(0⩽n<\tau ).\tau$ is the total annealing time (see figure 11) measured as the total number of MCS performed by the algorithm. We used $M=20$ walkers, and performed branching at each MCS, because the low
region is affected by severe weight instabilities which would otherwise make the algorithm unstable (for the initial, large
, part of the annealing one could consider branching less often, as weights are well under control). However, even with this very conservative choice, the weights go sometimes completely wild if a
of order $10^{-7}$ is reached (i.e., for long annealing times τ). As a consequence, we decided to cut off the ${\Gamma }_n$ annealing schedule in such a way that the final
is $\approx 10^{-4}$ and not 0. This finally guarantees a good weight stability. (Whenever it was possible to perform annealings with smaller cut-offs on
, and approximately the same slope $1/\tau$ , we checked that the results obtained are not very sensitive to the cut-off chosen.)

For each value of ${\Gamma }_n$ the trial wavefunction used is the Boltzmann-like one, defined in equation (55), with a variational parameter ${\beta }_{\mathrm{opt}}({\Gamma }_n)$ , which corresponds to the instantaneous optimal value. (Practically, we used for $\beta (\Gamma )$ the fitting function shown in figure 12, left panel.) Figure 13 shows the best residual energy per spin ever reached during the annealing simulation, for several values of $\tau ,$ averaged over ten independent repetitions of the whole annealing process (due to computer-time limitations, a single run is shown for the largest, $\tau >10^8$ , annealings). For comparison, the CA and PIMC–QA data obtained in [35] are also shown. Note first that the τ axes of the three calculations are completely unrelated. The GFMC τ is measured in units in which a MCS is just a single spin $\begin{array}{r}\#ip,\end{array}$ while MCS for CA and PIMC–QA are intended as sweeps of the entire lattice of N spins (including all the $P=20$ Trotter slices, for the PIMC case). For this reason, we also present the CA and PIMC–QA data in a shifted time axis where τ is multiplied by a factor $N=80^2$ (rightmost curves). Although the GFMC–QA data are strictly below both CA and PIMC–QA data, on the same perspin timeunit (i.e., compared to the shifted CA and PIMC–QA data), it is clear that the GFMCslope is still worse than that of PIMC–QA, and indeed surprisingly similar to CA. Moreover, the CPUtime needed for a single spin flip in GFMC is much larger than the corresponding single-spin move in CA or PIMC–QA (each GFMC move costs of order N operations).

Let us pause to consider the similarity between the CA and the GFMC–QA slopes that figure 13 suggests. This similarity must be somehow related to the fact that we have used, as importance wavefunction for the GFMC, a Boltzmann-like wavefunction, ${\psi }_T(x)\propto {\mathrm{e}}^{-(\beta /2)H_{\mathrm{cl}}(x)}$ . At the level of a plain variational Monte Carlo (VMC), we already pointed out that such a choice of wavefunction amounts to sampling ${\left|{\psi }_T^{(\beta )}(x)\right|}^2={\mathrm{e}}^{-\beta H_{\mathrm{cl}}(x)}$ and is therefore totally equivalent to a classical Metropolis MC at temperature $T=1/\beta$ . If, during the GFMC simulation, we neglect the weights associated with the walkers (as well as the associated branching process), we will be carrying over a completely classical simulation where the generated configurations are distributed (see equation (50)) according to

$$p_{x^{\prime },x}=\frac{{\bar{G}}_{x^{\prime },x}}{{\bar{b}}_x}\propto \Gamma \frac{{\psi }_T^{(\beta )}(x^{\prime })}{{\psi }_T^{(\beta )}(x)}=\Gamma \exp (-({\beta }_{\mathrm{opt}}(\Gamma )/2)[H_{\mathrm{cl}}(x^{\prime })-H_{\mathrm{cl}}(x)]).$$

As a consequence, the Markov process in x-space will obey a classical detailed balance condition

$$p_{x^{\prime },x}{\mathrm{e}}^{-{\beta }_{\mathrm{opt}}{H}_{\mathrm{cl}}(x)}=p_{x,x^{\prime }}{\mathrm{e}}^{-{\beta }_{\mathrm{opt}}{H}_{\mathrm{cl}}(x^{\prime })}.\text{\hspace{1em}(56)}$$

In other words, a GFMC–QA without weights would be just a computationally heavy way of doing a CA with a peculiar form of annealing of the effective temperature ${\beta }_{\mathrm{opt}}(\Gamma )$ (note, in passing, that such an effective temperature never gets too low, since ${\beta }_{\mathrm{opt}}$ saturates to around ${\beta }_{\mathrm{opt}}$ ≈ 2 for low
). Quantum mechanics enters, therefore, only through the weights that the GFMC carries over (and the unavoidable branching process which makes the multiplicative process of weight updating numerically stable). Evidently, such a weight updating is, in the present disordered case, not sufficiently strong and effective as to make the resulting averages really different from the underlying classical physics governing the Markov chain in x-space, and the resulting GFMC–QA data are rather close to the CA ones (although they are much more computer-time demanding).

Figure 13. The average best residual energy obtained by Green’s function Monte Carlo (GFMC) $\mathrm{QA}$ for the 80 × 80 instance of the random Ising model studied in section 4.2.1, versus the total annealing-time τ . Previous results obtained by classical simulated annealing (CA) and by pathintegral Monte Carlo quantum annealing $(PIMC–QA)$ with $P=20$ Trotter slices [35, 36] are shown for comparison. The GFMC time-unit is a single spin flip, while CA and PIMC–QA Monte Carlo time units are sweeps of the entire lattice (see [36]). Importance sampling is performed by using the optimal trial wavefunction ${\psi }_T^{(\mathrm{Boltz})}$ . The transverse field is linearly reduced down to $10^{-4}$ in a total annealing time τ, starting from ${\Gamma }_0=\mathrm{2.5.}$ . We used here $M=20$ walkers and performed branching at every Monte Carlo step (MCS) $(t_B=1)$

6. Summary and conclusions

We have illustrated several applications of quantum annealing strategies to a range of problems going from textbook toy models—displaying in a clear way the crucial differences between classical and quantum annealing—all the way to challenging optimization problems (random Ising models, TSP, 3-SAT). The techniques used to implement QA are either deterministic Schrodinger’s evolutions, for the toy models, or path-integral Monte Carlo (PIMC), and¨ Green’s function Monte Carlo (GFMC) approaches, for the hard optimization problems.

Worth mentioning, here, are a few other approaches or applications which we have had no occasion to discuss so far. An alternative $T=0$ projection QMC scheme (similar in spirit to the GFMC discussed above, but not using importance sampling) is presented in [84], and applied to an infinite range ±J spin glass model. Quantum annealing strategies have also been applied to a one-dimensional kinetically constrained model [84, 85], based on Ising spins ${\sigma }_i^z$ subject to a downward longitudinal field $h,$ and experiencing a barrier for flipping, whenever the neighbouring spin ${\sigma }_{i-1}^z$ is ↓; based on a semiclassical treatment of the problem, it has been shown [84, 85] that a quantum tunnelling of the barrier is more effective than the thermal barrier jump in annealing the system towards the ultimate ground state, with all spins ↓. Finally, let us also mention the application of PIMC–QA to problems of image reconstruction and information processing [86].

As a way of summary, we would like to mention here some of the major points touched upon.

Is quantum really better? Although, in the examples illustrated, QA often wins over CA, sometimes it does not (the negative result for the 3-SAT case is particularly instructive in this respect). The fact that QA wins over its classical counterpart is a priori not guaranteed, the outcome of the battle being strongly related to the landscape of the problem one deals with. The role of the disorder is particularly crucial in the quantum case. Both the two major algorithmic achievements of quantum computation, i.e., Grover’s and Schor’s algorithm, deal with problems without disorder, where the intrinsic parallelism of a quantum algorithm can be fully exploited. We suspect that even simple optimization problems, like for instance the random Ising model in one-dimension, can be quite hard, if tackled by an adiabatic quantum computation, precisely due to the all-important role played by disorder.

Role of kinetic energy. The choice of the kinetic energy is clearly all important in QA. Section 4.3, illustrating the improvements in annealing a double-well potential upon using a relativistic kinetic energy, is particularly instructive.

Concerning more technical points.

Limitations of PIMC. Path-integral Monte Carlo QA, although often successful over its classical counterpart, suffers from several limitations (finite temperature T, possible sampling problems for the action, difficulties with the Trotter break-up), which suggest investigating alternative quantum Monte Carlo approaches to QA.

Limitations of GFMC. Green’s function Monte Carlo seems, in principle, a very good tool for QA. However, GFMC needs, as a crucial ingredient, a reasonably good guess of trial variational wavefunctions. Finding a good trial function ψ_T(x) for the problem at hand is an essential part of a GFMC–QA application, and constitutes the delicate point of the whole algorithm.

In conclusion, it is quite plausible that quantum annealing (adiabatic quantum computation), although potentially useful and sometimes more convenient than classical annealing, is not capable, in general, of finding solutions to NP-complete problems in polynomial time. Nevertheless, understanding when and how quantum mechanics can quantitatively improve on the solution of hard optimization problems is still an open and timely issue. Moreover, there are fundamental physics issues behind the dynamics of quantum disordered systems which certainly deserve further study.

Acknowledgments

The work described has been sponsored, over the years, by MIUR through FIRB RBAU017S8R004, FIRB RBAU01LX5H, COFIN2003 and COFIN2004, and by INFM (‘Iniziativa trasversale calcolo parallelo’). It is a real pleasure to thank Roman Martonˇak,´ Roberto Car, Demian Battaglia, Lorenzo Stella and Osvaldo Zagordi for their collaboration on various stages with the work described in the present review. Finally, instructive discussions with Rosario Fazio, Riccardo Zecchina, Sandro Sorella and Silvio Franz are gratefully acknowledged.

References

[1] Parisi G, Mezard M and Virasoro M 1987 ´ Spin Glass Theory and Beyond (Singapore: World Scientific)

[2] Papadimitriou C H and Steiglitz K 1998 Combinatorial Optimization: Algorithms and Complexity (New York: Dover)

[3] Garey M R and Johnson D S 1979 Computers and Intractability (San Francisco: Freeman)

[4] Papadimitriou C H 1993 Computational Complexity (Reading, MA: Addison-Wesley)

[5] Barahona F 1982 J. Phys. A: Math. Gen. 15 3241

[6] Kirkpatrick S, Gelatt C D Jr and Vecchi M P 1983 Science 220 671

[7] Dubois O, Monasson R, Selman B and Zecchina R (ed) 2001 Special Issue on NP-hardness and phase transition Theor. Comput. Sci. 265 1–307

[8] Press W H, Teukolsky S A, Vetterling W T and Flannery B P 1992 Numerical Recipes in C: The Art ofScientific Computing 2nd edn (Cambridge: Cambridge University Press)

[9] Amara P, Hsu D and Straub J E 1993 J. Phys. Chem. 97 6715

[10] Finnila A B, Gomez M A, Sebenik C, Stenson C and Doll J D 1994 Chem. Phys. Lett. 219 343

[11] Kadowaki T and Nishimori H 1998 Phys. Rev. E 58 5355

[12] Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1792

[13] Ray P, Chakrabarti B K and Chakrabarti A 1989 Phys. Rev. B 39 11828

[14] Thirumalai D, Li Q and Kirkpatrick T R 1989 J. Phys. A: Math. Gen. 22 3339

[15] Brooke J, Bitko D, Rosenbaum T F and Aeppli G 1999 Science 284 779

[16] Aeppli G and Rosenbaum T F 2005 Quantum Annealing and Related Optimization Methods ed A Das and B K Chakrabarti (Berlin: Springer) p 159

[17] Nielsen M and Chuang I L 2000 Quantum Computation and Quantum Information (Cambridge: Cambridge University Press)

[18] Shor P W 1997 SIAM J. Comput. 26 1484

[19] Fisher D S 1995 Phys. Rev. B 51 6411

[20] Farhi E, Goldstone J, Gutmann S and Sipser M 2000 Quantum computation by adiabatic evolution Preprint quant-ph/0001106

[21] Aharonov D, van Dam W, Kempe J, Landau Z, Lloyd S and Regev O 2004 Adiabatic quantum computation i equivalent to standard quantum computation Preprint quant-ph/0405098

[22] Gantmakher F R 1960 The Theory ofMatrices (New York: Chelsea)

[23] Messiah A 1962 Quantum Mechanics vol 2 (Amsterdam: North-Holland)

[24] Suzuki S and Okada M 2005 J. Phys. Soc. Japan 74 1649

[25] Landau L D 1932 Phys. Z. Sowjetunion 2 46

[27] Farhi E, Goldstone J, Gutmann S, Lapan J, Lundgren A and Preda D 2001 Science 292 472

[26] Zener C 1932 Proc. R. Soc. A 137 696

[28] Stella L, Santoro G E and Tosatti E 2005 Phys. Rev. B 72 014303

[29] Ceperley D M 1995 Rev. Mod. Phys. 67 279

[30] Kalos M H and Whitlock P A 1986 Monte Carlo Methods vol 1 (New York: Wiley)

[31] Lee Y H and Berne B J 2000 J. Phys. Chem. A 104 86

[32] Lee Y H and Berne B J 2001 J. Phys. Chem. A 105 459

[33] Liu P and Berne B J 2003 J. Chem. Phys. 118 2999

[34] Gregor T and Car R 2005 Chem. Phys. Lett. 412 125

[35] Santoro G E, Martonˇak R, Tosatti E and Car R 2002´ Science 295 2427

[36] Martonˇak R, Santoro G E and Tosatti E 2002´ Phys. Rev. B 66 094203

[37] Sarjala M, Petaj¨ a V and Alava M 2006 ¨ J. Stat. Mech. P01008

[38] Martonˇak R, Santoro G E and Tosatti E 2004 ´ Phys. Rev. E 70 057701

[39] Battaglia D A, Santoro G E and Tosatti E 2005 Phys. Rev. E 71 066707

[40] Sen P and Das P K 2005 Quantum Annealing and Related Optimization Methods ed A Das and B K Chakrabarti (Berlin: Springer) p 305

[41] Stadler P F and Schnabl W 1992 Phys. Lett. A 161 337

[42] Stella L, Santoro G E and Tosatti E 2006 Phys. Rev. B 73 144302 (Preprint cond-mat/0512064)

[43] Das A and Chakrabarti B K 2005 Quantum Annealing and Related Optimization Methods (Lecture Notes in Physics) (Berlin: Springer)

[44] van Kampen N G 1992 Stochastic Processes in Physics and Chemistry ed (Amsterdam: North-Holland) (Revised and enlarged)

[45] Shinomoto S and Kabashima Y 1991 J. Phys. A: Math. Gen. 24 L141

[46] Huse D A and Fisher D S 1986 Phys. Rev. Lett. 57 2203

[47] Kabashima Y and Shinomoto S 1991 J. Phys. Soc. Japan 60 3993

[48] Grover L K 1997 Phys. Rev. Lett. 79 325

[49] Childs A M and Goldstone J 2003 Spatial search by quantum walk Preprint quant-ph/0306054

[50] Anderson P W 1958 Phys. Rev. 109 1492

[51] Lee P A and Ramakrishnan T V 1985 Rev. Mod. Phys. 57 287

[52] Fisher D S and Young A P 1998 Phys. Rev. B 58 9131

[53] Caneva T, Fazio R and Santoro G E Quantum annealing of a random Ising chain in a transverse field (in preparation)

[54] Suzuki M 1976 Prog. Theor. Phys. 56 1454

[55] Nishimori H 2001 Statistical Physics of Spin Glasses and Information Processing (Oxford: Oxford University Press)

[56] Hopfield J J and Tank D W 1986 Science 233 625

[57] Mezard M and Parisi G 1986 ´ Europhys. Lett. 2 913

[58] Percus A G and Martin O C 1999 J. Stat. Phys. 94 739

[59] Dean D S, Lancaster D and Majumdar S N 2005 J. Stat. Mech. L01001

[60] Mezard M, Parisi G and Zecchina R 2002´ Science 297 812

[61] Mezard M and Zecchina R 2002´ Phys. Rev. E 66 56126

[62] Weigt M and Hartmann A K 2000 Phys. Rev. Lett. 84 6118

[63] Braunstein A, Mulet R, Pagnani A, Weigt M and Zecchina R 2003 Phys. Rev. E 68 036702

[64] Mitchell D, Selman B and Levesque H 1992 Proc. ofthe Tenth National Conf. on Artificial Intelligence (Menlo Park, CA: AAAI Press) p 459

[65] Mezard M, Ricci-Tersenghi F and Zecchina R 2003 ´ J. Stat. Phys 111 505

[66] De Simone C, Diehl M, Junger M, Mutzel P, Reinelt G and Rinaldi G 1995¨ J. Stat. Phys. 80 487

[67] Johnson D S and McGeoch L A 2003 Local Search in Combinatorial Optimization (Princeton, NJ: Princeton University Press) p 215

[68] Kadowaki T 2002 Study of optimization problems by quantum annealing Preprint quant-ph/0205020

[69] Stadler P F 1995 Complex Systems and Binary Networks ed R Lopez-Pena, R Capovilla, R Garca-Pelayo, H Waelbroeck and F Zertuche (Berlin:Springer) p 77

[70] See http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95

[71] See http://www.intellektik.informatik.tu-darmstadt.de/SATLIB/

[72] Battaglia D A, Kola´ˇr M and Zecchina R 2004 Phys. Rev. E 70 036107

[73] Montanari A, Parisi G and Ricci-Tersenghi F 2004 J. Phys. A: Math. Gen. 37 2073

[74] Michalewicz Z 1996 Genetic Algorithms + Data Structures = Evolution Programs 3rd edn (Berlin: Springer)

[75] Takahashi M and Imada M 1984 J. Phys. Soc. Japan 53 3765

[76] Negele J W and Orland H 1988 Quantum Many-Particle Systems (Reading, MA: Addison-Wesley)

[77] Trivedi N and Ceperley D M 1990 Phys. Rev. B 41 4552

[78] Foulkes W M C, Mitas L, Needs R J and Rajagopal G 2001 Rev. Mod. Phys. 73 33

[79] Calandra Buonaura M and Sorella S 1998 Phys. Rev. B 57 11446

[80] Holland J H 1971 Associative Information Processing (New York: Elsevier)

[81] Stella L and Santoro G E Quantum annealing of an Ising spin glass by Green’s function Monte Carlo (in preparation)

[82] Binder K and Young A P 1986 Rev. Mod. Phys. 58 801

[83] Castellani T and Cavagna A 2005 J. Stat. Mech. P05012

[84] Das A and Chakrabarti B K 2005 Quantum Annealing and Related Optimization Methods ed A Das and B K Chakrabarti (Berlin: Springer) p 217

[85] Das A, Chakrabarti B and Stinchcombe R B 2005 Phys. Rev. E 72 026701

[86] Inoue J-I 2005 Quantum Annealing and Related Optimization Methods ed A Das and B K Chakrabarti (Berlin:Springer) p 171

J. Phys. A: Math. Theor. 41 (2008) 209801 (2pp)

Corrigendum

Optimization using quantum mechanics: quantum annealing through adiabatic evolution

Giuseppe E Santoro and Erio Tosatti 2006 J. Phys. A: Math. Gen. 39 R393-R431

In our recent paper (Giuseppe E. Santoro and Erio Tosatti 2006 J. Phys. A: Math. Gen. 39, R393-R431) we reviewed some of the recent work in the field of quantum annealing, alias adiabatic quantum computation. Here we point out two early references, due to Apolloni, de Falco and collaborators and dating back to 1988, where the idea of quantum annealing was first put forward and tested on hard combinatorial optimization problems.

The idea of quantum annealing is an elegant and fascinating alternative to classical thermal simulated annealing; it consists in helping the system escape the local minima using quantum mechanics — by tunneling through the barriers rather than thermally overcoming them, with an artificial and appropriate source of quantum fluctuations (the counterpart of the temperature) initially present and slowly (adiabatically) switched off.

In our recent review on this subject [1] we erroneously indicated what appeared to us as being the earliest references [2–4] in which similar ideas were first explicitly formulated, and tested in numerical simulations

We have in the meantime learned that the idea of using stochastic processes based on quantum mechanics with the goal of minimizing classical complex functions was indeed formulated and tested much earlier, in 1988, by the group of Apolloni, de Falco and collaborators, who applied these ideas to combinatorial optimization problems like graphpartitioning [5, 6].

We also stress that the idea of quantum annealing is inherently related to the idea of adiabaticity. Computing by adiabatic evolution of a quantum system has become a quite popular idea in the Quantum Computing community, where it is commonly known as adiabatic quantum computation, and commonly traced back to Ref. [7]. Quantum annealing and adiabatic quantum computation are, however, two names given by two different scientific communities to the very same idea.

Acknowledgments

It is a pleasure to thank Professor de Falco for bringing to our attention the references presented here.

References

[1] Santoro G E and Tosatti E 2006 J. Phys. A: Math. Gen. 39, R393-R431

[2] Amara P, Hsu D and Straub J E 1993 J. Phys. Chem. 97 6715

[3] Finnila A B, Gomez M A, Sebenik C, Stenson C and Doll J D 1994 Chem. Phys. Lett. 219 343

[4] Kadowaki T and Nishimori H 1998 Phys. Rev. E 58 5355

[5] Apolloni A, Cesa-Bianchi N and de Falco D (1990) in Stochastic Processes, Physics and Geometry, Proceedings

1751-8113/08/209801+02$30.00 © 2008 IOP Publishing Ltd Printed in the UK

ofthe Ascona-Locarno Conference, 4-9 July 1988, S. Albeverio et al. editors, pp. 97-111 (World Scientific, 1990)

[6] Apolloni A, Carvalho C and de Falco D (1989) Stoc. Proc. Appl. 33, 233-244

[7] Farhi E, Goldstone J, Gutmann S and Sipser M 2000 Quantum computation by adiabatic evolution (Preprint quant-ph/0001106)