Quantum annealing of pure and random Ising chains coupled to a bosonic environment

We study quantum annealing of pure and random transverse Ising chains that are coupled to boson baths, using a novel numerical method based on the combination of the quasi-adiabatic propagation path integral (QUAPI) and the matrix product state (MPS) formalisms. We present numerical results on systems up to 64 spins and conclude that the baths disturb quantum annealing of pure and random Ising chains. The numerical method to compute the reduced density matrix is described in detail.


Introduction
Controlling the time evolution of a quantum many-body system has been an important issue for decades in condensed matter physics, quantum physics, statistical physics, information science, and engineering. Today's swell in the study of this issue stems partly from recent advances in quantum computing technology. In particular, the appearance of quantum annealing processor has stimulated research interests quite strongly. Although the quantum annealing processor has not been as useful as a conventional computer so far, it has attracted quite a lot of hopes that the first practical quantum computer using quantum annealing should appear soon.
One of the biggest problems associated with the quantum annealing processor is that no one has known the theoretical power of it so far. Quantum annealing was proposed originally by assuming an isolated interacting spin system. 1,2) The power of quantum annealing in an isolated system has been studied for the last two decades and clarified to some extent. 3) The simplest model of quantum annealing with the transverse field is apt to fail in solving typical problems such as 3-Satisfiability because of a first-order quantum phase transition and/or a localization phenomena of the wave function. 4) In realistic situations, however, the system which quantum annealing is applied to cannot be isolated from but must be inevitably open to an environment. This is true as for the current quantum annealing processor made by D-Wave System Inc. It adopts superconducting flux qubits as physical spins, which are coupled with an environment through a fluctuation in the potential energy of a flux caused by the normal current flowing across the Josephson junction. In fact, the current quantum annealing processor shows properties that are different from the ideal quantum annealing. 7) Therefore it is inevitably important to consider an open system coupled to an environment in order to clarify the power of quantum annealing processors.
The study of time evolution in an open quantum system is significant in the context different from quantum computation as well. The time evolution of a closed (i.e., isolated) quantum many-body system with driving has been studied intensively and a lot of theoretical and experimental progresses have been attained in the last few decades. On the other hand, as we lack analytical or numerical tools that are useful to solve manybody problems with an environment, an open quantum many- * sei01@saitama-med.ac.jp body system has not been understood as much as a closed system. An open system, however, is with no doubt as signifant as a closed system. In fact, time evolution of an open system has started receiving attentions independently of quantum computation. Thus the interests in the time evolution of an open quantum many-body system are growing in the wide range of communities.
In the present paper, we focus on the transverse Ising chain coupled to bosonic baths. This model is standard and has been studied in some earlier works. The single spin version of it was studied extensively by Leggett et al. 8) thirty years ago. An accurate numerical method, named the quasiadiabatic propagator path integral (QUAPI), was developed to study the time evolution of a single spin coupled to a bosonic bath by Makarov and Makri in 1994. 9) After that, QUAPI was applied to study the Landau-Zener model coupled to an environment by Nalbach and Thorwart 12) in 2009. The time evolution of a driven many-spin system coupled to a bosonic bath was first studied by Patanè et al., 14,15) where modified Kibble-Zurek scalings were discussed. The Kibble-Zurek scaling in the presence of baths has been also studied by Nalbach et al. 16) and Dutta et al. 17) The pure transverse Ising chain coupled to a bosonic chain has been studied in the context of quantum annealing by Smelyanskiy et al. 18) and Arceci et al. 19) recently. Quantum annealing of random transverse Ising models in an environment has been studied in Refs. 7, 20-22, for the purpose to clarify the performance of D-Wave's quantum annealing processors.
Most of these theoretical studies on many-spin systems coupled to bosonic baths have assumed the so-called Born-Markov (BM) approximation and/or an integrable spin system. The BM approximation is expected to be valid for weak coupling limit between the system and the environment. It however involves a drawback that the approximation is uncontrollable, namely, one cannot improve the accuacy of the approximation systematically. We resort to neither the BM approximation nor the integrability of a spin system. We instead develop a novel method that is based on the QUAPI and the matrix product state (MPS) formalism. Our method employs approximations in a controllable manner and one can attain the exact result in an appropriate limit. The accessble size of the system reaches N ∼ 10 2 . Our method is applicable to one-dimensional quantum systems with bosonic baths. The present paper is devoted to describe this method and to present some results obtained using it. This paper is organizaed as follows. We introduce the transverse Ising model coupled to bosonic baths in Sec. 2. Some of the numerical results are presented in Sec. 3, where we mention how an error after quantum annealing is influenced by a bosonic environment. The numerical method is descirbed in detail in Sec. 4. Performing the partial trace in the density operator with respect to the bosonic degree of freedom of the baths, we obtain the QUAPI formula for the reduced density matrix. We provide an expression for it in Sec. 4.1. In order to compute the time evolution of the reduced density matrix in systems with more than 10 spins, we needs a trick to avoid an exponential increase in the number of states. We introduce the MPS formalism to reduce the number of bases in Sec. 4.2. We conclude the present paper in Sec. 5

