Path-Integral Quantum Monte Carlo simulation with Open-Boundary Conditions
Zhang Jiang, Vadim N. Smelyanskiy, Sergio Boixo, Hartmut Neven
PPath-Integral Quantum Monte Carlo Simulation with Open-Boundary Conditions
Zhang Jiang,
1, 2
Vadim N. Smelyanskiy, Sergio Boixo, and Hartmut Neven QuAIL, NASA Ames Research Center, Moffett Field, California 94035, USA SGT Inc., 7701 Greenbelt Rd., Suite 400, Greenbelt, MD 20770 Google, Venice, CA 90291, USA (Dated: October 24, 2017)The tunneling decay event of a metastable state in a fully connected quantum spin model can besimulated efficiently by path-integral quantum Monte Carlo (QMC) [Isakov et al. , Phys. Rev. Lett. , 180402 (2016).]. This is because the exponential scaling with the number of spins of the ther-mally assisted quantum tunneling rate and the Kramers escape rate of QMC are identical [Jiang etal. , Phys. Rev. A , 012322 (2017).], a result of a dominant instantonic tunneling path. In Isakov etal. , it was also conjectured that the escape rate in open-boundary QMC is quadratically larger thanthat of conventional periodic-boundary QMC; therefore, open-boundary QMC might be used as apowerful tool to solve combinatorial optimization problems. The intuition behind this conjectureis that the action of the instanton in open-boundary QMC is a half of that in periodic-boundaryQMC. Here, we show that this simple intuition—although very useful in interpreting some numer-ical results—deviates from the actual situation in several ways. Using a fully connected quantumspin model, we derive a set of conditions on the positions and momenta of the end points of theinstanton, which remove the extra degrees of freedom due to open boundaries. In comparison, thehalf-instanton conjecture incorrectly sets the momenta at the end points to zero. We also foundthat the instantons in open-boundary QMC correspond to quantum tunneling events in the sym-metric subspace (maximum total angular momentum) at all temperatures, whereas the instantonsin periodic-boundary QMC typically lie in subspaces with lower total angular momenta at finitetemperatures. This leads to a lesser-than-quadratic speedup at finite temperatures. The resultsprovide useful insights in utilizing open-boundary QMC to solve hard optimization problems. Wealso outline the generalization of the instantonic tunneling method to many-qubit systems withoutpermutation symmetry using spin-coherent-state path integrals. I. INTRODUCTION
Computational hard combinatorial optimization prob-lems can be mapped to spin glass models in statisticalphysics [1]. The energy landscapes of the correspondingproblem Hamiltonians H P are rough and posses largenumbers of spurious local minima. Conventional opti-mization strategies, such as simulated annealing, exploitthermal over–the–barrier transitions to escape from thelocal minima. In quantum annealing (QA) [2–7], tunnel-ing offers additional paths for systems to go to low-energystates [8]. In an archetypical example of QA, the stateevolution is determined by a time-dependent Hamilto-nian H ( t ) = H P − Γ( t ) (cid:80) j σ xj , where σ xj is the Pauli- x operator of the j th spin, and Γ( t ) slowly interpolatesfrom a large value to 0. All low-energy eigenstates of H ( t ) are localized in the vicinity of the minima of H P atsufficiently small values of Γ [9]. The order of the ener-gies of two minima changes at an avoided crossing, andquantum tunneling between the two minima gives rise tothe energy gap ∆. The quantum state remains pure inthe absence of an environment, and it closely follows theinstantaneous ground state of H ( t ) when Γ( t ) is changedslowly enough. Such a process is called adiabatic quan-tum computation [3], where the state dynamics can bedescribed as a cascade of Landau–Zener transitions atthe avoided crossings [4].Interaction with an environment causes relaxation, andthis can suppress tunneling. On the other hand, the environment also gives rise to thermal excitations tohigher energy levels from where the system can tun-nel faster [10]. This is called thermally assisted tunnel-ing [11], and it was recently discussed in QA with fluxqubits [12, 13]. We assume that the system has a free-energy minimum associated with a thermodynamicallymetastable state and that the incoherent tunneling decayrate W of this state is much smaller than the smallest re-laxation rate towards the quasiequilibrium distributionin the domain associated with this state. Then the in-coherent tunneling decay rate is W ∝ ∆ / ( (cid:126) γ ), where γ (cid:29) ∆ / (cid:126) is a largest relaxation rate (typically, dephasingrate).Path-integral quantum Monte Carl (QMC) is a clas-sical Markov chain Monte Carlo algorithm to calculatethermodynamic properties of quantum systems; it is themost efficient algorithm to compute exponentially smallgaps at first order phase transitions [5, 14, 15] and is of-ten used to simulate QA in spin glasses [4, 5, 16, 17]. Theequivalence between the exponential scaling of the QMCtransition rate and the QA tunneling rate in a fully con-nected quantum spin model was established numericallyat the effective “zero-temperature” limit [18], where theeffect of thermal excitations on the tunneling rate canbe neglected. A detailed theoretical calculation on thisequivalence based on an instanton technique was givenin [19] at finite temperatures. Therein, the instantonmethod was used to demonstrate the identical exponen-tial scaling of the QA tunelling rate and the QMC es-cape rate for a bit-symmetric cost function with a thin, a r X i v : . [ qu a n t - ph ] O c t high energy barrier (“Hamming weight with a spike”).Crosson and Harrow also considered the same cost func-tion, and they derived a bound on the mixing time of theunderlying Markov chain [20]; their results also showedthat QMC takes polynomial time to find the solution inthe regime with no tunneling.In QMC, paths are being sampled instead of the ac-tual configurations of the system. A path consists of asequence of replicas of the system, where changes in con-figurations between adjacent replicas are penalized en-ergetically. The Kramers escape event in the stochasticprocess describes the transition from a local minimumto the global minimum, which is dominated by a single“transition state” (a saddle point) that the system needsto reach in order to make an escape from the metastablestate. As we show (see also Refs. [18, 19]), this transitionstate corresponds to an instanton, and the change in freeenergy needed to reach this state is the same as that inthe quantum tunneling case. This explains the equiva-lence in the exponential scaling of the QMC transitionrate and thermally assisted quantum tunneling rate.In Ref. [18], the authors found numerically that aquadratic speedup may be achieved in QMC escapeevents by using open-boundary condition (OBC) insteadof periodic-boundary condition (PBC); it was conjec-tured that the speedup is due to an escape path of ahalf instanton. Physical quantities usually cannot be pre-dicted correctly by using open-boundary QMC; however,it can be used as a physics-inspired classical algorithm forcombinatorial optimization problems. In a recent numer-ical study [21], it was found that the autocorrelation timein QMC simulations with OBC is significantly smallerthan that of conventional QMC methods.Because QMC samples paths instead of the physicalstates, it might have conserved quantities not present inthe physical system, such as the number of world lines(particles, magnetization), braiding, or winding num-bers [22, 23]. Topological obstructions can arise fromthese conserved quantities, which prevent QMC fromreaching equilibrium even though the quantum Hamil-tonian is gapped. Open-boundary QMC is not immuneto this, but it can escape many topological obstructionsthat trap periodic-boundary QMC. Recently, Andriyashand Amin [24] argued that topological obstructions canmake QMC less efficient than QA when there are mul-tiple tunneling paths between two minima. These ob-structions partially prevent the forming of instantons inperiodic-boundary QMC by creating a barrier in the mid-dle. However, open-boundary QMC is not encumberedby such obstructions.Here, we consider open-boundary QMC of a fully con-nected spin model with a bit-symmetric cost function.We show that the saddle point of the mean-field QMCfree-energy functional satisfies the same differential equa-tions as those in periodic-boundary QMC. We also de-rived a set of conditions on the positions and momentaof the end points of instantons in open-boundary QMC,which uniquely determine the instanton solution by re- moving the extra degrees of freedom due to OBC. In-terestingly, we found nonzero initial and final momentaof instanton in open-boundary QMC. This can be at-tributed to the extra entropic factors at the boundaries.The free-energy at the local minimum also correspondsto a nonstationary solution in the path-integral formal-ism due to the same entropic factors. Another some-what surprising fact about open-boundary QMC is thatthe optimal tunneling path lies in the symmetric sub-space (maximum total angular momentum) at all tem-peratures. In contrast, the optimal tunneling paths inperiodic-boundary QMC typically have total angular mo-mentum less than the maximum value at finite temper-atures. Our analytical results show that the conjecturein Ref. [18] is only approximately correct: The QMC es-cape rate is not always enhanced quadratically by usingopen-boundary conditions.In Sec. II we formally establish the free-energy func-tional for QMC with OBC and derive its saddle-pointequation. In Sec. III we compute the free-energy func-tional of the open-boundary QMC at the saddle pointusing the instanton approach. In Sec. IV we map theopen-boundary QMC free-energy functional to a quan-tum propagator and use the Wentzel-Kramers-Brillouin(WKB) approach to derive the corresponding “tunnelingrate” in the model. The result conforms with the QMCescape rate derived in preceding sections. II. SADDLE POINTS OF PATH INTEGRALQMC WITH OBC
We consider a fully connected quantum system of N spins [10, 25],ˆ H = − S x − N g (2 S z /N ) , S α = 12 N (cid:88) j =1 σ ( j ) α , (1)where Γ is the strength of the transverse field, σ ( j ) α isthe Pauli matrix of the j th spin with α = x, y, z , S α isthe α component of the total spin operator, and g is anarbitrary function.A classical method to calculate the partition func-tion Z = Tr e − β ˆ H is the path-integral QMC, where the“imaginary-time” interval [0 , β ] is sliced into R pieces.In the limit R → ∞ , we have Z = (cid:88) σ (0) ,...,σ ( R − (cid:104) σ (0) | e − ∆ ˆ H | σ (1) (cid:105)(cid:104) σ (1) | e − ∆ ˆ H · · · | σ ( R − (cid:105)(cid:104) σ ( R − | e − ∆ ˆ H | σ (0) (cid:105) , (2)where ∆ = β/R and σ ( τ ) = { σ j ( τ ) } Nj =1 . Equation (2)also takes the form Z = Tr e − βH QMC , (3)where H QMC is a Hamiltonian on R × N classical spins, H QMC = − J N (cid:88) j =1 R − (cid:88) τ =0 σ j ( τ ) σ j ( τ + 1) − NR R − (cid:88) τ =0 g [ m ( τ )] . (4)The quantity m ( τ ) in Eq. (4) is the total magnetization, m ( τ ) = 1 N N (cid:88) j =1 σ j ( τ ) , (5)and J is the effective coupling strength between adjacentreplicas, J = − β ln tanh (cid:0) Γ∆ (cid:1) ≥ . (6)When the coupling strength J = ∞ (Γ → m ( τ ) = m for τ = 0 , . . . , R −
1; the free-energy functional F then reduces to the classical case, F classical ( m ) = − g ( m ) − Q ( m ) β , (7)where Q is the binary entropy, Q ( m ) ≡ − m m − − m − m . (8)In Appendix A, we review the mean-field approachto path-integral QMC with periodic and open-boundaryconditions. Here, we discuss how to calculate the saddlepoint of the free-energy functional of the QMC Hamilto-nian (4) with open-boundary condition. The free-energyfunctional of open-boundary QMC takes the followingform in the continuous limit R → ∞ (see Appendix A): F OB = 1 β (cid:90) β (cid:16) m ( τ ) λ ( τ ) − g [ m ( τ )] (cid:17) dτ − β ln (cid:101) V ( λ ) , (9)where λ ( τ ) is a to-be-determined function of the imagi-nary time τ ∈ [0 , β ], and (cid:101) V ( λ ) = Tr (cid:0) ωK β, ( λ ) (cid:1) , ω = 1 + σ x = (cid:18) (cid:19) . (10)The propagator K β, corresponds to a spin-1/2 particleevolving in imaginary time under the action of the time-dependent magnetic field B ( τ ), K τ ,τ = T + e − (cid:82) τ τ dτH ( τ ) , (11) H ( τ ) = − B ( τ ) · σ , B ( τ ) = (cid:0) Γ , , λ ( τ ) (cid:1) , (12)where T + is the time-ordering operator, and σ =( σ x , σ y , σ z ) is the vector of Pauli matrices.The saddle-point conditions for the free-energy func-tional (9) are λ ( τ ) = g (cid:48) [ m ( τ )] , m ( τ ) = δ ln (cid:101) V ( λ ) δλ ( τ ) . (13) To solve these equations, we introduce a vector function m ( τ ) = (cid:0) m x ( τ ) , m y ( τ ) , m z ( τ ) (cid:1) corresponding to the fol-lowing expectation values of the Pauli operators, m ( τ ) = Tr (cid:0) ωK β,τ σ K τ, (cid:1) Tr (cid:0) ωK β, (cid:1) = δδ B ( τ ) ln Tr (cid:0) ωK β, (cid:1) , (14)where K τ ,τ is defined in Eq. (11) and B ( τ ) in Eq. (12).Differentiating Eq. (14) with respect to τ and using (12),we obtain d m dτ = Tr (cid:0) ωK β,τ [ H ( τ ) , σ ] K τ, (cid:1) Tr (cid:0) ωK β, (cid:1) , (15)which can be rewritten into the form d m dτ = 2 i B × m = − i ∂ H ( m ) ∂ m × m , (16)where H ( m ) = − Γ m x − g ( m z ) . (17)In component form, Eq. (16) becomes . m x = − λim y , (18) i . m y = − λm x + 2Γ m z , (19) . m z = 2Γ im y . (20)Equation (16) allows for two integrals of motion, H ( m ) = ε, (21) m ( τ ) · m ( τ ) = (cid:96) (22)where m ( τ ) · m ( τ ) ≡ m x ( τ ) + m y ( τ ) + m z ( τ ). Usingthese two integrals of motion, we have m x = − ε + g ( m z )Γ = (cid:112) (cid:96) − m z cosh p , (23) im y = (cid:112) (cid:96) − m z sinh p , (24)where p serves as the “conjugate momentum” to m z .Solving m y as a function of m z from Eqs. (23) and (24)and putting the result into Eq. (19), we have the instan-ton equation | . m z | = 2 (cid:112) [ ε + g ( m z )] − Γ ( (cid:96) − m z ) . (25)We still need to specify ε , (cid:96) , and the initial condition of m z to fully determine the instanton trajectory based onthe differential equation (25). Using the identity σ x ω = ωσ x = ω , we haveTr( ωK β σ x ) = Tr( ωK β ) , i Tr( ωK β σ y ) = Tr( ωK β σ z ) . (26)As a result, we have the following initial conditions: m x (0) = 1 , m z (0) = im y (0) = 12Γ . m z (0) . (27)Because (cid:96) is a constant of motion, its value can be deter-mined at τ = 0, (cid:96) = (cid:88) j = x,y,z m j (0) = 1 . (28)This means that the parameter (cid:96) is independent of β under OBC; the same condition can also be derived usingthe WKB approach [see Eq. (62)]. Similarly, we have m x ( β ) = 1 , m z ( β ) = − im y ( β ) = − . m z ( β ) . (29)From the energy conservation condition Eq. (21) and thecondition m x (0) = m x ( β ) = 1, we have g [ m z (0)] = g [ m z ( β )] . (30)This condition relates the position of the end points ofthe instanton trajectory under OBC. III. FREE ENERGY OF THE INSTANTON
To calculate the free energy from Eq. (9), it remainsto evaluate the quantity (cid:101) V = Tr( ωK β, ) = 2 (cid:104) + | K β, | + (cid:105) , (31)where | + (cid:105) = √ ( | (cid:105) + | (cid:105) ). We introduce the doublepropagator K τ, ⊗ K τ, acting on the original qubit anda replica qubit, which is generated by the Hamiltonian H (2)0 = H ⊗ I + I ⊗ H , (32)where H = − Γ σ x − g (cid:48) ( m z ) σ z is defined in Eq. (12). Withthe Bell basis | Φ + (cid:105) = √ ( | (cid:105) + | (cid:105) ) , | Φ − (cid:105) = √ ( | (cid:105) −| (cid:105) ) , | Ψ + (cid:105) = √ ( | (cid:105) + | (cid:105) ) , | Ψ − (cid:105) = √ ( | (cid:105) − | (cid:105) ),we have − H (2)0 | Φ + (cid:105) = 2Γ | Ψ + (cid:105) + 2 g (cid:48) ( m z ) | Φ − (cid:105) , (33) − H (2)0 | Φ − (cid:105) = 2 g (cid:48) ( m z ) | Φ + (cid:105) , (34) − H (2)0 | Ψ + (cid:105) = 2Γ | Φ + (cid:105) , (35) − H (2)0 | Ψ − (cid:105) = 0 . (36)Because the time evolution is closed in the symmetricsubspace of the two qubits, we have | ξ ( τ ) (cid:105) = K τ, ⊗ K τ, | ξ (0) (cid:105) = − ξ x ( τ ) | Φ − (cid:105) + iξ y ( τ ) | Φ + (cid:105) + ξ z ( τ ) | Ψ + (cid:105) , (37)where ξ x , iξ y , and ξ z take real values. According toEqs. (33)–(35), the state vector ξ = ( ξ x , ξ y , ξ z ) satisfiesthe linear differential equation d ξ dτ = − i ∂ H ( m ) ∂ m × ξ , (38) where m is determined by the instanton solutionEqs. (23)–(25). A known solution to Eq. (38) is the in-stanton solution ξ ( τ ) = m ( τ ). Denote ζ ( τ ) as the solu-tion to Eq. (38) with the initial conditions ζ x (0) = 0 , iζ y (0) = 1 √ , ζ z (0) = 1 √ . (39)The corresponding state to ζ x (0) defined in Eq. (37) is | ζ (0) (cid:105) = 1 √ (cid:0) | Φ + (cid:105) + | Ψ + (cid:105) (cid:1) = | + (cid:105) ⊗ | + (cid:105) . (40)Thus, the quantity in Eq. (31) can be expressed as (cid:104) + | K β | + (cid:105) = (cid:104) ζ (0) | K β ⊗ K β | ζ (0) (cid:105) = (cid:104) ζ (0) | ζ ( β ) (cid:105) , (41)where K β is a shorthand for K β, . To solve ζ ( β ), wedefine the real-valued symmetric bilinear form B (cid:0) ξ ( τ ) , η ( τ ) (cid:1) = ξ x ( τ ) η x ( τ ) + ξ y ( τ ) η y ( τ ) + ξ z ( τ ) η z ( τ ) , (42)where ξ ( τ ) and η ( τ ) are any two solutions to Eq. (38).Because the bilinear form is a constant of motion, wehave the identities from the conditions (27) and (39), B (cid:0) ζ ( τ ) , ζ ( τ ) (cid:1) = B (cid:0) ζ ( τ ) , m ( τ ) (cid:1) = 0 , (43)for any τ ∈ [0 , β ]. Consequently, we have − (cid:0) ζ y + ζ z (cid:1) m x = ζ y m y + 2 ζ y ζ z m y m z + ζ z m z , (44)Introducing (cid:37) = iζ y /ζ z , we have( (cid:37) − m x = − (cid:37) m y − i(cid:37) m y m z + m z . (45)The solution to the above quadratic equation is (cid:37) = m x − im y m z − m z , (46)where we use Eq. (22) to simplify things and the con-ditions (39) [which implies (cid:37) (0) = 1] to rule out theother solution of Eq. (45). Putting the condition (29)into Eq. (46), we have (cid:37) ( β ) = 1 + m z ( β ) − m z ( β ) . (47)Combining Eqs. (43) and (47), we have ζ x ( β ) = −√ κm z ( β ) , (48) iζ y ( β ) = κ (cid:2) m z ( β ) (cid:3)(cid:14) √ , (49) ζ z ( β ) = κ (cid:2) − m z ( β ) (cid:3)(cid:14) √ , (50)where κ is to be determined. Putting the definition of (cid:37) into Eq. (38), we have . ζ z ( τ ) = 2Γ (cid:37) ( τ ) ζ z ( τ ) , (51)which can be solved exactly, ζ z ( τ ) = ζ z (0) e (cid:82) τ (cid:37) ( τ ) dτ . (52)Combining the above expression and Eq. (50), we haveln κ = 2Γ (cid:90) β (cid:37) ( τ ) dτ − ln (cid:2) − m z ( β ) (cid:3) . (53)Putting Eq. (46) into the above integral, we haveΓ (cid:90) β (cid:37) ( τ ) dτ = (cid:90) β (cid:12)(cid:12) ε + g ( m z ) (cid:12)(cid:12) − m z . m z / − m z dτ = (cid:90) β (cid:12)(cid:12) ε + g ( m z ) (cid:12)(cid:12) − m z dτ + 14 ln (cid:0) − m z (cid:1)(cid:12)(cid:12)(cid:12) β , (54)where the relations m x = (cid:12)(cid:12) ε + g ( m z ) (cid:12)(cid:12) / Γ and im y = . m z /
2Γ are used. The quantity κ can thus be determinedexplicitly,ln κ = 2 I − (cid:16) ln (cid:2) − m z (0) (cid:3) + ln (cid:2) − m z ( β ) (cid:3)(cid:17) , (55)where the integral I is defined as I = (cid:90) β | ε + g ( m ) | − m dτ . (56)Noticing that κ = (cid:104) ζ (0) | ζ ( β ) (cid:105) and using Eqs. (31) and(41), we have (cid:101) V = Tr( ωK β ) = 2 √ κ . (57)The free energy Eq. (9) can thus be evaluated, βF OB = (cid:90) β (cid:2) m z g (cid:48) ( m z ) − g ( m z ) (cid:3) dτ − ln 2 − I + 14 (cid:16) ln (cid:2) − m z (0) (cid:3) + ln (cid:2) − m z ( β ) (cid:3)(cid:17) . (58) IV. WENTZEL-KRAMERS-BRILLOUIN (WKB)APPROACH
The partition function of the QMC Hamiltonian (3)with OBC takes the form in the limit of large number ofreplicas Z OB = Tr(Ω e − β ˆ H ) , (59)where ˆ H is defined in Eq. (1), and Ω is the matrix ofones; i.e., Ω jk = 1 for all j, k = 1 , , . . . , N . The rank ofthe matrix Ω is one, an it can be expressed asΩ = | s (cid:105)(cid:104) s | , (60) where | s (cid:105) = (cid:0) | (cid:105) + | (cid:105) (cid:1) ⊗ N is in the symmetric subspace.The partition function (59) can also be written as Z OB = (cid:104) s | e − β ˆ H | s (cid:105) . (61)Letting | m (cid:105) ≡ | N/ , mN/ (cid:105) denote the normalized statein the symmetric subspace of N spin-1/2 particles withtotal magnetization mN/
2, we have Z OB = (cid:88) m ,m (cid:104) s | m (cid:105)(cid:104) m | e − β ˆ H | m (cid:105)(cid:104) m | s (cid:105) . (62)One distinction of QMC with OBC as opposed to PBCfor the Hamiltonian (1) is that the discussion can be car-ried out in the symmetric subspace even at finite tem-peratures. The inner products (cid:104) m | s (cid:105) and (cid:104) s | m (cid:105) can beinterpreted as additional “entropic factors” to the parti-tion function. Indeed, we haveln (cid:104) m | s (cid:105) N = 1 N ln (cid:115) N ! N + ! N − ! (cid:39) Q ( m ) , (63)where N + = N (1 + m ) / N − = N (1 − m ) /
2, and Q ( m ) is the binary entropy defined in Eq. (8). We willalso need the derivative of the entropic factor, Q (cid:48) ( m ) ≡ dQ ( m ) dm = 12 ln 1 − m m . (64)The WKB “free energy” for OBC thus takes the form β F OB = βε + A − N (cid:16) ln (cid:104) m | s (cid:105) + ln (cid:104) s | m (cid:105) (cid:17) (65) (cid:39) βε + 12 (cid:0) A − Q ( m ) − Q ( m ) (cid:1) , (66)where ε is the energy. The WKB action A is given bythe integral A = (cid:90) m m p dm (67)= m p − m p − (cid:90) m m mp (cid:48) dm , (68)where p is the conjugate momentum of the magnetiza-tion m . The first two terms in Eq. (65) are due to thepropagator (cid:104) m | e − β ˆ H | m (cid:105) .A small change in energy of the path takes the form δε = − β (cid:90) β (cid:0) . m δp − . p δm (cid:1) dτ (69)= − δ A β + 12 β ( p δm − p δm ) . (70)Using Eqs. (66) and (70), we have δ F OB (cid:39) δε + 12 β (cid:0) δ A − δQ ( m ) − δQ ( m ) (cid:1) = 12 β (cid:16) [ p − Q (cid:48) ( m )] δm − [ p + Q (cid:48) ( m )] δm (cid:17) . (71)The boundary conditions of any saddle point of F OB canthen be derived, p + Q (cid:48) ( m ) = p − Q (cid:48) ( m ) = 0 . (72)We now show that these conditions are identical toEqs. (27) and (29). The WKB instanton solution is givenby [same as Eq. (23)] | p | = cosh − (cid:12)(cid:12) ε + g ( m ) (cid:12)(cid:12) Γ √ − m . (73)Using conditions (27) and (29), we have (cid:12)(cid:12) ε + g ( m ) (cid:12)(cid:12) / Γ = m x = 1 for the end points. Putting this into Eq. (73),we recover the results in Eq. (72), | p | = cosh − √ − m = 12 ln 1 + | m | − | m | = (cid:12)(cid:12) Q (cid:48) ( m ) (cid:12)(cid:12) . (74)For the instanton solution, the integral in Eq. (68)takes the form (cid:90) m m mp (cid:48) dm = (cid:90) m m − mg (cid:48) ( m ) + m − m (cid:12)(cid:12) ε + g ( m ) (cid:12)(cid:12)(cid:113) [ ε + g ( m )] − Γ (cid:0) − m (cid:1) | dm | (75)= 2 (cid:90) β (cid:16) − mg (cid:48) ( m ) + m − m (cid:12)(cid:12) ε + g ( m ) (cid:12)(cid:12)(cid:17) dτ (76)= 2 (cid:0) βε + I (cid:1) − (cid:90) β (cid:2) mg (cid:48) ( m ) − g ( m ) (cid:3) dτ , (77)where we use the instanton solution given in Eq. (25),and I is defined in Eq. (56). For m < p <
0, wehave m p = m m − m − m − Q ( m ) −
12 ln(1 − m ) + ln 2 . (79)For m > p <
0, we have m p = m − m − m m Q ( m ) + 12 ln(1 − m ) − ln 2 . (81)Putting the above results together, we have A = m p − m p − (cid:90) m m mp (cid:48) dm = Q ( m ) + Q ( m ) + 12 (cid:2) ln(1 − m ) + ln(1 − m ) (cid:3) − (cid:0) ln 2 + βε + I (cid:1) + 2 (cid:90) β (cid:2) mg (cid:48) ( m ) − g ( m ) (cid:3) dτ . (82) Putting Eq. (82) into (66), we have the effective WKB“free energy” β F OB = (cid:90) β (cid:2) mg (cid:48) ( m ) − g ( m ) (cid:3) dτ − ln 2 − I + 14 (cid:16) ln (cid:2) − m (cid:3) + ln (cid:2) − m (cid:3)(cid:17) , (83)which is identical to the QMC free energy F OB in Eq. (58).Another surprising fact about the QMC with OBC isthat the minimum of the free-energy functional is notachieved by a static solution independent of τ . This canbe understood by the definition of F OB in Eq. (83); while A and ε can be minimized simultaneously by an optimalstatic solution, the boundary terms are not minimizedby that solution. As a consequence, the free-energy func-tional is minimized by a nontrivial solution caused by theentropic factors at the boundaries. The end points of thisinstanton satisfy m = m and p = − p . To get rid offthis effect, one can add boundary terms (pinning fields)to cancel the extra entropic factors. V. SPIN-COHERENT-STATE QUANTUMMONTE CARLO
The instantonic analysis presented here allows thestudy of spin tunneling for a class of Hamiltonians, sym-metric with respect to permutations of individual spin-1/2 particles (qubits). It can potentially be generalizedto the case without permutation symmetry using pathintegrals over spin-coherent states. The action along theinstanton trajectory in imaginary time is A ( n ) = i (cid:126) N (cid:88) i =1 χ ( n j ) + (cid:90) β dτ H [ n ( τ ) , . . . , n N ( τ )] , (84)where n j ( τ ) is the Bloch unit vector of the j th spin andthe Berry phases of individual spins are defined as χ ( n j ) = (cid:90) β dτ (cid:2) − cos θ j ( τ ) (cid:3) ˙ φ j ( τ ) , (85)where θ j ( τ ) and φ j ( τ ) are the polar angle and azimuthalangle of the Bloch unit vector n j ( τ ), respectively. Thefunction H [ n ( τ ) , . . . , n N ( τ )] in (84) is obtained from thesystem Hamiltonian H by replacing the Pauli matricesof the j th spin by the corresponding components of theBloch unit vector n j ( τ ).Within the incoherent tunneling framework, the in-stanton trajectory with open-boundary condition con-nects the spin coherent state { n j (0) } Nj =1 correspondingto the maximum of the wavefunction near the local min-ima of the energy landscape H [ n ( τ ) , . . . , n N ( τ )] to thecoherent state { n j ( β ) } Nj =1 corresponding to the remotetail of the wavefunction on the other side of the bar-rier. Such a tunneling transition corresponds to corre-lated motions of individual spins in imaginary time satis-fying δ A /δ n j ( τ ) = 0 (to minimize the action). Therefore,the Bloch unit vector evolves according to the instantonequation, (cid:126) d n j ( τ ) dτ = n j ( τ ) × ∂H∂ n j ( τ ) . (86)We note that the first (Berry phase) term in (84) containsadditional factor i compared to the second term. There-fore, the y component of n j ( τ ) is complex for the instan-ton trajectory [cf. Eqs. (15)–(20) above]. The Bloch unitvector can be written in the form n j = (sin θ j cosh ϕ j , − i sin θ j sinh ϕ j , cos θ j ) , (87)corresponding to a purely imaginary azimuthal angle φ j ( τ ) = − iϕ j ( τ ). This substitution makes the Berryphase terms in Eq. (84) real along the instanton path.The values of H [ n ( τ ) , . . . , n N ( τ )] are also real due tothe fact that the Hamiltonian H is Hermitian. There-fore, despite the presence of the imaginary Berry phase inEq. (84), the instanton trajectory equations (86) involveonly real quantities after the substitution of Eq. (87).A similar situation also happens in the spin tunnelingstudied above via Eqs. (15)–(20). Finally, the tunnel-ing matrix element ∆ tunn is determined with logarithmicequivalence as ∆ tunn = B exp[ −A ( n ∗ )] , (88)where n ∗ is an instanton trajectory and the prefactor Bcan be obtained in terms of the functional determinantof the kernel δ A ( n ) at n ∗ .The variational method outlined above can be used tostudy the tunneling matrix elements in transverse-fieldspin glasses between the computational basis states sep-arated by large Hamming distances. These matrix ele-ments are usually exponentially small in the number ofspins N , and alternative methods involving exact diago-nalization are not feasible for n (cid:38) { θ j ( τ ) , ϕ j ( τ ) } Nj =1 in imag-inary time according to the probability functional ∝ exp[ −A ( n )]. The explicit form of the prefactor can beobtained using the results in Ref. [26]. The trajectoriesof two variables for each spin need to be simulated insuch Monte Carlo methods which avoids topological ob-structions such as those considered in Ref. [24]. VI. DISCUSSION AND SUMMARY
In QMC, quasi-equilibrium distributions of paths aredetermined by classical free-energy functionals. Transi-tions from a local minimum to the global minimum aredescribed by Kramers escape events, where the systemreaches a single “transition state” (a saddle point) be-fore making an escape from the metastable state. We found the transition state in open-boundary QMC simu-lations analytically for a fully connected spin model withbit-symmetric cost functions. The transition state cor-responds to an instanton, governed by the same differ-ential equation as in the PBC case [19]. We derive theinstanton equation using two different approaches. One isthe mean-field approach, where the free-energy functionalis expressed using a two-dimensional propagator. Thispropagator corresponds to a particle evolving in imag-inary time under the action of a time-dependent mag-netic field. The instanton equation can be derived byobserving two constants of motion. Another approachis to map the free energy of the open-boundary QMCto an imaginary-time quantum propagator (cid:104) s | e − β ˆ H | s (cid:105) ,where | s (cid:105) = (cid:0) | (cid:105) + | (cid:105) (cid:1) ⊗ N ; this contrasts to the PBCcase, where the free energy is related to Tr e − β ˆ H . Theinstanton solution that dominates the escape event canbe derived using the WKB approach in the large-spinpicture, which coincides with our mean-field results. Theinstanton in periodic-boundary QMC is determined bythe instanton equation and the period (inverse tempera-ture); however, one also has to optimize the positions andmomenta at both ends of an instanton in open-boundaryQMC. We found that the initial and final momenta ofinstantons in open-boundary QMC are nonzero; this isbecause the extra entropic factors introduced at the openboundaries. To cancel this side effect, one can add extrapotential terms at the end replicas to suppress configu-rations with higher entropy. At finite temperatures, theperiodic-boundary QMC escape event corresponds to aquantum tunneling path in a subspace with total angu-lar momentum less than the maximum value [19]. WithOBCs, however, this result no longer holds; we foundthat the QMC escape event always corresponds to tun-neling events in the symmetric subspace (maximum totalangular momentum). We also outlined the generalizationof the instantonic tunneling method to systems withoutpermutation symmetry using spin-coherent-state path in-tegrals. VII. ACKNOWLEDGMENTS
Z.J. would like to acknowledge enlightening and usefuldiscussions with Salvatore Mandr`a and Andre Petukhov.This work is supported by the NASA Advanced Explo-ration Systems program and NASA Ames Research Cen-ter. The research is based in part upon work supportedby the Office of the Director of National Intelligence(ODNI), Intelligence Advanced Research Projects Activ-ity (IARPA), via Interagency Umbrella Agreement No.IA1-1198. The views and conclusions contained hereinare those of the authors and should not be interpretedas necessarily representing the official policies or endorse-ments, either expressed or implied, of the ODNI, IARPA,or the U.S. Government. The U.S. Government is autho-rized to reproduce and distribute reprints for Governmen-tal purposes notwithstanding any copyright annotationthereon. [1] Y. Fu and P. W. Anderson, “Application of statisticalmechanics to NP-complete problems in combinatorial op-timisation,” J. Phys. A: Math. Gen. , 1605 (1986).[2] T. Kadowaki and H. Nishimori, “Quantum annealingin the transverse ising model,” Phys. Rev. E , 5355(1998); J. Brooke, D. Bitko, T. F. Rosenbaum, andG. Aeppli, “Quantum Annealing of a Disordered Mag-net,” Science , 779 (1999).[3] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lund-gren, and D. Preda, “A quantum adiabatic evolution al-gorithm applied to random instances of an np-completeproblem,” Science , 472 (2001).[4] G. E. Santoro, R. Martonak, E. Tosatti, and R. Car,“Theory of quantum annealing of an ising spin glass,”Science , 2427 (2002).[5] S. Boixo, T. F. Ronnow, S. V. Isakov, Z. Wang,D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer,“Evidence for quantum annealing with more than onehundred qubits,” Nature Physics , 218 (2014).[6] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov,D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer,“Defining and detecting quantum speedup,” Science ,420 (2014).[7] J. King, S. Yarkoni, M. M. Nevisi, J. P. Hilton, and C. C.McGeoch, “Benchmarking a quantum annealing proces-sor with the time-to-target metric,” arXiv:1508.05087.[8] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov,M. Dykman, V. S. Denchev, M. Amin, A. Smirnov,M. Mohseni, and H. Neven, “Computational Roleof Collective Tunneling in a Quantum Annealer,”arXiv:1411.4036 .[9] B. Altshuler, H. Krovi, and J. Roland, “Anderson lo-calization makes adiabatic quantum optimization fail,”Proceedings of the National Academy of Sciences ,12446 (2010).[10] K. Kechedzhi and V. N. Smelyanskiy, “Open-SystemQuantum Annealing in Mean-Field Models with Ex-ponential Degeneracy,” Physical Review X , 021028(2016).[11] A. I. Larkin and Y. N. Ovchinnikov, “Quantum tun-neling with dissipation,” J. Exp. Theor. Phys. , 382(1983) [Pis’ma Zh. Eksp. Teor. Fiz. 37, No. 7, 322 (1983)];D. A. Garanin and E. M. Chudnovsky, “Thermally ac-tivated resonant magnetization tunneling in molecularmagnets: Mn Ac and others,” Physical Review B ,11102 (1997); I. Affleck, “Quantum-Statistical Metasta-bility,” Phys. Rev. Lett. , 388 (1981).[12] M. H. S. Amin and D. V. Averin, “Macroscopic Reso-nant Tunneling in the Presence of Low Frequency Noise,”Phys. Rev. Lett. , 197001 (2008).[13] N. G. Dickson, M. W. Johnson, Amin, et al. , “Thermallyassisted quantum annealing of a 16-qubit problem,” NatCommun , 1903 (2013).[14] A. Young, S. Knysh, and V. Smelyanskiy, “First-orderphase transition in the quantum adiabatic algorithm,”Physical review letters , 020502 (2010).[15] I. Hen and A. Young, “Exponential complexity of thequantum adiabatic algorithm for certain satisfiabilityproblems,” Physical Review E , 061152 (2011).[16] R. Martoˇn´ak, G. E. Santoro, and E. Tosatti, “Quantumannealing by the path-integral monte carlo method: The two-dimensional random ising model,” Phys. Rev. B ,094203 (2002); “Quantum annealing of the traveling-salesman problem,” Physical Review E , 057701(2004); D. A. Battaglia, G. E. Santoro, and E. Tosatti,“Optimization by quantum annealing: Lessons from hardsatisfiability problems,” Phys. Rev. E , 066707 (2005);G. E. Santoro and E. Tosatti, “Optimization using quan-tum mechanics: quantum annealing through adiabaticevolution,” Journal of Physics A: Mathematical and Gen-eral , R393 (2006).[17] B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer,“Quantum versus classical annealing of ising spinglasses,” Science , 215 (2015).[18] S. V. Isakov, G. Mazzola, V. N. Smelyanskiy, Z. Jiang,S. Boixo, H. Neven, and M. Troyer, “UnderstandingQuantum Tunneling through Quantum Monte Carlo Sim-ulations,” Physical Review Letters , 180402 (2016).[19] Z. Jiang, V. N. Smelyanskiy, S. V. Isakov, S. Boixo,G. Mazzola, M. Troyer, and H. Neven, “Scaling anal-ysis and instantons for thermally assisted tunneling andquantum Monte Carlo simulations,” Physical Review A , 012322 (2017).[20] E. Crosson and A. W. Harrow, “Simulated QuantumAnnealing Can Be Exponentially Faster Than ClassicalSimulated Annealing,” in (2016) pp. 714–723.[21] G. Mazzola and M. Troyer, “Accelerated nuclearquantum effects sampling with open path integrals,”arXiv:1608.01262 (2016).[22] H. G. Evertz, “The Loop Algorithm,” Advances inPhysics , 1 (2003), arXiv: cond-mat/9707221.[23] M. B. Hastings, “Obstructions to classically simulatingthe quantum adiabatic algorithm.” Quantum Informa-tion & Computation , 1038 (2013).[24] E. Andriyash and M. H. Amin, “Can quantum MonteCarlo simulate quantum annealing?” arXiv:1703.09277(2017).[25] V. Bapst and G. Semerjian, “On quantum mean-fieldmodels and their quantum annealing,” Journal of Statis-tical Mechanics: Theory and Experiment , P06007(2012).[26] A. Garg, E. Kochetov, K.-S. Park, and M. Stone, “Spincoherent-state path integrals and the instanton calculus,”Journal of mathematical physics , 48 (2003).[27] T. J¨org, F. Krzakala, J. Kurchan, A. Maggs, and J. Pu-jos, “Energy gaps in quantum first-order mean-field-liketransitions: The problems that quantum annealing can-not solve,” EPL (Europhysics Letters) , 40004 (2010). Appendix A: Meanfield description of path-integralquantum Monte Carlo
The partition function of H QMC defined in Eq. (3) canbe calculated with a mean-field approach. First, we con-sider N noninteracting spins under the mean-field Hamil-tonian H MF = − NR R − (cid:88) τ =0 λ ( τ ) m ( τ ) − J N (cid:88) j =1 R − (cid:88) τ =0 σ j ( τ ) σ j ( τ + 1) , (A1)where λ ( τ ) is an “average” local magnetic field (mean-field) on the τ th replica slice and m ( τ ) is defined inEq. (5). The partition functional corresponds to H MF is Z MF ( λ ) = [ V ( λ )] N , (A2)where λ = (cid:0) λ (0) , . . . , λ ( R − (cid:1) . The function V ( λ ) isdefined as V ( λ ) = (cid:88) σ ,...,σ R − e β (cid:80) R − τ =0 (cid:0) λ ( τ ) σ ( τ ) /R + Jσ ( τ ) σ ( τ +1) (cid:1) = Tr (cid:16) L [ λ ( R − · · · L [ λ ( τ )] · · · L [ λ (0)] (cid:17) , (A3)where the transition matrix takes the form L ( λ ) = (cid:18) e βJ e − βJ e − βJ e βJ (cid:19) (cid:18) e βλ/R e − βλ/R (cid:19) . (A4)The trace in Eq. (A3) is a consequence of the periodic-boundary condition σ j (0) = σ j ( R −
1) typically used inQMC. By defining the propagator K τ ,τ = L [ λ ( τ )] · · · L [ λ ( τ )] , τ > τ , (A5)we have V ( λ ) = Tr K R − , ( λ ) . (A6)The expectation value of m ( τ ) can be calculated fromthe partition function m ( τ ) = RβN ∂ ln Z MF ( λ ) ∂λ ( τ )= Rβ ∂ ln V ( λ ) ∂λ ( τ ) = Tr (cid:2) K R − , τ σ z K τ − , (cid:3) V ( λ ) , (A7)where σ z is the Pauli- z matrix. The expectation valuesof the magnetization m (0) , . . . , m ( R −
1) are functions ofthe parameters λ (0) , . . . , λ ( R − m ( τ ) vanishes in the large N limit,we neglect the difference between m ( τ ) and m ( τ ), andEq. (A7) becomes m ( τ ) = Rβ ∂ ln V ( λ ) ∂λ ( τ ) . (A8)The free-energy functional corresponding to the Hamil-tonian (A1) is F MF = − βN ln Z MF ( λ ) = − β ln V ( λ ) . (A9) Consequently, the free-energy functional of the QMCHamiltonian (4) reads (see also Ref. [27]), F = 1 R R − (cid:88) τ =0 (cid:16) λ ( τ ) m ( τ ) − g [ m ( τ )] (cid:17) − β ln V ( λ ) , (A10)where the term λ ( τ ) m ( τ ) is introduced to cancel the cor-responding term in Eq. (A1). When the QMC free-energyfunctional F is minimized, we have R ∂ F ∂m ( τ ) = λ ( τ ) − g (cid:48) [ m ( τ )]+ R − (cid:88) τ =0 (cid:18) m ( τ ) ∂λ ( τ ) ∂m ( τ ) − Rβ ∂λ ( τ ) ∂m ( τ ) ∂ ln V ( λ ) ∂λ ( τ ) (cid:19) = λ ( τ ) − g (cid:48) [ m ( τ )] = 0 , (A11)where Eq. (A8) is used to show that the sum over τ equals to zero. Equations (A8) and (A11) determine thesaddle point of the free-energy functional (A10).To calculate the exponential scaling of the QMC escaperate using the instanton method, we will need to go thecontinuous limit R → ∞ . In this limit, the propagatordefined in Eq. (A5) corresponds to a spin-1/2 particleevolving in imaginary time under the action of the time-dependent magnetic field B ( τ ), K τ ,τ = T + e − (cid:82) τ τ dτH ( τ ) , (A12) H ( τ ) = − B ( τ ) · σ , B ( τ ) = (cid:0) Γ , , λ ( τ ) (cid:1) , (A13)where τ ∈ [0 , β ] denotes the imaginary time, T + is thetime-ordering operator, and σ = ( σ x , σ y , σ z ) is the vec-tor of Pauli matrices. We will replace the notations K τ ,τ ( λ ) and V ( λ ) with K τ ,τ ( λ ) and V ( λ ), because λ ( τ ) is a function in the continuous limit.Similar to Eqs. (A2) and (A6), the meanfield partitionfunctional with open-boundary condition takes the formin the continuous limit Z OBMF = [ (cid:101) V ( λ )] N , (cid:101) V ( λ ) = Tr (cid:0) ωK β, ( λ ) (cid:1) , (A14)where the matrix ω = 1 + σ x = (cid:18) (cid:19) . (A15)The matrix ω makes it possible to sum over configu-rations prohibited by the periodic-boundary condition,which is unique to open-boundary QMC. Similar toEq. (A10), the free-energy functional of open-boundaryQMC takes the form, F OB = 1 β (cid:90) β (cid:16) m ( τ ) λ ( τ ) − g [ m ( τ )] (cid:17) dτ − β ln (cid:101) V ( λ ) . (A16)The saddle-point conditions for open-boundary QMC aresimilar to those for periodic-boundary QMC; Eq. (A11)0remains the same, and one just need to replace V with (cid:101) V in Eq. (A8), λ ( τ ) = g (cid:48) [ m ( τ )] , m ( τ ) = δ ln (cid:101) V ( λ ) δλ ( τ ) ..