Kibble-Zurek mechanism in a Dissipative Transverse Ising Chain

We study the Kibble-Zurek mechanism in the transverse Ising chain coupled to a dissipative boson bath, making use of a new numerical method with the infinite time evolving block decimation combined with the discrete-time path integral. We first show the ground-state phase diagram and confirm that a quantum phase transition takes place in the presence of the system-bath coupling. Then we present the time dependence of the energy expectation value of the spin Hamiltonian and the scaling of the kink density with respect to the time period over which the spin Hamiltonian crosses a quantum phase transition. The energy of spins starts to grow from the energy at the ground state of the full system near a quantum phase transition. The kink density decays as a power law with respect to the time period. These results confirm that the Kibble-Zurek mechanism happens. We discuss the exponent for the decay of the kink density in comparison with a theoretical result with the quantum Monte-Carlo simulation. A comparison to an experimental study is also briefly mentioned.


I. INTRODUCTION
When a parameter of an isolated macroscopic system is varied with a finite speed near a quantum phase transition, the system, departing from its ground state, develops into a non-trivial out-of-equilibrium state with topological defects. This spontaneous creation of inhomegeneity during the time evolution near a quantum phase transition is called the Kibble-Zurek mechanism (KZM) [1,2]. The most remarkable phenomenon of KZM is the universal scaling, named as the Kibble-Zurek scaling (KZS), of the density of defects with respect to the changing speed of the parameter. Intensive studies in past decades have established that, as far as a standard quantum phase transition of the second order is concerned, the density of defects obeys a power law with an exponent determined by critial exponents of the associated quantum phase transition [3][4][5].
Compared to an isolated system, dissipative systems have not been studied thoroughly so far. Previous studies on dissipative systems have been limited to ideal freefermion systems in one dimension with no mixing between different momentum sectors [6][7][8][9][10][11]. However, it is no doubt that dissipation is crucially important to out-ofequilibrium states in more real physical systems as well. Therefore it is significant to investigate KZM in dissipative quantum systems from more realistic point of view.
Recently, a commercialized quantum annealing machine manufactured by D-Wave Systems Inc. has received great attention. Quantum annealing was proposed * hiroki.oshiyama.e6@tohoku.ac.jp; Present address: Graduate School of Information Sciences, Tohoku University, Sendai 980-8578, Japan † shibata@cmpt.phys.tohoku.ac.jp ‡ sei01@saitama-med.ac.jp to solve combinatorial optimization problems represented by the energy minimization of an interacting Ising-spin model [12,13]. Following the standard method of quantum annealing, one applies a strong transverse field to interacting Ising spins at first and weaken it slowly with time. Then the spin state which is initialized as the trivial ground state of the transverse field evolves adiabatically, and eventually reaches the target ground state of the interacting Ising model when the transverse field vanishes. The adiabatic time evolution is guaranteed when the square of the instantaneous energy gap above the ground state is sufficiently large compared to the changing speed of the Hamiltonian [14,15]. However, it has been known that a quantum phase transtition usually takes place during quantum annealing and the instantaneous energy gap vanishes at a transition [16]. Therefore, quantum annealing is seen as a time evolution across a quantum phase transition and its performance is closely related to KZM. D-Wave's machine implements quantum annealing physically in its device. Although a lot of tests have been done so far to characterize the performance of the machine, its mechanism has not been fully clarified. Several studies have shown that the system embedded in D-Wave's machine is not closed but open to a dissipative environment. D-Wave's machine is now believed to implement quantum annealing in a dissipative system. Thus KZM in a dissipative transverse Ising model is worth studying from the viewpoint of not only out-ofequilibrium statistical mechanics but also revealing the mechanism of D-Wave's machine.
The difficulty in theoretically studying a dissipative quantum many-body system is that there have been few good method. In order to study the time evolution of a dissipative quantum system, the Born-Markov (BM) approximation is often used [17]. BM approximation reduces the complexity of computation, but it involves an approximation in an uncontrollable manner and one cannot improve accuracy systematically. The Lindblad formalism, which is derived with BM approximation, enables us to perform numerical computation for large systems, i.e., tens of quantum spins for instance. The approximation in the Lindblad formalism replaces the coupling to an environment with an abstract one and hence some of the features of the original system are lost completely. In the present work, we apply the infinite time evolving block decimation (iTEBD) method to the dissipative transverse Ising chain for the first time. The dissipative transverse Ising chain is mapped to a doublelayered two-dimensional Ising model after tracing out the bosonic degrees of freedom. The difference from the dissipationless model is that there are long-range interactions along the time axis. The use of matrix product representation for the spin chain along the time axis enables us to utilize the iTEBD method to the dissipative model. This is essantially an extention of the method for finite system that the present authors developed recently [18]. The present method does not rely on the BM approximation and the approximation involved is controlled by parameters. With the present method, one can study the time evolution as well as the thermal equilibrium of the translationally invariant infinite system without any finite-size effect.
Using our method, we study the equilibrium and outof-equilibrium states of the dissipative transverse Ising chain with the Caldeira-Leggett type of the system-bath coupling. As shown in Sec. II, this model cannot be mapped to a free-fermion model and has not been studied so far in the context of KZM. We first present the groundstate phase diagram, confirming that the quantum phase transition survives in a certain range of the system-bath coupling strength. Then we investigate the scaling of the density of defects associated with the change of the Hamiltonian linear in time across a quantum phase transition. It is shown that the density of defects scales as a power law of the changing speed of the transverse field, and that the exponent diminishes gradually with increasing the strength of the system-bath coupling. The result is discussed in comparison to the previous work. We finally comment on experiments of KZM in the present model using D-Wave's machines. This paper is organized as follows. We specify the dissipative transverse Ising chain in the next section. In Sec. III, we describe our new method for equilibrium and out-of-equlibrium states of the dissipative transverse Ising chain based on the discrete-time path integral and iTEBD. Section IV is devoted to the presentation of our results. The ground-state phase diagram of the present model is given and KZS is discussed in detail there. This paper conludes with Sec. V.