Model
We consider the transverse Ising chain coupled to bosonic baths. The Hamiltonian is composed of H S , H B and H int , representing the spin system, the baths, and the interaction between them, respectively: (1) The system is the transverse Ising chain that is written as where σ α j (α = x, z) are the Pauli's matrices at site j and N is the number of sites. The spin-spin coupling constant and the transverse field are denoted by J j (t) and h(t), respectively. We focus on spatially uniform h(t) for simplicity, though it is possible to consider nonuniform h(t). It is noted that J j (t) and h(t) may depend on time. We assume the open boundary condition. In standard quantum annealing, one considers H S (t) such that h(0) J j (0) and h(τ) J j (τ) for a given runtime τ, namely, the transverse field initially dominates the Hamiltonian while the Hamiltonian coincides with the Ising Hamiltonian finally. The solution of the optimization to be solved is encoded into the ground state of this Ising Hamiltonian. If the system is isolated and the change speed of the Hamiltonian is sufficiently slow, an adiabatic time evolution from the simple ground state of H S (0) to H S (τ) brings us the solution.
The baths are collections of harmonic oscillators whose Hamiltonian is written as where a stands for the mode of harmonic oscillators. b ja and b † ja are the bosonic creation and annihilation operators for the site j and the mode a. ω a is the enregy of the harmonic oscillator in the unit of = 1. We assume that each site has own bath independently.
The interaction between the spin system and the environment is given by where λ ja determines the strength of the interaction. This in-teraction Hamiltnoian implies that each spin couples with own bath independently. We note that one may in general consider coupling between σ α j with α = x as well as z and the bath. In fact, several previous studies [14][15][16][17]23) have considered coupling through σ x j , since this assumption makes the model a free fermionic model through the Jordan-Wigner transformation. However, it has been know that the coupling through σ z j is dominant over σ x j 24, 25) in the system of superconducting flux qubits of which the D-Wave's quantum annealing processor is composed. Hence, in the present paper, we focus our attention on the coupling through σ z j and consider Eq. (4). For the spectrum of the bosonic baths, we define the spectral density as where δ(ω − ω a ) stands for the Dirac's delta function. We then assume a continuous spectral density as and s = 1, namely, the Ohmic spectral density thoughout this paper for the sake of simplicity. Here g 2 and ω c are the coupling constant, between the spin system and the bath, and the cutoff frequency of the spectrum of the bath, respectively. θ(ω) is the Heaviside's step function. Now, we introduce the time-evolution operator U(t) of the composite system that obeys the Schrödinger equation and U S (t) of the system and U B (t) of the bath as Using these operators, we define One can easily see that this obeys where H I int (t) is the interaction picture of H int defined by The density operator is defined by We assume the initial condition as follows.
where |Ψ in is the (given) ground state of H S (0), β is the inverse temperature of the bath, and Z B is the partition function with respect to the bath, which is explicitly written as Z B = Tr B e −βH B = j,a (1 − e −βω a ) −1 . Note that Tr B stands for the trace with respect to the bosonic degree of freedom.
The reduced density operator is defined by Using Eqs. (10) and (13), this writes as To summarize the setup of the time evolution, our composite system is initially in the product state of the ground state of H S (0) and the thermal equilibrium state of H B at inverse temperature β. The spin system and the baths starts to interact at t = 0, and the composite system evolves according to H(t) = H S (t) + H B + H int . The physical quantity of the spin system at t > 0 is determined through the reduced density operator ρ S (t).

Results
In this section, we present results on the kink density for pure and random Ising chains computed by our method. Note that the results on the closed system (namely, g = 0) were obtained by solving the time-dependent Bogoliubov-de Gennes equation for the equivalent free fermion model. 26) The kink density is defined by where Tr S represents the trace with respect to the spin degree of freedom, τ is the final time, N is the number of spins in the system. This quantity vanishes when the spins are perfectly aligned along the z axis in the spin space. Hence the kink density is a measure that quantifies the deviation from the perfect ferromagnetic state. We recall here that we assume that in our model the bosonic operators couple to the longitudinal spin σ z j . This assumption breaks the integrability of the transverse Ising chain and has not been considered in the previous studies for onedimensional system. Therefore the results we present below are entirely new.

