Construction of the steady state density matrix and quasilocal charges for the spin-1/2 XXZ chain with boundary magnetic fields

# Construction of the steady state density matrix and quasilocal charges for the spin-1/2 XXZ chain with boundary magnetic fields

Chihiro Matsui and Tomaž Prosen
Graduate School of Mathematical Sciences, The University of Tokyo
3-8-1, Komaba, Meguro-ku, 153-8914 Tokyo, Japan
Department of Physics, Faculty of Mathematics and Physics,
University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia
July 27, 2019

Abstract

We construct the nonequilibrium steady state (NESS) density operator of the spin- XXZ chain with non-diagonal boundary magnetic fields coupled to boundary dissipators. The Markovian boundary dissipation is found with which the NESS density operator is expressed in terms of the product of the Lax operators by relating the dissipation parameters to the boundary parameters of the spin chain. The NESS density operator can be expressed in terms of a non-Hermitian transfer operator (NHTO) which forms a commuting family of quasilocal charges. The optimization of the Mazur bound for the high temperature Drude weight is discussed by using the quasilocal charges and the conventional local charges constructed through the Bethe ansatz.

## 1 Introduction

The Heisenberg spin chain with anisotropy, the so-called XXZ chain, is one of the most widely studied quantum systems. Its integrability allows us to diagonalize the Hamiltonian [1, 22], which leads to the derivation of exact physical quantities such as correlation functions [9, 20, 7, 6]. The integrability comes from the decomposability of many-body scatterings into a sequence of two-body scatterings. The decomposability is guaranteed by the Yang-Baxter equation satisfied by the scattering matrices.

However, it is still a challenging problem to unveil the model’s nonequilibrium properties. Peculiar nonequilibrium behaviors of integrable systems result from the existence of sufficiently many conserved quantities so that they are in one-to-one correspondence with the degrees of freedom. For instance, it was predicted that observables after relaxation are described by the generalized Gibbs ensemble (GGE) [18, 17] by maximizing the entropy, instead of the grand canonical ensemble. Subsequently, it has been demonstrated that integrable systems exhibit generalized thermalization [18, 23]. Another interesting question is, whether an integrable system exhibits ballistic transport at finite or high temperatures. The optimized lower bound on the ballistic transport coefficient – the so-called Drude weight – has been introduced [24] by using a commuting family of local charges, although this bound generically vanishes when the system possesses the -symmetry with respect to which the transporting current is odd.

A new and fruitful approach to the above questions came recently with exact solutions of boundary driven open quantum systems [13, 12]. The idea to derive the explicit nonequilibrium steady state (NESS) density operator lies in the matrix product ansatz (MPA), which was originally introduced for constructing the ground state of the spin chain [8] and widely used to solve classical boundary driven diffusive many-body systems in one-dimension [3]. The matrix product states serving as the NESS density operator for the quantum boundary driven system is given by the non-Hermitian transfer operator (NHTO) constructed from the product of the Lax operators but with the complex spin representation for the auxiliary space [13, 12]. Due to the Yang-Baxter equation, the NHTO satisfies the divergence condition [21, 15], which implies bulk cancellation in the steady-state Lindblad master equation with the remainder terms localized at the boundaries. These terms are in turn compensated by the boundary dissipation, which turns the NHTO into the NESS density operator. Amazingly, the NESS density operator constructed in this way can serve, in the limit of small dissipation, as a novel quasilocal conserved quantity, which can be used to evaluate the lower bound for the high temperature Drude weight on the corresponding non-dissipative quantum system [16, 5]. Indeed, the NESS satisfies both commutativity and quasilocality conditions [16, 14]. What made the long-standing problem, i.e. to find the NESS of the quantum boundary-driven diffusive many-body system, solved is the complex spin representation of the auxiliary space, which has not been considered before. This extension turned out to be important also in the context of integrable systems, since the NHTO provides a new commuting family of conserved quantities which contain parts that are orthogonal to the known local conserved quantities.