II. MODEL
We consider the dissipative transverse Ising chain described by the Hamiltonian where H S , H B and H int represent the transverse Ising chain, the boson bath, and the interaction between them, respectively, and given by Here σ c j (c = x, z) are the Pauli operators at site j, J and Γ stand for the coupling strength of Ising spins and the transverse field, respectively, and b † ja (b ja ) is the bosonic creation (annihilation) operator at site j with mode a. We make = 1 throughout the paper. We assume that the boson bath has the Ohmic spectral density where α is the coupling strength between the spin system and the bath, and ω c is the cutoff energy of the bath spectrum which is chosen to be larger than the other energy scales. The spin-boson interaction assumed in Eq. (4) is the many spin version of the Caldeira-Leggett model for the superconducting flux qubit [19,20]. The present model is the simplest and the most realistic one to study the dissipative system in D-Wave's machine. The presence of the system-bath coupling breaks the integrability of the transverse Ising chain, H S . The nonintegrability of the present model makes it difficult to study out-of-equilibrium states. We utilize the discretetime path integral and infinite time evolving block decimation (iTEBD) to study the infinitely large system.

III. METHOD
In this section, we describe the numerical method to compute the reduced density matrix and physical quantities.
A. Discrete-time path integral

Thermal equilibrium
We consider first the thermal equilibrium state of the full system. The reduced density operator for the spin subsystem is given by where β is the inverse temperature, Z is the partition function, and Tr B denotes the trace with respective to the bosonic degree of freedom. The exponential operator in RHS is arranged, using the symmetric Trotter formula [21,22], into where M denotes the Trotter number, ∆τ = β/M , and Using Eqs. (7) and (8) in Eq. (6), and inserting the completeness relations between exponential operators, one obtains a path integral representation for the reduced density matrix. After performing the Gaussian integrals with respect to bosonic degree of freedom, the resultant path-integral formula of the reduced density matrix is written as where |σ l denotes the eigenbasis of σ z j satifying σ z j |σ l = σ j,l |σ l , and C is a normalization factor. B eq and F eq are defined by and F eq (σ j,0 , · · · , σ j,M ) and Note that Eq. (9) represents the partition function of a (1 + 1) dimensional Ising model with the long-range interaction in imaginary time.

