Molecular Imprinted Sites Translate into ... - ACS Publications

Technology, 130 Meilong Road, Shanghai, 200237, China. KEYWORDS: sensors, rheometry, switch stiffness, donor-acceptor, host-guest chemistry. Abstract...
0 downloads 0 Views 2MB Size
Journal of Physics: Condensed Matter

PAPER

Quenching dissipative quantum Ising chain: an exact result for nonequilibrium dynamics To cite this article: Xin Li et al 2019 J. Phys.: Condens. Matter 31 235402

View the article online for updates and enhancements.

This content was downloaded from IP address 154.59.124.171 on 22/08/2019 at 10:40

Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 31 (2019) 235402 (8pp)

https://doi.org/10.1088/1361-648X/ab0e47

Quenching dissipative quantum Ising chain: an exact result for nonequilibrium dynamics Xin Li1,2 , YouYang Xu1 and ShunCai Zhao1 1

  Faculty of science, Kunming University of Science and Technology, Kunming 650093, People’s Republic of China 2   State Key Laboratory of Surface Physics and Department of Physics, Fudan University, Shanghai 200433, People’s Republic of China E-mail: [email protected] Received 13 December 2018, revised 26 February 2019 Accepted for publication 8 March 2019 Published 1 April 2019 Abstract

To create a more comprehensive understanding of nonequilibrium dynamics of open quantum many-body systems, we visit an exactly solvable example, that is a quenched transverse-field Ising chain coupled to Markovian baths, which act locally on the Jordan–Wigner fermionic space. Performing explicit calculations on the heat transfer, transverse magnetization, and kink density, we find that the imbalance of two opposite damping mechanisms play a crucial role in constructing a nontrivial nonequilibrium steady state accompanied with a dissipative quantum phase transition, also that the competition between unitary drive and decoherence does not necessarily lead to a quasi-stationary state or prethermalization under certain ordinary relaxation. Keywords: Ising chain, dissipative quantum phase transition, nonequilibrium dynamics (Some figures may appear in colour only in the online journal)

1. Introduction

shell in which some physical quantities are almost conserved [12]. Due to the inevitable interaction with an external environ­ ment, the dissipation eventually drives the system to a nonequilibrium steady state (NESS) ρ∞ (g). When the steady properties of the system abruptly changes with the variation of the parameter g, we say a dissipative phase transition (DPT) occurs [13–16], relevant to the low excitation spectrum of the Liouvillian superoperator. The interplay between the coupling within many body systems and the dissipation determines a new class of the non-equilibrium critical phenomena [17, 18], which can be demonstrated on a variety of highlycontrollable platforms such as Lattices of superconducting resonators [19, 20], Rydberg atoms in optical lattices [21, 22], and nano-optomechanical oscillators [23, 24]. At the same time, these quantum simulators provide a rare chance for further probing the driven-dissipative many-body dynamics out of equilibrium, accompanied with the occurrence of a series of new phenomena. For instance, a test on an open Ising chain coupling to bosonic baths shown that the quantum diffusion near the quantum critical point could speed up the quantum

Recent years there has been increasing attention paid to the nonequilibrium dynamics of quantum many-body systems both in theoretical and experimental aspects [1–5]. For the long-time dynamics of the isolated systems after a sudden quench a common belief is that certain observables at late time will relax to a stationary state described by either an effective thermal distribution or a generalized Gibbs ensemble (GGE) depending on the integrability of the Hamiltonian. In particular, suppose we first prepare a quantum many-body system in the ground state of an integrable Hamiltonian H0 and then quench the system into a new Hamiltonian H = H0 + V , where V is a small integrability breaking term. It is found that in the weakly non-integrable system some observables might pass through a prethermal stage predicted by a deformed GGE [6–9] before relaxing to a real thermal state. The relevant experimental evidences for prethermalization have been presented via manipulating trapped 1D Bose gases [10, 11]. Such a time-scale separation can also emerge if there exists a subspace of the energy

1361-648X/19/235402+8$33.00

1

© 2019 IOP Publishing Ltd  Printed in the UK

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402