We aim in this paper to investigate how boundary magnetic fields imposed on the spin chain affect its nonequilibrium behavior. We derive the NESS density operator of the boundary-driven quantum spin chain with arbitrary boundary magnetic fields. Interestingly, there always exists the corresponding boundary dissipation for which the NESS density operator is given in terms of the NHTO by properly choosing the dissipation rates as functions of boundary fields. We have also constructed the quasilocal charges for the corresponding non-dissipative spin chain. We showed that, even under the existence of arbitrary boundary magnetic fields, the NHTO forms a commuting family by keeping quasilocality. These two properties of the NHTO allow us to evaluate the lower bound for the high temperature Drude weight. We found optimization of the optimized Mazur bound due to the -symmetry breaking in the spin chain with non-diagonal boundary magnetic fields, which leads to a finite contribution of the conventional local charges to the lower bound.

This paper is organized as follows. In the next section, we give the basics of the open spin- XXZ chain and the boundary dissipator of the Lindblad type. We derive the NESS density operator in Section 3. The quasilocal charges are constructed in Section 4. The lower bound for the high temperature Drude weight is also evaluated. The last section is devoted to the concluding remarks.

## 2 The open spin-1/2 XXZ chain with boundary dissipation

### 2.1 The spin-1/2 XXZ chain with non-diagonal boundaries

Let us consider the spin- XXZ chain with arbitrary boundary magnetic fields:

 H=n−1∑x=11l2x−1⊗h⊗% 1l2n−x−1+hB,L⊗1%l2n−1+1l2n−1⊗hB,R, (1)

where represents the -fold tensor product of the identity matrix. The Hamiltonian density for the bulk part is expressed by the Pauli matrices () as

 h=2σ+⊗σ−+2σ−⊗σ++σz⊗σzcosη, (2)

while the boundary Hamiltonian density is expressed as

 hB,L=12σzsinηcotξL+σ+κLeθLsinηsinξL+σ−κLe−θLsinηsinξL,hB,R=12σzsinηcotξR+σ+κReθRsinηsinξR+σ−κRe−θRsinηsinξR, (3)

containing six free parameters and which uniquely parametrise arbitrary boundary magnetic fields.

The model (1) is known to be integrable in the sense that its transfer matrix forms a commuting family of infinitely many local operators. The local charges are obtained by expanding the logarithm of the transfer matrix around the permutation point. The leading term gives the momentum operator, while the next-to-leading term gives the Hamiltonian:

 H=ddφ(sinη2sinξLK1(φ,ξL)+n−1∑x=12ˇRx,x+1(φ))∣∣∣φ=η2+tr0K0(φ+η,ξR)hn,0trK0(φ+η,ξR)∣∣∣φ=η2. (4)

The -matrix and the -matrix satisfy the so-called relation and the reflection relation, respectively:

 ˇR1,2(φ1−φ2)L1(φ1)L2(φ2)=L1(φ2)L2(φ1)ˇR1,2(φ1−φ2), (5) ˇR2,1(φ1−φ2)K2(φ1)ˇR1,2(φ1+φ2)K2(φ2)=K2(φ2)ˇR2,1(φ1+φ2)K2(φ1)ˇR1,2(φ1−φ2), (6)

whose solutions are expressed in terms of the Pauli matrices :

 ˇR(φ)=sinφ2(h+1lcosη)−1+cosφ21lsinη+1−cosφ2σz⊗σzsinη, (7) K(φ;ξ,κ,θ)=1lsinξcosφ+σzcosξsinφ+σ+κeθsin(2φ)+σ−κe−θsin(2φ). (8)

On the other hand, the Lax operator consists of the physical space, which has the same representation as that of the -matrix, and the auxiliary space, which admits any representation including the complex spin representation:

 L(φ,s)=(sin(φ+ηsza)sinηs−asinηs+asin(φ−ηsza))=∑α∈{+,−,0,z}Lαa(φ,s)⊗σαp, (9)