Pure Ising chain
We first consider quantum annealing of the pure Ising chain described by the Hamiltonian where τ denotes the runtime of quantum annealing and the time t evolves from 0 to τ. Since this Hamiltonian coincides with the pure ferromanetic Ising chain at t = τ, the kink density measures an error of quantum annealing. We show numerical results on this model with N = 64 in the present subsection. Let us first consider the dependence of the kink density on the coupling strength between the system and the bath. We present results at T = 0 in Fig. 1. The kink denisty at T = 0 increases with inceasing g, implying that coupling to the baths never reduces an error of quantum annealing compared to the closed system even if the temperature of the baths is zero. On the other hand, the kink density decays monotonically with increasing τ for a fixed g. The slope of this decay becomes smaller for larger g. It is not clear from our results whether The kink density for g > 0 is never below the value of isolated system (namely, g = 0) and grows with increasing g at fixed τ. The monotonically decaying behavior with increasing τ is found for all g's we studied. The other parameters are the same as Fig. 1. The kink density grows with g, but has a minimum as a function of τ. The value of τ where the kink density is minimized decreses with increasing g. The results for g ≥ 0.05 suggest that the kink densities for different g's converge with increasing τ.
the kink density vanishes for τ → ∞ when T = 0 and g 0. Figure 2 shows the dependence of the kink density on g at T = 2. One finds that, at this temperature, the kink density turns to an increase with increasing τ for large τ. We guess that it should be true even for g = 0.001 if one has more data for much larger τ's. The minimum position of the kink density moves to a smaller τ with increasing g. This is because the coherence time in which the system keeps unaffected by the bath is shorter for larger g. Interestingly, the result for g = 0.2 suggests that the kink density decreases again with increasing τ for very long τ. A similar phenomenon has been pointed by Amin in Ref. 21 and explained by the equilibration with the bath. We have not yet confirmed whether the equilibration is true or not. This point is an open issue to be studied. Figure 3 shows results on the temperature dependence of the kink density for the system-bath coupling fixed at g =  Fig. 1. The kink density decays with increasing τ for small τ, while it grows with τ for long τ at T ≥ 2. At T = 1, it is not clear that the kink density grows for long τ. The lower the temperature is, the smaller the kink density is. The kink density is never below the value for the closed system (g = 0). 0.01. One finds that the kink density starts to grow for long τ at T ≥ 2, involving a minimum as a function of τ. This minimum shifts to larger τ for lower temperature. It is not clear from the figure whether there is a minimum at T = 1. Anyway, as shown in Fig. 1, the minimum disappears when T = 0. Quite recently, Arceci et al. showed that the minimum at a certain finite τ is a global minimum when T is sufficiently high, but it is in turn a local one with decreasing the temperature, and disappears finally for sufficiently low T . 19) Although we have not confirmed such transitions, there is no conflict between our results and previous ones.

Disordered Ising chain
We here consider quatnum annealing of the disordered Ising chain described by the Hamiltonian where the coupling constants J j 's are assumed to be chosen randomly from the uniform distribution between 0 and 2. Same as the pure Ising chain, the disordered Ising chain considered here has the ferromagnetic ground state and hence the kink density serves to measure an error of quantum annealing. Although the present disordered Ising chain has a trivial ground state, excited states or dynamics are nontrivial and hence this model has drawn a lot of attentions from the quantum annealing community as well as statistical physics community. [27][28][29][30][31][32] The excited states of the present system are characterized by the Anderson localization when the model is switched to a free fermoin model through the Jordan-Wigner transformation. It is a quite interesting problem to inquire into how the localized state is influenced by the baths. As for quantum annealing, the localization nature of excited states might change quantum annealing dynamics from the pure case. Hence, it is quit important to consider a random model in the context of quantum annealing. Here we present a few results on a single instance of this model with N = 64 spins. Figure 4 shows the results on the kink density for temper- Kink density atures T = 0, 1, and 2. The system-bath coupling is fixed at g = 0.01. The behavior of the kink density is qualitatively the same as the pure Ising chain. The kink density is never smaller than that for the closed system even at T = 0. With increasing the temperature, it grows and has a minimum as a function of τ for T = 1 and 2. The influence of environment looks qualitatively the same in pure and random Ising chains on this result. To summarize the results on the pure and random Ising chains, coupling to baths never reduces an error after quantum annealing from the isolated system, even if the temperature of the baths is zero. The negative influence on quantum annealing is more pronounced with increasing the coupling strength g between the system and the baths as well as with increasing the temperature of the baths. The kink density decays with increasing the annealing runtime in the limit of low temperature, while it has a minimum when T is sufficiently high.