annealing [25]; the prethermalization seemed to appear in a noisy Ising model after a sudden quench of the transverse field [26]; the competition of Markovian noise and disorder also displayed a prethermal Anderson localisation as a result of quenching the Ising chain [27]. However, at present a variety of intriguing dynamic characteristics in the dissipative manybody systems still need to be mined. On the theoretical treatment, it should be noted that the majority of problems about open quantum many-body sytem have to be discussed case-by-case. To date there have been very few general methods to deal with this within a unified framework except for the Keldysh formalism of non-equilibrium Green’s functions [28] and time-dependent density matrix renormalization group techniques [29] etc. The coupling to the external environment is so intractable as to mark the weak coupling, Markov, and secular approximations in general. Recent work [30, 31] sought out a notable class of exactly solvable quantum Liouville equation with a general quadratic Hamiltonian of N fermions and linear fermionic Lindblad operators, immediately extending to more general cases with thermal bath [32]. Although the actual microscopic architectures of the baths are completely unclear in general, through utilizing the phenomenological Lindblad operators locally on Jordan–Wigner fermionic space some physically meaningful results have yet been obtained. For example, applying the fermionic bath Keck et  al [33] pointed out that the competition between the nonadiabatic effects following a Kibble–Zurek mechanism [34, 35] and the dissipation leads to an optimal working time for adiabatic quantum computation; suddenly quenching an open Ising chain coupled to this bath may has detected a dynamic quantum phase transition [36]; in the presence of this bath the stability of Majorana edge modes at the boundaries of Kitaev chains has also been investigated numerically [37]. However, so far the exact analytical study on open quantum many-body system are still very lacking. In this paper, we focus on a solvable case, i.e. an open Ising chain subject to this type of bath local in Jordan–Wigner fermion space, aiming at the post-quench dissipative dynamics itself.

Hλ = P + Hλ+ P + + P − Hλ− P − , (2) N here P ± = [1 ± j=1 (1 − 2c†j cj )]/2 are the projectors onto the even (+) and odd (−) parity subspaces. While the parity subspace Hamiltonians Hλ± own the same form N 1 Hλ± = − [λ(1 − 2c†j cj ) + (c†j − cj )(c†j+1 + cj+1 )], (3) 2 j=1

different boundary conditions cN+1 = −c1 and cN+1 = c1 are set for Hλ+ and Hλ−, respectively. Hereafter we restrict our attention to the even parity subspace only, enabling the substitution Hλ ≡ Hλ+ for brevity. To retain the analytic solvability for the problem of the open quantum many-body system, we consider the fermion on every site j  is coupled to its own structured reservoir locally. A current topic in quantum control also involves the design of this artificial reservoir to obtain a desired, often exotic, steady state. The time evolution of the open Ising chain can be described as a quantum Markov process obeying the Lindblad master equation ˆ = −i[Hλ , ρ], ˆ + D[ρ], ˆ ρ˙ = H[ρ] H[ρ]  1 (4) ˆ = D[ρ] Γ [Lj ρLj† − {Lj† Lj , ρ}], 2 j,=p,d

ˆ governs the unitary evolution, and the dissipator D ˆ here H describes a nonunitary dissipative process that drives the system into a fixed point, involving two types of the damping channels with dissipation rates Γp,d > 0:

(i) Ljp = c†j pumping mechanism, (5) (ii) Ljd = cj dacaying mechanism.

Even though being surprisingly simple, such Lindblad formalism equation (4) is reasonable for the system whose typical time scale is much slower than the reservoir correlation time. It should be stressed that besides the simulation for the exper­ imental realities, such as in circuit-QED, cold-atom settings or trapped ions, choosing these dissipators is also due to their solvability, but the fermionic decaying (or pumping) operator is not local in spin representations. Using the Fourier trans−iπ/4  i2πjk/N formation cj = e √N (k ∈ ± 12 , ± 32 , ..., ± (N−1) k ck e 2 ), the system Hamiltonian can be written as a direct sum over  different subspaces: Hλ = k>0 Hλk , while the density matrix can be factorized into ρ = ⊗k>0 ρk for limiting the initial state to the ground state of the system throughout the work. Moreover, the dissipative part via Fourier transformation still remain local, hence the quantum evolution is only restricted to the quasimomentum subspace for every k  >  0:  1 ρ˙ k = Sˆk [ρk ] ≡ −i[Hλk , ρk ] + Γd [ci ρk c†i − 2 i=k,−k   1 {c†i ci , ρk }] + Γp [c†i ρk ci − {ci c†i , ρk }], (6) 2

2.  Treating master equation of open quantum Ising chain Begin with the Hamiltonian of transverse-field Ising chain which is of the form N