where

 L0a(φ,s)=sinφcos(ηsza),Lza(φ,s)=cosφsin(ηsza),L±a(φ,s)=(sinη)s∓a. (10)

We take the complex spin representations in the way introduced in Ref. [15]:

 sza=∞∑k=0(s−k)|k⟩⟨k|,s+a=∞∑k=0sin(k+1)ηsinη|k⟩⟨k+1|,s−a=∞∑k=0sin(2s−k)ηsinη|k+1⟩⟨k| (11)

by which the algebraic relations are satisfied:

 [sz,s±]=±s±,[s+,s−]=[2sz]q. (12)

### 2.2 Boundary dissipation of the Lindblad type

Within the theory of open quantum systems [2], incoherent Markovian quantum dissipation is completely described by a set of Lindblad operators . Such a system’s many-body density operator obeys the time evolution described by the Lindblad-Gorini-Kossakowski-Sudarshan master equation [4, 10]:

 (13)

Throughout this paper, we consider the ultra-local Lindblad operators

 Lμ=ℓμ⊗1l2n−1orLμ=1l2n−1⊗ℓμ,ℓμ∈End(Hp), (14)

especially of the following forms associated with three different dissipation rates and two additional dissipaiton parameters :

 L1=√ε(σ+1+ασ01),L2=√ε(σ−n+α′σ0n),L3=√ε′σz1,L4=√ε′′σzn. (15)

Notice that the dissipators and are coupled to the left boundary of the spin chain, while and to the right boundary. We let and free for the moment and determine their relations to the boundary parameters later.

The Liouvillian is then written as

 ^L=−i[H,ρ(t)]+ε^Dσ+1+ασ01+ε′^Dσz1+ε^Dσ−n+α′σ0n+ε′′^Dσzn, (16)

where the dissipator map is defined by

 ^DL(ρ)=2LρL†−{L†L,ρ}. (17)

From the master equation (13), the density operator at time is expressed as leading to, if the limit exists, the NESS density operator . Since the NESS is invariant under the time development, the NESS density operator gives the fixed point of the propagator:

 ^Lρ∞=0. (18)

### 3.1 Construction of the NESS density operator

Let us first write the NESS density operator in terms of the product of non-Hermitian amplitude operators:

 ρ∞=R∞trR∞,R∞=ΩnΩ†n. (19)

Contrary to the trivial open boundary case where a particular factorization occurs in terms of  [15], some modification is required under the presence of boundary magnetic fields.

The fixed point condition (18) for the NESS density operator under the choice of the Liouvillian (16) leads to the following condition on and :

 iε−1[H,ΩnΩ†n]=^Dσ+1+ασ01(ΩnΩ†n)+ε′ε^Dσz1(ΩnΩ†n)+^Dσ−n+α′σ0n(ΩnΩ†n)+ε′′ε^Dσzn(ΩnΩ†n). (20)

In order to deal with the product of the amplitude operators, it is useful to introduce the double Lax operator acting over a tensor product of a pair of complex spin representations [15]. Note that the complex spin operators have the highest weight representations given in (11). Besides, we introduce another complex spin operators of the transposed lowest weight representations:

 tz=∞∑k=0(t−k)|k⟩⟨k|,t+=∞∑k=0sin(k+1)ηsinη|k+1⟩⟨k|,t−=∞∑k=0sin(2t−k)ηsinη|k⟩⟨k+1|. (21)

The double Lax operators is then defined in the product representation :

 Lx(φ,θ,s,t)=LTa,x(φ,s)Lb,x(θ,t), (22)

where is the transposed lowest weight representation, while is the highest weight representation.

The double Lax operator obeys the Yang-Baxter-like equation:

 (23)

from which we obtain the commutativity condition for the bulk part:

 [Hbulk,L1⋯Ln]=2sinη(∂L1L2⋯Ln−L1⋯Ln−1∂Ln), (24)