Method
In the present section, we describe the numerical method to compute the reduced density matrix of Eq. (17) in detail. The equation of motion of ρ S (t) is too complicated to solve even for a single spin in general. When the coupling constant g 2 is sufficietly small compared to the other energy scales in the model, the Born-Markov (BM) approximation is often applied to reduce the difficulty of the problem. Briefly mentioning, the BM approximation neglects the order of g 4 and higher order terms in the equation of motion, and replaces the reduced density operator of the past (ρ S (s) with s < t) with the current one (ρ S (t)). 33) The resulting equation, called as the Redfield equation, is less complicated than the original one. However, there are a few problems on this equation or approximation. (i) The reduced density operator obtained with the BM approximation has no guarantee that it preserves the unitarity Tr S ρ S (t) = 1, where Tr S means the trace with respect to the system's degree of freedom. (ii) The approximation cannot be improved in the framework of the BM approximation. (iii) The Redfield equation reduces to a linear differential equation for 2 2N unknown complex functions of t. Therefore the computational cost increases exponentially with the number N of spins and this method is restricted to small N up to N ∼ 10. The first and the third problems might be cured by the socalled Lindblad approximation. This prescription, however, makes the second problem more serious.
In the present study, we develop a novel numerical method in order to avoid the problem of the BM approximation. Our approach is based on the path integral and utilizes the MPS formalism, in which the unitarity of the trace is preserved. It includes a sort of approximation, but the accuracy of the result can be improved systematically. Moreover, systems with N ∼ 10 2 spins are accessible. One limitation of our method is that it is suitable only for one dimensional system, as is the case with the DMRG of a closed system. We describe this method in the following subsections.

QUAPI
We are going to compute the matrix elements of ρ S (t), starting from Eq. (17). At first we focus on the trace with respect to the bosonic degree of freedom. This trace can be done by performing the path integral of the bosonic degree of freedom. In order to perform the path integral, we apply the Trotter decomposition to U int (t). Letting M be the Trotter number such that t = M∆t, the Trotter decomposition yields Each exponential operator can be written as where l = 1, · · · , M − 1. Note an identity e A+B = e A e B e −[A,B]/2 for a couple of operators A and B satifying [A, B] ∈ R, and the relation e iH B l∆t b ja e −iH B l∆t = b ja e −iω a l∆t . Inserting the completeness relation made from the product of an eigenstate of σ z j and the coherent state of bosons between the right-sided square bracket and U S (t) in Eq. (21), one obtains the path integral representation for Eq. (17). Fortunately, the path integral over the bosonic coherent-state parameter is Gaussian and can be performed. See Appendix A for details. The result, after taking the continuous limit of the spectral density, is given as follows: where |σ M and |τ M are eigenstates of σ z j with eigenvalues σ j,M and τ j,M ( j = 1, 2, · · · , N), respectively. The effective Hamiltonian H is divided into two parts as H = H S + H int , where H S comes from the isolated spin system and H int is from the interaction between the spin system and the baths. They are defined by respectively. We remark that σ j,0 , σ j,M , τ j,0 , and τ j,M with j = 1, 2, · · · , N in H int must be multiplied by the factor 1 2 . H z l ( σ) in Eq. (24) denotes the Ising-model part in H S (t) that is written as and γ l is defined by and The effecitve wave function Ψ 0 comes from the initial spin state, which is written as The normalization factor N is given by which is necessary to have Tr S ρ S (M∆t) = 1.
One can see from Eqs. (24) and (25) that the effective Hamiltonian H is equivalent to a double-layered (1+1)dimensional Ising model. If the system-bath coupling is absent, each layer is decoupled from the other and the interactions between different times work only between nearest neighbors within the same spatial site, see Eq. (24). In the presence of the system-bath coupling, however, H int brings interlayer as well as intralayer long-range interactions along the time axis within the same spatial site. Figure 5 shows interactions in H schematically.
We note that the long-range interactions induced by the system-bath coupling do work within a spatial site. They never work between spatially sparated spins. Moreover, the long-range interactions as well as the interlayer interactions are uniform in space. This is obvious from Eq. (25) because their interaction constants do not depend on the site index for space. These properties are used when we construct a numerical computation method.
The behavior of the long-range interactions given by K(t) along the time direction is shown in Fig. 6. The real part decays with t from a large positive value. It goes down below zero for low temperatures, while it decreases monotonically towards zero for high temperatures. For long times, it decays as 1/t 2 for T = 0 and 1/t for T > 0. The imaginary part, on the other hand, is independent of the tempratures. It decreases from zero at first with t and then turns into an increase up to zero. One can show that the decay for long times is as 1/t 2 . Now, in order to obtain the matrix elements of the reduced density matrix, it is necessary to perform the trace of exp(H)Ψ 0 with respect to σ j,l and τ j,l with j = 1, 2, · · · , N and l = 0, 1, · · · , M − 1. The brute-force method is useless for this computation, since e N M terms must be summed up. This computation is reminiscent of the transfer matrix