1 x x Hλ = − (σj σj+1 + λσjz ), (1) 2 j=1

where σjx,z are Pauli operators on site j , the transverse field α = σ1α with strength λ and periodic boundary conditions σN+1 are imposed. The quantum paramagnetic (|λ| > 1) and ferromagnetic (|λ| < 1) phases are separated at the critical point |λ| = 1. Actually, using the Jordan–Wigner transformation  σj† = l  0 can be written as 

ρ11 k (t)  0  ρk (t) =  0  xk14 (t) − iy14 k (t)

0 ρ22 k (t) 0 0

0 0 ρ33 k (t) 0

Figure 1.  The contour plots of the absorbed heat per spin Q A (∞) (see equation (14)) as a function of λ and λ0 with the different proportions of decaying rate to pumping rate: (A) Γd /Γp = 0, (B) Γp /Γd = 10, (C) Γd /Γp = 1, (D) Γp /Γd = 0. The green line expressing Q A (∞) = 0 forms the boundary between endothermic and exothermic zones. Set N  =  1000.

mechanisms are equal, i.e. Γm = 0; otherwise the NESS is a nontrivial mixture state implicating a quantum phase trans­ ition. Besides, if λ = ∞ we have   2 Γp 0 0 0 2  Γt    Γd Γp 0 0 0  2 Γt , ρk (t = ∞, λ = ∞) =  (11) Γd Γp 0 0 0   Γ2t   Γ2d 0 0 0 2 Γ

 xk14 (t) + iy14 k (t)  0  , 0  ρ44 k (t)

(9) of which the concrete expressions will be given in the appendix. In the long time limit, we have xk14 (∞) =

sin 2θk 2k Γm , Γ3t + 4Γt 2k

33 ρ22 k (∞) = ρk (∞) =

y14 k (∞) =

t

which is a same diagonal density matrix for any k, being reduced to a dark state |1k , 1−k  (or |0k , 0−k ) in case Γd = 0 (or Γp = 0 ).

sin θk k Γm , Γ2t + 42k

Γd Γp Γ2t + (4Γd Γp + sin2 θk Γ2m )2k , Γ4t + 4Γ2t 2k

3.  Nonequilibrium dynamics of observables

44 ρ11 k (∞) = αk (d, p), ρk (∞) = αk ( p, d),

Γ2d Γ2t + (sin2 θk (3Γd + Γp )Γm + 4Γ2d )2k , Γ4t + 4Γ2t 2k (10)

3.1.  Work done, heat transfer and residual energy

αk ( p, d) =

In this section we will visit some observables obtained from the analytic solution. Firstly, the whole process is divided into two steps: transient work done on system by suddenly quench at time t  =  0; then continually absorbing or releasing heat from or to the reservoir. We suppose ρ0 = |0λ0 0λ0 | is the initial state with vacuum energy E0λ0 , for the nonequilibrium processes governed by equation (4). The fluctuation of the work done on the system requires a description with a probability distribution, which according to a two-measurement protocol [38] can be written as

here Γm ≡ Γp − Γd denotes the intensity difference between the two damping mechanisms, Γt ≡ Γp + Γd is the Liouvillian gap, also called asymptotic decay rate [14–16], and ρk (∞) is the unique eigenmatrix of the superoperator Sˆk corresponding to the zero eigenvalue for each k mode, fulfilling the equality ∂t ρk (∞) = Sˆk [ρk (∞)] = 0. The equation  (10) indicates that the NESS is independent of the initial state, and will be the infinite temperature state if the strengths of the two dissipation

3

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402 θ0

cos ∆θk ) + 20k sin2 2k ] with ∆θk ≡ θk − θk0, no different from the one of the closed Ising chain. By considering the  λ injected total energy EI (t) = Tr[Hλ ρ(t)] − E0 0 = k>0 θ0 [k ((1 − e−tΓt ) cos θk Γm /Γt − e−tΓt cos ∆θk + cos θk ) + 20k sin2 2k ], we concentrate on the absorbed heat from the reservoir: QA (t) = Tr[Hλ ρ(t)] − 0λ0 |Hλ |0λ0   Γm (13) k (cos θk + cos θk ). = (1 − e−tΓt ) Γt k>0

In the thermodynamic limit, the sum over k transforms into an integral over the half of the Brillouin zone k ∈ [0, π]. The absorbed heat per spin Q A ≡ QA /N can be further obtained QA (t) 1 − e−tΓt ( = 2

π 0

dk 1 + λλ0 − (λ + λ0 ) cos k λΓm  + ) π Γt 1 − 2λ cos k + λ2 0