where

 ∂Lx(φ,θ,s,t)=∂δ(LTa,x(φ+δ,s)Lb,x(θ−δ,t))δ=0=∂φLTa,x(φ,s)Lb,x(θ,t)−LTa,x(φ,s)∂θLb,x(θ,t). (25)

The bulk commutativity condition implies that the dissipation should be localized at the boundaries of the spin chain. Using the bulk commutativity condition (24) and the boundary Hamiltonian density (3), we obtain the commutativity condition for the full spin chain:

 [H,L⊗n]=[Hbulk,L⊗n]+[HL,B,L⊗n]+[HR,B,L⊗n]=2(sinη)∂L⊗L⊗n−1−2(sinη)L⊗n−1⊗∂L+[12σzsinηcothξL+σ+κLeθLsinηsinξL+σ−κLe−θLsinηsinξL,L]⊗L⊗n−1+L⊗n−1⊗[12σzsinηcothξR+σ+κReθRsinηsinξR+σ−κRe−θRsinηsinξR,L]. (26)

We assume that the NESS density operator is given by the product of the Lax operators, similarly as for the trivial open boundary case [14]:

 Ωn=1sinn(φ+sη)a⟨0|LT1LT2…LTn|0⟩a, (27)

where is the highest weight vector. Note that we have used the partially transposed Lax operator:

 LT(φ,s)=∑α∈{+,−,0,z}Lα(φ,s)⊗(σα)T. (28)

From the definition (22), the product admits the expression in terms of the double Lax operator:

 ΩnΩ†n=1sinn(φ+sη)sinn(φ+¯sη)a⟨0|b⟨0|L1⋯Ln|0⟩a|0⟩b. (29)

By applying the double highest weight vector to the fixed point condition (20), we obtain the left and right boundary conditions:

 a⟨0|b⟨0|(−iε−12(sinη)∂L−iε−1[12σzsinηcothξL+σ+κLeθLsinηsinξL+σ−κLe−θLsinηsinξL,L]+^Dσ+1+ασ01(L)+ε′ε^Dσz1(L))=0,(iε−12(sinη)∂L−iε−1[12σzsinηcothξR+σ+κReθRsinηsinξR+σ−κRe−θRsinηsinξR,L]+^Dσ−n+α′σ0n(L)+ε′′ε^Dσzn(L))|0⟩a|0⟩b=0. (30)

Note that the dissipation operators act on a -by- matrix as

 ^Dσ++ασ0(abcd)=(2d+α(b+c)−b−α(a−d)−c−α(a−d)−2d−α(b+c)), (31) ^Dσ−+ασ0(abcd)=(−2a−α(b+c)−b+α(a−d)−c+α(a−d)2a+α(b+c)), (32) ^Dσz(abcd)=(0−4b−4c0). (33)

The nontrivial solution to (30) exists for the following case. We first need the spectral parameter restricted by

 φ=π2. (34)

This condition coincides with the trivial open boundary case [15]. The dissipation rates must be related to the anisotropy and the boundary parameters of the diagonal parts:

 ε=−2isinηtan(sη),ε′=−i4sinηcotξL,ε′′=−i4sinηcotξR. (35)

This result includes the trivial open boundary case for  [14]. The parameters and are determined by the boundary parameters of the non-diagonal parts:

 α=iε−1sinηsinξLκLeθL=−iε−1sinηsinξLκLe−θL,α′=−iε−1sinηsinξRκReθR=iε−1sinηsinξRκRe−θR, (36)

which require . As a result, the boundary Hamiltonian density is allowed to take only the following form:

 hB,L=12σzsinηcotξL−2σyκLsinηsinξL, (37) hB,R=12σzsinηcotξR−2σyκRsinηsinξR. (38)

### 3.2 Rotation in the xy-plane