MPS representation 4.2.1 The time direction
Let us focus on a ladder of Ising spins along the time direction with a fixed spatial site, which consists of σ j,l and τ j,l with a fixed j and l = 0, 1, · · · M. The interaction originated from the system-bath coupling works between any pair of Ising spins in this ladder and is restricted inside the ladder, see Fig. 7(a). In addition, there are nearest neighbor interactions due to the transverse field. Now we group σ j,l and τ j,l with the same l (and j as well) and define a pseudospin by this composite. Remark that this pseudospin has four states. Then the ladder of Ising spins is seen as a chain of pseudospins as shown in Fig. 7(b).
Let S j,l denote the state of a pseudospin, such as S j,l = (1 − σ j,l )/2 + 2{(1 − τ j,l )/2} for instance. Using this variable, the factor of exp(H int ) for the spatial site j is written as This can be seen as a wave function for a one-dimensional array of M pseudospins with the basis S j,0 , S j,1 , · · · , S j,M , and written in the form of a MPS as where u (l) S j,l ζ l+1 ;ζ l (l = 1, 2, · · · , M − 1) and ψ S j,0 ζ 1 are defined as follows.
Omitting the site index j the singular value decomposition for ψ S 0 S 1 ···S M defined the right-sided unitary matrix u (M−1) S M−1 ζ M ;ζ M−1 : where with l = M − 1, M − 2, · · · , 2. Finally, we define ψ S 0 ζ 1 as We assume that the singular values λ (l) ζ (l = 1, 2, · · · ) are ordered as λ (l) 1 ≥ λ (l) 2 ≥ · · · . As is usual in one-dimensional systems, the singular value decays rapidly with l. In a practical situation, we introduce a maximum number D t for the matrix dimension of u and discard u (l) S l ζ l+1 ;l l with l l > D t , so that the error due to tiny singular values becomes negligible without having the matrix size exponetially large. As a result, the matrix product state representation by Eq. (33) reduces the size of ψ S 0 S 1 ···S M from 4 M to 4D 2 t (M − 1) + 4D t ≈ 4D 2 t M. This exponential reduction is the greatest merit of the MPS representation.
In order to implement the singular value decompositions in Eqs. (34) and (35), one needs to take into account all elements of ψ S 0 S 1 ···S M which causes a cost exponential in M. To avoid this cost, we introduce the cutoff l c in the range of the interaction K(t) such that K(l∆t) = 0 for l > l c . Since K(t) decays with t, this approximation does not give rise to serious error if one choose as large l c as l c ∆t 1. Thanks to the interaction cutoff, one can neglect the variables S m with m < l − l c in the computation of u (l) .