0

λΓm ) iff λ0  0, Γt 4λ0 (λ + λ0 )(λ0 + 1) with A[λ, λ0 ] = EllipticE[ ] (1 + λ0 )2 πλ0 4λ0 (λ − λ0 )(λ0 − 1) + EllipticK[ ] , (14) (1 + λ0 )2 πλ0 =



1−e 2

−tΓt



(A[λ, λ0 ] +

here EllipticE[x] (or EllipticK[x]) gives a complete int­ egral (or the first kind of complete elliptic integrals) for x. Equations  (13) or (14) displays a very trivial dynamics of the absorbed heat, of which the asymptotic value could be discussed in detail. Figure  1 exhibits the final absorbed heat per spin Q A (∞) with different dissipative param­ eters, reminding us that the proportion of the decaying rate to the pumping rate has complicated effects on the division of the endothermic and exothermic zones in (λ, λ0 ) plane. Since Q A (λ, λ0 , Γm , t = ∞) = QA (−λ, −λ0 , −Γm , t = ∞), changing the pure pumping mechanism to the pure decaying one will lead to a mirror transformation of the two zones, as shown in figures  1(A) and (D); while figures  1(B) and (C) tell us that exothermic zone may be separated by the endothermic zone when the two mechanisms coexist. Moreover, as A[λ, λ0 ] is always greater than zero if λ and λ0  0, in order to release heat the pumping rate need to be significantly lesser than the dacaying one among the certain region of parameters. Otherwise, since the energy eigen spectrum changed to k from 0k after the quench, the factors k in the W, EI and QA have been cancelled totally, leading to that all of them are linear with λ. Hence, these three process volumes cannot contain the quantum criticality about λ, but can do the one about λ0 . We plot the final absorbed heat per spin Q A (∞) and its first derivative ∂λ0 Q A (∞) as functions of λ0 in figure  2. The latter highlights the nonequilibrium criticality (see figures  2(B) and (C)). In the thermodynamic limit, we have  π dk (λ−λ0 ) sin2 k ∂λ0 Q A (∞) = 0 2π 3 , indicating that, due to 2

Figure 2.  (A): The absorbed heat per spin Q A (∞) as a function of λ0 with different λ. (B): The first derivative ∂λ0 Q A (∞) is a fine indicator for exploring the quantum criticality of λ0 . (C): The variation of ∂λ0 Q A (∞) with λ0 for different N at λ = 2. We assume Γd = Γp = 0 in (A), N  =  1000 in (A) and (B). Actually, limN→∞ ∂λ0 Q A (∞) = ∞ at λ0 = 1, and irrelevant to Γd and Γp.

 P(W) = δ(W − Emλ + E0λ0 )|mλ |0λ0 |2 , (12) m0

where |mλ  is the eigenvector of the post-quench Hamiltonian Hλ with energy Emλ . The generating func +∞ tion takes the form G(t) = −∞ dWe−iWt P(W) λ0 G(t) = eitE0 0λ0 |e−itHλ |0λ0 . The first cumulant of work distribution provides the average work done  W = i∂t G(t)|t=0 = 0λ0 |Hλ |0λ0  − E0λ0 = k>0 [k (cos θk −

(1−2λ0 cos k+λ0 ) 2

the closed energy gap at k  =  0, there is a divergence at λ0 = 1

4

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402

Figure 4.  The steady value of the transverse magnetization per spin in the thermodynamic limit as a function of λ with different decaying rates at Γp = 0 . Exchanging the values of the two rates Γd and Γp only Inverts the sign of Mz (∞), which is independent of the initial value λ0 (see equation (18)).

Figure 3.  The time-evolution of the transverse magnetization per

spin Mz (t) (red line) under dissipation can be divided into two parts: the exponential decay term e−tΓt MzQ (t) (blue line) and the drifting correction Γm F(t) (green line), MzQ (t) (black line) as an ordinary one gives a contrast. N  =  1000.

except λ = λ0 = 1 (see the blue line in figure  2(B)). It also should be noted that, as Γd or Γp = 0 , this criticality is belong to DPT in essence, different from the equilibrium case where the quantum phase transitions are only determined by the low energy spectrum of the system Hamiltonian rather than disspation. Finally, the steady value of the injected energy per spin EI (∞) ≡ EI (∞)/N = λΓp /Γt − E0λ0 /N implies an interesting feature that as long as the pumping exists, the quenching can always adjust the steady energy of the sytem, otherwise the initial vacuum energy will be leaked out completely.

2 sin2 θ (sin2 (t )Γ2 +sin(2t )Γ  +22 )

t k k k k t k bk (t) = . With certain Γ3t +4Γt 2k parameters that would not cover up the basic features, we plotted the time-evolutions of Mz (t), e−tΓt MzQ (t), Γm F(t), and MzQ (t) in figure 3, respectively, for a contrast. Figure 3 implies that the coherent fluctuation comes mainly from the ordinary part, i.e. MzQ (t) and will be exponentially suppressed by the relaxation. Otherwise, under the dissipation, since Mz (∞) = Γm F(∞), keeping the discrepancy between the two mechanisms would bring the magnetization towards a nonzero result, which by coincidence could be equal to the GGE value, meaning Γm F(∞) = MzQ (∞). However, maybe the same relaxation rate Γt for every momentum mode k makes it hard to mark out a quasi-stationary state or a prethermal stage, this point will be reaffirmed in discussion on the kink density below. If Γd or Γp = 0 , we can obtain the dissipative steady value of the transverse magnetization per spin in the thermodynamic limit

3.2.  Transverse magnetization

Another noteworthy observable is the transverse magnetization per spin which reads: Mz (t) ≡ Tr[σ z ρ(t)] = e−tΓt MzQ (t) + Γm F(t), (15)

here

k>0

describes the drifting of the system towards an asym4 sin2 θk 2k Γ3t +4Γt 2k

−Γm Γt

π

(18) which captures a DPT in the NESS at λ = 1 (see figure  4). Analogous to the one of [15], in our case the criticality will also be suppressed by increasing the asymptotic decay rate Γt , but which degenerates into a constant here. Actually, setting Γp = 0 and expanding Mz (∞) with Γd: Mz (∞) = Mz0 (∞) + O2z (∞)Γ2d + · · ·, we find Mz0 (∞) = 12 for λ < 1 and 1 − 2λ1 2 for λ > 1, meaning that the criticality is indeed embed in the zero order of Γd.

e−tΓt − 1 2  (ak − e−tΓt bk (t)) F(t) = + (17) Γt N

potic steady state if Γd = Γp , here ak =



dk 4(λ − cos k)2 + Γ2t π 4(1 + λ2 − 2λ cos k) + Γ2t 0  2 2 Γm Γt − 12λ + 4 − (4(λ − 1)2 + Γ2t )(4(λ + 1)2 + Γ2t ) , = Γt 16λ2

Mz (∞) =

2  MzQ (t) = (cos ∆θk cos θk + cos(2tk ) sin ∆θk sin θk ) N k>0 (16) is the ordinary one only due to the sudden quench in closed Ising chain, and subject to a exponential decay. The additional correction

and

5

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402

determines relaxation scale, but no any distinct non-thermal quasi-stationary states exist. From equation (19) we can read off the steady value of the kink density per spin in the thermodynamic limit as  π 1 dk 2Γm λ sin2 k NK (∞) = + 2 π Γt (4 + 4λ2 − 8λ cos k + Γ2t ) 0  Γm 1  ( (4(λ − 1)2 + Γ2t )(4(λ + 1)2 + Γ2t ) = − 2 32λΓt (20) − 4 − 4λ2 − Γ2t ),

which signals the quantum criticality smeared out by increasing the total damping rate Γt like the magnetization (see figure  6). In the case Γp = 0 , expanding NK (∞) with Γd as NK (∞) = NK0 (∞) + O(∞)2K Γ2d + · · ·, we have 1 NK0 (∞) = 12 − λ4 for λ < 1 and = 12 − 4λ for λ > 1. Finally, we details the population of the kink density in the thermodynamic limit for every k mode:

Figure 5.  The dissipative evolution of the kink density following the sudden quench under different dissipative mechanisms. The black line without dissipation gives a contrast. We set λ0 = 0.5, λ = 4, and N  =  1000.

nk (∞) = Tr[γk† (λ = 0)γk (λ = 0)ρk (∞)] 1 Γm (2λ(3 + cos 2k) − (4 + 4λ2 + Γ2t ) cos k) k + sin2 , 2 2 π 2 2πΓt (4 + 4λ − 8λ cos k + Γt )  (21) by the way pointing out that the criticality originates from the edge zore mode k  =  0 (see the inset of figure 6). =

4. Summary In this work we presented a completely solvable scenario for the nonequilibrium dynamics of the open quantum manybody system, that is a quenched Ising chain subject to dissipation via a properly designed coupling to reservoir. The major roles of the decoherence include exponentially suppressing the coherent fluctuation of the observables and driving the system toward a steady state. Two types of damping channels have been considered simultaneously. The competition between them leads to a nontrivial NESS holding a dissipative quantum phase transition. The dependences of the steady observables on the intensity difference between the two damping rates is always linear. The energy is injected into the system by sudden quenching, then increasingly or decreasingly, depending on the amount of the post-quenching and the dissipation rates, approaches its saturation value. The amount of the transient work done on system in our case is equal to the one of the quenched closed Ising chain, while for the noisy Ising chain the work done is persistent, without heat transfer [26]. The quasi-stationary state or prethermalization herein is totally absent perhaps because the trivial asymptotic decay rate Γt [15] only leads to a monotonic relaxation dynamics regardless of the short-time coherent fluctuation, distinguishing from the noisy Ising chain that owns different decay rate for every k mode [26]. In fact, how to define the prethermalization and identify its time scales in open quantum systems is still an open issue.

Figure 6.  The steady value of the kink density per spin in the thermodynamic limit as a function of λ with different decaying rates at Γp = 0 (see equation (20)). The inset exhibits the relevant populations nk (∞) = Tr[γk† (λ = 0)γk (λ = 0)ρk (∞)] versus wave vector k at λ = 1.

3.3.  Density of kinks

We now move to the kink (defect) density pre spin  NK (t)(≡ N1 k>0 Tr[γk† (λ = 0)γk (λ = 0)ρk (t)]), which is also a primary observable of interest: 1 1  2kπ † † (c c + c−k ck ) NK (t) = + Tr[(sin 2 N N k −k k>0

1 2kπ † − cos (ck ck + c†−k c−k ))ρk (t)] = N 2 2kπ 14 2kπ 11 2 x (t) − cos (ρk (t) + ρ22 (sin + k (t))). N N k N

(19)

k>0

The information about the dissipation has been included in 22 such variables: xk14 (t), ρ11 k (t), and ρk (t), which all involve complicated expressions (see appendix). Figure 5 provides a visual image for the complete dynamic evolution of the kink density per spin, which shows that total decoherence rate Γt

6

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402

Optionally, we decided on an antiperiodic boundary condition in that the edge effect can be ignored in the thermodynamic limit, but exploring the finite-size relaxation in certain systems with peculiar edges, such as the twisted Kitaev chains [39], may be an intriguing work. In the end, we stress that the analytical display of the whole nonequilibrium dynamics process, despite through a special example, undoubtedly provides a useful reference to more complex circumstances, including nonintegrability or more elusive effects from the environment. A recent experimental exhibition of dissipative phase trans­ition in circuit QED lattices [20] also suggests that some of the constructive results mentioned here may be tested under quantum many-body simulation plus reservoir engineering [40, 41].

xk14 =

e−tΓt Γm sin 2θk (2(etΓt − 1)2k − sin2 (tk )Γ2t − sin(2tk )Γt k ) ( 2 Γ3t + 4Γt 2k

− sin2 (tk ) sin(2θk − θk0 ) − cos2 (tk ) sin θk0 ), y14 k =

e−tΓt Γm sin θk (2etΓt k − sin(2tk )Γt − 2 cos(2tk )k ) ( 2 Γ2t + 42k

− sin(2tk ) sin ∆θk ), ρ11 k =

e−2tΓt ((1 + e2tΓt − 2etΓt (cos ∆θk cos θk + 42k )

4Γ2t (Γ2t

+ cos(2tk ) sin ∆θk sin θk ))Γ2t (Γ2t + 42k ) + Γ2m ((1 + e2tΓt − 2etΓt (cos2 θk + cos(2tk ) sin2 θk ))Γ2t + 4(etΓt − 1)2 cos2 θk 2k )

+ 2Γm Γt ((e2tΓt + cos θk0 − etΓt (cos θk

(cos ∆θk + cos θk ) + cos(2tk ) sin θk (sin ∆θk + sin θk )))Γ2t − 2etΓt sin(2tk )(sin ∆θk − sin θk ) sin θk Γt k + 4

(etΓt − 1) cos θk (etΓt cos θk − cos ∆θk )2k )),

Acknowledgments

33 ρ22 k = ρk =

This work was supported by the National Science Foundation of China (Grant No. 11664021, 61565008, 11365013).

e−2tΓt ((e2tΓt − 1)(Γ2t + 42k )Γ2t 4Γ2t (Γ2t + 42k )

+ Γ2m ((2etΓt cos2 θk + 2etΓt cos(2tk ) sin2 θk − e2tΓt − 1)

Γ2t − 4(etΓt − 1)2 cos2 θk 2k ) + 2Γm Γt (etΓt (sin ∆θk sin θk Γt (cos(2tk )Γt + 2 sin(2tk )k ) + cos ∆θk cos θk (Γ2t

Appendix.  Solving the Lindblad equation

+ 42k )) − cos θk0 Γ2t − 4 cos θk cos ∆θk 2k )),

e−2tΓt ((1 + e2tΓt + 2etΓt (cos ∆θk cos θk + 42k )

Reminding that the ground state assumption at the initial time t  =  0 has restricted the only non-zero components of ρk to xk14 , 11 22 33 44 y14 k , ρk , ρk , ρk , and ρk as shown in equation (9). So in order to solve the Lindblad equation (6) for every k mode space, we can recast it into the set of linear coupled first order differ­ ential equations of dimension six as follows

ρ44 k =

∂t Ak = −Γt Ak , (A.1a)

− 2etΓt sin(2tk ) sin θk (sin ∆θk + sin θk )Γt k

4Γ2t (Γ2t

+ cos(2tk ) sin ∆θk sin θk ))Γ2t (Γ2t + 42k ) + Γ2m ((1 + e2tΓt − 2etΓt (cos2 (tk ) + cos 2θk sin2 (tk )))Γ2t + 4(etΓt − 1)2 cos2 θk 2k ) + 2Γm Γt ((cos θk0 − e2tΓt + etΓt (cos2 (tk )

− cos ∆θk cos θk + cos 2θk sin2 (tk ) − cos(2tk ) sin ∆θk sin θk ))Γ2t

(A.4) + 4 cos θk (etΓt ((1 − etΓt ) cos θk − cos ∆θk ) + cos ∆θk )2k )),  based on which the whole picture of the nonequilibrium dynamics evolution would be gotten in hand.

∂t Bk = Γt (Ck − Bk ) − Γm Dk , (A.1b) ∂t Ck = Γt (1 − 2Ck ) + Γm Dk , (A.1c) ∂t |Pk  = R|Pk  + |Γ, (A.1d)

ORCID iDs

33 33 22 11 44 here Ak ≡ ρ22 k − ρk , Bk ≡ ρk + ρk , Ck ≡ ρk + ρk , and 11 44 Dk ≡ ρk − ρk . Also, the vector |Pk  is defined as T T |Pk  ≡ [xk14 , y14 k , Dk ] , |Γ ≡ [0, 0, Γm ] , and the 3 × 3 matrix R is given by   −Γt 2k cos θk 0  −Γt k sin θk  R = −2k cos θk (A.2) . 0 −4k sin θk −Γt

Xin Li

https://orcid.org/0000-0002-9657-5508

References [1] Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885 [2] Polkovnikov A, Sengupta K, Silva A and Vengalattore M 2011 Rev. Mod. Phys. 83 863 [3] Cheneau M, Barmettler P, Poletti D, Endres M, Schauß P, Fukuhara T, Gross C, Bloch I, Kollath C and Kuhr S 2012 Nature 481 484 [4] Langen T, Erne S, Geiger R, Rauer B, Schweigler T, Kuhnert M, Rohringer W, Mazets I E, Gasenzer T and Schmiedmayer J 2015 Science 348 207 [5] Eisert J, Friesdorf M and Gogolin C 2015 Nat. Phys. 11 124

The formal solution of equation (A.1d) can be easily written as  t |Pk (t) = eRt ( e−Rτ dτ |Γ + |Pk (0)). (A.3) 0

Equation (22a) combined with the initial condition 33 33 22 ρ22 k (0) = ρk (0) means ρk (t) = ρk (t). So from equations  (A.1a)–(A.1d) it is not difficult to obtain the explicit expressions of the solution of equation (6):

7

X Li et al

J. Phys.: Condens. Matter 31 (2019) 235402

[6] Kollar M, Wolf F A and Eckstein M 2011 Phys. Rev. B 84 054304 [7] Marcuzzi M, Marino J, Gambassi A and Silva A 2013 Phys. Rev. Lett. 111 197203 [8] Worm M, Sawyer B C, Bollinger J J and Kastner M 2013 New J. Phys. 15 083007 [9] Essler F H L, Kehrein S, Manmana S R and Robinson N J 2014 Phys. Rev. B 89 165104 [10] Gring M, Kuhnert M, Langen T, Kitagawa T, Rauer B, Schreitl M, Mazets I, Smith D A, Demler E and Schmiedmayer J 2012 Science 337 1318 [11] Smith D A, Gring M, Langen T, Kuhnert M, Rauer B, Geiger R, Kitagawa T, Mazets I, Demler E and Schmiedmayer J 2013 New J. Phys. 15 075011 [12] Mori T, Ikeda T N, Kaminishi E and Ueda M 2018 J. Phys. B: At. Mol. Opt. Phys. 51 112001 [13] Bhaseen M J, Mayoh J, Simons B D and Keeling J 2012 Phys. Rev. A 85 013817 [14] Kessler E M, Giedke G, Imamoglu A, Yelin S F, Lukin M D and Cirac J I 2012 Phys. Rev. A 86 012116 [15] Horstmann B, Cirac J I and Giedke G 2013 Phys. Rev. A 87 012108 [16] Minganti F, Biella A, Bartolo N and Ciuti C 2018 Phys. Rev. A 98 042118 [17] Werner P, Völker K, Troyer M and Chakravarty S 2005 Phys. Rev. Lett. 94 047201 [18] Overbeck V R, Maghrebi M F, Gorshkov A V and Weimer H 2017 Phys. Rev. A 95 042133 [19] Houck A A, Tureci H E and Koch J 2012 Nat. Phys. 8 292 [20] Fitzpatrick M, Sundaresan N M, Li A C Y, Koch J and Houck A A 2017 Phys. Rev. X 7 011016 [21] Müller M, Diehl S, Pupillo G and Zoller P 2012 Adv. At. Mol. Opt. Phys. 61 1

[22] Bernien H et al 2017 Nature 551 579 [23] Aspelmeyer M, Kippenberg T J and Marquardt F 2014 Rev. Mod. Phys. 86 1391 [24] Gil-Santos E, Labousse M, Baker C, Goetschy A, Hease W, Gomez C, Lemaître A, Leo G, Ciuti C and Favero I 2017 Phys. Rev. Lett. 118 063605 [25] Smelyanskiy V N, Venturelli D, Perdomo-Ortiz A, Knysh S and Dykman M I 2017 Phys. Rev. Lett. 118 066802 [26] Marino J and Silva A 2014 Phys. Rev. B 89 024303 [27] Lorenzo S, Apollaro T, Palma G M, Nandkishore R, Silva A and Marino J 2018 Phys. Rev. B 98 054302 [28] Haug H and Jauho A-P 1996 Quantum Kinetics in Transport and Optics of Semiconductors (Berlin: Springer) [29] Schollwöck U 2005 Rev. Mod. Phys. 118 259 [30] Žunkovič B and Prosen T 2010 J. Stat. Mech. P08016 [31] Prosen T 2008 New J. Phys. 10 043026 [32] Prosen T and Žunkovič B 2010 New J. Phys. 12 025016 [33] Keck M, Montangero S, Santoro G E, Fazio R and Rossini D 2017 New J. Phys. 19 113029 [34] Kibble T W B 1976 J. Phys. A: Math. Gen. 9 1387 [35] Zurek W H 1985 Nature 317 505 [36] Bandyopadhyay S, Laha S, Bhattacharya U and Dutta A 2018 Sci. Rep. 8 11921 [37] Carmele A, Heyl M, Kraus C and Dalmonte M 2015 Phys. Rev. B 92 195107 [38] Talkner P, Lutz E and Hänggi P 2007 Phys. Rev. E 75 050102 [39] Kawabata K, Kobayashi R, Wu N and Katsura H 2017 Phys. Rev. B 95 195140 [40] Verstraete F, Wolf M M and Cirac J I 2009 Nat. Phys. 5 633 [41] Murch K W, Vool U, Zhou D, Weber S J, Girvin S M and Siddiqi I 2012 Phys. Rev. Lett. 109 183602

8