Out of equilibrium
We next consider the time-dependent dissipative transverse Ising chain with H S (t) = −J(t) j , and assume the initial state being a product state between a pure spin state denoted by |Ψ 0 and the thermal equilibrium of the boson bath: where β is the inverse temperature of the boson bath and Let U (t) be the time-evolution operator that satisfies i d dt U (t) = H(t)U (t). The reduced density operator is written as a In spite of the time dependence in the Hamiltonian, an algebra with the symmetric Trotter decomposition and Gaussian integrals similar to those used in the previous subsection yield the path-integral formula [18,23]: = exp i∆tJ l σ j,l σ j+1,l − σ j,l σ j+1,l , and F(σ j,0 , σ j,0 , · · · , σ j,M , σ j,M ) where D is a normalization factor. The coupling parameters between Ising spins are given by and We remark that two types Ising variables, σ j,l and σ j,l , are included in Eq. (16). This is because the reduced density operator in Eq. (15) involves the time-forward operator U (t) and the time-backward one U † (t). The former results in σ j,l and the latter in σ j,l in Eq. (16). Hence the effective Hamiltonian H represents coupled double layers of two-dimensional Ising models with nearest neighbor interactions in space and long-range interactions in time.
Since the strength of the long-range interaction K(t) decays with t, we set a cutoff t c in the numerical simulation so that K(t) = 0 for t > t c .

B. iTEBD
In this subsection, we briefly describe how to compute the reduced density matrix, Eq. (16), for an outof-equilibrium state. The method mentioned here can be easily applied to the thermal equilibrium state. We hereafter use a shortened notation S j,l for (σ j,l , σ j,l ) in such a way that S j,l = 1 2 (1 − σ j,l ) + (1 − σ j,l ). Using the singular value decomposition (SVD), one can write Eq. (18) in the form of a matrix product as F(S j,0 , · · · , S j,M ) The factor B l can be represented, in a similar manner, as Finally, we write σ 0 |Ψ 0 Ψ 0 |σ 0 as Strictly speaking, the largest matrix dimensions of φ (l)S and Λ (j,l)S increase exponentially with M and N , respectively. However, as is the case with the ground-state wave function in a one-dimensional quantum system, one can choose the basis of matrices so that the necessary information is preserved with keeping the largest matrix dimension constant in M or N [24,25].
Let us now consider the infinite system (N → ∞) with the translational invariance. In order to perform the trace over the spin degrees of freedom in Eq. (16), we exploit the iTEBD algorithm [26,27]. This algorithm performs the trace with respect to {S j,l } for all j from l = 0 to l = M − 1 iteratively. After the M times iteration, we obtain the infinite matrix product state (iMPS) representation of ρ S (t), from which we can calculate physical quantities.
Hereafter, we write the l-th level of iMPS representation generated from the l − 1 level as ψ (l)η j,l µ l ,µ l , where η j,l denotes the matrix index for φ (l)S j,l in Eq. (22), and µ l and µ l denote matrix indices. We denote the largest matrix dimension of φ (l)S j,l and ψ (l)η j,l by χ t and χ s , respectively.
where µ l+1 = (ζ j,l , µ l ) and µ l+1 = (ζ j,l , µ l ). For each step, one can keep the matrix dimension up to χ s by means of the truncation algorithm in iTEBD. The fianl iMPS is calculated as Thus the reduced density matrix is given in terms of iMPS as (28) where tr denotes the matrix trace. The expectation values of physical quantities are obtained in terms of ψ S as follows. Redefining ψ S andψ ≡ ψ S=0 + ψ S=3 in such a way that the largest eigenvalue ofψ is unity, the magnetization is given by where v t L and v R are the left and right eigenvector of ψ for the unit eigenvalue that satisfy the normalization condition v t L v R = 1. The correlation function is given by Equations (29) and (30) are used to compute expectations of physical quantities.

IV. RESULTS
In order to study properties of the dissipative transverse Ising chain in and out of equilibrium, we introduce a parameter s for the system Hamiltonian, so that J = s and Γ = 1 − s, as This Hamiltonian reduces to the decoupled spins in a transverse field when s = 0, while it corresponds to the Ising model when s = 1. In the absence of coupling to the bath, the ground state of H S (s) is an ordered one for 1 2 < s < 1, a disoredered one for 0 ≤ s < 1 2 , and quantum critical at s = 1 2 . The Hamiltonian of the full system is given by H = H S (s) + H B + H int . We will focus on the zero-temperature phase diagram of the full system in the next subsection and then move to the time evolution in Sec. IV B.