Transfer matrix method
Having obtained the MPS representation for a chain of pseudospins along the time direction, one can perform the trace for the spin degee of freedom. This computation can be accomplished using the numerical transfer matrix method combined with the MPS representation.
where v ( j,2) and ϕ ( j,2) with j = 2, · · · , N −1 are defined through the singular value decompositions as follows: S 1,1 ,S 2,1 ζ 1,1 ,ζ 2,1 ,ζ 3,1 q 2,1 ,q 3,1φ (1,1) ζ 1,1 q 2,1w (2,1) ζ 2,1 q 3,1 ;q 2,1w (3,1) S 3,1 ζ 3,2 q 4,1 ;p 2,2 =: S j,1 l j+1,1 q j+1,1 ϕ p j−1,1 S j,1 ζ j,2 q j+1,1w ( j+1,1) p j−1,2 ζ j,2 ;p j,2 κ ( j,2) p j,2 w ( j+1,1) S j+1,1 ζ j+1,2 q j+2,1 ;p j,2 =: p j,2 v ( j,2) p j−1,2 ζ j,2 ;p j,2 ϕ ( j,2) p j,2 S j+1,1 ζ j+1,2 q j+2,1 ( j = 3, · · · , N − 2), (56) Equation (54) is arranged into another MPS form in a similar manner to obtain Eq. (50) from Eq. (49) as ×w (2,2) ζ 2,2 q 3,2 ;q 2,2 · · ·w (N−2,2) Repeating the same procedure from Eq. (53) to (59), one can accomplish the trace with respect to S jl with j = 1, 2, · · · , N and l = 2, · · · , M − 1 to have where v ( j) and ϕ (N) are defined through the singular value decompositions as follows:  The energy expectation value at t = M∆t is given by Finally, the ground-state probability that the ground state of the system at time t = M∆t is found in the state after time evolution is given by where |Ψ G denotes the ground state of H S (t = M∆t). Note that, if there are degenerated ground states at t = M∆t, it is necessary to add the above quantities computed for each ground state. For instance, when the ground states are the fully polarized states along the σ z axis, the ground-state probability is written simply as We comment on the matrix dimension of v ( jl) ,ṽ ( jl) , w ( jl) andw ( jl) , or the range of indices p jl and q jl in other words. Speaking regorously, the maximum of matrix size grows exponentially with the number of spin N. However, in a practical situation, the matrix dimension is restricted to a number D s , omitting the bases corresponding vanishingly small singular values. The faster the singular value decays with increasing its index, the smaller D s can be. Note that the singular values are ordered in a descending way. Therefore the decaying behavior of the singular values is crucial to the present method. In fact, it has been well known that, due to the area law of the entanglement entropy, D s can be made so small in the computation of the ground state in one dimensional system. As for our problem of a time-dependent open system in one dimension, although there is no theoretical ground by now, we expect that an acceptable D s to the numerical computation suffices to obtain an accurate result.

Test on a single spin
Let us first look through a single spin for the sake of a test of our method. We consider the Landau-Zener model coupled to a bosonic bath. The Hamiltonian for the spin system is given by where v denotes the velocity of driving. The Hamiltonians for the bath and the system-bath coupling are given by Eq. (3) and (4) with omission of j.
The well-known solution of the Landau-Zener model without the bath is described briefly as follows. Let | ↑ and | ↓ be the eigenstates of σ z with the eigenvalues +1 and −1, respectively. Assuming the initial condition |ψ(t = −∞) = | ↓ which is the ground state of H LZ (t = −∞), the ground state probability denoted by P G that the system remains in the ground state at t = +∞ is given by P LZ G = 1 − exp(−π/2v). In the presence of the bath, however, the analytic solution is not known in general, except for the zero temperature and the high temperature limit. In the high temperature limit, P G is modified as P T →∞ G = 1 2 (1 − exp(−π/v)). 34) Figure 8 shows numerical results using our QUAPI-MPS method. For performing the simulation, we fixed parameters as g = 0.0282, ω c = 10, and ∆t = 0.1. The initial time and final time are set as t in = −40v and t fin = +40v, respectively. The initial state is set at the ground state of H LZ (t in ). The cutoff in the range of interaction along the time direction and the maximum number of states kept in the MPS representation are fixed at l c = 40 and D t = 12, respectively. The ground state probability decreases with increasing the temperature from P LZ G down to P T →∞ G . A similar calculation has been done by Nalbach and Thorwart, 12) who implemented QUAPI without MPS. Our results are quantitatively consistent with their results.