One may wonder why the solution to the boundary conditions requires the non-diagonal parts of the left and right boundary Hamiltonian to consist of only -terms, although the Hamiltonian itself possesses the symmetry with respect to the rotation in the -plane. Let us consider the transformation performed by the matrix defined by

 U:=⎛⎜⎝eiϕ200e−iϕ2⎞⎟⎠. (39)

The transformation makes invariant but rotates the other components of the Pauli matrices:

 UσxU−1=σxcosϕ+σysinϕ, (40) UσyU−1=σxsinϕ−σycosϕ, (41)

and, subsequently,

 Uσ+U−1=eiϕσ+,Uσ−U−1=e−iϕσ−. (42)

Therefore, it is easy to show that the bulk part of the Hamiltonian is invariant under the transformation . Also, the -terms of the left and the right boundary Hamiltonian are invariant under the transformation. Thus, rotates the Hamiltonian in the -plane by letting the non-diagonal boundary terms contain both - and -terms.

Let us denote the transformed operator by . The boundary Hamiltonian density is then deformed as

 hB,L→˜hB,L=12σzsinηcotξL−2(σxsinϕ−σycosϕ)κLsinηsinξL, (43) hB,R→˜hB,R=12σzsinηcotξR−2(σxsinϕ−σycosϕ)κRsinηsinξR. (44)

The transformation angle is related to the non-diagonal boundary parameters as

 iθL,R=ϕ+π2, (45)

which implies that the only condition for obtaining the NESS density operator is . Indeed, the fixed point condition (20) holds for the transformed Hamiltonian and double Lax operator if we deform the Lindblad operators as

 ^Dσ++ασ0→^D˜σ++α˜σ0=^Deiϕσ++ασ0,^Dσ−+ασ0→^D˜σ−+α˜σ0=^De−iϕσ−+ασ0,^Dσz→^D˜σz=^Dσz. (46)

## 4 Quasilocal charges

### 4.1 Optimized Mazur bound on the Drude weight

The ballistic transport is characterized by the finite Drude weight . The Drude weight is defined as the coefficient for the diverging part of d.c. conductivity in the context of linear-response transport. Since the Drude weight has the temperature dependence, its expansion at high temperature leads to

 D=βD∞+O(β2),D∞=limt→∞limn→∞12tn∫t0dt′(J(t′),J), (47)

where is the extensive current , is the local current, and is the Hilbert-Schmidt inner product defined by .

By using a set of extensive local conserved operators , the Mazur bound [11] allows rigorous estimation of the high-temperature Drude weight from below [24]:

 D∞≥limn→∞12n∑r,r′(J,Qr)(K−1)r,r′(Qr′,J):=DQ,Kr,r′:=(Qr,Qr′). (48)

However, in the system with the -symmetry such as the spin-reversal symmetry, with respect to which the conserved charges are even and the current is odd, the lower bound always becomes zero, which tells nothing about the ballistic behavior.

The situation drastically changes by the introduction of quasilocal charges [16, 14]. The quasilocal charges constructed by differentiating the NESS density operator for the corresponding boundary driven spin chain with respect to the representation parameter consists of both even and odd parity parts, which makes the lower bound finite. The optimized Mazur bound [14, 15] was obtained as

 D∞≥12Re∫Dmd2φZJ(φ)f(φ):=DZ, (49)

which is bounded by the conventional Drude weight . Here we used the notation , and is the family of quasilocal charges which can be generated from NESS density operator of the boundary driven chain in the limit of weak driving. The function is obtained as the solution of a complex Fredholm equation of the first kind:

 ∫Dmd2φ′limn→∞1n(Z(φ),Z(φ′))f(φ′)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯ZJ(φ),φ∈Dm⊂C. (50)

### 4.2 Construction of quasilocal charges

Thus, what we need for evaluating the lower bound of the high temperature Drude weight is to find a commuting family of quasilocal charges. Analogously to the periodic and trivial open boundary cases [14, 15], let us introduce the NHTO defined by the product of the Lax operators:

 Wn(φ,s)=⟨0|L(φ,s)⊗pn|0⟩. (51)

