Classical simulation and theory of quantum annealing in a thermal environment
CClassical simulation and theory of quantum annealing in a thermal environment
Hiroki Oshiyama, ∗ Sei Suzuki, † and Naokazu Shibata ‡ Department of Physics, Tohoku University, Sendai 980-8578, Japan § Department of Liberal Arts, Saitama Medical University, Moroyama, Saitama 350-0495, Japan
We study quantum annealing in the quantum Ising model coupled to a thermal environment. When the speedof quantum annealing is sufficiently slow, the system evolves following the instantaneous thermal equilibrium.This quasistatic and isothermal evolution, however, fails near the end of annealing because the relaxation timegrows infinitely, therefore yielding excess energy from the thermal equilibrium. We develop a phenomenologicaltheory based on this picture and derive a scaling relation of the excess energy after annealing. The theoreticalresults are numerically confirmed using a novel non-Markovian method that we recently proposed based on apath-integral representation of the reduced density matrix and the infinite time evolving block decimation. Inaddition, we discuss the crossover between the adiabatic and quasistatic regimes and propose experiments onthe D-Wave quantum annealer.
Introduction. —The quantum annealing (QA) device manufactured by D-Wave Systems has made an immense impact not only in thephysics community but also in the industrial community witha hope of developing quantum computers and simulators [1–9]. It is known that this device carries out QA imperfectly inthe sense that the system embedded in this device is affectedby its environment [1, 10, 11]. This fact raises issues regard-ing QA dynamics in a thermal environment [2, 12–18].QA was proposed as a quantum mechanical algorithm tosolve combinatorial optimization problems [19–24]. Theproblem to be solved is encoded in an Ising Hamiltonian suchthat the solution is given by its ground state. The original algo-rithm is based on the quantum adiabatic time evolution froma known trivial ground state of an initial Hamiltonian to theunknown ground state of the Ising Hamiltonian [25].However, a system in a quantum device cannot be free fromenvironmental effects. Studying QA in a thermal environmentis beneficial not only for QA devices but also to understandnonequilibrium statistical mechanics [16, 26–37].A plausible picture of QA in the presence of a thermal en-vironment is quasistatic and isothermal evolution, in which asystem evolves maintaining a thermal equilibrium state at thetemperature of its environment. This picture should be validwhen the QA duration is much longer than the relaxation timeof the system. Previous studies , based a system-bath couplingrealistic in the D-Wave quantum annealer [11, 15, 38, 39],have suggested that the relaxation time increases dramaticallyas the transverse field is reduced. Because of this increase, thequasistatic and isothermal evolution should fail near the endof annealing, and result in a final state with an effective tem-perature higher than that of the environment. Even though thispicture has only been studied in small-sized systems, the scal-ings of the physical quantities expected in the thermodynamiclimit have not yet been studied. In this letter, we develop aphenomenological theory and derive a novel scaling relationof the energy after QA. This scaling is valid when the tem-perature of the environment is sufficiently high and/or QA issufficiently slow.To study the QA of a system with an experimentally realis-tic system-bath coupling, we employ a novel numerical non- Markovian method proposed by the present authors in Ref.[40]. This method makes use of the discrete-time path integralfor a dissipative system [41] and the infinite time-evolvingblock decimation (iTEBD) algorithm [42], which enable thecomputation of the reduced density matrix (RDM) in and outof equilibrium of the translationally invariant quantum Isingchain in the thermodynamic limit. We verify the theoreti-cal consequences on QA in a thermal environment using thismethod.
Model. — We consider the dissipative quantum Ising chain(DQIC) described by the Hamiltonian H ( s ) = H S ( s )+ H B + H SB , where H S ( s ) represents the system Hamiltonian givenby the quantum Ising chain, H S ( s ) = A ( s ) H TF + B ( s ) H I , (1)with H TF = − P Nj =1 ˆ σ xj and H I = − P N − j =1 ˆ σ zj ˆ σ zj +1 . Here ˆ σ xj and ˆ σ zj denote the Pauli matrices for the site j , N is thenumber of sites, and s is a parameter ranging from 0 to 1. Theschedule functions, A ( s ) and B ( s ) , are assumed to be A ( s ) = (1 − s ) α , B ( s ) = s, (2)where an exponent α > in A ( s ) represents how the trans-verse field goes to zero at the end of annealing. The bathHamiltonian is represented by the collection of harmonic os-cillators, H B = P Nj =1 P k ω k ˆ b † j,k ˆ b j,k , where ˆ b j,k and ˆ b † j,k arethe annihilation and creation operators, respectively, of the bo-son for the site j and mode k with the frequency ω k . We usethe unit ~ = 1 throughout this letter. As for the interactionbetween the system and the bath, we assume the Caldeira–Leggett model [43, 44] for dissipative superconductor fluxqubits given by H SB = N X j =1 ˆ σ zj X k λ k (cid:16) ˆ b † j,k + ˆ b j,k (cid:17) , (3)where λ k is a coupling constant for the mode k . We assumethe Ohmic spectral density for the bath modes, J ( ω ) = X k λ k δ ( ω − ω k ) = η ωe − ω/ω c , (4) a r X i v : . [ qu a n t - ph ] F e b where η is the dimensionless coupling constant and ω c is thecut-off frequency of the spectrum of the bath, which is chosento be larger than the bath temperature. When we mention QA,we consider the time evolution with time t from t = 0 to t a by the Hamiltonian H ( t/t a ) .The dynamics of the spin system is specified by the RDMdefined by tracing out the bosonic degrees of freedom fromthe density matrix ρ ( t ) of the full system, ρ S ( t ) ≡ Tr B ρ ( t ) = Tr B (cid:2) U ( t ) ρ (0) U † ( t ) (cid:3) , (5)where Tr B(S) stands for the trace with respect to the boson(spin) degrees of freedom, U ( t ) is the time evolution operatorof the full system and ρ (0) is an initial density matrix. Weassume that ρ (0) is the direct product of the ground state of H S (0) denoted by | ψ (0) i and the thermal equilibrium stateof H B at the temperature T B : ρ (0) = | ψ (0) i h ψ (0) | ⊗ e − H B /T B /Z B , where Z B is the partition function of H B . Werefer to T B as the bath temperature. We chose the Boltzmannconstant k B to be the temperature unit throughout this letter.The spin state in the instantaneous thermal equilibrium at s and temperature T is given by ρ eqS ( s, T ) ≡ Tr B [ e − H ( s ) /T ] /Z ( s, T ) , (6)where Z ( s, T ) is the partition function of the full system. Wedefine the Gibbs state of H I as ρ eqI ( T ) ≡ e − H I /T / Tr S [ e − H I /T ] . (7)Note that ρ eqS ( s, T ) at s = 1 reduces to ρ eqI ( T ) because H S (1) commutes with H SB and Tr B e − ( H S + H SB ) /T is independentof σ zj , which is shown by introducing new boson operators ˜ b j,k ≡ b j,k + λ k σ zj /ω k for all j and k . Non-Markovian iTEBD. — We focus on a time-dependentstate and outline the numerical method [40] used to computeEq. (5). The application to the equilibrium RDM in Eq. (6) isstraightforward.Let us apply the Trotter decomposition [45, 46] with a stepsize ∆ t = t/M and the Trotter number M to U ( t ) in Eq. (5),and perform the Gaussian integral with respect to the bosonicdegrees of freedom. The resulting discrete-time path integralformula of the RDM is given by h σ ( M ) | ρ S ( t ) | σ ( M +1) i = X { σ ( l ) j = ± } l = M,M +1 e i S + S infl , (8)where σ ( l ) j is the Ising-spin variable at the site j and the time t l is defined as t l = (cid:26) l ∆ t (0 ≤ l ≤ M )(2 M + 1 − l )∆ t ( M + 1 ≤ l ≤ M + 1) , (9)and | σ ( l ) i denotes the eigenstate of σ zj with the eigenvalue σ ( l ) j [43, 47]. S denotes the action of the isolated spin system. The influence action S infl induced by coupling to the bathis given by S infl = N X j =1 | t l − t m | <τ c X l>m κ l,m σ ( l ) j σ ( m ) j , (10)where κ l,m = ∆ t Z ∞ dωJ ( ω ) cosh[ ω/ (2 T B ) − iω ( t l − t m )]sinh[ ω/ (2 T B )] . (11)Note that τ c in Eq. (10) is the memory time cut-off introducedto reduce the computational cost.The key idea of our method is to represent the part of exp S infl associated with a site j in terms of a matrix prod-uct state (MPS) as follows: exp | t l − t m | <τ c X l>m κ l,m σ ( l ) j σ ( m ) j ! ≈ χ t X { µ j,l } φ ( j, S (0) j µ j, φ ( j, S (1) j µ j, ,µ j, φ ( j, S (2) j µ j, ,µ j, · · · φ ( j,M ) S ( M ) j µ j,M − , (12)where S ( l ) j ≡ ( σ ( l ) j , σ (2 M +1 − l ) j ) denotes the composite vari-able and χ t is the bond dimension which controls the preci-sion of the approximation in this MPS representation. UsingEq. (12) in Eq. (8), we obtain a tensor network representationfor the RDM: h σ ( M ) | ρ S ( t ) | σ ( M +1) i ≈ (13) X { S ( l ) j ,µ j,l } l = M e i S N Y j =1 φ ( j, S (0) µ j, h M − Y l =1 φ ( j,l ) S ( l ) j µ j,l − ,µ j,l i φ ( j,M ) S ( M ) j µ j,M − . Having obtained this tensor network representation, theiTEBD algorithm can be applied to implement the sum withrespect to { S ( l ) j , µ j,l } l = M and compute local quantities, tak-ing N → ∞ and using the translational invariance in space[42]. Phenomenological theory. — Let us assume a finite bathtemperature T B > . In the limiting case of t a → ∞ , QA inthis thermal environment leads to the quasistatic and isother-mal process. Accordingly, the final state of the spin systemis described by ρ eqS (1 , T B ) = ρ eqI ( T B ) . When t a is finite, thespin system approximately maintains thermal equilibrium aslong as the relaxation time of the spin system is shorter thanthe annealing time scale. However, as we elaborate below, therelaxation time grows infinitely with s → . Therefore, thequasistatic and isothermal evolution must fail before QA ends,and the spin state is expected to be frozen at a time t ∗ = s ∗ t a .We refer to t ∗ or s ∗ as the freezing time. To develop a scal-ing theory for the freezing time, we employ the quasistatic-freezing approximation as follows. The quasistatic-freezingapproximation assumes that the spin state is frozen when thechanging rate with t of the instantaneous relaxation time of . . . . s = t / t a − − . − . − . E n e r g y p e r s p i n s ∗ (a) T B = 1 η = 0 . t a = 200 ⟨ H S ( t / t a ) ⟩ t ⟨ H S ( s ) ⟩ eq s,T B . T − − − D K L ( T ) T ∗ (b) T B = 1 η = 0 . t a = 200 t a . E e x c ( t a ) p e r s p i n (c) T B = 1 η = 0 . α = 310 . . α . . . . b (d) T B = 1 η = 0 . η = 0 . Eq . (20) FIG. 1. (a) Energy expectation values, h H S ( t/t a ) i t and h H S ( s ) i eq s,T B , per spin of the time-dependent state ρ S ( t ) and the instantaneous thermalequilibrium state ρ eqS ( s, T B ) , respectively, as functions of the rescaled time s for t a = 200 and η = 0 . at T B = 1 . We fixed α = 1 . Thedashed line and the solid vertical line indicate N h H I i t a t/t a and s ∗ ≡ T B /T ∗ , respectively, where T ∗ is determined by the minimization ofthe Kullback–Leibler (KL) divergence. (b) KL divergence D KL ( T ) between the final state after QA and the Gibbs state of H I with temperature T . See the main text for a detailed definition. (c) Excess energy E exc per spin from the thermal expectation value after QA as a function of t a for various α . Lines indicate the best power-law fits E exc = at − ba to the data for t a > with the fitting parameters a and b . (d) Exponent b as a function of α . The numerical results (symbols) are compared to the theoretical prediction shown by the solid line. The parameters used inthe numerical simulations are ω c = 5 , τ c = 10 , ∆ t = 0 . and N → ∞ . The bond dimensions are up to 128. the spin system exceeds unity. Writing the instantaneous re-laxation time at s = t/t a as τ rel ( s ) , the freezing time is thendetermined by ˙ τ rel ( s ∗ ) = 1 , (14)where the dot denotes differentiation by t . τ rel ( s ) is now es-timated from the transition rate γ ( s ) . Using Fermi’s goldenrule, the latter is given, up to an s -independent factor, as γ l,m ( s ) ∝ η | h ψ l ( s ) | X i ˆ σ zi | ψ m ( s ) i | (15)where | ψ l ( s ) i and | ψ m ( s ) i are the l -th and m -th eigenstatesof H S ( s ) , respectively. For sufficiently large t a , the freezingtime s ∗ must be close to 1. Accordingly, the perturbation ex-pansion for the eigenstates can be evoked around s = 1 , whichyields τ rel ( s ) ≈ γ l,m ( s ) − ∼ η − A ( s ) − ∼ η − (1 − s ) − α . (16)Using this in Eq. (14), the scaling relation of s ∗ is obtained asfollows: (1 − s ∗ ) ∼ ( ηt a ) − / (2 α +1) . (17)Now, the quasistatic-freezing approximation implies that theRDM after the freezing time is approximately replaced by thatof the instantaneous thermal equilibrium at s = s ∗ , namely, ρ S ( t ) ≈ ρ eqS ( s ∗ , T B ) for t > s ∗ t a . Moreover, ρ eqS ( s ∗ , T B ) is approximated as ρ eqS ( s ∗ , T B ) = ρ eqI ( T B /s ) + O (1 − s ∗ ) α .Therefore, neglecting the O (1 − s ∗ ) α and O ((1 − s ) α (1 − s ∗ ) α ) terms for α > , the energy of the spin system for t > s ∗ t a is estimated to be h H S ( t/t a ) i t ≈ B ( t/t a ) h H I i eqI ,T B /s ∗ , (18)where h·i t and h·i eqI ,T represent the expectation values withrespect to ρ S ( t ) and ρ eqI ( T ) , respectively. Therefore, the energy of the spin system approaches the thermal expecta-tion value of H I at the temperature T B /s ∗ in a linear fash-ion with s as s → . For generic α > , h H S (1) i t a ≈h H S (1) i eqI ,T B + c ( T B )(1 − s ∗ ) min { , α } , where c ( T ) is a func-tion of T which is independent of s ∗ . Applying Eq. (17), theexcess energy of the final state is obtained as E exc ≡ h H S (1) i t a − h H S (1) i eqI ,T B ∼ ( η t a ) − b , (19)with b = min { , α } / ( α + 1) . (20)Note that the excess energy decays the fastest when α = .Equations (19) and (20) are valid for DQIC in any dimensionand any lattice. Numerical results. — Figure 1(a) shows the energy expec-tation value per site of the time-dependent state during QAand that of the instantaneous thermal equilibrium as func-tions of the rescaled time s . After the initial relaxation, thesystem maintains thermal equilibrium until a certain time s ∗ ,when the quasistatic and isotheral evolution fails and the en-ergy deviates upwards from that of the instantaneous equi-librium state. This behavior is perfectly consistent with thequasistatic-freezing picture mentioned above. To evaluate thefreezing time s ∗ and identify the final energy h H I i t a , we fo-cus on the Kullback-Leibler (KL) divergence D KL of the finalstate and the Boltzmann distribution of H I as a measure ofthe distance between the two. Because this quantity is notaccessible for the RDMs of the entire spin system when us-ing our method, we instead consider the RDMs of eight spinsgiven by ρ ≡ Tr ¯8 ρ S ( t a ) and ρ eq8 ( T ) ≡ Tr ¯8 ρ eqI ( T ) to define D KL ( T ) ≡ Tr [ ρ (log ρ − log ρ eq8 ( T ))] , where Tr and Tr ¯8 denote the trace with respect to the eight adjacent spins andthe other spins, respectively. We show D KL ( T ) in Fig. 1(b). D KL ( T ) has a sharp minimum at a certain T labeled T ∗ . Thisimplies that the RDM after QA is approximated by the Gibbs . . . . η . t a E e x c p e r s p i n t a = 50 t a = 100 t a = 200 t a = 500 t a = 1000 0 1 η . ( η t a ) E e x c FIG. 2. Excess energy per spin after QA scaled by t − / as a func-tion of η for T B = 1 , α = 1 , and various t a ranging from 50 to1000. With increasing t a , the data collapse into a single nonmono-tonic curve, which implies E exc ∼ t − / for large t a . The insetshows the data rescaled by ( ηt a ) − / . The constancy of the datawith large t a near η = 0 corresponds to Eq. (19). The parameters inthe numerical simulations are the same as those in Fig. 1. state of H I with the temperature T ∗ . In addition, as shown inFig. 1(a), the curve of h H S ( t/t a ) i t is indistinguishable fromthe line of s h H I i I ,T ∗ near s = 1 . Assuming T ∗ = T /s ∗ ,this result implies Eq. (18) and that s ∗ determined by T /T ∗ isconsistent with the freezing time when the quasistatic evolu-tion fails (see the vertical line in Fig. 1(a)). Figure 1(c) showsthe excess energy as a function of the annealing time t a for η = 0 . and T B = 1 . It can be seen that the excess energydecays as a power law for large t a with an exponent denotedby b that depends on α . Figure 1(d) shows the α dependenceof the exponent b . There is excellent agreement between thenumerical results and the theoretical prediction.Figure 2 shows the η -dependence of E exc scaled by t / for T B = 1 , α = 1 , and various t a . It can be seen that E exc is nonmonotonic with respect to η . The decreasing be-havior of E exc with increasing η in the weak coupling regimeis consistent with Eq. (19), while its increasing behavior inthe strong-coupling regime ( η & . ) is not described bythe phenomenological theory mentioned above. This failureof the theory arises from the perturbative argument used forthe relaxation time in Eq. (15). The existence of the optimalstrength in the system-bath coupling to reduce E exc is first re-vealed by our numerical method based on a non-perturbativeformulation. Note that the scaling of E exc by t a is valid evenin the strong-coupling situation.So far, we have focused on slow QA in a thermal envi-ronment with a finite temperature and have discussed con-sequences of freezing near the end of annealing. Here, wecomment on two situations where the dynamics is governedby a quantum phase transition at zero temperature, assumingthe absence of a thermal phase transition at finite temperature.The first is a weak-coupling and short-annealing-time situa-tion. When the system-bath coupling is sufficiently weak, i.e. . . . . s = t / t a − − . − . − . − . E n e r g y p e r s p i n T = 1 T = 0 T B = 1 t a = 200 t a = 1 η = 0 . η = 0 . η = 0 FIG. 3. Energy expectation values per spin as functions of therescaled time s . The solid line shows the energy of the instanta-neous ground state of H S ( s ) . The filled symbols show the energy ofthe time-dependent state for t a = 1 of the closed system ( η = 0 ).The empty symbols denote energies of the time-dependent states ofthe dissipative system with η = 0 . and T B = 1 for t a = 1 and t a = 200 . We fixed α = 1 . The parameters in the numerical simula-tions are the same as those in Fig. 1. η (cid:28) , QA drives the spin system in the same way as a closedsystem as long as t a is not large, as demonstrated by the cir-cles in Fig. 3. This insensibility to the environment gives riseto the Kibble–Zurek scaling (KZS) [48–50] of the residual en-ergy to the ground state after QA, which is a manifestation ofa quantum phase transition. For larger t a , the system is ther-malized and the final state gains an excess energy as shownby the squares in Fig. 3. This crossover is accompanied by anincrease in the residual energy over a certain parameter range[51]. The second is a medium-coupling and low-temperaturesituation, where the dynamics is governed by a quantum phasetransition of the dissipative system at T = 0 . In this case, KZSwith a modified exponent [40] is observed. When the temper-ature is not low and/or t a is much larger, however, the spinsystem is not influenced by a quantum phase transition andthe quasistatic-freezing picture is valid because the time scaleof QA is beyond the characteristic time in the quantum criticalregion [51]. A recent experimental study suggested that sys-tems realized in the D-Wave device should be in a situationwith a medium η and a low T [36]. Therefore, if one per-forms experiments with still longer t a or higher T , the scalingof the excess energy given by Eqs. (19) and (20) should beobserved. Summary. — We developed a phenomenological theory forQA in a thermal environment based on the quasistatic-freezingapproximation. The theoretical results were numerically con-firmed using a novel non-Markovian iTEBD method. The the-oretical and numerical findings presented here will be benefi-cial in designing and evaluating QA devices.The authors acknowledge H. Nishimori and Y. Susa forvaluable discussions, and Y. Bando M. Ohzeki, F. J. G´omez-Ruiz, A. del Campo, and D. A. Lidar for collaboration on arelated experimental project. ∗ [email protected] † [email protected] ‡ [email protected] § Present address: Graduate School of Information Sciences, To-hoku University, Sendai 980-8578, Japan[1] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze,N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk,E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizin-sky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom,E. Tolkacheva, C. J. Truncik, S. Uchaikin, J. Wang, B. Wilson,and G. Rose, Nature , 194 (2011).[2] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A.Lidar, J. M. Martinis, and M. Troyer, Nature Physics , 218(2014).[3] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush,V. Smelyanskiy, J. Martinis, and H. Neven, Phys. Rev. X ,031015 (2016).[4] T. Albash and D. A. Lidar, Phys. Rev. X , 031016 (2018).[5] A. D. King, J. Carrasquilla, J. Raymond, I. Ozfidan, E. An-driyash, A. Berkley, M. Reis, T. Lanting, R. Harris, F. Altomare,K. Boothby, P. I. Bunyk, C. Enderud, A. Fr´echette, E. Hoskin-son, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich, Y. Sato,A. Y. Smirnov, L. J. Swenson, M. H. Volkmann, J. Whittaker,J. Yao, E. Ladizinsky, M. W. Johnson, J. Hilton, and M. H.Amin, Nature , 456 (2018).[6] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H.Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang,E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky,T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh,I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov,L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao,Science , 162 (2018).[7] A. D. King, J. Raymond, T. Lanting, S. V. Isakov, M. Mohseni,G. Poulin-Lamarre, S. Ejtemaee, W. Bernoudy, I. Ozfidan,A. Y. Smirnov, M. Reis, F. Altomare, M. Babcock, C. Baron,A. J. Berkley, K. Boothby, P. I. Bunyk, H. Christiani, C. En-derud, B. Evert, R. Harris, E. Hoskinson, S. Huang, K. Jooya,A. Khodabandelou, N. Ladizinsky, R. Li, P. A. Lott, A. J. R.MacDonald, D. Marsden, G. Marsden, T. Medina, R. Molavi,R. Neufeld, M. Norouzpour, T. Oh, I. Pavlov, I. Perminov,T. Prescott, C. Rich, Y. Sato, B. Sheldan, G. Sterling, L. J.Swenson, N. Tsai, M. H. Volkmann, J. D. Whittaker, W. Wilkin-son, J. Yao, H. Neven, J. P. Hilton, E. Ladizinsky, M. W.Johnson, and M. H. Amin, “Scaling advantage in quan-tum simulation of geometrically frustrated magnets,” (2019),arXiv:1911.03446 [quant-ph].[8] H. Ushijima-Mwesigwa, C. F. A. Negre, and S. M.Mniszewski, in Proceedings of the Second International Work-shop on Post Moores Era Supercomputing , PMES’17 (Associa-tion for Computing Machinery, New York, NY, USA, 2017) pp.22–29, event-place: Denver, CO, USA.[9] A. Teplukhin, B. K. Kendrick, S. Tretiak, and P. A. Dub, Sci-entific Reports , 20753 (2020).[10] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A.Lidar, Nature Communications , 2067 (2013).[11] J. Marshall, E. G. Rieffel, and I. Hen, Phys. Rev. Applied ,064025 (2017).[12] M. S. Sarandy and D. A. Lidar, Phys. Rev. Lett. , 250503 (2005).[13] M. H. S. Amin, P. J. Love, and C. J. S. Truncik, Phys. Rev. Lett. , 060503 (2008).[14] N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris,F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple,P. Chavez, F. Cioata, T. Cirip, P. deBuen, M. Drew-Brook,C. Enderud, S. Gildert, F. Hamze, J. P. Hilton, E. Hoskinson,K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Lanting, T. Ma-hon, R. Neufeld, T. Oh, I. Perminov, C. Petroff, A. Przybysz,C. Rich, P. Spear, A. Tcaciuc, M. C. Thom, E. Tolkacheva,S. Uchaikin, J. Wang, A. B. Wilson, Z. Merali, and G. Rose,Nature Communications , 1903 (2013).[15] M. H. Amin, Phys. Rev. A , 052323 (2015).[16] L. Arceci, S. Barbarino, R. Fazio, and G. E. Santoro, Phys.Rev. B , 054301 (2017).[17] V. N. Smelyanskiy, D. Venturelli, A. Perdomo-Ortiz, S. Knysh,and M. I. Dykman, Phys. Rev. Lett. , 066802 (2017).[18] A. Y. Smirnov and M. H. Amin, New Journal of Physics ,103037 (2018).[19] T. Kadowaki and H. Nishimori, Phys. Rev. E , 5355 (1998).[20] A. Finnila, M. Gomez, C. Sebenik, C. Stenson, and J. Doll,Chemical Physics Letters , 343 (1994).[21] A. Das and B. K. Chakrabarti, eds., Quantum annealing andrelated optimization methods (Springer, Berlin, Heidelberg,2005) lecture Notes in Physics, Vol. 679.[22] A. Das and B. K. Chakrabarti, Rev. Mod. Phys. , 1061 (2008).[23] T. Albash and D. A. Lidar, Rev. Mod. Phys. , 015002 (2018).[24] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, andW. D. Oliver, Reports on Progress in Physics , 054401(2020).[25] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, andD. Preda, Science , 472 (2001).[26] D. Patan`e, A. Silva, L. Amico, R. Fazio, and G. E. Santoro,Phys. Rev. Lett. , 175701 (2008).[27] D. Patan`e, L. Amico, A. Silva, R. Fazio, and G. E. Santoro,Phys. Rev. B , 024302 (2009).[28] S. Yin, P. Mai, and F. Zhong, Phys. Rev. B , 094108 (2014).[29] M.-J. Hwang, R. Puebla, and M. B. Plenio, Phys. Rev. Lett. , 180404 (2015).[30] P. Nalbach, S. Vishveshwara, and A. A. Clerk, Phys. Rev. B ,014306 (2015).[31] A. Dutta, A. Rahmani, and A. del Campo, Phys. Rev. Lett. ,080402 (2016).[32] M. Keck, S. Montangero, G. E. Santoro, R. Fazio, andD. Rossini, New Journal of Physics , 113029 (2017).[33] P. Weinberg, M. Tylutki, J. M. R¨onkk¨o, J. Westerholm, J. A.˚Astr¨om, P. Manninen, P. T¨orm¨a, and A. W. Sandvik, Phys. Rev.Lett. , 090502 (2020).[34] R. Puebla, A. Smirne, S. F. Huelga, and M. B. Plenio, Phys.Rev. Lett. , 230602 (2020).[35] D. Rossini and E. Vicari, Phys. Rev. Research , 023211 (2020).[36] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J.G´omez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, andH. Nishimori, Phys. Rev. Research , 033369 (2020).[37] S. Bandyopadhyay, S. Bhattacharjee, and A. Dutta, Phys. Rev.B , 104307 (2020).[38] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, New Journalof Physics , 123016 (2012).[39] J. Marshall, D. Venturelli, I. Hen, and E. G. Rieffel, Phys. Rev.Applied , 044083 (2019).[40] H. Oshiyama, N. Shibata, and S. Suzuki, J. Phys. Soc. Japan , 1 (2020).[41] N. Makri, Chemical Physics Letters , 435 (1992).[42] R. Or´us and G. Vidal, Phys. Rev. B , 155117 (2008). [43] A. O. Caldeira and A. J. Leggett, Ann. Phys. (N. Y). , 374(1983).[44] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher,A. Garg, and W. Zwerger, Rev. Mod. Phys. , 1 (1987).[45] H. F. Trotter, Proceedings of the American Mathematical Soci-ety , 545 (1959).[46] M. Suzuki, Progress of Theoretical Physics , 1454 (1976). [47] D. E. Makarov and N. Makri, Chemical Physics Letters ,482 (1994).[48] T. W. B. Kibble, Journal of Physics A: Mathematical and Gen-eral , 1387 (1976).[49] W. H. Zurek, Nature , 505 (1985).[50] J. Dziarmaga, Phys. Rev. Lett. , 245701 (2005).[51] See Supplementary Materials. upplementary Meterials: Classical simulation and theory of quantum annealing in a thermalenvironment Hiroki Oshiyama, Sei Suzuki, and Naokazu Shibata Department of Physics, Tohoku University, Sendai 980-8578, Japan Department of Liberal Arts, Saitama Medical University, Moroyama, Saitama 350-0495, Japan
In this Supplementary Materials, we supplement numerical results on the residual energy to the ground state after QA, whichis defined by E res ( t a ) ≡ h H S (1) i t a − E g (1) where the ground energy E g (1) of H S (1) is given by E g (1) = − .At first, we briefly review the Kibble-Zurek scaling (KZS) of the residual energy in a closed system. When η = 0 , the spinsystem encounters a quantum phase transition during QA with the time-dependent Hamiltonian H S ( t ) changing from H TF to H I . The system loses adiabaticity near a quantum critical point and thereby ends up with a residual energy and topologicaldefects. The residual energy is smaller for slower QA, i.e., larger t a . This mechanism for the residual energy due to timeevolution across a continuous phase transition is called as the Kibble-Zurek mechanism (KZM) [ ? ? ]. The scaling of theresidual energy with respect to t a is universal and given by E res ∼ t − dν/ ( nν +1) a [ ? ? ], where ν and z are the correlation-lengthand the dynamical critical exponents of the associated quantum critical point. This scaling of the residual energy is the KZS. Incase of the transverse Ising chain, one has ν = z = 1 and hence KZS leads to E res ∼ t − / a [ ? ].Before proceeding to numerical results, we comment on KZS in the dissipative transverse Ising chain. Consider the case with η > and T = 0 . This dissipative spin system also involves a quantum phase transition. The quantum critical point forms aquantum critical line s c ( η ) in the s − η plane. An earlier study using the quantum Monte-Carlo simulation have estimated criticalexponents on this critical line as ν ≈ . and z ≈ . for η > [ ? ]. These exponents suggest a modification of the KZSexponent b of E res ∼ t − ba as b QMC ≈ . from / of the closed system. The present authors have studied KZS of the samemodel by means of iTEBD in the previous work [ ? ]. The results imply that the exponent of KZS decreases gradually from 0.5to 0.25 with increasing η from 0 to 0.7. The qualitative tendency that the exponent declines in the presence of the environment isshared by the theoretical prediction with quantum Monte-Carlo data and numerical experiments. The quantitative inconsistencyis attributed to the smallness of t a in numerical experiments.We now move on to numerical results for finite temperatures. Figure S1(a) shows results for η = 0 . as a weak couplingsituation. When the temperature is sufficiently low, the residual energy decreases monotonically with t a . We expect that itshould converge to the thermal expectation value, h H S (1) i eqI ,T B + 1 , for t a → ∞ . For higher T B , the residual energy variesnonmonotonically with t a . This nonmonotonic curve can be divided in three regions: (i) the KZS for sufficiently small t a where E res decays following the power law as that of the closed system, (ii) the anti-KZS for the middle range of t a where E res deviatesfrom the power law and increases with t a , and (iii) the quasistatic-freezing region for sufficiently large t a where E res decreasestoward h H S (1) i eqI ,T B + 1 . Note that the excess energy, defined by E exc ≡ h H S (1) i t a − h H S (1) i eq I,T B follows Eq. (19) in (iii) aswe mentioned in the main text. This nonmonotonic anti-KZS behavior of the residual energy has been discussed in literatures[ ? ? ]. It arises as a consequence of a large KZS region such that the residual energy at the end of the KZS region is below thethermal expectation value.We next focus on η = 0 . as a medium coupling situation. Figure S2 compares the energy of a time-dependent state drivenby QA with that of the instantaneous equilibrium state at T = 0 . The spin system equilibrates with the bath in an early time ofannealing, keeps equilibrium approximately for a while, and then drop out of equilibrium. The result for η = 0 . is in contrastto the one for η = 0 . where the energy of the time-dependent state is mostly indistinguishable from that in a closed system.The upward deviation of the curve of the time-dependent state at T = 0 from that of the equilibrium state comes from KZM,namely, the growth of the relaxation time near a quantum critical point and the failure of a quasistatic evolution. In Fig. S3,the residual energy at T B = 0 is observed to decay, following a power law with the exponent b ≈ . . As mentioned above,although this exponent is larger than b QMC ≈ . predicted from the quantum Monte-Carlo simulation, a power law confirmsKZS of the residual energy in a dissipative system. Raising temperature of the bath in Fig. S3, one finds that the range of t a forKZS shrinks and eventually vanishes. The anti-KZS is not observed here. These results are in sharp contrast to the results forthe weak coupling. At finite temperatures, QA does not cross a quantum critical point. The separation from a quantum criticalpoint is larger and the relaxation time is shorter for higher temperature. When the relaxation time maximized near a quantumcritical point at finite temperature is shorter than the time scale of QA, KZM does not work and a quasistatic evolution survivesuntil freezing near the end of QA. This is why KZS disappears for higher temperatures. The absence of the anti-KZS is becausethe residual energy does not drop enough with KZS. a r X i v : . [ qu a n t - ph ] F e b t a . E r e s p e r s p i n η = 0 . Thermal value( T = 1 ) T = 1 . T = 0 . T = 0 . T = 0 . T = 0 . η = 0 FIG. S1. Residual energy after QA as a function of annealing time t a for a weak coupling with η = 0 . and various temperatures. The resultof the closed system is shown by the solid line for comparison. The horizontal line indicates the thermal expectation value at T = 1 . . . . . s = t / t a − − . − . − . E n e r g y p e r s p i n η = 0 . t a = 10 T B = 0 η = 0 . η = 0 . η = 0 FIG. S2. Energies per spin as function of the rescaled time s of time-dependent states with t a = 10 at T B = 0 and of equilibrium states at T = 0 . Squares and triangles represent results of time-dependent states for η = 0 . and η = 0 . , respectively, while circles are for theclosed system. The solid and dashed lines show the instantaneous energies per spin of the equilibrium states of the closed system at T = 0 and the dissipative system at T = 0 with η = 0 . . t a . E r e s p e r s p i n η = 0 . Thermal value( T = 1 ) T = 1 . T = 0 . T = 0 . T = 0 . T = 0 . η = 0 FIG. S3. Residual energy after QA as a function of annealing time t a for a medium coupling with η = 0 . and various temperatures. Theresult of the closed system is shown by the solid line for comparison. The horizontal line indicates the thermal expectation value at T = 1= 1