Test on eight spins
We next consider quantum annealing of the pure Ising chain with eight spins. The Hamiltonian is given by Eq. (19). Figure 9 shows the ground-state probability P G that the ground states are found in the final state at t = τ. We fixed ω c = 5, ∆t = 0.025, g = 0.01, and l c = 100. The numbers of states kept fot the MPS representation are set as D t = 20 and D s = 80. When T = 0, P G is almost identical with the values for the closed system (g = 0). When T > 0, P G is lower than T = 0, and increases with increasing τ for small τ but turns into a decrease for large τ. This behavior of P G with respect to τ for finite temperatures is universal as shown above for larger systems. This is understood as follows. When τ is small, quantum annealing ends before the system is influenced by the Results for the closed system were obtained by solving the time-dependent Bogoliubov-de Gennes equation for the equivalent free fermion model. The result at T = 0 is identical to that for the closed system (i.e., g = 0). With increasing the temperature, the ground-state probability decreases. When T > 0, the ground-state probability increases with increasing τ for small τ, and decreases with τ for large τ.
bath. However, when τ is large, the influence by the bath is strong so that P G is lowered. As mentioned in Sec. 1, the exact computation of the timedependent reduced density matrix is too difficult to achieve even in N = 8 spin systems. Hence it is difficult to compare results by our method to exact ones. Here we focus on the perturbative expansion with respect to the system-bath coupling g 2 . As shown in Appendix B, the ground-state probability at final time t = τ is written up to the order of g 2 as where σ zI j (t) := U † S (t)σ z j U S (t) is the interaction picture of σ z j , and the summation with respect to Ψ G implies the summation over the two fully polarized ground states of H S (τ). When N = 8, the quantity P 2 can be accurately computed by solving the Schrödinger equation directly. For instance, we estimated P 2 = −670.30±0.02 for T = 2 and ω c = 5, by discretizing the integrals over t 1 and t 2 and extrapolating the data to the continuous limit. Figure 10 shows the comparison between QUAPI-MPS and the perturbative expansion. For QUAPI-MPS, we chose ∆t = 0.1, l c = 100, D t = 12 and D s = 64. One can see that the extrapolated value of (P G (g) − P G (0))/g 2 by QUAPI- MPS to g = 0 agrees excellently with P 2 obtained by the perturbative calculation.

Concluding remarks
We showed numerical results on quantum annealing of pure and random Ising chains coupled to bosonic baths. When the system-bath coupling is weak (g = 0.01), the baths with zero temperature hardly influences quantum annealing. However, even if the temperature of the bath is zero, the bath hinders quantum annealing with strengthening the system-bath coupling at least in the pure system. This result is nontrivial bacause the dissipation due to a sufficiently low temperature may suppress the nonadiabatic excitation during the time evolution. Our result implies that this does not happen. Although the kink density is larger than the closed situation, it decreases monotonically with increasing the annealing time τ even though the system-bath coupling is not weak whenever the temperature is zero. This does not mean, however, that one can get the solution with probability one by infinitely slow quantum annealing. This is because, when the systembath coupling is not weak, the system does not necessarily equilibriate at its ground state even if the bath is at the ground state initially. Such a picture is different when the temperature is finite. When T > 0, the kink density first decays and then turns to an increase with increasing the annealing time τ. This is a universal feature for pure and random systems.
In the present paper, we gave an explanation of our QUAPI-MPS method for the numerical computation of the reduced density matrix. Advantages of the present method are as follows. (i) It does not rely on the Born-Markov approximation. Therefore one can apply it to a strong system-bath coupling.
(ii) The approximations used in the method, the Trotter decomposition, the limited matrix dimension, and the cutoff in the interaction range along the time direction, is controllable in the sense that one can improve these approximations by changing parameters. (iii) The accessible system size is up to N ∼ 10 2 . The complexity of the present method scales as a polynomial in N and M, where M stands for the Trotter num-ber. (iv) One can engage the infinite size scheme in present method. Then one can apply the QUAPI-MPS method to the infinite size, if the system is homogeneous. (v) The extention of the present method to other type of spin-spin interactions as well as spin-boson interactions is possible. The disadvantage of the present method, on the other hand, is that it should not work if a kind of the entanglement entropy computed by the (square of) singular values obeys so-called the volume law. We are not sure so far when this happens, though we have not encountered this in studying quantum annealing of Ising chains. It is an important open issue to determine the limitation of the present method.
The study of the time evolution of an open quantum manybody system is important with no doubt in the development of near-term quantum computers. In order to confirm the correctness of the operation, the comparison of experimental results with numerical calculation is necessary. Moreover, numerical study may provide suggestions on designing a novel quantum computation, such as a dissipation assisted quantum computation. 18) The present method will be useful for many such studies. Apart from the issue of quantum computer, a lot of interesting problems remain to be solved regarding the time evolution of an open quantum many-body system. In the present paper, we have not discussed scaling properties of the kink density, that are associated with the Kibble-Zurek mechanism. The modification of the Kibble-Zurek scaling in an open system is an urgent issue. Another problem is the relaxation or thermalization of a quantum many-body system coupled to an environment. It is interesting to ask what the steady state or the equilibrium state of an open quantum system is in the presence of driving, such as a quantum quench. We believe that our method should mark the beginning of a new era in the study of open quantum many-body systems.
One of the authors (S.S.) acknowledges fruitful discussions with L. Arceci, S. Barbarino, and G. E. Santoro. The present work was partly supported by JSPS KAKENHI, Japan through Grant No. 26400402.