Note that the commutativity holds for the NHTO with any pair of spectral and representation parameters in spite of the existence of non-diagonal boundary magnetic fields:

 Wn(φ,s)Wn(θ,t)=⟨0|a⟨0|bˇRa,b(φ−θ,s,t)(n∏x=1La,x(φ,s)Lb,x(θ,t))|0⟩a|0⟩b=⟨0|a⟨0|b(n∏x=1La,x(φ,s)Lb,x(θ,t))ˇRa,b(φ−θ,s,t)|0⟩a|0⟩b=Wn(θ,t)Wn(φ,s). (52)

The NHTO satisfies the following commutativity condition:

 (53)

Here we introduced

 W+n(φ,s)=⟨1|L1(φ,s)⋯Ln(φ,s)|0⟩,W−n(φ,s)=⟨0|L1(φ,s)⋯Ln(φ,s)|1⟩. (54)

From now on, we show that the differentiation of the NHTO with respect to the representation parameter forms a family of quasilocal charges. Precisely, the quasilocal charge is given by

 Zn(φ)=12(sinφ)n−2ηsinη∂sWn(φ,s)∣∣s=0−sinφcosφ2sinηMzn, (55)

where . By using the time derivative of NHTO (53), we have

 (56)

The third and fourth lines are evaluated through

 ∂sW+n(φ,s)∣∣s=0=n∑k=1⟨1|L1(φ,s)⋯∂sLk(φ,s)⋯Ln(φ,s)|1⟩∣∣s=0=n∑k=1(sinφ)n−k∑αj∈{+,−,0,z}⟨1|Lα11(φ,s)…(∂sLk(φ,s))αk|0⟩∣∣s=0σα1⊗⋯⊗σαk⊗1l2n−k=n∑k=12η(sinφ)n−k∑αj∈{+,−,0,z}⟨1|Lα11(φ,0)…Lαk−1k−1(φ,0)|1⟩σα1⊗⋯⊗σαk−1⊗σ+⊗1l2n−k, (57)
 W−n(φ,0)=⟨0|L1(φ,0)⋯Ln(φ,0)|1⟩=∑αj∈Jsinφ⟨0|Lα22(φ,0)⋯Lαnn(φ,0)|1⟩σ0⊗σα2⊗⋯⊗σαn+∑αj∈{+,−,0,z}sinη⟨1|Lα22(φ,0)⋯Lαnn(φ,0)|1⟩σ−1⊗σα22⊗⋯⊗σαnn=n∑k=1sinη(sinφ)k−1∑αj∈{+,−,0,z}⟨1|Lαkk(φ,0)⋯Lαnn(φ,0)|1⟩1l2k−1⊗σ−⊗σαk+1⊗⋯⊗σαn, (58)

by using

 Lj(φ,0)|0⟩=σ0jsinφ|0⟩,∂sLj(φ,s)|s=0|0⟩=2ησ+j|1⟩,⟨0|Lj(φ,0)=σ0jsinφ⟨0|+σ−jsinη⟨1|, (59)

and subsequently,

 Wn(φ,0)=(sinφ)n1l2n,W+n(φ,0)=0. (60)

Setting

 qr(φ)=(sinφ)−r+2∑αj∈J⟨1|Lα2(φ,0)⋯Lαr−1(φ,0)|1⟩σ−⊗σα2⊗⋯⊗σαr−1⊗σ+, (61) p+r(φ)=(sinφ)−r+2∑αj∈{+,−,0,z}⟨1|Lα2(φ,0)⋯Lαr−1(φ,0)|1⟩σz⊗σα2⊗⋯⊗σαr−1⊗σ+, (62) p−r(φ)=(sinφ)−r+2∑αj∈{+,−,0,z}⟨1|Lα2(φ,0)⋯Lαr−1(φ,0)|1⟩σ−⊗σα2⊗⋯⊗σαr−1⊗σz (63)