A. Phase diagram
Consider that the full system is in the thermal equilibrium with the inverse temperature β. The reduced density operator is given by Eq. (6) and the correlation function with respect to the longitudinal spin is given by It is expected that, since the longitudinal spin is disordered in the equilibrium at a finite temperature, the correlation function decays exponentially for sufficiently large r as where ξ represents the correlation length. Figure 1(a) demonstrates an exponential decay of C z (r) obtained by the iTEBD simulation. Based on the ansatz of Eq. (33), we estimate the correlation length from the correlation function by With lowering the temperature, the correlation length obeys the scaling relation to β depending upon the property of the ground state. When the ground state is disordered, the correlation length becomes constant with respect to the temperature. When the ground state is ordered, on the other hand, the correlation length diverges exponentially as ξ ∼ e aβ as β → ∞ with a constant a. Finally, a quantum critical ground state bears a powerlaw scaling, ξ ∼ β 1/z , where z is the dynamical critical exponent. We show results on ξ for several s in Fig. 1(b).
One can see that ξ grows faster than a power law with β for sufficiently large s's, while it looks converging for smaller s's. At a certain critical value of s in between, ξ behaves as a power law of β, from which one can determine the location of the quantum critical point s c for a fixed system-bath coupling α. Thus we show the ground-state phase diagram of the dissipative Transverse Ising chain in Fig. 2. Note that, at α = 0, the transverse Ising chain is isolated from the bath and the critical point s = 1 2 is exactly known. One finds that the critical point declines with increasing the system-bath coupling and the ordered phase extends down to a smaller s. Although our numerical method becomes unreliable for α > 0.64, previous studies have revealed that when s = 0 an order-disorder transition takes place at α ≈ 0.625 for ω c = 10 and at larger α for smaller ω c [20,[28][29][30]. This implies that the critical line should terminate at α larger than 0.625 and s = 0. Our results for ω c = 5 are perfectly consistent to the previous studies.
Quantum critical properties of the dissipative transverse Ising chain have been studied earlier by Werner et Since the difference between these two lines is invisible, we consider that the solid line represents the ground state. As for the time evolution, the initial state is assumed to be the product state of the ground state of the spin system and that of the boson system. The dashed vertical line shows the location of the quantum phase transition. After an initial relaxation, the state follows the ground state but starts to deviate from the ground state near the quantum phase transition. The parameters for numerical simulations were chosen as ∆t = 0.05, χt = 128, χs = 128, and ωc = 5.
al. using the quantum Monte-Carlo (QMC) simulation for an effective model [29]. They have found that the critical exponents (ν, z) associated with the quantum phase transition change from (1, 1) to (0.63, 2.0) by introducing an infinitesimal system-bath coupling α and do not vary on the critical line. As for our simulation, the exponent z can be evaluated from the power law of ξ when s is fixed at a quantum critical point. The result is in fact 1 z ≈ 0.46 for α = 0.09, which is in reasonable agreement with 1 2 from the QMC simulation. It should be noted that the parameters involved in the QMC simulation are defined for an effective model with the unit Trotter imaginary-time step, i.e., ∆τ = 1. Therefore there is not necessarily a quantitative consistency among the phase diagrams given in ref. [29] and Fig. 2.

B. Time evolution across a quantum phase transition
Having obtained the ground state phase diagram, we next focus on the time evolution across a quantum phase transition. Let us consider the time-dependent Hamiltonian, where time t moves from 0 to t 0 . The parameter t 0 controls the changing speed of H S . The larger t 0 is, the slower H S varies. As seen in the previous subsection, the present model involves a quantum phase transition.
When the system-bath coupling α is not so large, The system undergoes a quantum phase transition at a certain normalized time s = t/t 0 for sufficiently large t 0 . See Fig. 2.
Suppose that the full system is initially in the ground state. On approaching a quantum phase transition, we expect that the full system is excited because the characteristic time of the instantaneous ground state exceeds the time scale t 0 of Hamiltonian's variation. This is illustrated by Fig. 3, where the energy expectation value of the spin Hamiltonian H S (t/t 0 ) is plotted as a function of the normalized time t/t 0 for various t 0 . We have assumed in our simulation that the full system is initialized as the product state of the ground states of the spin system and the boson system. Although this is not the true ground state of the full system, we find that the state relaxes into the ground state of the full system in an initial period. After that, the full system keeps the ground state up to a time near the quantum phase transition, where the state starts to deviate from the ground state.
It has been known that the excitation during time evolution near a quantum phase transition is well explained by KZM. According to the phenomenological theory of KZM, the state after passing a quantum phase transition involves topologial defects with a characteristic lengthξ which is scaled by t 0 asξ ∼ t ν/(zν+1) 0 , where ν and z are the critical exponents [3,4]. This universal scaling relation is KZS. In our model, the topological defects is identified by the kinks between neighboring spins in the final state. The kink density is then defined by Using KZS, this should scale as (37) Figure 4 shows results of our simulation on the kink density n as a function of t 0 . The kink density decays monotonically with t 0 and obeys the power law n ∼ t −b 0 approximately for large t 0 . The exponent b drawn by fitting the data decreases from 0.5 to around 0.25 with strengthening the system-bath coupling α from 0 up to 0.49, beyond which it is uncertain from our data whether a power law is valid or not.
As mentioned above, the exponents ν and z are known to be ν = z = 1 for α = 0 by the exact solution of the transverse Ising chain, whereas they are estimated as ν ≈ 0.63 and z ≈ 2.0 for α > 0 by the QMC simulation [29]. These values predict that b = 0.5 for α = 0 and b ≈ 0.28 for α > 0. Comparing them to our result, the fact that b is reduced by coupling spins to the bath is shared by the theory and our simulation. However, although the theory predicts a discontinuous change of b from 0.5 to 0.28 by an infinitesimal α, our results implies a continuos change of b. Moreover, the smallest value of b by our simulation, which is obtained for α = 0.49, is slightly smaller than 0.28. These discrepancies might be attributed to the range of t 0 . Although our method of simulation is not able to access t 0 > 200, the slope of n might be closer to that of t −0.28 0 if one has much larger t 0 .
Recently, KZM of the transverse Ising chain was studied experimentally using D-Wave's machines. As mentioned in Sec. I, the physical system in D-Wave's machine can be regarded as the dissipative transverse Ising model and hence D-Wave's machine should serve as a quantum simulator of the present model. In fact, it has been reported in ref. [31] that the kink density decays as t −b 0 with b ≈ 0.20 by the device in NASA and ≈ 0.34 by the low-noise one in Burnaby. We note that the range of t 0 investigated in these experiments is roughly 4-digit larger than ours. Anyway, the fact that the exponent decreases from 0.5 is common to the experiment, our simulation, and the theory with the QMC simulation. As for the numerical value, the comparison between the experiment and our simulation is difficult, since the strength of the system-bath coupling is unknown for D-Wave's machines. Nontheless, b = 0.34 obtained experimentally lies in the range from 0.5 for α = 0 to 0.25 for α = 0.49 by our simulation. Therefore our results are consistent with the experimental result in ref. [31].

V. CONCLUSION
We investigated KZM of a dissipative transverse Ising chain. First we presented the ground-state phase dia-gram and confirmed that a quantum phase transition survives even at a finite coupling strength between the spin system and the boson bath. We next showed the energy expectation value of the spin subsystem as a function of time when the spin Hamiltonian varies linearly in time from that of the decoupled spins in the transverse field to that of the Ising chain and when the full system is initialized as the product state of the ground states for spins and bosons. We found that the energy of the spin subsystem relaxes into the value at the ground state of the full system soon and starts to deviate from it at around the quantum critical point. This result implies that KZM takes place in our model. We finally showed the kink density as a function of the time period with which the spin Hamiltonian varies. We observed that the kink density decays as a power law with an exponent b and that b decreases with increasing the system-bath coupling. This fact partly agrees with a theoretical deduction with a QMC result that b changes from 0.5 to 0.28 by an intro-ductin of the system-bath coupling, though the values of b do not always agree. Our results on the kink density, on the other hand, agrees well with an experimental result by D-Wave's machine. Geting better understanding on b remains to study. To this end, quantum simulators like D-Wave's machine would be useful.
We used a new numerical algorithm in the present work, which is a combination of the discrete-time path integral representation for the reduced density matrix and iTEBD. This method was for the first time applied to a dissipative quantum many-body model. So far there had been no method to study out-of-equilibrium states of such a model with large system size and without the Born-Markov approximation. The present study will pave the new way to the study of out-of-equilibrium physics in a dissipative quantum many-body system.