Appendix A: Path integral
In this appendix, we derive the QUAPI formula, Eq. (23), for the reduced density matrix.
The reduced density operator ρ S (t) is written using the time-evolution operators as Applying the Trotter decomposition, U int (t = M∆t) can be written up to the order of (∆t) 2 as and each exponential operator is written as (A·5) Using the notation for the basis |σ of the spin state such that the completeness relation is given by Dz ja e −z * ja z ja |z}{z|, where the summation is taken over σ 1 , · · · , σ N and Dz ja := (2π) −1 dz ja dz * ja . Inserting the completeness relation between the right square bracket and U S (l∆t) in Eq. (A·3), one obtains the path-integral representation of ρ S (M∆t) as follows: where we defined i∆tλ ja e iω ja (M−1)∆t σ jM−1 . . .
i∆tλ ja e iω ja ∆t τ j1 i ∆t Since the integrals over z's and w's are Gaussian, they are performed analytically. The result is given by where L and K(t) are defined by Eqs. (28) and (29), respectively. We remark that σ j,0 , σ j,M , τ j,0 , and τ j,M with j = 1, · · · , N must be multiplied by the factor 1 2 in the above equation.
We move on to the spin degree of freedom. One can notice that the product U S (l∆t)U † S ((l−1)∆t) is the time evolution operator from t = (l − 1)∆t to t = l∆t. Applying the symmmetric decomposition for the exponential operator, one obtains U S (l∆t)U † S ((l − 1)∆t) = e −iH S (l∆t) ∆t 2 e −iH S ((l−1)∆t) ∆t 2 + O(∆t 3 ). (A·15) Using the notation H z l = − N−1 j=1 J j (l∆t)σ z j σ z j+1 and H x l = − N j=1 h(l∆t)σ x j , an exponential operator can be further decomposed as where H z l (σ m ) = − N−1 j=1 J j (l∆t)σ j,m σ j+1,m and γ l is defined by Eq. (27). Substituting Eq. (A·19) for the matrix elements in Eq. (A·14), one obtains Eq. (23).

Appendix B: Perturbation expansion of the reduced density matrix
In this appendix, we describe the perturbation expantion of the reduced density matrix up to the second order with respect to the system-bath coupling.
As mentioned in Sec. 2, we decompose the time evolution operator U(t) as U(t) = U S (t)U B (t)U int (t), where U S (t) and U B (t) are the time-evolution operatores for the isolated system and the bath, respectively. Then U int (t) obeys an equation where H I int (t) = U † B (t)U † S (t)H int U S (t)U B (t) is the interaction picture of H int . We hereafter assume U S (0) = U B (0) = U int (0) = 1. Now we consider the perturbation expansion of U int up to the second order in H int as U int (t) ≈ 1 + U (1) int (t) + U The density matrix is defined by ρ(t) = U(t)ρ in U † (t). We assume the initial condition: ρ(0) = ρ in and ρ in = |Ψ in Ψ in | ⊗ e −βH B /Z B , where |Ψ in denotes the initial state of the system, β is the inverse temperature of the bath, H B is the Hamiltonian of the bath, and Z B = Tr B e −βH B . Using Eqs. (B·2)-(B·4), the density matrix is expanded into the perturbation series up to the second order as ρ(t) ≈ U S (t)U B (t)ρ in U † B (t)U † S (t) +U S (t)U B (t) U (1) int (t)ρ in + ρ int U (1) † int (t) U † S (t)U † B (t) The perturbation expansion of the reduced density matrix ρ S (t) := Tr B ρ(t) is obtained by taking the trace with respect to the degree of freedom of the bath. Now we assume the Hamiltonians H B and H int in Eqs. (3) and (4), and the Ohmic spectral density defined by Eqs. (5) and (6). Noting Tr B (U (1) int ρ in ) = Tr B (ρ in U (1) int ) = 0 and the kernel function K(t) defined by Eq. (29), we obtain ρ S (t) ≈ U S (t)|Ψ in Ψ in |U † S (t) where σ zI j (t) := U † S (t)σ z j U S (t) denotes the interaction picture of σ z j . We note that the matrix elements U S (t) can be numerically computed for small systems with N up to about N = 10 by solving the Schrödinger equation. Using them, one can evaluate the matrix elements of the right-hand side of Eq. (B·6) through a discrete approximation on integrals.