1204.2907

## QUANTUM ANNEALING: FROM VIEWPOINTS OF STATISTICAL PHYSICS, CONDENSED MATTER PHYSICS, AND COMPUTATIONAL PHYSICS ## SHU TANAKA Department of Chemistry, University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan E-mail: shu-t@chem.s.u-tokyo.ac.jp ## RYO TAMURA Institute for Solid State Physics, University of Tokyo, 5-1-5, Kashiwanoha, Kashiwa-shi, Chiba, 277-8501, Japan International Center for Young Scientists, National Institute for Materials Science, 1-2-1, Sengen, Tsukuba-shi, Ibaraki, 305-0047, Japan E-mail: tamura.ryo@nims.go.jp In this paper, we review some features of quantum annealing and related topics from viewpoints of statistical physics, condensed matter physics, and computational physics. We can obtain a better solution of optimization problems in many cases by using the quantum annealing. Actually the efficiency of the quantum annealing has been demonstrated for problems based on statistical physics. Then the quantum annealing has been expected to be an efficient and generic solver of optimization problems. Since many implementation methods of the quantum annealing have been developed and will be proposed in the future, theoretical frameworks of wide area of science and experimental technologies will be evolved through studies of the quantum annealing. Keywords : Quantum annealing; Quantum information; Ising model; Optimization problem ## 1. Introduction Optimization problems are present almost everywhere, for example, designing of integrated circuit, staff assignment, and selection of a mode of transportation. To find the best solution of optimization problems is difficult in general. Then, it is a significant issue to propose and to develop a method for obtaining the best solution (or a better solution) of optimiza- tion problems in information science. In order to obtain the best solution, a couple of algorithms according to type of optimization problems have been formulated in information science and these methods have yielded practical applications. Furthermore, since optimization problem is to find the state where a real-valued function takes the minimum value, it can be regarded as problem to obtain the ground state of the corresponding Hamiltonian. Thus, if we can map optimization problem to well-defined Hamiltonian, we can use knowledge and methodologies of physics. Actually, in computational physics, generic and powerful algorithms which can be adopted for wide application have been proposed. One of famous methods is simulated annealing which was proposed by Kirkpatrick et al. 1,2 In the simulated annealing, we introduce a temperature (thermal fluctuation) in the considered optimization problems. We can obtain a better solution of the optimization problem by decreasing temperature gradually since thermal fluctuation effect facilitates transition between states. It is guaranteed that we can obtain the best solution definitely if we decrease temperature slow enough. 3 Then, the simulated annealing has been used in many cases because of easy implementation and guaranty. The quantum annealing was proposed as an alternative method of the simulated annealing. 4-11 In the quantum annealing, we introduce a quantum field which is appropriate for the considered Hamiltonian. For instance, if the considered optimization problem can be mapped onto the Ising model, the simplest form of the quantum fluctuation is transverse field. In the quantum annealing, we gradually decrease quantum field (quantum fluctuation) instead of temperature (thermal fluctuation). The efficiency of the quantum annealing has been demonstrated by a number of researchers, and it has been reported that a better solution can be obtained by the quantum annealing comparison with the simulated annealing in many cases. Figure 1 shows schematic picture of the simulated annealing and the quantum annealing. In optimization problems, our target is to obtain the stable state at zero temperature and zero quantum field, which is indicated by the solid circle in Fig. 1. Recently, methods in which we decrease temperature and quantum field simultaneously have been proposed and as a result, we can obtain a better solution than the simulated annealing and the simple quantum annealing. 12-14 Moreover, as an another example of methods in which we use both thermal fluctuation and quantum fluctuation, novel quantum annealing method with the Jarzynski equality 15,16 was also proposed, 17 which is based on nonequilibrium statistical physics. Fig. 1. Schematic picture of the simulated annealing and the quantum annealing. Our purpose is to obtain the ground state at the point indicated by the solid circle. <details> <summary>Image 1 Details</summary>  ### Visual Description ## Diagram: Quantum vs. Simulated Annealing ### Overview The image is a diagram illustrating the relationship between "Quantum field" and "Simulated annealing" (Temperature). It uses two perpendicular arrows to represent these concepts, with the intersection point indicating a common origin or starting point. ### Components/Axes * **Vertical Axis:** Labeled "Quantum field" with an arrow pointing upwards. The text "Quantum annealing" is written vertically alongside the arrow, indicating the direction of increasing quantum field strength. * **Horizontal Axis:** Labeled "Temperature" with an arrow pointing to the right. The text "Simulated annealing" is written horizontally alongside the arrow, indicating the direction of increasing temperature. * **Origin:** A black circle marks the intersection of the two arrows, representing the starting point for both quantum and simulated annealing processes. ### Detailed Analysis * The vertical arrow, labeled "Quantum field," points upwards, suggesting an increase in the quantum field strength. The text "Quantum annealing" is written vertically alongside the arrow. * The horizontal arrow, labeled "Temperature," points to the right, suggesting an increase in temperature. The text "Simulated annealing" is written horizontally alongside the arrow. * The intersection of the two arrows is marked by a black circle, indicating a common origin or starting point for both processes. ### Key Observations * The diagram visually represents the relationship between quantum annealing and simulated annealing, suggesting they are orthogonal or independent processes. * The arrows indicate the direction of increasing strength or temperature for each process. * The common origin suggests a shared starting point or initial state. ### Interpretation The diagram illustrates the conceptual difference between quantum annealing and simulated annealing. Quantum annealing relies on manipulating the "Quantum field," while simulated annealing relies on manipulating "Temperature." The orthogonal arrows suggest that these two approaches are distinct and potentially complementary methods for optimization or problem-solving. The common origin implies that both processes might start from the same initial state or problem formulation. </details> In this paper, we review the quantum annealing method which is the generic and powerful tool for obtaining the best solution of optimization problems from viewpoints of statistical physics, condensed matter physics, and computational physics. The organization of this paper is as follows. In Sec. 2, we review the Ising model which is a fundamental model of magnetic systems. The realization method of the Ising model by nuclear magnetic resonance is also explained. In Sec. 3, we show a couple of implementation methods of the quantum annealing. In Sec. 4, we explain two optimization problems - traveling salesman problem and clustering problem. The quantum annealing based on the Monte Carlo method for the traveling salesman problem is also demonstrated. In Sec. 5, we review related topics of the quantum annealing - Kibble-Zurek mechanism of the Ising spin chain and order by disorder in frustrated systems. In Sec. 6, we summarize this paper briefly and give some future perspectives of the quantum annealing. ## 2. Ising Model In this section we introduce the Ising model which is a fundamental model in statistical physics. A century ago, the Ising model was proposed to explain cooperative nature in strongly correlated magnetic systems from a microscopic viewpoint. 18 The Hamiltonian of the Ising model is given by $$\mathcal { H } _ { I s i n g } = - \sum _ { i , j } J _ { i j } \sigma _ { i } ^ { z } \sigma _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } \sigma _ { i } ^ { z } , \quad \sigma _ { i } ^ { z } = \pm 1 ,$$ where the summation of the first term runs over all interactions on the defined graph and N represents the number of spins. If the sign of J ij is positive/negative, the interaction is called ferromagnetic/antiferromagnetic interaction. Spins which are connected by ferromagnetic/antiferromagnetic interaction tend to be the same/opposite direction. The second term of the Hamiltonian denotes the site-dependent longitudinal magnetic fields. Although the Ising model is quite simple, this model exhibits inherent rich properties e.g. phase transition and dynamical behavior such as melting process and slow relaxation. For instance, the ferromagnetic Ising model with homogeneous interaction ( J ij = J for ∀ i, j ) and no external magnetic fields ( h i = 0 for ∀ i ) on square lattice exhibits the second-order phase transition, whereas no phase transition occurs in the Ising model on onedimensional lattice. Onsager first succeeded to obtain explicitly free energy of the Ising model without external magnetic field on square lattice. 19 After that, a couple of calculation methods were proposed. Furthermore, these calculation methods have been improved day by day, and the new techniques which were developed in these methods have been applied for other more complicated problems. Since the Ising model is quite simple, we can easily generalize the Ising model in diverse ways such as the Blume-Capel model, 20,21 the clock model, 22,23 and the Potts model. 24,25 By analyzing these models, relation between nature of phase transition and the symmetry which breaks at the transition point has been investigated. Then, it is not too much to say that the Ising model has opened up a new horizon for statistical physics. The Ising model can be adopted for not only magnetic systems but also systems in wide area of science such as information science. Optimization problem is one of important topics in information science. As we mention in Sec. 4, optimization problem can be mapped onto the Ising model and its generalized models in many cases. Then some methods which were developed in statistical physics often have been used for optimization problem. In Sec. 2.1, we show a couple of magnetic systems which can be well represented by the Ising model. In Sec. 2.2, we review how to create the Ising model by Nuclear Magnetic Resonance (NMR) technique as an example of experimental realization of the Ising model. ## 2.1. Magnetic Systems In many cases, the Hamiltonian of magnetic systems without external magnetic field is given by $$\begin{array} { r l } { \hat { \mathcal { H } } = - \sum _ { i , j } J _ { i j } \hat { \sigma } _ { i } \cdot \hat { \sigma } _ { j } } \\ { = - \sum _ { i , j } J _ { i j } \left ( \hat { \sigma } _ { i } ^ { x } \cdot \hat { \sigma } _ { j } ^ { x } + \hat { \sigma } _ { i } ^ { y } \cdot \hat { \sigma } _ { j } ^ { y } + \hat { \sigma } _ { i } ^ { z } \cdot \hat { \sigma } _ { j } ^ { z } \right ) , } \end{array}$$ where ˆ σ α i denotes the α -component of the Pauli matrix at the i -th site. The form of this interaction is called Heisenberg interaction. The definitions of Pauli matrices are $$\hat { \sigma } ^ { x } \colon = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} , \quad \hat { \sigma } ^ { y } \colon = \begin{pmatrix} 0 & - i \\ i & 0 \end{pmatrix} , \quad \hat { \sigma } ^ { z } \colon = \begin{pmatrix} 1 & 0 \\ 0 & - 1 \end{pmatrix} , \quad ( 3 )$$ where the bases are defined by $$| \uparrow \rangle \colon = { 1 \choose 0 } \, , \quad | \downarrow \rangle \colon = { 0 \choose 1 } \, .$$ In this case, magnetic interactions are isotropic. However, they become anisotropic depending on the surrounded ions in real magnetic materials. In general, the Hamiltonian of magnetic systems should be replaced by $$\hat { \mathcal { H } } = - \sum _ { i , j } J _ { i j } \left ( c _ { x } \hat { \sigma } _ { i } ^ { x } \hat { \sigma } _ { j } ^ { x } + c _ { y } \hat { \sigma } _ { i } ^ { y } \hat { \sigma } _ { j } ^ { y } + c _ { z } \hat { \sigma } _ { i } ^ { z } \hat { \sigma } _ { j } ^ { z } \right ) .$$ When | c x | , | c y | > | c z | , the xy -plane is easy-plane and the Hamiltonian becomes XY-like Hamiltonian. On the contrary, when | c z | > | c x | , | c y | , the z -axis is easy-axis and the Hamiltonian becomes Ising-like Hamiltonian. Such anisotropy comes from crystal structure, spin-orbit coupling, and dipoledipole coupling. Moreover, even if there is almost no anisotropy in magnetic interactions, magnetic systems can be regarded as the Ising model when the number of electrons in the magnetic ion is odd and the total spin is halfinteger. In this case, doubly degenerated states exist because of the Kramers theorem. These states are called the Kramers doublet. When the energy difference between the ground states and the first-excited states ∆ E is large enough, these doubly-degenerated ground states can be well represented by the S = 1 / 2 Ising spins. Table 1 shows examples of the magnetic materials which can be well represented by the Ising model on one-dimensional chain, two-dimensional square lattice, and three-dimensional cubic lattice. Table 1. Examples of magnetic materials which can be represented by the Ising model on chain (one-dimension), square lattice (two-dimension), and cubic lattice (three-dimension). | Material | Spatial dimension | Total spin | Type of interaction | J/k B | References | |-------------------------------|---------------------|--------------|-----------------------|-------------|--------------| | K 3 Fe(CN) 6 | One (chain) | 1 2 | Antiferromagnetic | - 0 . 23 K | 26-28 | | CsCoCl 3 | One (chain) | 1 2 | Antiferromagnetic | - 100 K | 29,30 | | Dy(C 2 H 5 SO 4 ) 2 · 9 H 2 O | One (chain) | 1 2 | Ferromagnetic | 0 . 2 K | 31-33 | | CoCl 2 · 2NC 5 H 5 | One (chain) | 1 2 | Ferromagnetic | 9 . 5 K | 34,35 | | CoCs 3 Br 5 | Two (square) | 1 2 | Antiferromagnetic | - 0 . 23 K | 36-38 | | Co(HCOO) 2 · 2 H 2 O | Two (square) | 1 2 | Antiferromagnetic | - 4 . 3 K | 39-42 | | Rb 2 CoF 4 | Two (square) | 1 2 | Antiferromagnetic | - 91 K | 43,44 | | FeCl 2 | Two (square) | 1 | Ferromagnetic | 3 . 4 K | 45,46 | | DyPO 4 | Three (cubic) | 1 2 | Antiferromagnetic | - 2 . 5 K | 47-50 | | Dy 3 Al 5 O 12 | Three (cubic) | 1 2 | Antiferromagnetic | - 1 . 85 K | 51-53 | | CoRb 3 Cl 5 | Three (cubic) | 1 2 | Antiferromagnetic | - 0 . 511 K | 54,55 | | FeF 2 | Three (cubic) | 2 | Antiferromagnetic | - 2 . 69 K | 56-59 | ## 2.2. Nuclear Magnetic Resonance In condensed matter physics, Nuclear Magnetic Resonance (NMR) has been used for decision of the structure of organic compounds and for analysis of the state in materials by using resonance induced by electromagnetic wave. The NMR can create the Ising model with transverse fields, which is expected to become an element of quantum information processing. In this processing, we use molecules where the coherence times are long compared with typical gate operations. Actually a couple of molecules which have nuclear spins were used for demonstration of quantum computing. 60-75 In this section we explain how to create the Ising model by NMR. The setup of the NMR spectrometer as a tool of quantum computing is as follows. We first put molecules which contain nuclear spins under the strong magnetic field B 0 . Next we apply radio frequency ω (rf) magnetic field which is perpendicular to the strong magnetic field B 0 . For simplicity, we here consider a molecule which contains two spins. We also assume that the considered molecule can be well described by the Heisenberg Hamiltonian. Then the Hamiltonian of this system is given by $$\begin{array} { r l } & { \hat { \mathcal { H } } = \hat { \mathcal { H } } _ { m o l } + \hat { \mathcal { H } } _ { 1 } ^ { ( r f ) } + \hat { \mathcal { H } } _ { 2 } ^ { ( r f ) } , } & { ( 6 ) } \end{array}$$ where ˆ H mol , ˆ H (rf) 1 , and ˆ H (rf) 2 are defined by $$\begin{array} { r l r } & { \hat { \mathcal { H } } _ { m o l } \colon = - h _ { 1 } \hat { \sigma } _ { 1 } ^ { z } - h _ { 2 } \hat { \sigma } _ { 2 } ^ { z } - J ( \hat { \sigma } _ { 1 } ^ { x } \cdot \hat { \sigma } _ { 2 } ^ { x } + \hat { \sigma } _ { 1 } ^ { y } \cdot \hat { \sigma } _ { 2 } ^ { y } + \hat { \sigma } _ { 1 } ^ { z } \cdot \hat { \sigma } _ { 2 } ^ { z } ) , } & { ( 7 ) } \\ & { \hat { \mathcal { U } } ^ { ( r f ) } \cdot } & { \Gamma _ { - + } ( \mu ^ { x } t _ { + } ( \hat { \sigma } _ { 1 } ^ { x } + \hat { \sigma } _ { 2 } ^ { y } ) , } & { ( 9 ) } \end{array}$$ $$\hat { \mathcal { H } } _ { 2 } ^ { ( r f ) } \colon = - \Gamma _ { 2 } \cos ( \omega ^ { ( r f ) } t - \phi _ { 2 } ) ( \gamma ^ { \prime - 1 } \hat { \sigma } _ { 1 } ^ { x } + \hat { \sigma } _ { 2 } ^ { x } ) , & & ( 9 )$$ $$\begin{array} { r l } & { \hat { \mathcal { H } } _ { 1 } ^ { ( r f ) } \colon = - \Gamma _ { 1 } \cos ( \omega ^ { ( r f ) } t - \phi _ { 1 } ) ( \hat { \sigma } _ { 1 } ^ { x } + \gamma ^ { \prime } \hat { \sigma } _ { 2 } ^ { x } ) , } \\ & { \hat { \mathcal { U } } ^ { ( r f ) } \colon \Gamma _ { 1 } \cos ( \omega ^ { ( r f ) } t - \phi _ { 1 } ) ( \gamma ^ { \prime - 1 } \hat { \mu } _ { 2 } ^ { x } + \hat { \mu } _ { 3 } ^ { x } ) . } \end{array}$$ respectively. We take the natural unit in which /planckover2pi1 = 1. The values of φ 1 and φ 2 are the phases at the time t = 0 of the first spin and that of the second spin, respectively. The quantities of h i are defined by h i := γ i B 0 , where γ i denotes the gyromagnetic ratio of the i -th spin ( i = 1 , 2). The values of h 1 and h 2 represent energy differences between |↑〉 and |↓〉 of the first spin and the second spin, respectively. The coefficients Γ 1 and Γ 2 in ˆ H (rf) 1 and ˆ H (rf) 2 are the effective amplitudes of the ac magnetic field, whose definitions are Γ i := γ i B ac , where B ac is amplitude of the ac magnetic field. The value of γ ′ is defined by the ratio of the gyromagnetic ratios γ ′ := γ 2 /γ 1 . We define the following unitary transformation: $$\hat { U } ^ { ( R ) } \colon = e ^ { - i h _ { 1 } \hat { \sigma } _ { 1 } ^ { z } t } \cdot e ^ { - i h _ { 2 } \hat { \sigma } _ { 2 } ^ { z } t } . & & ( 1 0 )$$ We can change from the laboratory frame to a frame rotating with h i around the z -axis by using the above unitary transformation. The dynamics of a density matrix can be calculated by $$i \frac { d \hat { \rho } } { d t } = [ \hat { \mathcal { H } } , \hat { \rho } ] . <text><loc_465><loc_48><loc_498><loc_100>(11)</text>$$ The density matrix on the rotating frame is given by $$\hat { \rho } ^ { ( R ) } \colon = \hat { U } ^ { ( R ) } \hat { \rho } \hat { U } ^ { ( R ) \dagger } . \quad ( 1 2 )$$ To be the same form as Eq. (11) on the rotating frame, the Hamiltonian on the rotating frame should be $$\hat { \mathcal { H } } ^ { ( R ) } = \hat { U } ^ { ( R ) } \hat { \mathcal { H } } \hat { U } ^ { ( R ) \dagger } - i \hat { U } ^ { ( R ) } \frac { d \hat { U } ^ { ( R ) \dagger } } { d t } .$$ Here we decompose the Hamiltonian on the rotating frame as $$\begin{array} { r l } & { \hat { \mathcal { H } } ^ { ( R ) } = \hat { \mathcal { H } } _ { m o l } ^ { ( R ) } + \hat { \mathcal { H } } _ { 1 } ^ { ( R ) ( r f ) } + \hat { \mathcal { H } } _ { 2 } ^ { ( R ) ( r f ) } , } & { ( 1 4 ) } \end{array}$$ where the three terms are defined by $$\begin{array} { r l } & { \hat { \mathcal { H } } _ { m o l } ^ { ( R ) } \colon = \hat { U } ^ { ( R ) } \hat { \mathcal { H } } _ { m o l } \hat { U } ^ { ( R ) \dagger } - i \hat { U } ^ { ( R ) } \frac { d \hat { U } ^ { ( R ) \dagger } } { d t } , } & { ( 1 5 ) } \\ & { \hat { \mathcal { H } } _ { m o l } ^ { ( R ) } \colon = \hat { U } ^ { ( R ) } \hat { \mathcal { H } } _ { m o l } \hat { U } ^ { ( R ) \dagger } - i \hat { U } ^ { ( R ) } \frac { d \hat { U } ^ { ( R ) \dagger } } { d t } , } & { ( 2 5 ) } \end{array}$$ $$\begin{array} { r l } & { \hat { \mathcal { H } } _ { 1 } ^ { ( R ) ( r f ) } \colon = \hat { U } ^ { ( R ) } \hat { \mathcal { H } } _ { 1 } ^ { ( r f ) } \hat { U } ^ { ( R ) \dagger } , } & { ( 1 6 ) } \\ & { \hat { \mathcal { U } } ^ { ( R ) ( r f ) } \colon = \hat { U } ^ { ( R ) } \hat { \mathcal { U } } ^ { ( r f ) } \hat { U } ^ { ( R ) \dagger } . } & { ( 1 7 ) } \end{array}$$ $$\begin{array} { r l } & { \hat { \mathcal { H } } _ { 2 } ^ { ( R ) ( r f ) } \colon = \hat { U } ^ { ( R ) } \hat { \mathcal { H } } _ { 2 } ^ { ( r f ) } \hat { U } ^ { ( R ) \dagger } . } & { ( 1 7 ) } \end{array}$$ The intramolecular magnetic interaction Hamiltonian on the rotating frame ˆ H (R) mol can be calculated as $$\hat { \mathcal { H } } _ { m o l } ^ { ( R ) } = J \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & e ^ { i ( h _ { 2 } - h _ { 1 } ) t } & 0 \\ 0 & e ^ { - i ( h _ { 2 } - h _ { 1 } ) t } & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} - J \hat { \sigma } _ { 1 } ^ { z } \hat { \sigma } _ { 2 } ^ { z } \simeq - J \hat { \sigma } _ { 1 } ^ { z } \hat { \sigma } _ { 2 } ^ { z } .$$ The approximation is valid when | h 2 -h 1 | τ /greatermuch 1, where τ is a characteristic time scale since the exponential terms are averaged to vanish. The radio frequency magnetic field Hamiltonian on the rotating frame ˆ H (R)(rf) 1 under the resonance condition ω (rf) = h i can be calculated as $$\begin{array} { r l } & { t h e r o n a l i c e c o n d u m \omega ^ { \prime } = h _ { i } \, c a l l u p e d a s } \\ & { \quad \hat { H } _ { 1 } ^ { ( R ) ( f r ) } = - \Gamma _ { 1 } \left [ \begin{array} { l l l } { 0 } & { 0 } & { e ^ { - i \phi _ { 1 } } } & { 0 } \\ { 0 } & { 0 } & { 0 } & { e ^ { - i \phi _ { 1 } } } \\ { e ^ { i \phi _ { 1 } } } & { 0 } & { 0 } & { 0 } \\ { 0 } & { e ^ { i \phi _ { 1 } } } & { 0 } & { 0 } \end{array} \right ] + \gamma ^ { \prime } \left ( \begin{array} { l l l } { 0 } & { a _ { - - } } & { 0 } & { 0 } \\ { a _ { + + } } & { 0 } & { 0 } & { 0 } \\ { 0 } & { 0 } & { 0 } & { a _ { - - } } \end{array} \right ) \right ] , } \end{array}$$ $$\hat { \mathcal { H } } _ { 1 } ^ { ( R ) ( r f ) } = - \Gamma _ { 1 } ( \cos \phi _ { 1 } \hat { \sigma } _ { 1 } ^ { x } + \sin \phi _ { 1 } \hat { \sigma } _ { 1 } ^ { y } ) .$$ where a --:= e -i ( h 2 -h 1 ) t + φ 1 +e -i ( h 1 + h 2 ) t -φ 1 and a ++ := e i ( h 2 -h 1 ) t + φ 1 + e i ( h 1 + h 2 ) t -φ 1 . The second term of ˆ H (R)(rf) 1 vanishes when | h 1 + h 2 | τ, | h 2 -h 1 | τ /greatermuch 1. Then under these conditions, the Hamiltonian becomes In the same way, the Hamiltonian ˆ H (R)(rf) 2 can be calculated as $$\begin{array} { r l } & { \hat { \mathcal { H } } _ { 2 } ^ { ( R ) ( r f ) } = - \Gamma _ { 2 } ( \cos \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { x } + \sin \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { y } ) . \quad ( 2 0 ) } \end{array}$$ By taking the rotation operators on the individual sites, we can rewrite the Hamiltonians ˆ H (R)(rf) 1 and ˆ H (R)(rf) 2 by only the x -component of the Pauli matrix: $$\begin{array} { r l } & { e ^ { i \phi _ { 1 } \hat { \sigma } _ { 1 } ^ { z } } \hat { H } _ { 1 } ^ { ( R ) ( r f ) } e ^ { - i \phi _ { 1 } \hat { \sigma } _ { 1 } ^ { z } } = - \Gamma _ { 1 } \hat { \sigma } _ { 1 } ^ { x } , \quad ( 2 1 ) } \\ & { i \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { z } \hat { u } ^ { ( R ) ( r f ) } e ^ { - i \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { z } } = \Gamma _ { 1 } \hat { \sigma } _ { 2 } ^ { x } , \quad ( 2 3 ) } \end{array}$$ $$e ^ { i \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { z } } \hat { \mathcal { H } } _ { 2 } ^ { ( R ) ( r f ) } e ^ { - i \phi _ { 2 } \hat { \sigma } _ { 2 } ^ { z } } = - \Gamma _ { 2 } \hat { \sigma } _ { 2 } ^ { x } . \quad & & ( 2 2 )$$ Then, the total Hamiltonian can be represented by the Ising model with site-dependent transverse fields: $$\hat { \mathcal { H } } ^ { ( R ) } = - J \hat { \sigma } _ { 1 } ^ { z } \hat { \sigma } _ { 2 } ^ { z } - \Gamma _ { 1 } \hat { \sigma } _ { 1 } ^ { x } - \Gamma _ { 2 } \hat { \sigma } _ { 2 } ^ { x } . & & ( 2 3 )$$ It should be noted that the above procedure is not restricted for two spin system. Then, the NMR technique can be create the Ising model with sitedependent transverse fields in general. ## 3. Implementation Methods of Quantum Annealing As stated in Sec. 1, the quantum annealing method is expected to be a powerful tool to obtain the best solution of optimization problems in a generic way. The quantum annealing methods can be categorized according to how to treat time-development. One is a stochastic method such as the Monte Carlo method which will be shown in Sec. 3.1. Other is a deterministic method such as mean-field type method and real-time dynamics. We will explain the mean-field type method and the method based on real-time dynamics in Secs. 3.2 and 3.3. Although in the Monte Carlo method and the mean-field type method, we introduce time-development in an artificial way, the merit of these methods is to be able to treat large-scale systems. The methods based on the Schr¨ odinger equation can follow up real-time dynamics which occurs in real experimental systems. However, these methods can be used for very small systems and/or limited lattice geometries because of limited computer resources and characters of algorithms. Each method has strengths and limitations based on its individuality. Then when we use the quantum annealing, we have to choose implementation methods according to what we want to know. In this section, we explain three types of theoretical methods for the quantum annealing and some experimental results which relate to the quantum annealing. ## 3.1. Monte Carlo Method In this section we review the Monte Carlo method as an implementation method of the quantum annealing. In physics, the Monte Carlo method is widely adopted for analysis of equilibrium properties of strongly correlated systems such as spin systems, electric systems, and bosonic systems. Originally the Monte Carlo method is used in order to calculate integrated value of given function. The simplest example is 'calculation of π '. Suppose we consider a square in which -1 ≤ x, y ≤ 1 and a circle whose radius is unity and center is ( x, y ) = (0 , 0). We generate pair of uniform random numbers ( -1 ≤ x i , y i ≤ 1) many times and calculate the following quantity: $$\frac { \text {number of steps when $\sqrt{x_{i}^{2}}+y_{i}^{2}\leq 1$ is satisfied} } { \text {number of steps} } .$$ Hereafter we refer to the denominator as Monte Carlo step. The quantity should converge to π/ 4 in the limit of infinite Monte Carlo step. This is a pedagogical example of the Monte Carlo method. We first explain how to implement and theoretical background of the Monte Carlo method which is used in physics. In equilibrium statistical physics, we would like to know the equilibrium value at given temperature T . The equilibrium value of the physical quantity which is represented by the operator O is defined as $$\langle O \rangle _ { T } ^ { ( e q ) } \colon = \frac { T r \, O e ^ { - \beta \mathcal { H } } } { T r \, e ^ { - \beta \mathcal { H } } } , \quad \quad ( 2 5 )$$ where Tr means the trace of matrix and β denotes the inverse temperature β = ( k B T ) -1 . Hereafter we set the Boltzmann constant k B to be unity. For small systems, we can obtain the equilibrium value by taking sum analytically, on the contrary, it is difficult to obtain the equilibrium value for large systems except few solvable models. Then in order to evaluate equilibrium value of the physical quantity, we often use the Monte Carlo method. We consider the Ising model given by $$\mathcal { H } _ { I s i n g } = - \sum _ { \langle i , j \rangle } J _ { i j } \sigma _ { i } ^ { z } \sigma _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } \sigma _ { i } ^ { z } , \quad \sigma _ { i } ^ { z } = \pm 1 .$$ The Ising model without transverse field can be expressed as a diagonal matrix by using 'trivial' bit representation |↑〉 and |↓〉 which were introduced in Sec. 2. Then, in this case, we can easily calculate the eigenenergy once the eigenstate is specified. We can use the Monte Carlo method for obtaining the equilibrium value defined by Eq. (25) as well as the calculation of π : $$\frac { \sum _ { \Sigma } O ( \Sigma ) e ^ { - \beta E ( \Sigma ) } } { \sum _ { \Sigma } e ^ { - \beta E ( \Sigma ) } } \to \langle O \rangle _ { T } ^ { ( e q ) } ,$$ where O (Σ) and E (Σ) denote the physical value of O and the eigenenergy of the eigenstate Σ, respectively. Here the eigenstate Σ is generated by uniform random number and ∑ Σ 1 is equal to Monte Carlo step. In the limit of infinite Monte Carlo step, LHS of Eq. (27) should be converge to the equilibrium value. Equilibrium statistical physics says that the probability distribution at equilibrium state can be described by the Boltzmann distribution which is proportional to e -βE (Σ) . In this case, since we know the form of the probability distribution, it is better to use the distribution function to generate a state according to the Boltzmann distribution instead of uniform random number. This scheme is called importance sampling. When we use the importance sampling, we can obtain the equilibrium value as follows: $$\frac { \sum _ { \Sigma } O ( \Sigma ) } { \sum _ { \Sigma } 1 } \rightarrow \langle O \rangle _ { T } ^ { ( e q ) } . </doctag>$$ In order to generate a state according to the Boltzmann distribution, we use the Markov chain Monte Carlo method. Let P (Σ a , t ) be the probability of the a -th state at time t . In this method, time-evolution of probability distribution is given by the master equation: /negationslash $$P ( \Sigma _ { a } , t + \Delta t ) = \left [ \sum _ { b \neq a } P ( \Sigma _ { b } , t ) w ( \Sigma _ { a } | \Sigma _ { b } ) + P ( \Sigma _ { a } , t ) w ( \Sigma _ { a } | \Sigma _ { a } ) \right ] \Delta t ,$$ /negationslash where w (Σ a | Σ b ) represents the transition probability from the b -th state to the a -th state in unit time. The transition probability w (Σ a | Σ b ) obeys $$\sum _ { \Sigma _ { a } } w ( \Sigma _ { a } | \Sigma _ { b } ) = 1 \quad ( \forall \, \Sigma _ { b } ) .$$ For convenience, let P ( t ) be a vector-representation of probability distribution { P (Σ a , t ) } . Then the master equation can be represented by $$P ( t + \Delta t ) = \mathcal { L } P ( t ) , & & ( 3 1 )$$ where L is the transition matrix whose elements are defined as $$\begin{array} { r l } & { \mathcal { L } _ { b a } \colon = w ( \Sigma _ { b } | \Sigma _ { a } ) \Delta t , } \\ & { \quad \mathcal { C } _ { b a } \colon = 1 \sum \mathcal { C } _ { b } ^ { 1 } \sum \Sigma _ { b } | \Sigma _ { a } ) \Lambda ( 3 2 ) } \end{array}$$ /negationslash $$\mathcal { L } _ { a a } \colon = 1 - \sum _ { b \neq a } \mathcal { L } _ { b a } = 1 - \sum _ { b \neq a } w ( \Sigma _ { b } | \Sigma _ { a } ) \Delta t .$$ /negationslash Here the matrix L is a non-negative matrix and does not depend on time. Then this time-evolution is the Markovian. If the transition matrix L is prepared appropriately, which satisfies the detailed balance condition and the ergordicity, we can obtain the equilibrium probability distribution in the limit of infinite Monte Carlo step regardless of choice of the initial state because of the Perron-Frobenius theorem. We can perform the Monte Carlo method easily as following process. - Step 1 We prepare a initial state arbitrary. - Step 2 We choose a spin randomly. - Step 3 We calculate the molecular field at the chosen site in Step 2. The molecular field at the chosen site i is defined as $$h _ { i } ^ { ( e f f ) } \colon = \sum _ { j } { ^ { \prime } J _ { i j } \sigma _ { j } ^ { z } + h _ { i } } ,$$ where the summation takes over the nearest neighbor sites of the i -th site. - Step 4 We flip the chosen spin in Step 2 according to a probability defined by some way. - Step 5 We continue from Step 2 to Step 4 until physical quantities such as magnetization converge. In this Monte Carlo method, we only update the chosen single spin, and thus we refer to this method as single-spin-flip method. There is an ambiguity how to define w (Σ a | Σ b ) in Step 4. Here we explain two famous choices of w (Σ a | Σ b ) as follows. Transition probability in the heat-bath method is given by $$w _ { H B } ( \sigma _ { i } ^ { z } \rightarrow - \sigma _ { i } ^ { z } ) = \frac { e ^ { - \beta h _ { i } ^ { ( e f f ) } \sigma _ { i } ^ { z } } } { 2 \cosh ( \beta h _ { i } ^ { ( e f f ) } ) } .$$ Transition probability in the Metropolis method is given by $$w _ { \text {MP} } ( \sigma ^ { z } _ { i } \to - \sigma ^ { z } _ { i } ) = \begin{cases} 1 & ( h ^ { ( \text {eff)} } _ { i } \sigma ^ { z } _ { i } < 0 ) \\ e ^ { - 2 \beta h ^ { ( \text {eff)} } _ { i } \sigma ^ { z } _ { i } } & ( h ^ { ( \text {eff} ) } _ { i } \sigma ^ { z } _ { i } \geq 0 ) \end{cases} .$$ Since both two transition probabilities satisfy the detailed balance condition, the equilibrium state can be obtained definitely in the limit of infinite Monte Carlo step a . It is important to select how to choice the transition probability since it is known that a couple of methods can sample states in an efficient fashion. 76-83 So far we considered the Monte Carlo method for systems where there is no off-diagonal matrix element. To perform the Monte Carlo method, in a precise mathematical sense, we only have to know how to choice the basis or appropriate transformation so as to diagonalize the given Hamiltonian. However, it is difficult to obtain equilibrium values of physical quantities of quantum systems, since we have to calculate the exponential of the given Hamiltonian e -β ˆ H in general. If we know all eigenvalues and the corresponding eigenvectors of the given Hamiltonian, we can easily calculate e -β ˆ H by the unitary transformation which diagonalizes the Hamiltonian ˆ H . In contrast, if we do not know all eigenvalues and eigenvectors, we have to calculate any power of the Hamiltonian ˆ H m since the matrix exponential is given by $$e ^ { \hat { A } } = \sum _ { m = 0 } ^ { \infty } \frac { 1 } { m ! } \hat { A } ^ { m } .$$ It is difficult to calculate the matrix exponential in general. Then we have to consider the following procedure in order to use the framework of the Monte Carlo method for quantum systems. In many cases, the Hamiltonian of quantum systems can be represented as $$\hat { \mathcal { H } } & = \hat { \mathcal { H } } _ { c } + \hat { \mathcal { H } } _ { q } . & ( 3 8 )$$ Hereafter we refer to ˆ H c and ˆ H q as classical Hamiltonian and quantum Hamiltonian, respectively. The classical Hamiltonian ˆ H c is a diagonal matrix. Here we assume that ˆ H q can be easily diagonalized b . This is a key of the quantum Monte Carlo method as will be shown later. Since ˆ H c and ˆ H q cannot commute in general: [ ˆ H c , ˆ H q ] = 0, then e -β ˆ H = e -β ˆ H c e -β ˆ H q . We /negationslash /negationslash a Recently, the algorithm which does not use the detailed balance condition was proposed. 76,77 It should be noted that the detailed balance condition is just a necessary condition. This novel algorithm is efficient for general spin systems. b This fact does not seem to be general. However we can prepare the matrices which can be easily diagonalized by the decomposition as ˆ H q = ∑ /lscript ˆ H ( /lscript ) q in many cases. decompose the matrix exponential by introducing large integer m , $$\begin{array} { r l } { \exp \left ( - \frac { \beta } { m } \hat { \mathcal { H } } \right ) = \exp \left [ - \frac { \beta } { m } ( \hat { \mathcal { H } } _ { c } + \hat { \mathcal { H } } _ { q } ) \right ] } \\ { = \exp \left ( - \frac { \beta } { m } \hat { \mathcal { H } } _ { c } \right ) \exp \left ( - \frac { \beta } { m } \hat { \mathcal { H } } _ { q } \right ) + \mathcal { O } \left ( \left ( \frac { \beta } { m } \right ) ^ { 2 } \right ) . \, ( 3 9 ) } \end{array}$$ This is a concrete representation of the Trotter formula. 84 From now on, we refer to m as Trotter number. By using this relation, we can perform the Monte Carlo method for quantum systems. To illustrate it, we consider the Ising model with longitudinal and transverse magnetic fields. The considered Hamiltonian is given as $$\hat { \mathcal { H } } = - \sum _ { \langle i , j \rangle } J _ { i j } \hat { \sigma } _ { i } ^ { z } \hat { \sigma } _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } ^ { z } \hat { \sigma } _ { i } ^ { z } - \Gamma \sum _ { i = 1 } ^ { N } \hat { \sigma } _ { i } ^ { x } = \hat { \mathcal { H } } _ { c } + \hat { \mathcal { H } } _ { q } ,$$ $$\hat { \mathcal { H } } _ { c } \colon = - \sum _ { ( i , j ) } J _ { i j } \hat { \sigma } _ { i } ^ { z } \hat { \sigma } _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } ^ { z } \hat { \sigma } _ { i } ^ { z } , \quad \hat { \mathcal { H } } _ { q } \colon = - \Gamma \sum _ { i = 1 } ^ { N } \hat { \sigma } _ { i } ^ { x } , \quad ( 4 1 )$$ where optimization problems often can be expressed by this classical Hamiltonian ˆ H c . The partition function of the Hamiltonian at temperature T (= β -1 ) is given by $$Z = \text {Tr} e ^ { - \beta \mathcal { H } } = \sum _ { \Sigma } \left \langle \Sigma \left | e ^ { - \beta ( \mathcal { H } _ { c } + \mathcal { H } _ { q } ) } \right | \Sigma \right \rangle . <text><loc_36><loc_499><loc_37><loc_500>###</text>$$ Using Eq. (39) we obtain $$\text {Using Eq. (39) we obtain} \\ Z \, = \, \lim _ { m \to \infty } \sum _ { \{ \Sigma _ { k } \} , \{ \Sigma ^ { \prime } _ { k } \} } \left \langle \Sigma _ { 1 } \Big | e ^ { - \beta \mathcal { H } _ { c } / m } \, \Big | \, \Sigma ^ { \prime } _ { 1 } \right \rangle \left \langle \Sigma ^ { \prime } _ { 1 } \, \Big | e ^ { - \beta \mathcal { H } _ { q } / m } \, \Big | \, \Sigma _ { 2 } \right \rangle \\ \times \left \langle \Sigma _ { 2 } \Big | e ^ { - \beta \mathcal { H } _ { c } / m } \, \Big | \, \Sigma ^ { \prime } _ { 2 } \right \rangle \left \langle \Sigma ^ { \prime } _ { 2 } \, \Big | e ^ { - \beta \mathcal { H } _ { q } / m } \, \Big | \, \Sigma _ { 3 } \right \rangle \\ \times \dots \\ \times \left \langle \Sigma _ { m } \, \Big | e ^ { - \beta \mathcal { H } _ { c } / m } \, \Big | \, \Sigma ^ { \prime } _ { m } \, \right \rangle \left \langle \Sigma ^ { \prime } _ { m } \, \Big | e ^ { - \beta \mathcal { H } _ { q } / m } \, \Big | \, \Sigma _ { 1 } \right \rangle ,$$ ∣ ∣ ∣ ∣ where | Σ k 〉 represents the direct-product space of N spins: $$\left | \Sigma _ { k } \right \rangle \colon = \left | \sigma _ { 1 , k } ^ { z } \right \rangle \otimes \left | \sigma _ { 2 , k } ^ { z } \right \rangle \otimes \cdots \left | \sigma _ { N , k } ^ { z } \right \rangle ,$$ where the first and the second subscripts of | σ z i,k 〉 indicate coordinates of the real space and the Trotter axis, respectively. Here | σ z i,k 〉 = |↑〉 or |↓〉 . Equation (42) consists of two elements 〈 Σ k | e -β ˆ H c /m | Σ ′ k 〉 and 〈 Σ ′ k | e -β ˆ H q /m | Σ k +1 〉 . Since the classical Hamiltonian ˆ H c is a diagonal matrix, the former is easily calculated: $$& \left \langle \Sigma _ { k } \Big | e ^ { - \beta \mathcal { H } _ { c } / m } \Big | \Sigma _ { k } \right \rangle \\ & = \exp \left [ \frac { \beta } { m } \left ( \sum _ { \langle i , j \rangle } J _ { i j } \sigma _ { i , k } ^ { z } \sigma _ { j , k } ^ { z } + \sum _ { i = 1 } ^ { N } h _ { i } \sigma _ { i , k } ^ { z } \right ) \right ] \prod _ { i = 1 } ^ { N } \delta ( \sigma _ { i , k } ^ { z } , \sigma _ { i , k } ^ { \prime z } ) ,$$ $$& \left \langle \Sigma _ { k } ^ { \prime } \Big | e ^ { - \beta \hat { H } _ { q } / m } \Big | \Sigma _ { k + 1 } \right \rangle \\ & = \left [ \frac { 1 } { 2 } \sinh \left ( \frac { 2 \beta \Gamma } { m } \right ) \right ] ^ { N / 2 } \exp \left [ \frac { 1 } { 2 } \ln \coth \left ( \frac { \beta \Gamma } { m } \sum _ { i = 1 } ^ { N } \sigma _ { i , k } ^ { \prime z } \sigma _ { i , k + 1 } ^ { z } \right ) \right ] .$$ where σ z i,k = ± 1. On the other hand, the latter 〈 Σ ′ k ∣ ∣ ∣ e -β ˆ H q /m ∣ ∣ ∣ Σ k +1 〉 is calculated as Then the partition function given by Eq. (43) can be represented as $$Z = \lim _ { m \to \infty } A \sum _ { \{ \sigma _ { i , k } ^ { z } = \pm 1 \} } \exp \, \left \{ \sum _ { k = 1 } ^ { m } \left [ \sum _ { \langle i , j \rangle } \left ( \frac { \beta J _ { i j } } { m } \sigma _ { i , k } ^ { z } \sigma _ { j , k } ^ { z } \right ) + \sum _ { i = 1 } ^ { N } \frac { \beta h _ { i } } { m } \sigma _ { i , k } ^ { z } \right ] \right \} ,$$ where A is just a parameter which does not affect physical quantities. It should be noted that the partition function of the d -dimensional Ising model with transverse field ˆ H is equivalent to that of the ( d +1)-dimensional Ising model without transverse field H eff which is given by $$\mathcal { H } _ { \text{eff} } = - \sum _ { ( i, j ) } \sum _ { k = 1 } ^ { m } \frac { J _ { i j } } { m } \sigma _ { i, k } ^ { z } \sigma _ { j, k } ^ { z } - \sum _ { i = 1 } ^ { N } \sum _ { k = 1 } ^ { m } \frac { h _ { i } } { m } \sigma _ { i, k } ^ { z } - \text{ coefficient of the third term of RHS is always negative} , and thus the \text{$\mathcal{C}$} = - \sum _ { ( i, j ) } \sum _ { k = 1 } ^ { m } \frac { J _ { i j } } { m } \sigma _ { i, k } ^ { z } \sigma _ { j, k } ^ { z }.$$ The coefficient of the third term of RHS is always negative, and thus the interaction along the Trotter axis is always ferromagnetic. This ferromagnetic interaction becomes strong as the value of Γ decreases. This is called the Suzuki-Trotter decomposition. 84,85 So far we explained the Monte Carlo method as a tool for obtaining the equilibrium state. However we can also use this method to investigate stochastic dynamics of strongly correlated systems, since the Monte Carlo method is originally based on the master equation. In terms of optimization problem, our purpose is to obtain the ground state of the given Hamiltonian. Then we decrease transverse field gradually and obtain a solution. There are many Monte Carlo studies in which the quantum annealing succeeds to obtain a better solution than that by the simulated annealing. 5,8-10,12,14,86 ## 3.2. Deterministic Method Based on Mean-Field Approximation In the previous section, we considered the Monte Carlo method in which time-evolution is treated as stochastic dynamics. In this section, on the other hand, we explain a deterministic method based on mean-field approximation according to Refs. [87,88]. Before we consider the quantum annealing based on the mean-field approximation, we treat the Ising model with random interactions and site-dependent longitudinal fields given by $$\mathcal { H } _ { \text {Ising} } = - \sum _ { \langle i , j \rangle } J _ { i j } \sigma _ { i } ^ { z } \sigma _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } \sigma _ { i } ^ { z } .$$ When the transverse field is absent, the molecular field of the i -th spin is given by Eq. (34). Then an equation which determines expectation value of the i -th spin at temperature T (= β -1 ) is given by $$m _ { i } ^ { z } = \frac { e ^ { \beta h _ { i } ^ { ( e f f ) } } - e ^ { - \beta h _ { i } ^ { ( e f f ) } } } { e ^ { \beta h _ { i } ^ { ( e f f ) } } + e ^ { - \beta h _ { i } ^ { ( e f f ) } } } = \tanh ( \beta h _ { i } ^ { ( e f f ) } ) .$$ In the mean-field level, we approximate that the state σ z j is equal to the expectation value m z j in Eq. (34), and we obtain $$m _ { i } ^ { z } = \tanh \left [ \beta \left ( \sum _ { j } ^ { \prime } J _ { i j } m _ { j } ^ { z } + h _ { i } \right ) \right ] ,$$ which is often called self-consistent equation. We can obtain equilibrium value in the mean-field level by iterating the following equation until convergence: $$m _ { i } ^ { z } ( t + 1 ) = \tanh ( \beta h _ { i } ^ { ( e f f ) } ( t ) ) , \quad h _ { i } ^ { ( e f f ) } ( t ) = \sum _ { j } { ^ { \prime } J _ { i j } m _ { j } ^ { z } ( t ) + h _ { i } . } \ \ ( 5 2 )$$ In order to judge the convergence, we introduce a distance which represents difference between the state at t -th step and that at ( t +1)-th step as follows: $$d ( t ) \colon = \frac { 1 } { N } \sum _ { i = 1 } ^ { N } | m _ { i } ^ { z } ( t + 1 ) - m _ { i } ^ { z } ( t ) | \, .$$ When the quantity d ( t ) is less than a given small value (typically ∼ 10 -8 or more smaller value), we judge that the calculation is converged. We summarize this method: - Step 1 We prepare a initial state arbitrary. - Step 2 We choose a spin randomly. - Step 3 We calculate the molecular field given by Eq. (34) at the chosen site in Step 2. - Step 4 We change the value of the chosen spin in Step 2 according to the obtained molecular field in Step 3. - Step 5 We continue from Step 2 to Step 4 until the distance d ( t ) converges to small value. The differences between the Monte Carlo method and this method are Step 4 and Step 5. We can perform the simulated annealing by decreasing temperature and using the state obtained in Step 5 as the initial state in Step 1 at the time changing temperature c . Next we explain a quantum version of this method. Here we apply transverse field as a quantum field. We consider the Hamiltonian given by $$\hat { \mathcal { H } } = - \sum _ { \langle i , j \rangle } J _ { i j } \hat { \sigma } ^ { z } _ { i } \hat { \sigma } ^ { z } _ { j } - \sum _ { i = 1 } ^ { N } h _ { i } \hat { \sigma } ^ { z } _ { i } - \Gamma \sum _ { i = 1 } ^ { N } \hat { \sigma } ^ { x } _ { i } . \hat { \mathcal { H } } & = - \sum _ { \langle i , j \rangle } J _ { i j } \hat { \sigma } ^ { z } _ { i } \hat { \sigma } ^ { z } _ { j } - \sum _ { i = 1 } ^ { N } h _ { i } \hat { \sigma } ^ { z } _ { i } - \Gamma \sum _ { i = 1 } ^ { N } \hat { \sigma } ^ { x } _ { i } . \\$$ The density matrix of the equilibrium state is $$\hat { \rho } = \frac { \exp ( - \beta \hat { H } ) } { T r \exp ( - \beta \hat { H } ) } = \frac { \sum _ { n = 1 } ^ { 2 ^ { N } } \exp ( - \beta \epsilon _ { n } ) \, | \lambda _ { n } \rangle \, \langle \lambda _ { n } | } { \sum _ { n = 1 } ^ { 2 ^ { N } } \exp ( - \beta \epsilon _ { n } ) } ,$$ where /epsilon1 n and | λ n 〉 denote the n -th eigenenergy and the corresponding eigenvector. The density matrix satisfies the variational principle that minimizes free energy: $$F = \min _ { \hat { \rho } } \left [ T r \left ( \hat { H } + \beta ^ { - 1 } \ln \hat { \rho } \right ) \hat { \rho } \right ] , <text><loc_465><loc_152><loc_500><loc_200>(56)</text>$$ where the logarithm of the matrix is defined by the series expansion as well as the definition of the matrix exponential (see Eq. (37)). Since it is difficult to obtain the density matrix, we have to consider alternative strategy as follows. c If we want to decrease temperature rapidly, we choose not so small value for judgement of convergence. A reduced density matrix is defined as $$\hat { \rho } _ { i } \colon = T r ^ { \prime } \, \hat { \rho } = \frac { 1 } { 2 } \left ( \hat { I } + m _ { i } ^ { z } \hat { \sigma } ^ { z } + m _ { i } ^ { x } \hat { \sigma } ^ { x } \right ) ,$$ where Tr ′ indicates trace over spin states except the i -th spin. The values m z i and m x i are calculated by $$m _ { i } ^ { z } = T r \left ( \hat { \sigma } _ { i } ^ { z } \hat { \rho } \right ) , \quad m _ { i } ^ { x } = T r \left ( \hat { \sigma } _ { i } ^ { x } \hat { \rho } \right ) . \quad \ \ ( 5 8 )$$ The reduced density matrix satisfies the following relations: $$T r \left ( \hat { \rho } _ { i } \right ) = 1 , \quad T r \left ( \hat { \sigma } _ { i } ^ { z } \hat { \rho } _ { i } \right ) = m _ { i } ^ { z } , \quad T r \left ( \hat { \sigma } _ { i } ^ { x } \hat { \rho } _ { i } \right ) = m _ { i } ^ { x } .$$ Here we assume that the density matrix can be represented by direct products of the reduced density matrices: $$\hat { \rho } \simeq \prod _ { i = 1 } ^ { N } \hat { \rho } _ { i } ,$$ which is mean-field approximation (in other words, decoupling approximation). Then, the free energy is expressed as $$F \simeq \min _ { \{ \hat { \rho } _ { i } \} } \mathcal { F } ( \{ \hat { \rho } _ { i } \} ) , \quad \ \ ( 6 1 )$$ $$\mathcal { F } ( \{ \hat { \rho } _ { i } \} ) = - \sum _ { ( i , j ) } J _ { i j } m _ { i } ^ { z } m _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } m _ { i } ^ { z } - \Gamma \sum _ { i = 1 } ^ { N } m _ { i } ^ { x } $$ From the variation of F ( { ˆ ρ i } ) under the normalization condition, we obtain the following relations: $$\begin{array} { r l } & { \hat { \rho } _ { i } = \frac { \exp ( - \beta \hat { \mathcal { H } } _ { i } ) } { T r \left [ \exp ( - \beta \hat { \mathcal { H } } _ { i } ) \right ] } , } \\ & { t = \left ( - h _ { i } = \sum _ { j = 1 } ^ { \prime } J _ { i j } m _ { j } ^ { z } + \Gamma \right ) } \end{array}$$ Then the reduced density matrix is represented by using the n -th ( n = 1 , 2) eigenvalues /epsilon1 ( i ) n and the corresponding eigenvectors | λ ( i ) n 〉 of ˆ H i : $$\hat { \mathcal { H } } _ { i } = \begin{pmatrix} - h _ { i } - \sum _ { j } ^ { \prime } J _ { i j } m _ { j } ^ { z } & - \Gamma \\ - \Gamma & + h _ { i } + \sum _ { j } ^ { \prime } J _ { i j } m _ { j } ^ { z } \end{pmatrix} .$$ $$\hat { \rho } _ { i } = \frac { \exp ( - \beta \epsilon ^ { ( i ) } _ { 1 } ) \, | \lambda ^ { ( i ) } _ { 1 } \rangle \, \langle \lambda ^ { ( i ) } _ { 1 } | + \exp ( - \beta \epsilon ^ { ( i ) } _ { 2 } ) \, | \lambda ^ { ( i ) } _ { 2 } \rangle \, \langle \lambda ^ { ( i ) } _ { 2 } | } { \exp ( - \beta \epsilon ^ { ( i ) } _ { 1 } ) + \exp ( - \beta \epsilon ^ { ( i ) } _ { 2 } ) } .$$ We can also obtain the equilibrium values of physical quantities as well as the case for Γ = 0: $$m _ { i } ^ { z } ( t + 1 ) = T r ( \hat { \sigma } _ { i } ^ { z } \hat { \rho } _ { i } ( t ) ) , \quad m _ { i } ^ { x } ( t + 1 ) = T r ( \hat { \sigma } _ { i } ^ { x } \hat { \rho } _ { i } ( t ) ) , \quad ( 6 6 )$$ $$\begin{array} { r l } & { \hat { \rho } _ { i } ( t ) = \frac { \exp ( - \beta \mathcal { H } _ { i } ( t ) ) } { T r \exp ( - \beta \mathcal { H } _ { i } ( t ) ) } , } \\ & { t = \left ( - h _ { i } - \sum _ { j = 1 } ^ { \prime } I _ { i j } m _ { z } ^ { z } ( t ) \right ) + \frac { - \Gamma } { } . } \end{array}$$ We continue the above self-consistent equation until the following distance converges: $$\hat { \mathcal { H } } _ { i } ( t ) = \left ( - h _ { i } - \sum _ { j } ^ { \prime } J _ { i j } m _ { j } ^ { z } ( t ) - \Gamma \\ - \Gamma + h _ { i } + \sum _ { j } ^ { \prime } J _ { i j } m _ { j } ^ { z } ( t ) \right ) .$$ $$d ( t ) \coloneqq \frac { 1 } { 2 N } \sum _ { i = 1 } ^ { N } \left ( | m _ { i } ^ { z } ( t + 1 ) - m _ { i } ^ { z } ( t ) | + | m _ { i } ^ { x } ( t + 1 ) - m _ { i } ^ { x } ( t ) | \right ) .$$ If the temperature is zero, the reduced density matrix should be $$\hat { \rho } _ { i } = | \lambda ^ { ( i ) } _ { 1 } \rangle \, \langle \lambda ^ { ( i ) } _ { 1 } | \, ,$$ where we consider the case for /epsilon1 ( i ) 1 < /epsilon1 ( i ) 2 . Note that if and only if -h i -∑ ′ j J ij m z j = Γ = 0, /epsilon1 ( i ) 1 = /epsilon1 ( i ) 2 is satisfied. Then if we perform the quantum annealing at T = 0, we have to know only the ground state of the local Hamiltonian ˆ H i . The procedure is the same as the case for finite temperature. By using the method, we can obtain a better solution than that obtained by the simulated annealing for some optimization problems. Recently, other type of implementation method based on mean-field approximation was proposed. 13 The method is a quantum version of the variational Bayes inference. 89 We can also obtain a better solution than the conventional variational Bayes inference. ## 3.3. Real-Time Dynamics In Sec. 3.1 and Sec. 3.2, we considered artificial time-development rules such as the Markov chain Monte Carlo method and mean-field dynamics. In this section, we explain real-time dynamics which is expressed by the time-dependent Schr¨ odinger equation: $$i \frac { \partial } { \partial t } \left | \psi ( t ) \right \rangle = \hat { \mathcal { H } } ( t ) \left | \psi ( t ) \right \rangle , \quad \, ( 7 1 )$$ where ˆ H ( t ) and | ψ ( t ) 〉 denote the time-dependent Hamiltonian and the wave function at time t , respectively. The solution of this equation is given by If we use the time-dependent Hamiltonian including time-dependent quantum field, we can perform the quantum annealing by decreasing the quantum field gradually. To obtain the solution, it is necessary to decide the initial state for Eq. (72). Since our purpose is to obtain the ground state of the given Hamiltonian which represents the optimization problem, we have no way to know the preferable initial state that leads to the ground state definitely in the adiabatic limit. However, in general, we often use a 'trivial state' as the initial state. Actually, it goes well in many cases. For instance, when we consider the Ising model with time-dependent transverse field which is given by $$| \psi ( t ) \rangle = \exp \left [ - i \int _ { 0 } ^ { t } \hat { \mathcal { H } } ( t ^ { \prime } ) d t ^ { \prime } \right ] | \psi ( t = 0 ) \rangle \, .$$ $$\hat { \mathcal { H } } ( t ) = - \sum _ { i , j } J _ { i j } \hat { \sigma } _ { i } ^ { z } \hat { \sigma } _ { j } ^ { z } - \Gamma ( t ) \sum _ { i = 1 } ^ { N } \hat { \sigma } _ { i } ^ { x } ,$$ we set the ground state for large Γ as the initial state, hence the initial state is set as $$\left | \psi ( t = 0 ) \right \rangle = \left | \rightarrow \rightarrow \cdots \rightarrow \right \rangle ,$$ where |→〉 denotes the eigenstate of ˆ σ x : $$| \rightarrow \rangle \colon = { \frac { 1 } { \sqrt { 2 } } } ( | \uparrow \rangle + | \downarrow \rangle ) .$$ In real-time dynamics, in order to obtain the ground state by using given initial condition, it is important whether there is level crossing. If there is no level crossing, the system can necessarily reach the ground state by the quantum annealing in the adiabatic limit. To show this fact, we first consider a single spin system under time-dependent longitudinal magnetic field. The Hamiltonian is given by $$\hat { \mathcal { H } } _ { \text {single} } ( t ) = - h ( t ) \hat { \sigma } ^ { z } = \left ( \begin{array} { c c } - h ( t ) & 0 \\ 0 & h ( t ) \end{array} \right ) .$$ Suppose we set | ψ (0) 〉 = |↓〉 as the initial state. For arbitrary sweeping schedules, the state at arbitrary positive t is obtained by $$| \psi ( t ) \rangle = \exp \left [ - i \int _ { 0 } ^ { t } \hat { \mathcal { H } } _ { \text {single} } ( t ^ { \prime } ) d t ^ { \prime } \right ] | \psi ( 0 ) \rangle = | \downarrow \rangle \, .$$ This is because the state |↓〉 is the eigenstate of the instantaneous Hamiltonian for arbitrary time t . In general, when there is a good quantum number Fig. 2. Eigenenergies of the single spin system under longitudinal and transverse magnetic fields for Γ = 0 . 5 (left panel) and Γ = 1 (right panel). The dotted lines represent eigenenergies for Γ = 0. <details> <summary>Image 2 Details</summary>  ### Visual Description ## Chart Type: Line Graphs ### Overview The image contains two line graphs, each displaying the relationship between a variable 'h' on the x-axis and 'ε+', 'ε-' on the y-axis. Both graphs show two solid lines and two dotted lines. The left graph shows the solid and dotted lines overlapping at h=0, while the right graph shows the solid lines diverging from the dotted lines near h=0. ### Components/Axes * **X-axis (Horizontal):** 'h', with scale markers at -5, 0, and 5. * **Y-axis (Vertical):** 'ε+', 'ε-', with scale markers at -5, 0, and 5. * **Lines:** Each graph contains two solid lines and two dotted lines. ### Detailed Analysis **Left Graph:** * **Solid Lines:** * One solid line slopes upwards from approximately (-5, -5) to (5, 5). * The other solid line slopes downwards from approximately (-5, 5) to (5, -5). * **Dotted Lines:** * One dotted line slopes upwards from approximately (-5, -5) to (5, 5), overlapping with the solid line. * The other dotted line slopes downwards from approximately (-5, 5) to (5, -5), overlapping with the solid line. * The dotted lines intersect at (0,0). **Right Graph:** * **Solid Lines:** * One solid line starts at approximately (-5, -5), curves upwards, flattens near (0, -2), then continues upwards to (5, 5). * The other solid line starts at approximately (-5, 5), curves downwards, flattens near (0, 2), then continues downwards to (5, -5). * **Dotted Lines:** * One dotted line slopes upwards from approximately (-5, -5) to (5, 5). * The other dotted line slopes downwards from approximately (-5, 5) to (5, -5). * The dotted lines intersect at (0,0). ### Key Observations * In the left graph, the solid and dotted lines are indistinguishable, suggesting they represent the same relationship under certain conditions. * In the right graph, the solid lines deviate from the dotted lines near h=0, indicating a change in the relationship between 'h' and 'ε+', 'ε-' under different conditions. * Both graphs show symmetry around the origin (0,0). ### Interpretation The graphs likely represent energy levels (ε+, ε-) as a function of some parameter 'h'. The left graph could represent a simplified model where the relationship is linear. The right graph shows a more complex relationship where the energy levels exhibit a gap or splitting near h=0. This could be due to some interaction or perturbation that is not present in the simplified model. The dotted lines likely represent a theoretical or unperturbed state, while the solid lines represent the actual or perturbed state. The divergence of the solid lines from the dotted lines in the right graph indicates the effect of this perturbation. </details> and the initial state is set to be the corresponding eigenstate, the good quantum number is conserved. Then when we perform the quantum annealing method based on the real-time dynamics, we should take care of the symmetries of the considered Hamiltonian. From this, we can obtain the ground state of the considered system in the adiabatic limit if there is no level crossing. In practice, however, since we change magnetic field with finite speed, a nonadiabatic transition is inevitable. To show this fact, we consider a single spin system under longitudinal and transverse magnetic fields. The Hamiltonian of this system is given by $$\begin{array} { r } { \hat { \mathcal { H } } _ { s i n g l e } = - h \hat { \sigma } ^ { z } - \Gamma \hat { \sigma } ^ { x } = \left ( \begin{matrix} - h & - \Gamma \\ - \Gamma & h \end{matrix} \right ) . } \end{array}$$ Since the eigenenergies are /epsilon1 ± = ± √ h 2 +Γ 2 , the smallest value of the energy difference between the ground state and the excited state is 2Γ at h = 0 as shown in Fig. 2. Suppose we consider the single spin system under time-dependent longitudinal magnetic field and fixed transverse magnetic field. The Hamiltonian is given by $$\begin{array} { r } { \hat { \mathcal { H } } _ { s i n g l e } ( t ) = - h ( t ) \hat { \sigma } ^ { z } - \Gamma \hat { \sigma } ^ { x } = \binom { - v t - \Gamma } { - \Gamma \ v t } , } \end{array}$$ where we adopt h ( t ) = vt as time-dependent longitudinal field. Here we set t = -∞ as the initial time. The initial state is set to be the ground state of the Hamiltonian at the initial time | ψ ( t = -∞ ) 〉 = |↓〉 . The ground state at t = + ∞ in the adiabatic limit is | ψ (ad) ( t = + ∞ ) 〉 = |↑〉 . Then a characteristic value which represents the nature of this dynamics is a probability of staying in the ground state at t = + ∞ which is defined by ∣ ∣ The probability of staying in the ground state should depend on the sweeping speed v and the characteristic energy gap and can be obtained by the Landau-Zener-St¨ uckelberg formula: 90-92 $$\begin{array} { r l } & { p b o b i t y o f s t a y n g i n t h e g r o u n d s t a t \, t = + \infty w h i c h i s d e n f i e n d b y } \\ & { P _ { s t a y } = \left \langle \psi ^ { ( a d ) } ( t = + \infty ) \right | \exp \left [ - i \int _ { - \infty } ^ { + \infty } \hat { H } _ { s i n g l e } ( t ^ { \prime } ) d t ^ { \prime } \right ] \left | \psi ( t = - \infty ) \right \rangle } \\ & { = \left \langle \uparrow \left | \exp \left [ - i \int _ { - \infty } ^ { + \infty } \hat { H } _ { s i n g l e } ( t ^ { \prime } ) d t ^ { \prime } \right ] \right | \downarrow \right \rangle . } \end{array}$$ $$P _ { s t a y } = 1 - \exp \left [ - \frac { \pi ( \Delta E ) ^ { 2 } } { 4 v \Delta m } \right ] ,$$ where ∆ E and ∆ m represent the energy gap at the avoided level-crossing point and the difference of the magnetizations in the adiabatic limit, respectively. In this case ∆ E = 2Γ and ∆ m = 2. In many cases, typical shape of energy structure can be approximated by simple systems such as the single spin system. Then the knowledge of the simple transitions such as the Landau-Zener-St¨ ukelberg transition and the Rosen-Zener transition 93 is useful to analyze the efficiency of the quantum annealing based on the real-time dynamics. ## 3.4. Experiments Transverse field response of the Ising model has been also established in experimentally. 94-103 A dipolar-coupled disordered magnet LiHo x Y 1 -x F 4 has easy-axis anisotropy and can be represented by the Ising model. 104,105 If we apply the longitudinal magnetic field (in other words, the magnetic field is parallel to the easy-axis), phase transition does not take place. 106,107 However, when we apply the transverse magnetic field (in other words, the magnetic field is perpendicular to the easy-axis), phase transitions occur and interesting dynamical properties shown in Ref.[ 6] were observed. In the phase diagram of this material, there are three phases. The ferromagnetic phase appears at intermediate temperature and low transverse magnetic field, whereas at low temperature and low transverse magnetic field, the glassy critical phase 108 appears. The paramagnetic phase exists at the other region. The glassy critical phase exhibits slow relaxation in general. It found that the characteristic relaxation time obtained by ac field susceptibility for quantum cooling in which we decrease transverse field after temperature is decreased is lower than that for temperature cooling case. 6 From this result, it has been expected that the effect of the quantum fluctuation helps us to obtain the best solution of the optimization problem. ## 4. Optimization Problems Optimization problems are defined by composition elements of the considered problem and real-valued cost/gain function. They are problems to obtain the best solution such that the cost/gain function takes the minimum/maximum value. In general, the number of candidate solutions increases exponentially with the number of composition elements in optimization problems. Although we can obtain the best solution by a brute force in principle, it is difficult to obtain the best solution by such a naive method in practice. Then we have to invent an innovative method for obtaining the best solution in a practical time and limited computational resource. Optimization problems can be expressed by the Ising model in many cases. Once optimization problems are mapped onto the Ising model, we can use methods that have been considered in statistical physics and computational physics such as the quantum annealing. In the anterior half of this section, we explain the correspondence between the Ising model and the traveling salesman problem which is one of famous optimization problems. We demonstrate the quantum annealing based on the quantum Monte Carlo simulation for this problem. In the posterior half, we explain the clustering problem as the example expressed by the Potts model which is a straightforward extension of the Ising model. ## 4.1. Traveling Salesman Problem In this section, we consider the traveling salesman problem which is one of famous optimization problems. The setup of the traveling salesman problem is as follows: - There are N cities. - We move from the i -th city to the j -th city where the distance between them is /lscript i,j . - We can pass through a city only once. - We return the initial city after we pass through all the cities. The traveling salesman problem is to find the minimum path under above conditions. The length of a path is given by $$L \coloneqq \sum _ { a = 1 } ^ { N } \ell _ { c _ { a } , c _ { a + 1 } } ,$$ where c a denotes the city where we pass through at the a -th step. In the traveling salesman problem, the length of a path is a cost function. From the fourth condition, the following relation should be satisfied: $$c _ { N + 1 } = c _ { 1 } . \quad ( 8 3 )$$ In terms of mathematics, the traveling salesman problem is to find { c a } N a =1 so as to minimize the path L under the above four conditions. If the number of cities N is small, it is easy to obtain the shortest path by a brute force. We can easily find the best solution of the traveling salesman problem for N = 6 shown in Fig. 3. Figure 3 (a) and (b) represent a bad solution and the best solution where the length of the path L is minimum, respectively. As the number of cities increases, the traveling salesman problem becomes seriously difficult since the number of candidate solutions is ( N -1)! / 2. Then if we want to deal with the traveling salesman problem with large N , we have to adopt smart and easy practical methods such as the simulated annealing instead of a brute force. To use the simulated annealing, we map the traveling salesman problem onto the Ising model with a couple of constraints as follows. We consider N × N two-dimensional lattice. Let n i,a be the microscopic state which represents the state at the i -th city at the a -th step. The value of n i,a can be taken either 0 or 1. If we pass through the i -th city at the Fig. 3. Traveling salesman problem for N = 6. Thin lines and thick lines denote the permitted paths and selected paths, respectively. (a) Bad solution. (b) The best solution in which the length of the path is minimum. <details> <summary>Image 3 Details</summary>  ### Visual Description ## Network Diagrams: Comparison of Two Network Configurations ### Overview The image presents two network diagrams, labeled (a) and (b), each consisting of six nodes arranged in a hexagonal configuration. The diagrams illustrate different connection patterns between these nodes, with some connections emphasized by thicker lines. ### Components/Axes * **Nodes:** Six nodes represented as white circles, arranged in a hexagonal shape in both diagrams. * **Edges:** Lines connecting the nodes, representing relationships or connections. Some edges are thin, while others are thick, indicating a distinction in the type or strength of the connection. * **Labels:** * **(a):** Located in the top-left corner of the left diagram. * **(b):** Located in the top-left corner of the right diagram. ### Detailed Analysis **Diagram (a):** * **Thick Edges:** * A thick edge connects the top-left node to the bottom-left node. * A thick edge connects the top-left node to the center node. * A thick edge connects the top-right node to the bottom-right node. * A thick edge connects the bottom-left node to the center node. * A thick edge connects the bottom-right node to the center node. * **Thin Edges:** All other possible connections between the nodes are represented by thin lines. **Diagram (b):** * **Thick Edges:** * The six nodes are connected by thick edges, forming a hexagon. * **Thin Edges:** All possible diagonal connections between the nodes are represented by thin lines. ### Key Observations * Diagram (a) shows a more clustered connection pattern, with a central node having multiple thick connections. * Diagram (b) shows a more uniform connection pattern, with all nodes having a thick connection to their immediate neighbors in the hexagon. ### Interpretation The diagrams likely represent different network topologies or configurations. Diagram (a) could represent a network with a central hub, while diagram (b) could represent a more distributed or ring-like network. The thickness of the edges could indicate the strength or importance of the connection. The diagrams could be used to compare the properties of these different network configurations, such as their resilience to node failure or their efficiency in transmitting information. </details> a -th step, n i,a is unity whereas n i,a = 0 if we do not pass through the i -th city at the a -th step. The third condition can be represented by $$locally, the city at which it is obvious that we can pass through only one city at time t = 1. (for ∀ i).$$ Furthermore, since it is obvious that we can pass through only one city at the a -th step, this constraint is expressed by $$\sum _ { i = 1 } ^ { N } n _ { i , a } = 1 \quad ( \text {for $\forall a$} ) .$$ Then the length of the path L can be rewritten as $$L = \sum _ { a = 1 } ^ { N } \sum _ { i , j } \ell _ { i , j } n _ { i , a } n _ { j , a + 1 } = \frac { 1 } { 4 } \sum _ { a = 1 } ^ { N } \sum _ { i , j } \ell _ { i , j } \sigma _ { i , a } ^ { z } \sigma _ { j , a + 1 } ^ { z } + c o n s t . ,$$ where the Ising spin variable σ z i,a = ± 1 is defined by $$\sigma _ { i , a } ^ { z } \colon = 2 n _ { i , a } - 1 . & & ( 8 7 )$$ Here we used the following relation derived by Eqs. (84) and (85): $$\sum _ { a = 1 } ^ { N } \sum _ { i , j } \ell _ { i , j } \sigma _ { i , a } ^ { z } = \text {const.}$$ Then the length of the path can be represented by the Ising spin Hamiltonian on N × N two-dimensional lattice. In general, it is difficult to obtain the stable state of the Ising model with some constraints regarded as some kind of frustration which will be shown in Sec. 5.2. ## 4.1.1. Monte Carlo Method We explain how to implement the Monte Carlo method in the traveling salesman problem. We cannot use the single-spin-flip method which was explained in Sec. 3.1 because of existence of two constraints given by Eqs. (84) and (85). The simplest way of transition between states is realized by flipping four spins simultaneously as shown in Fig. 4. Suppose we consider the case that we pass through at the i -th city at the a -th step and pass through at the j -th city at the a ′ -th step, which is described as $$\sigma _ { i , a } ^ { z } = + 1 , \, \sigma _ { j , a } ^ { z } = - 1 , \, \sigma _ { i , a ^ { \prime } } ^ { z } = - 1 , \, \sigma _ { j , a ^ { \prime } } ^ { z } = + 1 .$$ Fig. 4. The simplest way of flipping method in traveling salesman problem. Transition between the state depicted in (a) and that depicted in (b) occurs. In this case, i = 3, j = 6, a = 2, and a ′ = 5. <details> <summary>Image 4 Details</summary>  ### Visual Description ## Diagram: Graph Representations and Adjacency Matrices ### Overview The image presents two graph diagrams, labeled (a) and (b), each accompanied by an adjacency matrix representation. The graphs consist of six nodes (numbered 1 to 6) connected by edges of varying thickness, indicating different weights or strengths of connection. The adjacency matrices use filled circles to represent the presence of an edge between two nodes. The "step" label with a downward arrow suggests a sequential process or state transition. ### Components/Axes * **Nodes:** Numbered 1 through 6 in each graph. * **Edges:** Lines connecting the nodes, with varying thickness. * **Adjacency Matrices:** 6x6 grids representing the connections between nodes. Rows and columns are labeled 1 through 6. Filled circles indicate the presence of an edge. * **Step Indicator:** A downward arrow labeled "step" next to each matrix. ### Detailed Analysis **Diagram (a):** * **Graph:** * Nodes 1-6 are arranged in a hexagon. * Thick edges connect: * 1 and 2 * 2 and 3 * 3 and 4 * 4 and 5 * 5 and 6 * 6 and 1 * 2 and 5 * 1 and 4 * Thin edges connect: * 1 and 3 * 1 and 5 * 2 and 4 * 2 and 6 * 3 and 5 * 3 and 6 * 4 and 6 * **Adjacency Matrix:** * Row 1: Filled circle in column 1, column 2, and column 4. * Row 2: Filled circle in column 2, column 1, column 3, and column 5. * Row 3: Filled circle in column 3, column 2, and column 4. * Row 4: Filled circle in column 4, column 1, column 3, and column 5. * Row 5: Filled circle in column 5, column 2, column 4, and column 6. * Row 6: Filled circle in column 6, column 5. **Diagram (b):** * **Graph:** * Nodes 1-6 are arranged in a hexagon. * Thick edges connect: * 1 and 6 * 2 and 5 * 3 and 4 * Thin edges connect: * 1 and 2 * 1 and 3 * 1 and 4 * 1 and 5 * 2 and 3 * 2 and 4 * 2 and 6 * 3 and 5 * 3 and 6 * 4 and 5 * 4 and 6 * 5 and 6 * **Adjacency Matrix:** * Row 1: Filled circle in column 1 and column 6. * Row 2: Filled circle in column 2 and column 5. * Row 3: Filled circle in column 3 and column 4. * Row 4: Filled circle in column 4 and column 3. * Row 5: Filled circle in column 5 and column 2, and column 6. * Row 6: Filled circle in column 6, column 1, and column 5. ### Key Observations * The thick edges in the graphs correspond to filled circles in the adjacency matrices. * The "step" indicator suggests a transition from the graph in (a) to the graph in (b). * The adjacency matrices are symmetric, indicating undirected graphs. * The diagonal elements of the adjacency matrices are filled in row 1, row 2, row 3, row 4, and row 5 in diagram (a), and in row 1, row 2, row 3 in diagram (b), indicating self-loops. ### Interpretation The image illustrates the relationship between graph representations and adjacency matrices. The graphs visually depict connections between nodes, while the adjacency matrices provide a numerical representation of these connections. The varying thickness of the edges could represent different weights or strengths of the connections. The "step" indicator suggests a transformation or evolution of the graph structure over time or through a process. The change from graph (a) to graph (b) shows a shift in the primary connections between nodes, potentially representing a change in the system being modeled. The adjacency matrices provide a way to quantify and analyze these changes. The self-loops in diagram (a) and diagram (b) indicate that some nodes are connected to themselves. </details> The trial state generated by flipping four spins is as follows: $$\sigma _ { i , a } ^ { z } = - 1 , \, \sigma _ { j , a } ^ { z } = + 1 , \, \sigma _ { i , a ^ { \prime } } ^ { z } = + 1 , \, \sigma _ { j , a ^ { \prime } } ^ { z } = - 1 .$$ The heat-bath method and the Metropolis method can be adopted for the transition probability between the present state and the trial state. In Fig. 4, i = 3, j = 6, a = 2, and a ′ = 5. It should be noted that without loss of generality the initial condition can be set as $$\sigma _ { 1 , 1 } = + 1 , \, \sigma _ { i , 1 } = - 1 \quad ( i \neq 1 ) , & & ( 9 1 )$$ /negationslash and thus we can fix the states at the first step ( a = 1) during calculation. The number of interactions in which we try to flip all spins in each Monte Carlo step is ( N -1)( N -2) / 2. ## 4.1.2. Quantum Annealing In order to perform the quantum annealing, we introduce the transverse field as the quantum fluctuation effect as shown in Sec. 3. The quantum Hamiltonian is given by $$\hat { \mathcal { H } } = \frac { 1 } { 4 } \sum _ { a = 1 } ^ { N } \sum _ { i , j } \ell _ { i , j } \hat { \sigma } _ { i , a } ^ { z } \hat { \sigma } _ { j , a + 1 } ^ { z } - \Gamma \sum _ { a = 1 } ^ { N } \sum _ { i = 1 } ^ { N } \hat { \sigma } _ { i , a } ^ { x } ,$$ where the first-term corresponds to the length of path and the second-term denotes the transverse field. We can map this quantum Hamiltonian on N × N two-dimensional lattice onto N × N × m three-dimensional Ising model as well as the case which was considered in Sec. 3.1. The effective classical Hamiltonian derived by the Suzuki-Trotter decomposition is written as $$\mathcal { H } _ { \text {eff} } & = \frac { 1 } { 4 m } \sum _ { a = 1 } ^ { N } \sum _ { i , j } \sum _ { k = 1 } ^ { m } \ell _ { i , j } \sigma _ { i , a , k } ^ { z } \sigma _ { j , a + 1 , k } ^ { z } \\ & - \frac { 1 } { \beta } \sum _ { a = 1 } ^ { N } \sum _ { i = 1 } ^ { N } \sum _ { k = 1 } ^ { m } \frac { 1 } { 2 } \ln \coth \left ( \frac { \beta \Gamma } { m } \right ) \sigma _ { i , a , k } ^ { z } \sigma _ { i , a , k + 1 } ^ { z } , \quad \sigma _ { i , a , k } ^ { z } = \pm 1 .$$ In the quantum annealing procedure, we have to take care of the constraints given by Eqs. (84) and (85) as stated before. Then the simplest way of changing state is to flip simultaneously four spins on the same layer ( m is fixed) along the Trotter axis. ## 4.1.3. Comparison with Simulated Annealing and Quantum Annealing In order to demonstrate the comparison with the simulated annealing and the quantum annealing, we perform the Monte Carlo simulation for the traveling salesman problem. As an example, we consider N = 20 cities depicted in Fig. 5 (a). The positions of these cities were generated by pair of uniform random numbers (0 ≤ x i , y i ≤ 1). The time schedules of temperature T ( t ) for the simulated annealing and transverse field Γ( t ) for the Fig. 5. Traveling salesman problem for N = 20. (a) Positions of cities. (b) The best solution in which the length of the path is minimum. <details> <summary>Image 5 Details</summary>  ### Visual Description ## Scatter Plot and Traveling Salesman Path ### Overview The image presents two scatter plots side-by-side. The left plot, labeled "(a)", shows a random distribution of points. The right plot, labeled "(b)", displays the same points connected by lines, forming a path that appears to minimize the total distance traveled between them, resembling a solution to the Traveling Salesman Problem (TSP). ### Components/Axes Both plots share the same axes: * **x-axis:** Labeled "x", ranges from 0 to 1 in increments of 0.1. * **y-axis:** Labeled "y", ranges from 0 to 1 in increments of 0.1. * **Data Points:** Represented by gray circles. * **Path (Plot b):** Represented by black lines connecting the data points. ### Detailed Analysis **Plot (a): Scatter Plot** The points are scattered relatively randomly across the plot area. Here are approximate coordinates for some of the points: * (0, 0): One point * (0.1, 0): One point * (0.1, 0.73): One point * (0.2, 0.9): One point * (0.2, 0.6): One point * (0.3, 0.4): One point * (0.35, 0.35): One point * (0.4, 0.4): One point * (0.6, 0.25): One point * (0.7, 0.15): One point * (0.75, 0.6): One point * (0.8, 0.65): One point * (0.8, 0.5): One point * (0.9, 0.6): One point * (0.9, 0.4): One point **Plot (b): Traveling Salesman Path** The points from plot (a) are connected by a path. The path starts at approximately (0,0) and visits each point once before returning to the starting point. The path appears to be a reasonable approximation of the shortest possible route connecting all the points. Here's a breakdown of the path's segments and approximate coordinates of the points it connects: 1. (0, 0) to (0.1, 0) 2. (0.1, 0) to (0.1, 0.73) 3. (0.1, 0.73) to (0.2, 0.9) 4. (0.2, 0.9) to (0.2, 0.6) 5. (0.2, 0.6) to (0.35, 0.35) 6. (0.35, 0.35) to (0.4, 0.4) 7. (0.4, 0.4) to (0.6, 0.25) 8. (0.6, 0.25) to (0.7, 0.15) 9. (0.7, 0.15) to (0.9, 0.4) 10. (0.9, 0.4) to (0.9, 0) 11. (0.9, 0) to (0.9, 0.6) 12. (0.9, 0.6) to (0.8, 0.65) 13. (0.8, 0.65) to (0.75, 0.6) 14. (0.75, 0.6) to (0.8, 0.5) 15. (0.8, 0.5) to (0.2, 0.6) 16. (0.2, 0.6) to (0, 0) ### Key Observations * Plot (a) shows a random distribution of points. * Plot (b) demonstrates a possible solution to the Traveling Salesman Problem for the points in plot (a). * The path in plot (b) attempts to minimize the total distance traveled. ### Interpretation The image illustrates the Traveling Salesman Problem (TSP), a classic optimization problem. Plot (a) represents the set of locations (cities) that need to be visited. Plot (b) shows a possible solution to the TSP, where a path is constructed that visits each location exactly once and returns to the starting point, while attempting to minimize the total distance traveled. The TSP is NP-hard, meaning that finding the optimal solution becomes computationally expensive as the number of locations increases. The path shown in plot (b) is likely a heuristic solution, meaning it's a good approximation but not necessarily the absolute shortest path. </details> quantum annealing are defined as $$T ( t ) \colon = T _ { 0 } + T _ { 1 } \left ( 1 - \frac { t } { \tau } \right ) , & & ( 9 4 )$$ where T 0 and Γ 0 are temperature and transverse field at the final time ( t = τ ), and T 0 + T 1 and Γ 0 +Γ 1 are temperature and transverse field at the initial time ( t = 0). The value of τ -1 indicates the annealing speed, and the annealing speed becomes slow as the value of τ increases. In our simulations, we adopt T 0 = Γ 0 = 0 . 01 and T 1 = Γ 1 = 5. Furthermore, we fix the transverse field as Γ = 0 during the simulation in the simulated annealing and the temperature as T = 0 . 01 during the simulation in the quantum annealing. $$\Gamma ( t ) \coloneqq \Gamma _ { 0 } + \Gamma _ { 1 } \left ( 1 - \frac { t } { \tau } \right ) ,$$ We execute 100 independent simulations of simulated annealing based on the heat-bath type Monte Carlo method where each initial state generated by the uniform random number is different. To compare the efficiency of the simulated annealing and quantum annealing in an equitable manner, in the quantum annealing, the Trotter number is putted as m = 10, and we execute 10 independent simulations. We also calculate the minimum length of path L min ( t ) := min { L ( t ′ ) | 0 ≤ t ′ ≤ t } . It should be noted that L min ( t ) is a monotonic decreasing function. The upper panel of Fig. 6 shows the time dependence of minimum length of path L min ( t ) for various τ . From the upper panel of Fig. 6, we can see that the convergence of minimum length of path in the quantum annealing is faster than that in the simulated annealing. We also show the sweeping time τ dependence of the minimum length of path at the final state L min ( τ ) in the lower panel of Fig. 6. This figure indicates that the obtained solution in the quantum annealing is always better than that in the simulated annealing. Figure 5 (b) shows the obtained best solution in both the simulated annealing and the quantum annealing with slow schedule. In this way, we can obtain a better solution (in this case, the best solution) by both annealing methods with slow schedule. Moreover, in our calculation, the convergence of solution in the quantum annealing is faster than that in the simulated annealing, and the obtained solution in the quantum annealing is better than that in the simulated annealing regardless of sweeping time τ . Thus, we can say that the quantum annealing method is appropriate as the annealing method for the traveling salesman problem in comparison with the simulated annealing. This fact has been confirmed in some researches. 86,109 Fig. 6. (Upper panel) Time dependence of minimum length of path L min ( t ) for τ = 10, 100, and 1000 obtained by the simulated annealing (SA) and the quantum annealing (QA). (Lower panel) Sweeping-time τ dependence of minimum length of path at the final state L min ( τ ) obtained by the simulated annealing indicated by squares and the quantum annealing indicated by circles. <details> <summary>Image 6 Details</summary>  ### Visual Description ## Line Charts: Comparison of SA and QA Algorithms ### Overview The image presents a series of line charts comparing the performance of two algorithms, Simulated Annealing (SA) and Quantum Annealing (QA), across different time scales (τ = 10, 100, 1000). The charts display the minimum loss, Lmin(t), as a function of time, t, for each algorithm and time scale. A final chart compares Lmin(τ) for SA and QA across a range of τ values. ### Components/Axes **Top Six Charts (3x2 Grid):** * **Y-axis:** Lmin(t) - Minimum Loss, ranging from 4 to 8. * **X-axis:** t - Time, ranging from 1 to 1000, logarithmic scale. * **Titles:** Each chart has a title indicating the time scale (τ) and the algorithm (SA or QA). The titles are: * Top-left: τ = 10, SA * Top-right: τ = 10, QA * Middle-left: τ = 100, SA * Middle-right: τ = 100, QA * Bottom-left: τ = 1000, SA * Bottom-right: τ = 1000, QA **Bottom Chart:** * **Y-axis:** Lmin(τ) - Minimum Loss, ranging from 4 to 5. * **X-axis:** τ - Time Scale, ranging from 1 to 10000, logarithmic scale. * **Legend (Top-right):** * SA: Square symbol * QA: Circle symbol ### Detailed Analysis **Top Six Charts (SA):** * **τ = 10, SA:** Lmin(t) starts at approximately 8, decreases to around 7.5 by t=5, and then drops sharply to approximately 5 by t=10. It remains relatively constant thereafter. * **τ = 100, SA:** Lmin(t) starts at approximately 8, decreases in steps to approximately 6.5 by t=20, then drops to approximately 5.5 by t=50, and finally drops sharply to approximately 4.5 by t=100. It remains relatively constant thereafter. * **τ = 1000, SA:** Lmin(t) starts at approximately 8, decreases in steps to approximately 7 by t=20, then drops to approximately 6.5 by t=100, and finally drops sharply to approximately 4 by t=1000. It remains relatively constant thereafter. **Top Six Charts (QA):** * **τ = 10, QA:** Lmin(t) starts at approximately 5.5, decreases to approximately 4.5 by t=10, and remains relatively constant thereafter. * **τ = 100, QA:** Lmin(t) starts at approximately 5, decreases to approximately 4.2 by t=100, and remains relatively constant thereafter. * **τ = 1000, QA:** Lmin(t) starts at approximately 4.8, decreases to approximately 4.1 by t=1000, and remains relatively constant thereafter. **Bottom Chart:** * **SA (Squares):** * τ = 10: Lmin(τ) ≈ 4.8 * τ = 100: Lmin(τ) ≈ 4.3 * τ = 1000: Lmin(τ) ≈ 4.2 * τ = 10000: Lmin(τ) ≈ 4.2 * **QA (Circles):** * τ = 10: Lmin(τ) ≈ 4.4 * τ = 100: Lmin(τ) ≈ 4.2 * τ = 1000: Lmin(τ) ≈ 4.1 * τ = 10000: Lmin(τ) ≈ 4.1 ### Key Observations * For SA, the minimum loss Lmin(t) decreases more significantly as τ increases, but the time it takes to reach a stable value also increases. * For QA, the minimum loss Lmin(t) is generally lower than SA for smaller values of t. * In the bottom chart, both SA and QA show a decreasing trend in Lmin(τ) as τ increases, but the decrease becomes less pronounced at higher τ values. * QA consistently achieves a lower minimum loss Lmin(τ) than SA across all τ values in the bottom chart. ### Interpretation The data suggests that Quantum Annealing (QA) generally outperforms Simulated Annealing (SA) in terms of achieving a lower minimum loss, especially at smaller time scales (t). While SA can eventually reach a comparable minimum loss, it requires a longer time scale (τ) to do so. The bottom chart reinforces this observation, showing that QA consistently achieves a lower Lmin(τ) across a range of τ values. This indicates that QA is more efficient at finding the minimum loss within a given time frame. The diminishing returns observed at higher τ values suggest that there is a limit to how much further the minimum loss can be reduced by increasing the time scale. </details> ## 4.2. Clustering Problem In Sec. 4.1, we explained the traveling salesman problem which can be mapped onto the Ising model with some constraints. Many optimization problems can also be mapped onto the Ising model. However, there are a number of optimization problems that can be described by the other models which are straightforward extensions of the Ising model. In this section, we review the concept of clustering problem as such an example. Clustering problem is also one of important optimization problems in information science and engineering. 12-14 We need to categorize much data in the real world according to its contents in various situations. For instance, suppose we play stock market. In order to see the socioeconomic situation, we want to extract efficiently important information related to stock market from an enormous quantity of information in news sites and newspapers. In this case, it is better to categorize many articles in news sites and newspapers according to their contents. This is an example of clustering problem which is adopted for many applications in wide area of science such as cognitive science, social science, and psychology. The clustering problem is to divide the whole set into a couple of subsets. Here we refer to the subsets as 'cluster'. Figure 7 shows schematic picture of the clustering problem. Suppose we consider much data in the whole set which represents the square frame in Fig. 7 (a). The points in Fig. 7 denote individual data. In the clustering problem, our target is to find which the best division is. Figure 7 (b), (c), and (d) represent typical clustering states Σ 1 , Σ 2 , and Σ ∗ , respectively. The states Σ 1 and Σ 2 are an unstable solution and a metastable solution, respectively. The state Σ ∗ denotes the best solution of clustering problem. In order to consider how to implement the quantum annealing, the clustering problem can be described by the Potts model with random interac- Fig. 7. Schematic pictures of clustering problem. The points represent data and the square denote the whole set. (a) Data set. (b) Unstable solution Σ 1 . (c) Metastable solution Σ 2 . (d) The best solution Σ ∗ . <details> <summary>Image 7 Details</summary>  ### Visual Description ## Diagram: Point Distribution and Clustering ### Overview The image presents four diagrams (a, b, c, d) illustrating the distribution of points within a rectangular area. Diagrams (b), (c), and (d) show different ways of clustering these points using dashed-line rectangles. ### Components/Axes * **Diagrams:** (a), (b), (c), (d) * **Points:** Represented by small black circles. * **Clustering:** Indicated by dashed-line rectangles in diagrams (b), (c), and (d). ### Detailed Analysis **Diagram (a):** * Shows the initial distribution of points. * There are approximately 35 points in total. * The points are grouped into several clusters: * A cluster of 9 points in the top-left. * A cluster of 9 points in the top-right. * A cluster of 6 points in the bottom-left. * A cluster of 6 points in the bottom-right. * A linear cluster of 5 points in the center. **Diagram (b):** * Shows a clustering of the points from diagram (a). * The points are grouped into three clusters: * A cluster of 9 points in the top-left. * A cluster of 9 points in the top-right. * A cluster of 17 points in the bottom. **Diagram (c):** * Shows a different clustering of the points from diagram (a). * The points are grouped into three clusters: * A cluster of 9 points in the top-left. * A cluster of 6 points in the bottom-left. * A cluster of 11 points in the right. **Diagram (d):** * Shows a different clustering of the points from diagram (a). * The points are grouped into three clusters: * A cluster of 9 points in the top-left. * A cluster of 6 points in the bottom-left. * A cluster of 11 points in the right. ### Key Observations * Diagram (a) serves as the baseline, showing the initial point distribution. * Diagrams (b), (c), and (d) demonstrate different approaches to clustering the same set of points. * The number of points remains constant across all diagrams; only the clustering changes. ### Interpretation The diagrams illustrate different ways to group or cluster data points. Diagram (a) shows the raw data, while diagrams (b), (c), and (d) show different interpretations of the data through clustering. The choice of clustering method depends on the specific application and the desired outcome. The image demonstrates that the same data can be organized and interpreted in multiple ways. </details> tions d . The Hamiltonian of the Potts model is given by $$\mathcal { H } _ { P o t t s } = - \sum _ { i , j } J _ { i j } \delta _ { \sigma _ { i } , \sigma _ { j } } , \quad \sigma _ { i } = 1 , \cdots , Q ,$$ where the summation runs over all pairs of the i -th and j -th data. The spin variable σ i represents individual data. Here the value of Q represents the number of clusters. When σ i = σ j , the i -th and j -th data are in the same cluster. It is natural to adopt ferromagnetic/antiferromagnetic interaction between data in the same/different cluster. It should be noted that the Potts model is a straightforward extension of the Ising model since the Potts model is equivalent to the Ising model if Q = 2. Then the clustering problem is a problem to obtain the ground state of the Hamiltonian of the Potts model with given random interactions. Here we assume that the number of clusters is fixed. Next we explain how to introduce quantum field in order to perform the quantum annealing. In optimization problems which can be represented by the Ising model, we can use transverse field as the quantum fluctuation which is represented as -Γ ∑ i σ x i . However, we cannot use this transverse field -Γ ∑ i σ x i for the clustering problem directly, since the matrix which represents the state is Q × Q matrix. Thus, we generalize the x -component of the Pauli matrix of the Ising model as follows: $$\hat { \tau } ^ { x } \colon = \mathbb { E } _ { Q } - \mathbb { I } _ { Q } = \left ( \begin{array} { c c c c } 0 & - 1 - 1 & \cdots & - 1 \\ - 1 & 0 & - 1 & \cdots & - 1 \\ - 1 - 1 & 0 & \vdots & - 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ - 1 - 1 & - 1 & \cdots & 0 \end{array} \right ) ,$$ where E Q and I Q represent the Q × Q unit matrix and the Q × Q matrix whose all elements are unity. By using this generalized Pauli matrix, we can apply the quantum annealing for clustering problem. 12-14 Here we consider the following Hamiltonian: $$\hat { \mathcal { H } } = \hat { \mathcal { H } } _ { P o t t s } + \hat { \mathcal { H } } _ { q } ^ { ( P o t t s ) } , \quad \hat { \mathcal { H } } _ { q } ^ { ( P o t t s ) } \colon = - \Gamma \sum _ { i = 1 } ^ { N } \hat { \tau } _ { i } ^ { x } , \quad \quad ( 9 8 )$$ d In practice, we do not know { J ij } and have to estimate interactions when we consider the clustering problem. However, we assume the Hamiltonian for simple explanation. As shown in this section, the implementation method does not depend on the specific form of interactions. where N is the number of individual data. As well as the case for the Ising model, we can calculate the partition function of the Hamiltonian: $$\text {model, we can calculate the partition function of the Hamiltonian} \colon \\ Z _ { P o t s } \, = \, T r e ^ { - \beta \hat { H } } = \sum _ { \Sigma } \left \langle \Sigma \left | e ^ { - \beta ( \hat { H } _ { P o t s } + \hat { H } _ { q } ^ { ( P o t s ) } ) } \right | \Sigma \right \rangle \\ = \lim _ { m \to \infty } \sum _ { \{ \Sigma _ { k } \} , \{ \Sigma _ { k } ^ { \prime } \} } \left \langle \Sigma _ { 1 } \left | e ^ { - \beta \hat { H } _ { P o t s } / m } \right | \Sigma _ { 1 } ^ { \prime } \right \rangle \left \langle \Sigma _ { 1 } ^ { \prime } \left | e ^ { - \beta \hat { H } _ { q } ^ { ( P o t s ) } / m } \right | \Sigma _ { 2 } \right \rangle \\ \times \left \langle \Sigma _ { 2 } \left | e ^ { - \beta \hat { H } _ { P o t s } / m } \right | \Sigma _ { 2 } ^ { \prime } \right \rangle \left \langle \Sigma _ { 2 } ^ { \prime } \right | e ^ { - \beta \hat { H } _ { q } ^ { ( P o t s ) } / m } \right | \Sigma _ { 3 } \right \rangle \\ \times \left \langle \Sigma _ { m } \left | e ^ { - \beta \hat { H } _ { P o t s } / m } \right | \Sigma _ { m } ^ { \prime } \right \rangle \left \langle \Sigma _ { m } ^ { \prime } \right | e ^ { - \beta \hat { H } _ { q } ^ { ( P o t s ) } / m } \right | \Sigma _ { 1 } \right \rangle , \quad ( 9 9 ) \\$$ ∣ ∣ ∣ ∣ where | Σ k 〉 represents the direct-product space of N spins: $$\begin{array} { r l r } & { | \Sigma _ { k } \rangle = | \sigma _ { 1 , k } \rangle \otimes | \sigma _ { 2 , k } \rangle \otimes \cdots | \sigma _ { N , k } \rangle . } & { ( 1 0 0 ) } \\ & { \quad } & { ( n o t s ) } \end{array}$$ There are two elements 〈 Σ k | e -β ˆ H Potts /m | Σ ′ k 〉 and 〈 Σ ′ k | e -β ˆ H (Potts) q /m | Σ k +1 〉 . These factors are calculated as follows: $$\begin{array} { r l } & { \left \langle \Sigma _ { k } \left | e ^ { - \beta \mathcal { H } _ { P o t t s } / m } \right | \Sigma _ { k } ^ { \prime } \right \rangle = \exp \left ( \frac { \beta } { m } \sum _ { i , j } J _ { i j } \delta _ { \sigma _ { i , k } , \sigma _ { j , k } } \right ) \prod _ { i = 1 } ^ { N } \delta _ { \sigma _ { i , k } , \sigma _ { i , k } ^ { \prime } } , \quad ( 1 0 1 ) } \\ & { \left \langle \Sigma _ { k } ^ { \prime } \left | e ^ { - \beta \mathcal { H } _ { q } ^ { ( P o t t s ) / m } } \right | \Sigma _ { k + 1 } ^ { \prime } \right \rangle = \prod _ { i , j } ^ { N } \left [ e ^ { - \frac { \beta \Gamma } { m } } \delta _ { \sigma _ { i , k } ^ { \prime } , \sigma _ { i , k + 1 } } ^ { \prime } + \frac { 1 } { \Omega } \left ( e ^ { - \frac { \beta \Gamma } { m } ( 1 - Q ) } - 1 \right ) \right ] . } \end{array}$$ $$\left \langle \Sigma _ { k } ^ { \prime } \left | e ^ { - \beta \hat { H } _ { q } ^ { ( P o t t s ) } / m } \right | \Sigma _ { k + 1 } \right \rangle = \prod _ { i = 1 } ^ { N } \left [ e ^ { - \frac { \beta \Gamma } { m } \delta _ { \sigma _ { i , k } ^ { \prime } \sigma _ { i , k + 1 } } + \frac { 1 } { Q } \left ( e ^ { - \frac { \beta \Gamma } { m } ( 1 - Q ) } - 1 \right ) } \right ] .$$ By using the above expressions, we can perform the quantum Monte Carlo simulation as well as the Ising model with transverse field. If the spin variable is not S = 1 / 2 Ising spin as in the case just described, we can implement the quantum annealing by considering appropriate quantum field. There are some studies that the quantum annealing succeeds to obtain the better solution than the simulated annealing for clustering problems. 12-14 ## 5. Relationship between Quantum Annealing and Statistical Physics In the preceding sections we explained the Ising model, a couple of implementation methods of the quantum annealing, and the optimization problems. There are a couple of studies that clarify the efficiency and feature of the quantum annealing in terms of statistical physics. In this section we take two examples which display relationship between quantum annealing and statistical physics focusing on the thermal fluctuation effect and the quantum fluctuation effect for ordering phenomena. In the first half, we review the Kibble-Zurek mechanism which characterizes the efficiency of the quantum annealing for systems where a second-order phase transition occurs comparing with the efficiency of the simulated annealing. In the last half, we show similarities and differences between thermal fluctuation and quantum fluctuation for frustrated Ising spin systems. ## 5.1. Kibble-Zurek Mechanism In statistical physics, it has been an important topic to investigate the ordering process in systems where a phase transition takes place. 110-116 Especially, dynamical properties during changing control variables such as temperature and external fields are interesting. 111,113,115 Recently, the Kibble-Zurek mechanism has been drawing attention not only in statistical physics and condensed matter physics but also for the quantum annealing. In this section, we explain the Kibble-Zurek mechanism relating to a dynamics which passes across a second-order phase transition point. The Kibble-Zurek mechanism can make clear what happens in systems where the second-order phase transition occurs during the simulated annealing and the quantum annealing from a viewpoint of statistical physics. Before we consider the efficiency of the quantum annealing comparing with the simulated annealing by using the Kibble-Zurek mechanism, we show the general feature of the Kibble-Zurek mechanism. As an example, we consider the Kibble-Zurek mechanism in the ferromagnetic system where the second-order phase transition occurs at finite temperature. At the second-order phase transition point, the correlation length diverges in the equilibrium state, and thus the relaxation time should be infinite. Hence, the system cannot reach the equilibrium state, when we decrease temperature to the transition temperature with finite speed. Furthermore, since the relaxation time is long around the transition temperature, it is difficult to equilibrate the system. Here, we assume that growth of correlation length stops at the temperature where the system is less able to reach the equilibrium state. If we decrease temperature slow enough, the system can reach the equilibrium state even near the transition point. Thus, it is expected that the value of stopped correlation length because of the long relaxation time depends on the annealing speed. As we will see below, the value of stopped correlation length can be scaled by the annealing speed. To consider the second-order phase transition at finite temperature in the ferromagnetic systems, we define the dimensionless temperature g as $$g \colon = \frac { T - T _ { c } } { T _ { c } } , \quad ( 1 0 3 )$$ where T c is the phase transition temperature. When the absolute value of g is small, it is believed that the scaling ansatz is valid. By the scaling ansatz, the temperature-dependent correlation length ξ ( g ) is given as 117 $$\xi ( g ) \, \infty \, | g | ^ { - \nu } , & & ( 1 0 4 )$$ where ν is one of the critical exponents. Moreover, the relaxation time τ rel is scaled by the following relation: 117 $$\tau _ { r e l } ( g ) \, \infty \, [ \xi ( g ) ] ^ { z } \, \infty \, | g | ^ { - z \nu } , & & ( 1 0 5 )$$ where z is the dynamical critical exponent. Here, we decrease the temperature T ( t ) against the time t as following schedule: $$T ( t ) = T _ { c } \left ( 1 - \frac { t } { \tau _ { Q } } \right ) \quad ( - \infty < t \leq \tau _ { Q } ) .$$ The value of τ -1 Q corresponds to the annealing speed. When the value of τ Q is large/small, the system is annealed to low temperature slowly/quickly. At t = 0, the temperature is the phase transition temperature ( T (0) = T c ), and the temperature is zero ( T ( τ Q ) = 0) at t = τ Q . From Eq. (106), the dimensionless temperature g becomes the time-dependent function as follows: $$g ( t ) = \frac { T ( t ) - T _ { c } } { T _ { c } } = - \frac { t } { \tau _ { Q } } , \quad \ \ ( 1 0 7 )$$ In the Kibble-Zurek mechanism, we assume the following situation: $$\begin{cases} \tau _ { r e l } ( g ( t ) ) < | t | & \colon s y s t e m \ c a n \ r e a c h \ e q u i b l i r i u m \ s t a t e \\ \tau _ { r e l } ( g ( t ) ) > | t | & \colon s y s t e m \ c a n n o t \ r e a c h \ e q u i b l i r i u m \ s t a t e \end{cases} , \quad ( 1 0 8 )$$ where | t | is a remaining time to transition temperature. That is, when a remaining time | t | is longer/shorter than the relaxation time τ rel ( g ( t )), the system can/cannot reach the equilibrium state. Note that the value of considered t should be negative since the relaxation time diverges before the temperature reaches the transition temperature ( t = 0). From this assumption, the time ˜ t at which the system is less able to reach the equilibrium state is defined by following relation: $$\tau _ { r e l } ( g ( \tilde { t } ) ) = | \tilde { t } | . & & ( 1 0 9 )$$ Furthermore, since we have assumed that the growth of correlation length stops at t = ˜ t , the value of correlation length is always ξ ( g ( ˜ t )) below T ( ˜ t ) as shown in Fig. 8. Moreover, the dimensionless temperature at ˜ t is expressed as $$g ( \tilde { t } ) = \frac { | \tilde { t } | } { \tau _ { Q } } = \frac { \tau _ { r e l } ( g ( \tilde { t } ) ) } { \tau _ { Q } } \, \infty \, \frac { | g ( \tilde { t } ) | ^ { - z \nu } } { \tau _ { Q } } . \quad ( 1 1 0 )$$ From this relation, g ( ˜ t ) is scaled by the annealing speed, and from Eqs. (104) and (110), the correlation length at t = ˜ t is scaled as follows: $$g ( \tilde { t } ) \, \infty \, \tau _ { Q } ^ { - \frac { 1 } { 1 + z \nu } } , \quad \xi ( g ( \tilde { t } ) ) \, \infty \, \tau _ { Q } ^ { \frac { \nu } { 1 + z \nu } } . \quad ( 1 1 1 )$$ Furthermore, the density of domain wall n ( t ) is written as $$\begin{array} { r l r } & { n ( t ) \infty \xi ( g ( t ) ) ^ { - d } , } & { ( 1 1 2 ) } \\ & { \quad } & { ( 1 1 2 ) } \end{array}$$ where d is the spatial dimension, and n ( ˜ t ) at t = ˜ t is scaled as follows: $$\begin{array} { r l } & { n ( \tilde { t } ) \varpropto \tau _ { Q } ^ { - \frac { d \nu } { 1 + z \nu } } . } \\ & { ( 1 1 3 ) } \end{array}$$ For instance, in the ferromagnetic Ising model on two-dimensional lattice ( d = 2, ν = 1) when we adopt the Monte Carlo dynamics based on the single-spin-flip method ( z = 2 . 132), 118 the correlation length and the density of domain wall at t = ˜ t are naively obtained as $$\begin{array} { r l r } & { \xi ( g ( \tilde { t } ) ) \, \infty \, \tau _ { Q } ^ { 0 . 3 1 9 } , \quad n ( \tilde { t } ) \, \infty \, \tau _ { Q } ^ { - 0 . 6 3 9 } . } & { ( 1 1 4 ) } \\ \end{array}$$ In this way, in the dynamics which passes across the second-order phase transition point at finite temperature, the correlation length and the density of domain wall (topological defect) are scaled by the annealing speed. Fig. 8. Schematic of the annealing speed dependence of correlation length ξ ( g ( t )). τ -1 Q is annealing speed and τ Q 1 > τ Q 2 > τ Q 3 . We define ˜ T i := T c (1 + | ˜ t | /τ Q i ) and ˜ ξ i := ξ ( | ˜ t | /τ Q i ). The dotted curve represents correlation length in the equilibrium state. <details> <summary>Image 8 Details</summary>  ### Visual Description ## Line Graphs: Correlation Length vs. Temperature ### Overview The image presents three line graphs illustrating the relationship between correlation length (ξ) and temperature (T) for different scenarios denoted as τQ1, τQ2, and τQ3. Each graph shows a distinct behavior of the correlation length as temperature changes, particularly around a critical temperature (Tc). ### Components/Axes * **Vertical Axis (y-axis):** Represents the correlation length, labeled as "ξ" (xi). The y-axis markers are labeled as ξ̃1, ξ̃2, and ξ̃3 for the three graphs respectively. * **Horizontal Axis (x-axis):** Represents the temperature, labeled as "T". * **Curves:** Each graph contains a solid line representing the correlation length as a function of temperature. Additionally, each graph contains a dashed line. * **Vertical Lines:** Each graph has two vertical dotted lines, one at Tc (critical temperature) and another at a temperature denoted as T̃1, T̃2, and T̃3 for the three graphs respectively. * **Titles:** The graphs are titled as τQ1, τQ2, and τQ3. ### Detailed Analysis **Graph 1: τQ1** * **Trend:** The solid line starts at ξ̃1 and remains constant until T̃c. After T̃c, the line decreases sharply until T̃1, then continues to decrease at a slower rate as temperature increases. The dashed line increases sharply as it approaches T̃c from the right. * **Key Points:** * Correlation length is constant at ξ̃1 for T < T̃c. * Correlation length decreases rapidly for T > T̃c. **Graph 2: τQ2** * **Trend:** The solid line starts at ξ̃2 and remains constant until T̃c. After T̃c, the line decreases sharply until T̃2, then continues to decrease at a slower rate as temperature increases. The dashed line increases sharply as it approaches T̃c from the right. * **Key Points:** * Correlation length is constant at ξ̃2 for T < T̃c. * Correlation length decreases rapidly for T > T̃c. **Graph 3: τQ3** * **Trend:** The solid line starts at ξ̃3 and remains constant until T̃c. After T̃c, the line decreases sharply until T̃3, then continues to decrease at a slower rate as temperature increases. The dashed line increases sharply as it approaches T̃c from the right. * **Key Points:** * Correlation length is constant at ξ̃3 for T < T̃c. * Correlation length decreases rapidly for T > T̃c. ### Key Observations * All three graphs exhibit a similar pattern: a constant correlation length below a critical temperature (Tc) followed by a decrease in correlation length as temperature increases beyond Tc. * The dashed lines in each graph approach infinity as T approaches Tc from the right. * The temperatures T̃1, T̃2, and T̃3 are greater than Tc in each respective graph. ### Interpretation The graphs illustrate how the correlation length (ξ) of a system changes with temperature (T) near a critical point (Tc). The constant correlation length below Tc suggests a state of order or coherence. Above Tc, the rapid decrease in correlation length indicates a transition to a disordered state. The dashed lines likely represent a theoretical or limiting behavior of the correlation length, diverging at Tc, which is characteristic of second-order phase transitions. The different values of ξ̃1, ξ̃2, and ξ̃3 and the different temperatures T̃1, T̃2, and T̃3 suggest that the specific details of this transition depend on the parameter τQ, which could represent different physical conditions or system properties. </details> This argument is called the Kibble-Zurek mechanism. Since the KibbleZurek mechanism explains the creation of topological defects induced by cooling of the system which takes place the second-order phase transition, this relates to the evolution of cosmic strings by spontaneous symmetry breaking in the Big Bang theory. 119-121 The Kibble-Zurek mechanism can also describe the creation of topological defects in magnetic models, 122,123 superfluid helium systems, 124,125 and Bose-Einstein condensations. 126,127 Next we consider the efficiency of the simulated annealing and the quantum annealing using the Kibble-Zurek mechanism by taking examples which can be treated analytically. ## 5.1.1. Efficiency of Simulated Annealing and Quantum Annealing Next, we consider the efficiency of the simulated annealing and the quantum annealing according to the Kibble-Zurek mechanism. As an example, we treat the case where the non-domain wall state is the best solution. In this case, the value of n ( ˜ t ) approximately represents the difference between the obtained solution and the best solution. Thus, by using the Kibble-Zurek mechanism, we can compare the efficiency of annealing methods from the behavior of n ( ˜ t ) against the annealing speed. Suppose we solve optimization problems by using annealing methods, we would like to obtain a better solution as fast as possible, in other words, as small τ Q as possible. Then, the comparison obtained by the Kibble-Zurek mechanism is expected to become an useful information for the optimization problems. As an example, we consider the efficiency of the simulated annealing and the quantum annealing for the random ferromagnetic Ising chain in terms of the Kibble-Zurek mechanism according to Refs. [128,129]. ## 5.1.2. Simulated Annealing for Random Ferromagnetic Ising Chain The model Hamiltonian of the random ferromagnetic Ising chain is given as $$\mathcal { H } = - \sum _ { i } J _ { i } \sigma _ { i } ^ { z } \sigma _ { i + 1 } ^ { z } , \quad \sigma _ { i } ^ { z } = \pm 1 , <text><loc_465><loc_152><loc_500><loc_177>(115)</text> <text><loc_466><loc_240><loc_499><loc_273>The site</text> </doctag>$$ where J i is the interaction between the i -th site and the ( i +1)-th site. The value of J i is given by the uniform distribution between 0 < J i ≤ 1. The distribution function P (u) ( J i ) is given by $$P ^ { ( u ) } ( J _ { i } ) \colon = \begin{cases} 1 & f o r $ 0 < J _ { i } \leq 1 $ \\ 0 & o t h e r w i s e \end{cases} .$$ Since the interaction J i is always positive value, the ground state spin configuration is the all-up spin state or the all-down spin state. In this model, the ferromagnetic transition occurs at zero temperature. The correlation function between two sites where the distance is r is written as $$[ \langle \sigma _ { i } \sigma _ { i + r } \rangle ] _ { a v } = \left ( \frac { 1 } { \beta } \ln \cosh \beta \right ) ^ { r } ,$$ where 〈· · · 〉 and [ · · · ] av denote the thermal average and the random average. Physical quantities should depend on the specific spatial pattern of the random interactions { J i } . Then, these averages are defined by $$\langle O ( \{ J _ { i } \} ) \rangle \colon = \frac { T r \, O ( \{ J _ { i } \} ) e ^ { - \beta \mathcal { H } } } { T r \, e ^ { - \beta \mathcal { H } } } ,$$ $$[ O ( \{ J _ { i } \} ) ] _ { a v } \colon = \int \prod _ { i } d J _ { i } P ^ { ( u ) } ( J _ { i } ) O ( \{ J _ { i } \} ) , <text><loc_465><loc_192><loc_500><loc_200>(119)</text> <text><loc_466><loc_250><loc_499><loc_257>1119)</text> </doctag>$$ respectively. We omit the argument ( { J i } ) for simplicity. The relationship between the correlation function and the correlation length ξ is given by $$[ \langle \sigma _ { i } \sigma _ { i + r } \rangle ] _ { a v } = e ^ { - r / \xi } .$$ Here we mainly focus on the low-temperature limit, since the correlation length grows as temperature decreases. Then the correlation length is given as $$\xi = - \frac { 1 } { \ln ( \beta ^ { - 1 } \ln \cosh \beta ) } \simeq \frac { \beta } { \ln 2 } . \quad \ \ ( 1 2 1 )$$ Here, we adopt the Glauber dynamics 130 as the time development, and thus the relaxation time τ rel can be written as $$\tau _ { r e l } = \frac { 1 } { 1 - t a n h \, 2 \beta } \simeq \frac { 1 } { 2 } e ^ { 4 \beta } = \frac { 1 } { 2 } e ^ { 4 \xi \ln 2 } .$$ As we can see, in this model, the correlation length ξ and the relaxation time τ rel are not the power function of temperature unlike the case of the systems where the second-order phase transition occurs at finite temperature (Eqs. (104) and (105)). This is because properties are different between phase transition which exhibits at finite temperature and that occurs at zero temperature. We decrease temperature T ( t ) against the time t as following schedule: $$T ( t ) = - \frac { t } { \tau _ { Q } } \quad ( - \infty < t \leq 0 ) . \quad \quad ( 1 2 3 )$$ Here T c = 0 in this system. According to the Kibble-Zurek mechanism, we define ˜ t by following relation: $$\tau _ { r e l } ( T ( \tilde { t } ) ) = | \tilde { t } | , & & ( 1 2 4 )$$ $$T ( \tilde { t } ) = \frac { | \tilde { t } | } { \tau _ { Q } } = \frac { \tau _ { r e l } ( T ( \tilde { t } ) ) } { \tau _ { Q } } .$$ By using Eqs. (121) and (122), low-temperature limit of Eq. (125) is written as $$\frac { 1 } { \xi ( T ( \tilde { t } ) ) \ln 2 } \simeq \frac { 1 } { 2 \tau _ { Q } } e ^ { 4 \xi ( T ( \tilde { t } ) ) \ln 2 } , \quad ( 1 2 6 )$$ and, we obtain and, we obtain $$\xi ( T ( \tilde { t } ) ) = \frac { \ln \tau _ { Q } + \ln 2 - \ln ( \xi ( T ( \tilde { t } ) ) \ln 2 ) } { 4 \ln 2 } \, \infty \, \frac { \ln \tau _ { Q } } { 4 \ln 2 } . \quad ( 1 2 7 )$$ The approximation of RHS is valid in the case of τ Q /greatermuch 1 which indicates very slow annealing speed. Thus, we can estimate the density of domain wall n SA ( ˜ t ) at t = ˜ t as follows: $$n _ { S A } ( \tilde { t } ) \, \infty \, \frac { 4 \ln 2 } { \ln \tau _ { Q } } . \quad \ \ ( 1 2 8 )$$ ## 5.1.3. Quantum Annealing for Random Ferromagnetic Ising Chain We study the Kibble-Zurek mechanism for the random ferromagnetic Ising chain with transverse field Γ. The model Hamiltonian is given as $$\hat { \mathcal { H } } = - \sum _ { i } J _ { i } \hat { \sigma } _ { i } ^ { z } \hat { \sigma } _ { i + 1 } ^ { z } - \Gamma \sum _ { i } \hat { \sigma } _ { i } ^ { x } ,$$ where the value of J i is given by the uniform distribution between 0 < J i ≤ 1 as well as the case of simulated annealing. In this model, the quantum phase transition from the paramagnetic phase to the ferromagnetic phase occurs at Γ c = exp([ln J i ] av ). 131 Here, we define the dimensionless transverse field g as $$g \colon = \frac { \Gamma - \Gamma _ { c } } { \Gamma _ { c } } . \quad ( 1 3 0 )$$ When | g | /lessmuch 1, it has been known that the correlation length obtained by the renormalization group analysis 132 is scaled by the following relation: $$\xi ( g ) \varpropto | g | ^ { - \nu } \quad ( \nu = 2 ) . \quad & ( 1 3 1 )$$ Moreover, a coherence time τ coh is scaled by $$\tau _ { c o h } ( g ) \, \infty \, [ \xi ( g ) ] ^ { z } \, \infty \, | g | ^ { - \nu z } \quad ( \nu = 2 ) , \quad \quad ( 1 3 2 )$$ where the dynamical exponent z is scaled as $$z \, \infty \, \frac { 1 } { | g | } , <text><loc_435><loc_45><loc_464><loc_65>(133)</text>$$ which is also obtained by the renormalization group analysis. 132 This means that the dynamical exponent diverges at the transition point, and this behavior is a qualitative difference between the random system and the pure system ( z = 1). From this fact, τ coh cannot be expressed by the power function of g unlike the case of the second-order phase transition at finite temperature. We decrease transverse field Γ( t ) against the time t as following schedule: $$\Gamma ( t ) = \Gamma _ { c } \left ( 1 - \frac { t } { \tau _ { Q } } \right ) \quad ( - \infty < t \leq \tau _ { Q } ) .$$ According to the Kibble-Zurek mechanism, we define ˜ t by following relation: $$\tau _ { \text {coh} } ( g ( \tilde { t } ) ) = | \tilde { t } | , & & ( 1 3 5 )$$ $$g ( \tilde { t } ) = \frac { | \tilde { t } | } { \tau _ { Q } } = \frac { \tau _ { c o h } ( g ( \tilde { t } ) ) } { \tau _ { Q } } .$$ By using Eqs. (131), (132), and (133), Eq. (136) is written as $$\frac { 1 } { \sqrt { \xi ( g ( \tilde { t } ) ) } } \, \infty \, \frac { 1 } { \tau _ { Q } } | \xi ( g ( \tilde { t } ) ) | ^ { z } \, \infty \, \frac { 1 } { \tau _ { Q } } | \xi ( g ( \tilde { t } ) ) | ^ { \sqrt { \xi ( g ( \tilde { t } ) ) } } , \frac { 1 } { \sqrt { \xi ( g ( \tilde { t } ) ) } } \, \infty \, \frac { 1 } { \tau _ { Q } } | \xi ( g ( \tilde { t } ) ) | ^ { z } \, \infty \, \frac { 1 } { \tau _ { Q } } | \xi ( g ( \tilde { t } ) ) | ^ { \sqrt { \xi ( g ( \tilde { t } ) ) } } , \\ \text {d} , \, w e o b t a i n$$ and, we obtain $$\left ( \sqrt { \xi ( g ( \tilde { t } ) ) } + \frac { 1 } { 2 } \right ) \ln \xi ( g ( \tilde { t } ) ) \, \infty \, \ln \tau _ { Q } . ( 1 3 8 ) \\ \intertext { ( 1 3 8 ) } \left ( \sqrt { \xi ( g ( \tilde { t } ) ) } + \frac { 1 } { 2 } \right ) \ln \xi ( g ( \tilde { t } ) ) \, \infty \, \ln \tau _ { Q } .$$ In the limit of τ Q /greatermuch 1, since the value of ξ ( g ( ˜ t )) is very large, and we obtain 133 $$\sqrt { \xi ( g ( \tilde { t } ) ) } + \frac { 1 } { 2 } \simeq \sqrt { \xi ( g ( \tilde { t } ) ) } , 3 3 & & ( 1 3 9 ) \\ & & ( \xi ( g ( \tilde { t } ) ) ) \\ & & ( \xi ( g ( \tilde { t } ) ) )$$ $$\xi ( g ( \tilde { t } ) ) \, \infty \, \left ( \frac { \ln \tau _ { Q } } { \ln \xi ( g ( \tilde { t } ) ) } \right ) ^ { 2 } .$$ and we obtain Moreover, since the change of ln ξ ( g ( ˜ t )) is gradual in comparison with that of ξ ( g ( ˜ t )), we neglect ln ξ ( g ( ˜ t )) and obtain $$\xi ( g ( \tilde { t } ) ) \, \infty \, ( \ln \tau _ { Q } ) ^ { 2 } \, .$$ From this relation, we can estimate the density of domain wall n QA ( ˜ t ) at t = ˜ t as follows: $$n _ { Q A } ( \tilde { t } ) \, \infty \, ( \ln \tau _ { Q } ) ^ { - 2 } \, . & & ( 1 4 2 )$$ ## 5.1.4. Comparison between Simulated and Quantum Annealing Methods We have shown analysis of the domain wall density in the random ferromagnetic Ising chain during the simulated annealing and the quantum annealing by the Kibble-Zurek mechanism. The obtained densities of domain wall are $$\begin{array} { r l } { n _ { S A } ( \tilde { t } ) \, \infty \, ( \ln \tau _ { Q } ) ^ { - 1 } } & { \colon \, s i m u l a t e d a n e a l i n g , } \\ { ( \tilde { t } ) } & { ( 1 4 3 ) } \end{array}$$ From these relations, it is clear that the decay of n QA ( ˜ t ) is faster than that of n SA ( ˜ t ) against the value of τ Q . Thus, from the Kibble-Zurek mechanism, it is concluded that the quantum annealing method is appropriate as the annealing method for the random ferromagnetic Ising chain in comparison with the simulated annealing method. Suppose we consider the ferromagnetic Ising chain with homogeneous interaction ( J i = 1 for all i ). In this case, both the domain wall density in the simulated annealing and that in the quantum annealing are obtained as $$n _ { Q A } ( \tilde { t } ) \, \infty \, ( \ln \tau _ { Q } ) ^ { - 2 } \quad \colon \quad q u a n t u m \, a n e l i n g .$$ $$n ( \tilde { t } ) \, \infty \, \frac { 1 } { \sqrt { \tau _ { Q } } } , \quad ( 1 4 5 )$$ This relation for the simulated annealing can be obtained by a simple calculation as well as the case of the random Ising spin chain. On top of that, the relation for the quantum annealing can be derived by Eq. (113). Here the critical exponent ν of the transverse Ising chain with homogeneous interaction is ν = 1 and the dynamical exponent of this system is z = 1. Then there is no difference between the simulated annealing and the quantum annealing in the case of the homogeneous ferromagnetic Ising chain. However, since the optimization problem has some kind of randomness, the abovementioned result encourages that the quantum annealing is better than the simulated annealing for optimization problems. In general, the existence of the phase transition in optimization problems negatively influences performance of annealing methods. Here, we have introduced the Kibble-Zurek mechanism relating to the dynamics which passes across the second-order phase transition point. As the specific example, we have analyzed the efficiencies of the simulated annealing and the quantum annealing for the random ferromagnetic Ising chain according to the Kibble-Zurek mechanism. For this model, the efficiency of the quantum annealing is better than that of the simulated annealing. Of course, since the efficiency of annealing methods depends on the details of optimization problems, it is not to say that the quantum annealing is always appropriate as the annealing method for general optimization problems in comparison with the simulated annealing. Moreover, we have to develop a theory based on the Kibble-Zurek mechanism itself, 134 since we assume the growth of the correlation length stops at t > ˜ t . For example, if we adapt the Kibble-Zurek mechanism to two- or three-dimensional models and more complicated models, it is difficult to estimate the correlation length analytically, and thus we should execute numerical simulations such as the Monte Carlo simulation. For example, in the two-dimensional Ising model with random interactions, it has been shown that the efficiency of the quantum annealing is better than that of the simulated annealing by Monte Carlo simulation. 129 Although the efficiency of annealing methods for a number of optimization problems has been clarified by the Kibble-Zurek mechanism, it remains to be an open problem to investigate when to use the quantum annealing exhaustively. In the above-mentioned argument, the phase transition under consideration is of the second order. What happens if we adapt the same argument for the other type phase transitions such as first-order phase transition and Kosterlitz-Thouless (KT) transition? In these phase transitions, the behaviors of correlation length are different from that in systems where a secondorder phase transition occurs: the finite-correlation length at the first-order phase transition point and the quasi-long-range correlation length at the KTtransition point. Thus, it is an interesting problem to clarify relationship between behaviors of correlation length and the generalized Kibble-Zurek mechanism. By considering dynamical nature of the optimization problems in terms of non-equilibrium statistical physics in a deeper way, we believe that the quantum annealing method will become a central part of practical method for optimization problems. ## 5.2. Frustration Effects for Simulated Annealing and Quantum Annealing In many cases optimization problems can be represented by the Ising model with random interactions and magnetic fields as mentioned before. The Hamiltonian of this system is given by $$\mathcal { H } = - \sum _ { i , j } J _ { i j } \sigma _ { i } ^ { z } \sigma _ { j } ^ { z } - \sum _ { i = 1 } ^ { N } h _ { i } \sigma _ { i } ^ { z } , \quad \sigma _ { i } ^ { z } = \pm 1 .$$ When all interactions are ferromagnetic as the previous example in Sec. 5.1, the ground state is the all-up or the all-down states. However, if there are antiferromagnetic interactions in the system, the situation becomes different. In order to show the difference between ferromagnetic interaction and antiferromagnetic interaction, we first consider three spin system on triangle cluster as shown in Fig. 9. In this section, we treat the case for h i = 0 for all i . The dotted and solid lines in Fig. 9 represent ferromagnetic and antiferromagnetic interactions, respectively. The considered Hamiltonian is written as $$\mathcal { H } _ { t r i a n g l e } = - J ( \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } + \sigma _ { 2 } ^ { z } \sigma _ { 3 } ^ { z } + \sigma _ { 3 } ^ { z } \sigma _ { 1 } ^ { z } ) .$$ Here we set the all interactions are the same value for simplicity. The ground states for positive J (ferromagnetic interaction) are the all-up or the alldown states shown in Fig. 9 (a). In these states, all spins between all interactions are energetically favorable states. In the case of negative J (antiferromagnetic interaction), while on the other hand, six states shown in Fig. 9 (b) are ground states. These ground states have unfavorable interactions Fig. 9. Three spin system on triangle cluster. The dotted and solid lines represent ferromagnetic and antiferromagnetic interactions, respectively. The open and solid circles are the +1-state and the -1-state, respectively. The crosses indicate the positions of unfavorable interactions. (a) Ground states for ferromagnetic case. (b) Ground states for antiferromagnetic case. <details> <summary>Image 9 Details</summary>  ### Visual Description ## Diagram: Triangle Configurations ### Overview The image shows two sets of triangle configurations, labeled (a) and (b). Each triangle has vertices represented by circles, which are either filled (black) or unfilled (white). In set (a), there are two triangles, one with all white vertices and one with all black vertices. In set (b), there are six triangles, each with a mix of black and white vertices, and each containing an "X" mark near one of the vertices. ### Components/Axes * **Vertices:** Represented by circles, either filled (black) or unfilled (white). * **Edges:** Represented by solid lines connecting the vertices. * **Dashed Lines:** In diagram (a), dashed lines connect the vertices of each triangle. * **"X" Marks:** Present in all triangles in diagram (b), located near one of the vertices. * **Labels:** (a) and (b) identify the two sets of configurations. ### Detailed Analysis or ### Content Details **Diagram (a):** * **Top Triangle:** All three vertices are unfilled (white). The edges are dashed. * **Bottom Triangle:** All three vertices are filled (black). The edges are dashed. **Diagram (b):** There are six triangles, arranged in two rows of three. Each triangle has two vertices of one color and one vertex of the other color. Each triangle also has an "X" mark near one of the vertices. * **Top Row, Left Triangle:** Top vertex is white with an "X" mark. Bottom-left vertex is white. Bottom-right vertex is black. * **Top Row, Middle Triangle:** Top vertex is white. Bottom-left vertex is black. Bottom-right vertex is white with an "X" mark. * **Top Row, Right Triangle:** Top vertex is black. Bottom-left vertex is white. Bottom-right vertex is white with an "X" mark. * **Bottom Row, Left Triangle:** Top vertex is black with an "X" mark. Bottom-left vertex is white. Bottom-right vertex is black. * **Bottom Row, Middle Triangle:** Top vertex is black. Bottom-left vertex is black with an "X" mark. Bottom-right vertex is white. * **Bottom Row, Right Triangle:** Top vertex is white. Bottom-left vertex is black. Bottom-right vertex is black with an "X" mark. ### Key Observations * Diagram (a) shows two extreme configurations: all vertices white and all vertices black. * Diagram (b) shows all possible configurations with two vertices of one color and one vertex of the other color, along with an "X" mark near each vertex. * The "X" mark appears to be associated with a specific vertex in each triangle in diagram (b). ### Interpretation The image likely represents different states or configurations of a system with three components. The filled and unfilled circles could represent two different states of each component. Diagram (a) shows the two homogeneous states, while diagram (b) shows the heterogeneous states. The "X" mark might indicate a specific property or modification of the vertex it is near. The dashed lines in diagram (a) might indicate a weaker or different type of connection compared to the solid lines in diagram (b). Without further context, it's difficult to determine the exact meaning of these configurations. </details> indicated by the crosses in Fig. 9 (b). This situation is called frustration. In the homogeneous antiferromagnetic Ising spin systems on lattices based on triangle such as triangular lattice and kagom´ e lattice, frustration appears in all triangles. Since such frustration comes from lattice geometry, this is called geometrical frustration. It should be noted that the homogeneous antiferromagnetic Ising spin systems on square lattice and hexagonal lattice have no frustration. Since these systems are bipartite systems which can be decomposed by two sublattices, these systems can be transformed on the ferromagnetic systems by local gauge transformation of all spins belonging to one of the sublattices. Frustration appears in also inhomogeneous systems as shown in Fig. 10. The squares pointed by stars in Fig. 10 represent frustration plaquettes which are satisfied following relation: $$\kappa _ { k } \colon = \prod _ { i , j \in \square _ { k } } J _ { i j } < 0 ,$$ where /square k indicates the smallest square plaquette at the position k . If κ k for all k is positive, the system is not frustrated. In general, frustration prevents the system from conventional magnetic ordering such as ferromagnetic order and N´ eel order, since there is no state where all interactions are satisfied energetically in frustrated systems. Frustration makes peculiar density of states which induces unconventional phase transition and slow dynamics. 112,115,135-144 Although many optimization problems can be represented by the Ising model with random interactions and magnetic fields, here we focus on the frustration effect which comes Fig. 10. A ground state of the Ising spin system with random interactions. The dotted and solid lines represent ferromagnetic and antiferromagnetic interactions, respectively. The open and solid circles are the +1-state and the -1-state, respectively. The stars and crosses indicate frustration plaquettes and unfavorable interactions, respectively. <details> <summary>Image 10 Details</summary>  ### Visual Description ## Diagram: Lattice Network with Defects ### Overview The image depicts a 5x5 lattice network composed of nodes (circles) connected by lines. Some nodes are filled black, while others are white. The connections between nodes are either solid or dashed. Star symbols and "X" symbols are interspersed throughout the lattice, representing defects or other features. ### Components/Axes * **Nodes:** Represented by circles, either filled black or white. * **Connections:** Solid lines indicate one type of connection, while dashed lines indicate another. * **Stars:** Five-pointed star symbols are scattered throughout the lattice. * **"X" Symbols:** "X" symbols are also scattered throughout the lattice. ### Detailed Analysis The lattice structure is a 5x5 grid. The nodes alternate in color (black and white) in a checkerboard pattern, but this pattern is disrupted by the presence of defects. * **Node Colors:** * Black nodes are present in the top-left, second row first column, third row first column, fourth row first column, fifth row second column, third row third column, fourth row fourth column, fifth row fifth column, first row fifth column, second row fifth column. * White nodes are present in the remaining locations. * **Connections:** * Solid lines connect nodes vertically and horizontally. * Dashed lines connect nodes diagonally. * **Defect Placement:** * Stars are placed in the spaces between nodes, often near dashed lines. * "X" symbols are placed near nodes, often disrupting the checkerboard pattern. ### Key Observations * The lattice has a generally regular structure, but the presence of stars and "X" symbols introduces irregularities. * The black and white nodes form a disrupted checkerboard pattern. * The solid and dashed lines create a network of connections between the nodes. ### Interpretation The diagram likely represents a physical or abstract system with a lattice structure. The black and white nodes could represent different states or types of elements within the system. The solid and dashed lines could represent different types of interactions or connections between these elements. The stars and "X" symbols could represent defects, impurities, or other perturbations in the system. The diagram could be used to model various phenomena, such as crystal structures, social networks, or computer networks. The placement of the defects and the disruption of the checkerboard pattern suggest that these defects have a significant impact on the overall structure and properties of the system. </details> from non-random interactions. In terms of statistical physics, this is a firststep study to investigate similarities and differences between thermal fluctuation and quantum fluctuation for frustrated systems. Furthermore, it is important topic for the optimization problems to consider the thermal fluctuation and quantum fluctuation effects for frustrated systems. To obtain the ground state of frustrated systems is to find how to put the unsatisfied bonds represented by the crosses. Since the unsatisfied bonds are regarded as some kind of constraints, this situation is similar with the traveling salesman problem in which there are some constraints as mentioned before. We explain two topics in this section. In the first half, we consider the order by disorder effect in fully-frustrated systems. In the last half, we explain non-monotonic dynamics in decorated bond systems. ## 5.2.1. Thermal Fluctuation and Quantum Fluctuation Effect of Geometrical Frustrated Systems In general, there are many degenerated ground states in geometrical frustrated systems such as triangular antiferromagnetic Ising spin systems and kagom´ e antiferromagnetic Ising spin systems. In these cases, non-zero residual entropy which is entropy at zero temperature exists. Typical configurations of ground states of the triangular antiferromagnetic Ising spin systems are shown in Fig. 11. The residual entropy per spin of this system is S (tri) res /similarequal 0 . 323 k B , 145-148 where k B is the Boltzmann constant. Since the total entropy per spin is k B ln 2 /similarequal 0 . 693 k B , 46 . 6% of the total entropy remains even at zero temperature. In other words, there are macroscopic degenerated ground states in this system. In the antiferromagnetic Ising spin system on kagom´ e lattice, there are also macroscopic degenerated ground states. The residual entropy per spin of this system is S (kag) res /similarequal 0 . 502 k B , which is 72 . 4% of the total entropy. 149 Suppose we apply the simulated annealing or the quantum annealing with slow schedule for geometrical frustrated spin systems. Since there are macroscopically degenerated ground states in these systems, our purpose is to clarify whether all ground states are obtained with the same probabilities or biased probabilities. We first consider the obtained ground states in the case of the simulated annealing with slow schedule. If we decrease temperature slow enough, the obtained state should satisfy the equilibrium probability distribution. When the temperature is k B T /lessmuch | J | , the equilibrium probabilities of the ground states are dominant and that of any excited states can be neglected. The principle of equal weight which is the keystone in the equilibrium statistical physics says that if the eigenenergies of the Fig. 11. Typical configurations of ground states of antiferromagnetic Ising spin system on triangular lattice. The open and solid circles are the +1-state and the -1-state, respectively. The dotted circles indicate free spin where the molecular field is zero. <details> <summary>Image 11 Details</summary>  ### Visual Description ## Diagram: Triangular Lattice with Black and White Nodes ### Overview The image depicts two sections of a triangular lattice. The lattice consists of nodes connected by lines, forming equilateral triangles. The nodes are colored either black or white. Dashed circles surround the nodes, possibly indicating a specific region or selection. The two sections are separated by a gap. ### Components/Axes * **Nodes:** Represented by circles, colored either black or white. * **Edges:** Represented by straight lines connecting the nodes, forming equilateral triangles. * **Lattice Structure:** The nodes and edges form a repeating triangular pattern. * **Dashed Circles:** Encircle each node. * **Sections:** The diagram is split into two distinct sections, separated by a gap. ### Detailed Analysis **Left Section:** * The left section starts with a single black node at the bottom. * The next row has two nodes, one white and one black. * The subsequent rows alternate between black and white nodes, creating a pattern. * The top row has several nodes, with the rightmost node being black. **Right Section:** * The right section begins with two white nodes at the bottom. * The nodes alternate between black and white, similar to the left section. * The top row has several nodes, with the rightmost node being white. **Node Color Distribution:** * The distribution of black and white nodes appears somewhat random, but there might be underlying patterns. * Both sections exhibit a mix of black and white nodes. ### Key Observations * The diagram showcases a triangular lattice structure. * The nodes are colored black and white, potentially representing different states or properties. * The dashed circles around the nodes might indicate a selection or highlight. * The two sections are separated, suggesting a break or discontinuity in the lattice. ### Interpretation The diagram likely represents a system or model where nodes in a triangular lattice can exist in two states (black or white). The lattice structure could represent physical connections or interactions between the nodes. The dashed circles might indicate nodes that are being considered or analyzed. The separation into two sections could represent different regions or conditions within the system. Without further context, it's difficult to determine the specific meaning or purpose of the diagram, but it could be related to topics such as: * Statistical physics (e.g., Ising model on a triangular lattice) * Materials science (e.g., arrangement of atoms in a crystal) * Graph theory (e.g., coloring of nodes in a graph) </details> microscopic state Σ A and Σ B are the same, the equilibrium probability of Σ A and that of Σ B are also the same. Then we obtain all macroscopic degenerated ground states with the same probability after the simulated annealing with slow schedule. Next we consider the obtained ground states in the case of the quantum annealing where the transverse field decreases slow enough. Here we assume that the initial state is set to be the ground state of the Hamiltonian at the initial time. In order to capture the feature of the ground states in a graphical way, it is convenient to introduce the concept of free spin where the molecular field is zero. The molecular field at the i -th site is given by $$h _ { i } ^ { ( \text {eff} ) } \colon = J \sum _ { j } { ^ { \prime } \sigma _ { j } ^ { z } } ,$$ where the summation runs over the nearest-neighbor sites of the i -th site. For instance, in Fig. 11, spins indicated by dotted circles are free spins. Here, the transverse field is expressed as $$- \Gamma \sum _ { i } \hat { \sigma } ^ { x } _ { i } = - \Gamma \sum _ { i } ( \hat { \sigma } ^ { + } _ { i } + \hat { \sigma } ^ { - } _ { i } ) , <text><loc_465><loc_152><loc_500><loc_175>(150)</text> </doctag>$$ where ˆ σ + i and ˆ σ -i denote the raising and lowering operators at the i -th site, respectively. They are defined by $$\hat { \sigma } ^ { + } \colon = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} , \quad \hat { \sigma } ^ { - } \colon = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} .$$ The x -component of the Pauli matrix corresponds to the operator which flips the considered spin: $$\hat { \sigma } ^ { x } \left | \uparrow \right > = \left | \downarrow \right > , \quad \hat { \sigma } ^ { x } \left | \downarrow \right > = \left | \uparrow \right > . & & ( 1 5 2 )$$ From this, the states which have large number of free spins are expected to become stable at the limit of Γ → 0+ and T = 0. Actually, in the adiabatic limit, the amplitudes of the states which have the maximum number of free spins are larger than the others. 150-154 When we decrease the transverse field slow enough, the state at each time can be well approximated by the ground state of the instantaneous Hamiltonian. Then we obtain specific ground states with high probability after the quantum annealing with slow schedule. In this section, we considered the thermal fluctuation effect and the quantum fluctuation effect in the adiabatic limit. The simulated annealing can obtain all the ground states with the same probability, while on the other hand, the quantum annealing can obtain specific ground states in this limit. The biased probability distribution can be explained by the character of the quantum Hamiltonian. The selected states should depend on how to choice the quantum Hamiltonian. When we adopt the exchange type interaction as the quantum field, the states that have the maximum value of the 'free spin pair' should be selected. Moreover, it is an interesting topic to investigate differences between the simulated annealing and the quantum annealing with finite speed not only in terms of the quantum annealing but also in nonequilibrium statistical physics and condensed matter physics. At the present stage, to consider dynamic phenomena in strongly correlated systems is difficult, since a small number of theoretical methods for obtaining dynamic phenomena have been developed. If the technology of the artificial lattices develops more than ever, real-time dynamics and time-dependent phenomena of frustrated spin systems can be observed in real experiments. ## 5.2.2. Non-Monotonic Behavior of Correlation Function in Decorated Bond System In the ferromagnetic Ising spin systems, the correlation function behaves monotonic against the temperature and transverse field. However, the behavior of the correlation function is non-monotonic as a function of temperature in some frustrated spin systems. As an example of non-monotonic correlation function, we introduce equilibrium properties of the correlation function in decorated bond systems in which the frustration exists. The Hamiltonian of the decorated bond systems where the number of system spins is two shown in Fig. 12 is given by $$\mathcal { H } = - J _ { d i r } \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } - J \sum _ { i = 1 } ^ { N _ { d } } s _ { i } ^ { z } ( \sigma _ { 1 } ^ { z } + \sigma _ { 2 } ^ { z } ) ,$$ where σ z i = ± 1 and s z i = ± 1 are, respectively, called system spins and decorated spins, and N d is the number of decorated spins. The circles and the squares in Fig. 12 represent the system spins and the decorated spins, respectively. When the direct interaction between system spins J dir is zero and the decorated bond J is positive, the correlation function between system spins 〈 σ z 1 σ z 2 〉 is always positive and monotonic decaying function against the temperature. When the direct interaction between system spins J dir is negative and the decorated bond J is zero, on the other hand, the correlation function 〈 σ z 1 σ z 2 〉 is always negative and monotonic increasing function against the temperature. From this, the correlation function 〈 σ z 1 σ z 2 〉 is expected to behave non-monotonic in some cases for negative J dir and positive J or positive J dir and negative J . In order to obtain temperature dependence of the correlation function between system spins, we trace over spin states except the system spins: $$T r _ { \{ s _ { i } ^ { z } \} } e ^ { - \beta \mathcal { H } } = A e ^ { K _ { e f f } \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } } , \quad ( 1 5 4 )$$ where A is just a constant which does not affect any physical quantities and the effective coupling K eff is given by $$K _ { e f f } = \frac { N _ { d } } { 2 } \ln \cosh ( 2 \beta J ) + \beta J _ { d i r } . \quad \ \ ( 1 5 5 )$$ Temperature dependence of the correlation function between system spins Fig. 12. Decorated bond system where the number of system spins is two and the number of decorated spins is four ( N d = 4). The circles and squares represent system spins and decorated spins, respectively. The dotted and solid lines indicate the direct interaction between system spins and the decorated bonds, respectively. <details> <summary>Image 12 Details</summary>  ### Visual Description ## Diagram: Bipartite Graph ### Overview The image depicts a bipartite graph consisting of two sets of nodes: circles and squares. The circles are positioned on the left, and the squares are aligned horizontally on the right. Lines connect the circles to the squares, indicating relationships between them. One line between the two circles is dashed. ### Components/Axes * **Nodes:** * Circles: Two circles on the left side. * Squares: Four squares on the right side. * **Edges:** * Solid Lines: Connect each circle to each of the four squares. * Dashed Line: Connects the two circles. ### Detailed Analysis * **Circle Nodes:** There are two circle nodes positioned vertically on the left side of the diagram. * **Square Nodes:** There are four square nodes positioned horizontally on the right side of the diagram. * **Connections:** Each circle node is connected to each of the four square nodes by a solid line. Additionally, the two circle nodes are connected to each other by a dashed line. ### Key Observations * The graph is bipartite because it can be divided into two disjoint sets (circles and squares) such that every edge connects a node in one set to a node in the other set. * The dashed line between the two circle nodes indicates a possible relationship or connection within the circle node set itself. ### Interpretation The diagram represents a bipartite graph, a common structure in computer science and mathematics. The circles and squares represent different entities, and the lines represent relationships between them. The dashed line suggests a special relationship or connection between the two circle nodes, which might indicate a different type of interaction or a shared attribute. The graph could be used to model various scenarios, such as relationships between users and items in a recommendation system, or connections between variables and factors in a statistical model. </details> is represented by using K eff : $$C ^ { ( c ) } ( T ) \coloneqq \langle \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } \rangle = \frac { \text {Tr} \, \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } e ^ { - \beta \mathcal { H } } } { \text {Tr} \, e ^ { - \beta \mathcal { H } } } = \frac { \text {Tr} \, \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } e ^ { K _ { \text {eff} } \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } } } { \text {Tr} \, e ^ { K _ { \text {eff} } \sigma _ { 1 } ^ { z } \sigma _ { 2 } ^ { z } } } = \tanh K _ { \text {eff} } .$$ Hereafter we set J as the energy unit and J is positive. In order to compare the effect of the direct interaction J dir fairly, we assume the form such as J dir = -xN d J . This is because the effective coupling K eff is proportional to the number of decorated spins N d under the assumption. Figure 13 shows temperature dependence of correlation function between the system spins for N d = 1 and N d = 10 for several x . For small x and large x , the correlation function C (c) ( T ) is monotonic decreasing and increasing functions, respectively, against the temperature. However, the correlation function C (c) ( T ) behaves non-monotonic as a function of temperature for intermediate x . At the temperatures where the effective coupling K eff is larger than the critical value of the ferromagnetic Ising spin system on square lattice 19 K (square) c = 1 2 ln(1 + √ 2), ferromagnetic phase appears. On the other hand, at the temperature where K eff is less than -K (square) c , antiferromagnetic phase appears. In this case, successive Fig. 13. The correlation function between system spins C (c) ( T ) as a function of temperature for N d = 1 (left panel) and for N d = 10 (right panel) in the cases of x = 0 . 1, 0 . 2, 0 . 5, 1 . 0, and 2 . 0. <details> <summary>Image 13 Details</summary>  ### Visual Description ## Line Charts: Correlation vs. Temperature ### Overview The image presents two line charts side-by-side, displaying the relationship between a correlation function C^(c)(T) and temperature (T) for different values of a parameter 'x'. The left chart shows a gradual change in correlation with temperature, while the right chart exhibits a more abrupt transition. ### Components/Axes * **Y-axis (Left and Right Charts):** C^(c)(T), ranging from -1 to 1, with markers at -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, and 1. * **X-axis (Left and Right Charts):** T (Temperature), ranging from 0 to 8, with markers at intervals of 2. * **Legend (Top-Right):** Identifies the line styles corresponding to different values of 'x': * x = 0.1: Solid line * x = 0.2: Dashed line * x = 0.5: Dotted line * x = 1.0: Dotted-dashed line * x = 2.0: Dashed-dashed line ### Detailed Analysis **Left Chart:** * **x = 0.1 (Solid Line):** Starts at approximately 1 and decreases to approximately 0.05 as T increases from 0 to 8. * **x = 0.2 (Dashed Line):** Starts at approximately 0.8 and decreases to approximately 0.02 as T increases from 0 to 8. * **x = 0.5 (Dotted Line):** Starts at approximately 0.3 and decreases to approximately -0.25 as T increases from 0 to 8. * **x = 1.0 (Dotted-Dashed Line):** Starts at approximately -0.3 and increases to approximately -0.1 as T increases from 0 to 8. * **x = 2.0 (Dashed-Dashed Line):** Starts at approximately -1 and increases to approximately -0.1 as T increases from 0 to 8. **Right Chart:** * **x = 0.1 (Solid Line):** Remains at approximately 1 until T reaches approximately 1, then decreases to approximately 0.15 as T increases to 8. * **x = 0.2 (Dashed Line):** Remains at approximately 1 until T reaches approximately 1, then decreases to approximately 0.05 as T increases to 8. * **x = 0.5 (Dotted Line):** Remains at approximately 1 until T reaches approximately 0.5, then decreases to approximately -0.6 as T increases to 2, then increases to approximately 0.01 as T increases to 8. * **x = 1.0 (Dotted-Dashed Line):** Starts at approximately -1 and increases to approximately 0.01 as T increases from 0 to 8. * **x = 2.0 (Dashed-Dashed Line):** Starts at approximately -1 and increases to approximately 0.01 as T increases from 0 to 8. ### Key Observations * In the left chart, the correlation function C^(c)(T) generally decreases with increasing temperature for x = 0.1, 0.2, and 0.5, while it increases for x = 1.0 and 2.0. * In the right chart, the correlation function C^(c)(T) exhibits a more abrupt change around T = 1 for x = 0.1, 0.2, and 0.5. For x = 1.0 and 2.0, the correlation function increases with temperature. * The value of 'x' significantly influences the behavior of the correlation function with respect to temperature. ### Interpretation The charts illustrate how the correlation function C^(c)(T) changes with temperature for different values of the parameter 'x'. The left chart shows a more gradual response to temperature changes, while the right chart indicates a critical temperature around T = 1 where the correlation function undergoes a more rapid transition, particularly for lower values of 'x'. The parameter 'x' appears to control the sensitivity of the correlation function to temperature changes and the overall magnitude of the correlation. The data suggests that for lower 'x' values, the system maintains a high correlation until a certain temperature threshold is reached, after which the correlation rapidly diminishes. Conversely, for higher 'x' values, the correlation starts at a negative value and gradually increases with temperature. </details> phase transitions such as paramagnetic → antiferromagnetic → paramagnetic → ferromagnetic phases occur. Such phase transitions are called reentrant phase transitions which are sometimes appeared in frustrated systems. 115,139,155-160 We consider transverse field response of the decorated bond systems in the ground state. The Hamiltonian of the decorated bond system with transverse field is expressed as $$\hat { \mathcal { H } } = - J _ { d i r } \hat { \sigma } ^ { z } _ { 1 } \hat { \sigma } ^ { z } _ { 2 } - J \sum _ { i = 1 } ^ { N _ { d } } \hat { s } ^ { z } _ { i } ( \hat { \sigma } ^ { z } _ { 1 } + \hat { \sigma } ^ { z } _ { 2 } ) - \Gamma ( \hat { \sigma } ^ { x } _ { 1 } + \hat { \sigma } ^ { x } _ { 2 } + \sum _ { i = 1 } ^ { N _ { d } } \hat { s } ^ { x } _ { i } ) ,$$ where ˆ s α i denotes the α -component of the Pauli matrix of the i -th decorated spin. Here we consider transverse-field dependence of the correlation function in the ground state given by $$C ^ { ( q ) } ( \Gamma ) \colon = \langle \psi ^ { ( g s ) } ( \Gamma ) | \hat { \sigma } _ { 1 } ^ { z } \hat { \sigma } _ { 2 } ^ { z } | \psi ^ { ( g s ) } ( \Gamma ) \rangle \, ,$$ where | ψ (gs) (Γ) 〉 denotes the ground state at the transverse field Γ. Figure 14 shows transverse-field dependence of C (q) (Γ) for N d = 1 and N d = 10 for several x . Fig. 14. The correlation function between system spins C (q) (Γ) as a function of transverse field for N d = 1 (left panels) and for N d = 10 (right panels) in the cases of x = 0 . 1, 0 . 2, 0 . 5, 1 . 0, and 2 . 0. <details> <summary>Image 14 Details</summary>  ### Visual Description ## Line Chart: C^(q)(Γ) vs. Γ for different values of x ### Overview The image presents two line charts side-by-side, both plotting C^(q)(Γ) against Γ. The charts display the relationship between these variables for different constant values of 'x', ranging from 0.1 to 2.0. The left chart shows the behavior of the function from Γ = 0 to 8, while the right chart appears to focus on a specific region where the function exhibits more rapid changes. ### Components/Axes * **X-axis (Horizontal):** Γ (Gamma), ranging from 0 to 8 in both charts. * **Y-axis (Vertical):** C^(q)(Γ), ranging from -1 to 1 in both charts. * **Legend (Top-Right of each chart):** * x = 0.1 (Solid Line) * x = 0.2 (Dashed Line) * x = 0.5 (Dotted Line) * x = 1.0 (Dash-Dot Line) * x = 2.0 (Long Dash Line) ### Detailed Analysis **Left Chart:** * **x = 0.1 (Solid Line):** Starts at approximately 0.95 at Γ = 0 and decreases rapidly, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 0.2 (Dashed Line):** Starts at approximately 0.8 at Γ = 0 and decreases, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 0.5 (Dotted Line):** Starts at approximately 0.3 at Γ = 0 and decreases, becoming negative around Γ = 1, reaching approximately -0.25 at Γ = 2, and then slowly approaching 0 for Γ > 4. * **x = 1.0 (Dash-Dot Line):** Starts at approximately -0.3 at Γ = 0 and increases, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 2.0 (Long Dash Line):** Starts at approximately -1 at Γ = 0 and increases, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. **Right Chart:** * **x = 0.1 (Solid Line):** Starts at approximately 1 and decreases rapidly, approaching 0.1 at Γ = 4 and remaining close to 0.1 for Γ > 4. * **x = 0.2 (Dashed Line):** Starts at approximately 1 and decreases rapidly, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 0.5 (Dotted Line):** Starts at approximately -0.9 at Γ = 0 and increases rapidly, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 1.0 (Dash-Dot Line):** Starts at approximately -0.9 at Γ = 0 and increases rapidly, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. * **x = 2.0 (Long Dash Line):** Starts at approximately -0.9 at Γ = 0 and increases rapidly, approaching 0 at Γ = 4 and remaining close to 0 for Γ > 4. ### Key Observations * For smaller values of 'x' (0.1 and 0.2), C^(q)(Γ) starts positive and decreases towards zero as Γ increases. * For larger values of 'x' (0.5, 1.0, and 2.0), C^(q)(Γ) starts negative and increases towards zero as Γ increases. * The right chart shows a zoomed-in view of the initial behavior of the curves, highlighting the rapid changes near Γ = 0. * All curves tend to converge towards zero as Γ increases, suggesting a limiting behavior of the function. ### Interpretation The charts illustrate how the function C^(q)(Γ) behaves as Γ varies, for different fixed values of 'x'. The data suggests that 'x' influences the initial value and the direction of change of C^(q)(Γ). Smaller 'x' values lead to positive initial values and a decreasing trend, while larger 'x' values lead to negative initial values and an increasing trend. The convergence of all curves towards zero as Γ increases indicates that the influence of 'x' diminishes as Γ becomes large. The right chart emphasizes the sensitivity of C^(q)(Γ) to changes in Γ, especially when Γ is small. The relationship between 'x' and C^(q)(Γ) appears to be inversely proportional at Γ = 0. </details> For small x and large x , the correlation function C (q) (Γ) behaves monotonic decreasing and increasing, respectively as a function of transverse field, whereas for intermediate x , transverse-field dependence of the correlation function behaves nonmonotonic as well as the case of thermal fluctuation. Then, the reentrant phase transition also occurs by changing the transverse field. However there is a difference between the thermal fluctuation effect and the quantum fluctuation effect for decorated bond system. The temperature where C (c) ( T ) = 0 is satisfied is the same when we change the number of decorated spins N d , whereas the transverse field at C (q) (Γ) = 0 is different when N d is changed. The thermal fluctuation and the quantum fluctuation have similar properties for the phase transition phenomena in general. Indeed, the reentrant phase transitions occur by changing the thermal fluctuation and also the quantum fluctuation as shown in this section. However as described in Sec. 5.1, in order to obtain the best solution of optimization problems, it is better to erase phase transition. By dealing with thermal and quantum fluctuation effects for frustrated systems exhaustively, we can construct the best form of the adding fluctuation which erases phase transition e . ## 6. Conclusion In this paper, we described some aspects of the quantum annealing from viewpoints of statistical physics, condensed matter physics, and computational physics. Originally, the quantum annealing has been proposed as a method which can solve efficiently optimization problems in a generic way. Since many optimization problems can be mapped onto the Ising model or generalized Ising model such as the clock model and the Potts model, it has been considered that we can obtain a better solution by using methods which were developed in computational physics. For instance, we can obtain a better solution by decreasing temperature (thermal fluctuation) gradually in the simulated annealing which is one of the most famous practical methods. In the quantum annealing, we decrease an introduced quantum field (quantum fluctuation) instead of temperature (thermal fluctuation). In many studies, it was reported that a better solution can be obtained e It is not necessary that the adding fluctuation is restricted in quantum physics. From a viewpoint of optimization problems, we can arbitrary form adding term. Furthermore, it has studied that other novel fluctuation which may be able to erase phase transition as an alternative to thermal and quantum fluctuations. 14,116,161,162 Of course, if we want to realize experimentally, it is better that the added fluctuation term should be some kind of quantum fluctuation. by the quantum annealing efficiently in comparison with the simulated annealing as we explained in Sec. 4. Thus, the quantum annealing method is expected to be a generic and powerful solver of optimization problems as an alternative to the simulated annealing. The quantum annealing has become a milestone of some related fields under the situation in which the quantum annealing itself has been studied exhaustively. Since we use the quantum fluctuation in the quantum annealing with ingenuity, to obtain a better solution by using the quantum annealing is a kind of quantum information processing. Thus, many implementation methods of the quantum annealing in theoretical and experimental ways have been proposed by many researchers. A number of theoretical implementation methods are proposed based on knowledge of statistical physics. As we shown in Sec. 5, question of what are differences between the simulated annealing and the quantum annealing and question of which is efficient in the given optimization problem are catalysts to investigate differences between the thermal fluctuation and the quantum fluctuation in a deeper way. On top of that, studies on the quantum annealing are expected to open the door to consider equilibrium and nonequilibrium statistical physics. Recently, preparation methods of intended Hamiltonian have been established in some experimental systems such as artificial lattices and nuclear magnetic resonance because of recent development of experimental techniques. As long as we use classical computer and our present knowledge, there are a huge number of problems where to obtain the best solution is difficult without any and every approximation in theoretical methods. However if we prepare the Hamiltonian which expresses our intended problem, we can calculate experimentally the stable state of the prepared Hamiltonian in near future. The quantum annealing transcends just a method for obtaining the best solution of optimization problems and it will make a development in wide area of science. Although it seems that studies on the quantum annealing itself have been well established, we believe that the quantum annealing plays a role as a bridge with the abovementioned area of science and the quantum information. ## Acknowledgement The authors are grateful to Bernard Barbara, Bikas K. Charkrabarti, Naomichi Hatano, Masaki Hirano, Naoki Kawashima, Kenichi Kurihara, Yoshiki Matsuda, Seiji Miyashita, Hiroshi Nakagawa, Mikio Nakahara, Hidetoshi Nishimori, Masayuki Ohzeki, Hans de Raedt, Per Arne Rikvold, Issei Sato, Sei Suzuki, Eric Vincent, and Yoshihisa Yamamoto for their valuable comments. S.T. acknowledges Keisuke Fujii, Yoshifumi Nakada, and Takahiro Sagawa for their useful discussion during the lecture. S.T. is partly supported by Grand-in-Aid for JSPS Fellows (23-7601). R.T. is partly supported financially by National Institute for Materials Science (NIMS). The computation in the present work was performed on computers at the Suprecomputer Center, Institute for Solid State Physics, University of Tokyo. ## References 1. S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, Science 220 , 671 (1983). 2. S. Kirkpatrick, J. Stat. Phys. 34 , 975 (1984). 3. S. Geman and D. Geman, IEEE Transactions on Pattern Analysis and Machine Intelligence 6 , 721 (1984). 4. A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll, Chem. Phys. Lett. 219 , 343 (1994). 5. T. Kadowaki and H. Nishimori, Phys. Rev. E 58 , 5355 (1998). 6. J. Brooke, D. Bitko, T. F. Rosenbaum, and G. Aeppli, Science 284 , 779 (1999). 7. E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292 , 472 (2001). 8. G. E. Santoro, R. Martoˇ n´ ak, E. Tosatti, and R. Car, Science 295 , 2427 (2002). 9. A. Das and B. K. Chakrabarti, Quantum Annealing and Related Optimization Methods (Springer, Heidelberg, 2005). 10. A. Das and B. K. Chakrabarti, Rev. Mod. Phys. 80 , 1061 (2008). 11. M. Ohzeki and H. Nishimori, J. Comp. and Theor. Nanoscience 8 , 963 (2011). 12. K. Kurihara, S. Tanaka, and S. Miyashita, Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence (2009). 13. I. Sato, K. Kurihara, S. Tanaka, H. Nakagawa, and S. Miyashita, Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence (2009). 14. S. Tanaka, R. Tamura, I. Sato, and K. Kurihara, to appear in Kinki University Quantum Computing Series: 'Summer School on Diversities in Quantum Computation/Information' . 15. C. Jarzynski, Phys. Rev. Lett. 78 , 2690 (1997). 16. C. Jarzynski, Phys. Rev. E 56 , 5018 (1997). 17. M. Ohzeki, Phys. Rev. Lett. 105 , 050401 (2010). 18. E. Ising, Z. Phys. 31 , 253 (1925). 19. L. Onsager, Phys. Rev. 65 , 117 (1944). 20. M. Blume, Phys. Rev. 141 , 517 (1966). 21. H. W. Capel, Phys. Lett. 23 , 327 (1966). 22. J. Tobochnik, Phys. Rev. B 26 , 6201 (1982). 23. M. S. S. Challa and D. P. Landau, Phys. Rev. B 33 , 437 (1986). 24. R. B. Potts, Proc. Cambridge Philos. Soc. 48 , 106 (1952). 25. F. Y. Wu, Rev. Mod. Phys. 54 , 235 (1982). 26. T. Ohtsuka, J. Phys. Soc. Jpn. 16 , 1549 (1961). 27. M. Rayl, O. E. Vilches, and J. C. Wheatley, Phys. Rev. 165 , 698 (1968). 28. K. ˆ Ono, M. Shinohara, A. Ito, N. Sakai, and M. Suenaga, Phys. Rev. Lett. 24 , 770 (1970). 29. N. Achiwa, J. Phys. Soc. Jpn. 27 , 561 (1969). 30. M. Mekata and K. Adachi, J. Phys. Soc. Jpn. 44 , 806 (1978). 31. A. H. Cooke, D. T. Edmonds, F. R. McKim, and W. P. Wolf, Proc. Roy. Soc. London Ser. A 252 , 246 (1959). 32. A. H. Cooke, D. T. Edmonds, C. B. P. Finn, and W. P. Wolf, Proc. Roy. Soc. London Ser. A 306 , 313 (1968). 33. A. H. Cooke, D. T. Edmonds, C. B. P. Finn, and W. P. Wolf, Proc. Roy. Soc. London Ser. A 306 , 335 (1968). 34. K. Takeda, M. Matsuura, S. Matsukawa, Y. Ajiro, and T. Haseda, Proc. 12th Int. Conf. Low Temp. Phys., Kyoto 803 (1970). 35. K. Takeda, S. Matsukawa, and T. Haseda, J. Phys. Soc. Jpn. 30 , 1330 (1971). 36. B. N. Figgis, M. Gerloch, and R. Mason, Acta. Crystallogr. 17 , 506 (1964). 37. R. F. Wielinga, H. W. J. Blote, J. A. Roest, and W. J. Huiskamp, Physica 34 , 223 (1967). 38. K. W. Mess, E. Lagendijk, D. A. Curtis, and W. J. Huiskamp, Physica 34 , 126 (1967). 39. G. R. Hoy and F. de S. Barros, Phys. Rev. 139 , A929 (1965). 40. M. Matsuura, H. W. J. Blote, and W. J. Huiskamp, Physica 50 , 444 (1970). 41. R. D. Pierce and S. A. Friedberg, Phys. Rev. B 3 , 934 (1971). 42. K. Takeda and S. Matsukawa, J. Phys. Soc. Jpn. 30 , 887 (1971). 43. E. Stryjewski and N. Giordano, Adv. Phys. 26 , 487 (1977). 44. D. J. Breed, K. Gilijamse, and A. R. Miedema, Physica 45 , 205 (1969). 45. K. ˆ Ono, A. Ito, and T. Fujita, J. Phys. Soc. Jpn. 19 , 2119 (1964). 46. R. J. Birgeneau, W. B. Yelon, E. Cohen, and J. Makovsky, Phys. Rev. B 5 , 2607 (1972). 47. J. C. Wright, H. W. Moos, J. H. Colwell, B. W. Magnum, and D. D. Thornton, Phys. Rev. B 3 , 843 (1971). 48. G. T. Rado, Phys. Rev. Lett. 23 , 644 (1969). 49. W. Scharenberg and G. Will, Int. J. Magnetism 1 , 277 (1971). 50. H. Fuess, A. Kallel, and F. Tch´ eou, Solid State Commun. 9 , 1949 (1971). 51. M. Ball, M. J. M. Leask, W. P. Wolf, and A. F. G. Wyatt, J. Appl. Phys. 34 , 1104 (1963). 52. J. C. Norvell, W. P. Wolf, L. M. Corliss, J. M. Hastings, and R. Nathans, Phys. Rev. 186 , 557 (1969). 53. J. C. Norvell, W. P. Wolf, L. M. Corliss, J. M. Hastings, and R. Nathans, Phys. Rev. 186 , 567 (1969). 54. G. A. Baker, Jr., Phys. Rev. 129 , 99 (1963). 55. M. F. Sykes, D. L. Hunter, D. S. McKenzie, and B. R. Heap, J. Phys. A: Gen. Phys. 5 , 667 (1972). 56. J. W. Stout and E. Catalano, J. Chem. Phys. 23 , 2013 (1955). 57. C. Domb and A. R. Miedema, Progress in low Temperature Physics, Vol. 4, edited by C. J. Gorter (North-Holland, Amsterdam, 1964). 58. G. K. Wertheim and D. N. E. Buchanan, Phys. Rev. 161 , 478 (1967). 59. Y. Shapira, Phys. Rev. B 2 , 2725 (1970). 60. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, 2000). 61. M. Nakahara and T. Ohmi, Quantum Computing: From Linear Algebra to Physical Realizations (Taylor & Francis, London, 2008). 62. D. G. Cory, A. F. Fahmy, and T. F. Havel, Proc. Natl. Acad. Sci. USA 94 , 1634 (1997). 63. D. G. Cory, M. D. Price, W. Maas, E. Knill, R. Laflamme, W. H. Zurek, T. F. Havel, and S. S. Somaroo, Phys. Rev. Lett. 81 , 2152 (1998). 64. N. A. Gershenfeld and I. L. Chuang, Science 275 , 350 (1997). 65. I. L. Chuang, L. M. K. Vandersypen, X. Zhou, D. W. Leung, and S. Lloyd, Nature 393 , 143 (1998). 66. J. A. Jones and M. Mosca, J. Chem. Phys. 109 , 1648 (1998). 67. E. Knill, I. Chuang, and R. Laflamme, Phys. Rev. A 57 , 3348 (1998). 68. R. Laflamme, E. Knill, W. H. Zurek, P. Catasti, and S. V. S. Mariappan, Phil. Trans. R. Soc. Lond. A 356 , 1941 (1998). 69. J. A. Jones and M. Mosca, Phys. Rev. Lett. 83 , 1050 (1999). 70. M. D. Price, S. S. Somaroo, A. E. Dunlop, T. F. Havel, and D. G. Cory, Phys. Rev. A 60 , 2777 (1999). 71. L. M. K. Vandersypen, C. S. Yannoni, M. H. Sherwood, and I. L. Chuang, Phys. Rev. Lett. 83 , 3085 (1999). 72. L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, R. Cleve, and I. L. Chuang, Phys. Rev. Lett. 85 , 5452 (2000). 73. L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, M. H. Sherwood, and I. L. Chuang, Nature (London) 414 , 883 (2001). 74. M. Nakahara, Y. Kondo, K. Hata, and S. Tanimura, Phys. Rev. A 70 , 052319 (2004). 75. Y. Kondo, J. Phys. Soc. Jpn. 76 , 104004 (2007). 76. H. Suwa and S. Todo, Phys. Rev. Lett. 105 , 120603 (2010). 77. H. Suwa and S. Todo, arXiv :1106.3562. 78. R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58 , 86 (1987). 79. U. Wolff, Phys. Rev. Lett. 62 , 361 (1989). 80. K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65 , 1604 (1996). 81. N. Kawashima and K. Harada, J. Phys. Soc. Jpn. 73 , 1379 (2004). 82. T. Nakamura, Phys. Rev. Lett. 101 , 210602 (2008). 83. S. Morita, S. Suzuki, and T. Nakamura, Phys. Rev. E 79 , 065701(R) (2009). 84. H. F. Trotter, Proc. Am. Math. Soc. 10 , 545 (1959). 85. M. Suzuki, Prog. Theor. Phys. 56 , 1454 (1976). 86. T. Kadowaki, Ph. D thesis, Tokyo Institute of Technology (1998). 87. K. Tanaka and T. Horiguchi, Electronics and Communications in Japan, Part 3: Fundamental Electronic Science 83 , 84 (2000). 88. K. Tanaka and T. Horiguchi, Interdisciplinary Information Science 8 , 33 (2002). 89. H. Attias, Proceedings of the 15th Conference on Uncertainly in Artificial Intelligence 21 (1999). 90. L. Landau, Phys. Z. Sowjetunion 2 , 46 (1932). 91. C. Zener, Proc. R. Soc. London Ser. A 137 , 696 (1932). 92. E. C. G. St¨ uckelberg, Helv. Phys. Acta 5 , 369 (1932). 93. N. Rosen and C. Zener, Phys. Rev. 40 , 502 (1932). 94. B. K. Chakrabarti, A. Dutta, and P. Sen, Quantum Ising Phases and Transitions in Transverse Ising Models (Springer Verlag, Berlin, 1996). 95. G. T. Trammel, J. Appl. Phys. 31 , 362S (1960). 96. A. H. Cooke, D. T. Edmonds, C. B. P. Finn, and W. P. Wolf, J. Phys. Soc. Jpn. 17 , Suppl. B1 481 (1962). 97. J. W. Stout and R. C. Chisolm, J. Chem. Phys. 36 , 979 (1962). 98. V. L. Moruzzi and D. T. Teaney, Sol. State. Comm. 1 , 127 (1963). 99. A. Narath and J. E. Schriber, J. Appl. Phys. 37 , 1124 (1966). 100. R. F. Wielinga and W. J. Huiskamp, Physica 40 , 602 (1969). 101. W. P. Wolf, J. Phys. (Paris) 32 Suppl. C1 26 (1971). 102. W. Wu, B. Ellman, T. F. Rosenbaum, G. Aeppli, and D. H. Reich, Phys. Rev. Lett. 67 , 2076 (1991). 103. W. Wu, D. Bitko, T. F. Rosenbaum, and G. Aeppli, Phys. Rev. Lett. 71 , 1919 (1993). 104. D. H. Reich, B. Ellman, J. Yang, T. F. Rosenbaum, G. Aeppli, and D. P. Belanger, Phys. Rev. B 42 , 4631 (1990). 105. T. F. Rosenbaum, J. Phys.: Condens. Matter 8 , 9759 (1996). 106. D. H. Reich, T. F. Rosenbaum, G. Aeppli, and H. Guggenheim, Phys. Rev. B 34 , 4956 (1986). 107. J. A. Mydosh, Spin Glasses: An Experimental Introduction (Taylor & Francis, London, 1993). 108. P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59 , 381 (1987). 109. R. Martoˇ n´ ak, G. E. Santoro, and E. Tosatti, Phys. Rev. E 70 , 057701 (2004). 110. S. Tanaka and S. Miyashita, J. Phys.: Condens. Matter 19 , 145256 (2007). 111. H. Takayama and K. Hukushima, J. Phys. Soc. Jpn. 76 , 013702 (2007). 112. S. Tanaka and S. Miyashita, J. Phys. Soc. Jpn. 76 , 103001 (2007). 113. S. Miyashita, S. Tanaka, and M. Hirano, J. Phys. Soc. Jpn. 76 , 083001 (2007). 114. S. Tanaka and S. Miyashita, J. Phys. Soc. Jpn. 78 , 084002 (2009). 115. S. Tanaka and S. Miyashita, Phys. Rev. E 81 , 051138 (2010), Virtual Journal of Quantum Information 10 , (2010). 116. S. Tanaka and R. Tamura, J. Phys.: Conf. Ser. 320 , 012025 (2011). 117. H. Nishimori and G. Ortiz, Elements of Phase Transitions and Critical Phenomena (Oxford Univ Press, Oxford, 2010). 118. N. Ito, M. Taiji, and M. Suzuki, J. Phys. Soc. Jpn. 56 , 4218 (1987). 119. T. W. B. Kibble, J. Phys. A 9 , 1387 (1976). 120. T. W. B. Kibble, Phys. Rep. 67 , 183 (1980). 121. W. H. Zurek, Nature (London) 317 , 505 (1985). 122. B. Damski, Phys. Rev. Lett. 95 , 035701 (2005). 123. W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95 , 105701 (2005). 124. V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble, M. Krusius, Y. G. Makhlin, B. Placais, G. E. Volovik, and W. Xu, Nature (London) 382 , 334 (1996). 125. V. B. Eltsov, T. W. B. Kibble, M. Krusius, V. M. H. Ruutu, and G. E. Volovik, Phys. Rev. Lett. 85 , 4739 (2000). 126. H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. A 76 , 043613 (2007). 127. C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature (London) 455 , 948 (2008). 128. S. Suzuki, J. Stat. Mech. P03032 (2009). 129. S. Suzuki, J. Phys.: Conf. Ser. 302 , 012046 (2011). 130. R. J. Glauber, J. Math. Phys. 4 , 294 (1963). 131. R. Shankar and G. Murthy, Phys. Rev. B 36 , 536 (1987). 132. D. S. Fisher, Phys. Rev. B 51 , 6411 (1995). 133. J. Dziarmaga, Phys. Rev. B 74 , 064416 (2006). 134. G. Biroli, L. F. Cugliandolo, and A. Sicilia, Phys. Rev. E 81 , 050101(R) (2010). 135. G. Toulouse, Commun. Phys. (London) 2 , 115 (1977). 136. R. Liebmann, Statistical Mechanics of Periodic Frustrated Ising Systems (Springer-Verlag, Berlin/Heidelberg, GmbH, Heidelberg, 1986). 137. H. Kawamura, J. Phys.: Condens. Matter 10 , 4707 (1998). 138. H. T. Diep (ed.), Frustrated Spin Systems (World Scientific, Singapore, 2005). 139. S. Tanaka and S. Miyashita, Prog. Theor. Phys. Suppl. 157 , 34 (2005). 140. S. Tanaka and S. Miyashita, J. Phys. Soc. Jpn. 76 , 103001 (2007). 141. R. Tamura and N. Kawashima, J. Phys. Soc. Jpn. 77 , 103002 (2008). 142. S. Tanaka and S. Miyashita, J. Phys. Soc. Jpn. 78 , 084002 (2009). 143. R. Tamura and N. Kawashima, J. Phys. Soc. Jpn. 80 , 074008 (2011). 144. R. Tamura, N. Kawashima, T. Yamamoto, C. Tassel, and H. Kageyama, Phys. Rev. B 84 , 214408 (2011). 145. K. Husimi and I. Syozi, Prog. Theor. Phys. 5 , 177 (1950). 146. R. M. F. Houtappel, Physica 16 , 425 (1950). 147. G. H. Wannier, Phys. Rev. 79 , 357 (1950). 148. G. H. Wannier, Phys. Rev. B 7 , 5017 (1973). 149. K. Kano and S. Naya, Prog. Theor. Phys. 10 , 158 (1953). 150. Y. Matsuda, H. Nishimori, and H. G. Katzgraber, J. Phys.: Conf. Ser. 143 , 012003 (2009). 151. Y. Matsuda, H. Nishimori, and H. G. Katzgraber, New J. Phys. 11 , 073021 (2009). 152. S. Tanaka, M. Hirano, and S. Miyashita, Lecture Note in Physics 'Quantum Quenching, Annealing, and Computation' (Springer) 802 , 215 (2010). 153. S. Tanaka, to appear in proceedings of Kinki University Quantum Computing Series: 'Symposium on Quantum Information and Quantum Computing' (2011). 154. S. Tanaka and R. Tamura, in preparation . 155. E. H. Fradkin and T. P. Eggarter, Phys. Rev. A 14 , 495 (1976). 156. S. Miyashita, Prog. Theor. Phys. 69 , 714 (1983). 157. H. Kitatani, S. Miyashita, and M. Suzuki, Phys. Lett. 108A , 45 (1985). 158. H. Kitatani, S. Miyashita, and M. Suzuki, J. Phys. Soc. Jpn. 55 , 865 (1986). 159. P. Azaria, H. T. Diep, and H. Giacomini, Phys. Rev. Lett. 59 , 1629 (1987). 160. S. Miyashita and E. Vincent, Eur. Phys. J. B 22 , 203 (2001). 161. R. Tamura, S. Tanaka, and N. Kawashima, Prog. Theor. Phys. 124 , 381 (2010). 162. S. Tanaka and R. Tamura, and N. Kawashima, J. Phys.: Conf. Ser. 297 , 012022 (2011).