Dissipative Quantum Ising chain as a non-Hermitian Ashkin-Teller model
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J un Dissipative Quantum Ising chain as a non-Hermitian Ashkin-Teller model
Naoyuki Shibata ∗ and Hosho Katsura
1, 2 Department of Physics, Graduate School of Science,the University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan Institute for Physics of Intelligence, The University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan (Dated: July 2, 2019)We study a quantum Ising chain with tailored bulk dissipation, which can be mapped onto anon-Hermitian Ashkin-Teller model. By exploiting the Kohmoto-den Nijs-Kadanoff transformation,we further map it to a staggered XXZ spin chain with pure-imaginary anisotropy parameters. Thisallows us to study the eigenstates of the original Liouvillian in great detail. We show that thesteady state in each parity sector is a completely mixed state. The uniqueness of each is provedrigorously. We then study the decay modes on the self-dual line corresponding to the uniform XXZchain and obtain an exact formula for the Liouvillian gap g , the inverse relaxation time, in thethermodynamic limit. The gap g as a function of dissipation strength ∆ has a cusp, implying a kindof phase transition for the first decay mode. I. INTRODUCTION
Open quantum systems have recently attracted muchattention in a variety of fields including condensed mat-ter physics [1–5], quantum information [6, 7], quantumcomputing [8], and mathematical physics [9]. The Lind-blad equation [10] is a general quantum master equationdescribing such open quantum systems under Markovianand Completely Positive and Trace Preserving (CPTP)conditions. In general, the analysis of the Lindblad equa-tion is more challenging than that of the Schr¨odingerequations for closed systems, as one has to deal withthe space of opertors rather than the Hilbert space.One possible approach is to develop some approximatemethods such as the perturbative [11, 12] or numericalones [13, 14]. In particular, due to the recent develop-ment of machine learning approaches [15–18], the numberof these studies has been increasing. Another way whichis complementary to the above methods is to constructexactly solvable models. Although many such modelsare known in closed many-body systems, very few exactresults are available for open many-body systems [9, 19–23].The quantum Ising chain is a paradigmatic exampleof an exactly solvable model for closed systems and alsoserves as a textbook example of a quantum phase tran-sition [24, 25]. It is mapped to a free-fermion model,and hence integrable, which allows for explicit compu-tation of various quantities. With the recent surge ofinterest in open quantum systems, a number of studieson dissipative quantum Ising models have been reportedrecently [26–29]. However, to the best of our knowledge,exact solutions are available only for models subject todissipation at the end of the chain [30].In this paper, we present a dissipative dissipative quan-tum Ising model with bulk dissipation, which is inte-grable with judiciously chosen dissipators and param- ∗ [email protected] eters. The Hamiltonian and dissipators of the modelconserve parity, and hence it has two degenerate non-equilibrium steady states (NESSs). By vectorizing thedensity matrix, one can identify the Liouvillian of thesystem with a non-Hermitian analog of the Ashkin-Tellermodel, the Hermitian counterpart of which was studiedin Refs. [31–35]. Due to the Z × Z symmetry of themodel, the space of states splits into four sectors labeledby the eigenvalues of the parity operators. In each sector,the Hamiltonian of the model can be mapped to that ofa staggered XXZ chain with pure-imaginary anisotropyparameters. This enables us to prove that the two NESSsare unique. By further exploiting this mapping, we inves-tigate the Liouvillian gap g , the inverse relaxation time,and the corresponding first decay mode on the self-dualline at which the bulk Hamiltonian of the XXZ chain isspatially uniform. Due to the non-local nature of thetransformation, the boundary terms in the XXZ chaindiffers from sector to sector. With this in mind, weprove rigorously that the Liouvillian gap in two sectorsare exactly 4∆. Furthermore, we find that the gap inthe other two sectors in the thermodynamic limit can beobtained analytically from the Bethe ansatz solution ofthe model [36]. Combining these two results, we derivean explicit formula for the global Liouvillian gap g as afunction of the dissipation strength ∆, which has a cuspat ∆ = 1 / √ Z parafermion chain with non-Hermitian interactions. II. MODEL AND NESSs
Consider the Lindblad equationd ρ d t = L [ ρ ] := − i[ H, ρ ] + X i (cid:18) L i ρL † i − n L † i L i , ρ o(cid:19) (1)for a quantum Ising chain described by the Hamiltonian H = − h N X i =1 σ zi − J N X i =1 σ xi σ xi +1 (2)with the following Lindblad operators L i = p ∆ σ zi ( i = 1 , . . . , N ) , (3) L N + i = p ∆ σ xi σ xi +1 ( i = 1 , . . . N ) . (4)Here ρ is the density matrix, N is the number of sites, σ αi ( i = 1 , . . . , N, α = x, z ) are Pauli matrices at site i , and ∆ , ∆ ≥ σ xN +1 = σ x . Note thatthe Lindblad operators take the same form as the localHamiltonians of H . Since we can change the signs of h and J by an appropriate unitary transformation, we canassume that h, J ≥ parity operator P := N Y i =1 σ zi , (5)which satisfies[ H, P ] = 0 , [ L i , P ] = 0 ( i = 1 , . . . , N ) . (6)Due to the existence of this conserved charge, there aretwo non-equilibrium steady states (NESSs) ρ ± := 1 ± P N , (7)which are eigenstates of the Liouvillian with eigenvalue0: L [ ρ ± ] = 0 . (8)We briefly explain this construction of NESSs accord-ing to Ref. [23]. Since all Lindblad operators are Her-mitian, there is a trivial NESS, i.e., completely mixedstate ρ c := / N . By noting that P ρ c ∝ P also satis-fies L [ P ] = 0, we see that ρ ± obeys Eq. (8). We notethat the operators ρ ± are positive semi-definite, as theireigenvalues are nonnegative. (a)(b) FIG. 1. Schematic representations of (a) the set-up and(b) the Liouvillian (times i) on the Ket ⊗ Bra space [seeEqs. (9) and (10)]. Rectangles in (a) denotes reservoirs. Solidand wavy lines in (b) are referred to as Hermitian and non-Hermitian interactions, respectively.
III. MAPPING TO THE NON-HERMITIANSTAGGERED XXZ MODELA. Non-Hermitian Ashkin-Teller model andstaggered XXZ model
A 2 N × N density matrix ρ can be thought of as a 2 N -dimensional vector [12, 23, 37, 38]. In this sense, we canidentify the Liouvillian L (times i) as a non-HermitianHamiltonian on a “Ket ⊗ Bra space” asi
L ∼ = H ⊗ − ⊗ H T + i X i (cid:18) L i ⊗ L ∗ i − L † i L i ⊗ − ⊗ L T i L ∗ i (cid:19) , (9)where the Hilbert space of the RHS is the “Ket ⊗ Braspace”. For our model, the corresponding non-Hermitian
FIG. 2. Mapping procedure of our model. The left, middle, and right figures refer to Eqs. (11), (17), and (24), respectively. Inthe left and middle panels, the dashed lines indicate the ordering paths used for the Jordan-Wigner transformation.
Hamiltonian readsi L + const. ∼ = H = H bulk + H boundary , (10) H bulk = − h N X i =1 σ zi − J N − X i =1 σ xi σ xi +1 + h N X i =1 τ zi + J N − X i =1 τ xi τ xi +1 + i∆ N X i =1 σ zi τ zi + i∆ N − X i =1 σ xi σ xi +1 τ xi τ xi +1 , (11) H boundary = Jσ xN σ x + Jτ xN τ x + i∆ σ xN σ x τ xN τ x (12)where τ αi ( α = x, z ) are the Pauli matrices for the i th Bra site (see Fig. 1 (b)). First, let us concen-trate on the bulk Hamiltonian H bulk . It corresponds tothe quantum Ashkin-Teller model [31, 32] with imagi-nary anisotropy parameters by an appropriate unitarytransformation [39]. Furthermore, it is mapped to thestaggered XXZ model. It was already mentioned inRefs. [31, 32], but let us derive it in another way usingMajorana fermions. By a Jordan-Wigner transformation σ zi = − i c i − c i , (13) σ xi = i − Y j =1 − i c j − c j c i − , (14) τ zi = − i d i − d i , (15) τ xi = N Y j =1 − i c j − c j i − Y k =1 − i d k − d k ! d i − , (16) the model is mapped to the interacting Majorana fermionmodel with the Hamiltonian (Fig. 2) H bulk = i h N X i =1 c i − c i + i J N − X i =1 c i c i +1 − i h N X i =1 d i − d i − i J N − X i =1 d i d i +1 − N X i =1 i∆ c i − c i d i − d i − N − X i =1 i∆ c i c i +1 d i d i +1 . (17)Here, c i and d i ( i = 1 , . . . , N ) are Majorana operatorsobeying the following relations c † i = c i , d † i = d i , (18) { c i , c j } = { d i , d j } = 2 δ ij , { c i , d j } = 0 . (19)Now we apply another Jordan-Wigner transformationand rewrite the Majorana fermions as c i − = ( − i i − Y j =1 Z j X i − , (20) d i − = ( − i i − Y j =1 Z j Y i − , (21) c i = ( − i i − Y j =1 Z j Y i , (22) d i = ( − i i − Y j =1 Z j X i , (23)where X i , Y i and Z i ( i = 1 , . . . N ) are Pauli operatorsat site i . One can express these new Pauli operators asnon-local products of the original ones, i.e., σ j and τ j .Using the new Paulli operators, H bulk can be recast as H bulk = N X i =1 [ h ( X i − X i + Y i − Y i ) + i∆ Z i − Z i ]+ N − X i =1 [ J ( X i X i +1 + Y i Y i +1 ) + i∆ Z i Z i +1 ] , (24)which is nothing but the Hamiltonian of the staggeredXXZ model, with the caveat that the two anisotropyparameters are purely imaginary [40]. When h = J and ∆ = ∆ , the model reduces to the uniform XXZchain with pure imaginary anisotropy parameter, whichis further investigated in Sec. V. From another point ofview, the bulk Hamiltonian can be thought of as a Z parafermion chain with non-Hermitian interactions (seeAppendix A for more detail). B. Boundary terms
Here we consider the boundary terms in Eqs. (2) and(4), which lead to Eq. (12). It is useful to define thefollowing operators P σ := N Y i =1 σ zi , P τ := N Y i =1 τ zi . (25)They are conserved charges of H satisfying[ H , P σ/τ ] = 0 , [ P σ , P τ ] = 0 . (26)By the first Jordan-Wigner transformation, they can bewritten in terms of Majorana fermions as P σ = ( − i) N N Y j =1 c j − c j , P τ = ( − i) N N Y j =1 d j − d j , (27)and the boundary terms are mapped to H boundary = − i J ( P σ c N c − P τ d N d ) − i∆ P σ P τ c N c d N d . (28)By the second Jordan-Wigner transformation, we have P σ = ( − N N Y j =1 Y j , P τ = N Y j =1 X j , (29) H boundary = ( − N J ( P τ X N X + P σ Y N Y )+ i∆ P σ P τ Z N Z . (30)We then define a new set of Z charges as Q Z := P σ P τ = N Y j =1 Z j , Q X := P τ = N Y j =1 X j , (31) in terms of which the the boundary terms can be rewrit-ten as H boundary = ( − N JQ X ( X N X + Q Z Y N Y )+ i∆ Q Z Z N Z . (32)Equation (32) clearly shows that the Hilbert space ofthe Hamiltonian H = H bulk + H boundary splits into thefollowing sectors labeled by the eigenvalues of Q Z and( − N Q X :(i) Q Z = +1 , ( − N Q X = +1: periodic(ii) Q Z = +1 , ( − N Q X = −
1: anti-periodic(iii) Q Z = − , ( − N Q X = +1: anti-diagonaltwisted [41, 42](iv) Q Z = − , ( − N Q X = −
1: anti-periodic and anti-diagonal twistedWe now define H boundary ( a, b ) = bJ ( X N X + aY N Y )+ i∆ aZ N Z , (33) e H ( a, b ) = H bulk + H boundary ( a, b ) (34)where a and b are c-numbers, instead of q-numbers, tak-ing ±
1. One can construct the eigenstates of H fromsimultaneous eigenstates of e H ( a, ±
1) and Q Z . This canbe seen by noting that if | φ i is a simultaneous eigenstateof e H ( a, ±
1) and Q Z , then ( | φ i ± ( − N Q X | φ i ) / H with the same eigenvalue.It is known that e H ( a, b ) for any a = ± b = ± e H (+1 , b ) has a U(1) symmetry, i.e., the z com-ponent of the total spin Z tot = P i Z i commutes with e H (+1 , b ). Therefore, eigenvalues and eigenvectors in sec-tor (i) and (ii) can be obtained by means of the algebraicBethe ansatz (ABA). However, for e H ( − , b ), Z tot is nolonger a conserved quantity and the standard ABA failsto work in sector (iii) and (iv). Instead, there is a generalmethod for studying such models without U(1) symmetrycalled the off-diagonal Bethe ansatz (ODBA) [36, 42–44]. IV. PROOF OF THE UNIQUENESS OF THETWO NESSs
It can be proved that all eigenvalues of the Liouvillian L have non-positive real parts [10, 45]. Then, the eigen-value 0 of L , whose corresponding eigenstate is a NESS,has the largest real part. After the above mapping, thecorresponding eigenvalue of H has the largest imaginary part and its real part is 0. It turns out that such statesare ferromagnetic states : |⇑i := |↑ . . . ↑i , |⇓i := |↓ . . . ↓i , (35)where |↑i i (resp. |↓i i ) is the eigenstates of Z i with theeigenvalue +1 (resp. − LEMMA IV.1.
Let H = T + i K , with T † = T and K † = K , be a non-Hermitian matrix, λ i be eigenvaluesof H , and ǫ K be the largest eigenvalue of K . Then,max i (Im λ i ) ≤ ǫ K . (36) Proof.
Let | ψ i i be a right-eigenstate of H with eigenvalue λ i . We can assume that | ψ i i is normalized as h ψ i | ψ i i = 1.Then we have λ i = h ψ i | H | ψ i i = h ψ i | T | ψ i i + i h ψ i | K | ψ i i (37)Since T and K are Hermitian, h ψ i | T | ψ i i and h ψ i | K | ψ i i are real. Thus we find thatIm λ i = h ψ i | K | ψ i i ≤ ǫ K . (38)This holds for all i , so Eq. (36) also holds.For our Hamiltonian H = H bulk + H boundary , K = ( H − H † ) / (2i) (39)= ∆ N X i =1 Z i − Z i + ∆ N − X i =1 Z i Z i +1 + ∆ Q Z Z N Z (40)is already diagonalized. The largest eigenvalue is N (∆ +∆ ) and the corresponding (unique) eigenstates are |⇑i and |⇓i . Moreover, T |⇑i = T |⇓i = 0 , (41)where T = ( H + H † ) / N X i =1 h ( X i − X i + Y i − Y i )+ N − X i =1 J ( X i X i +1 + Y i Y i +1 )+ ( − N JQ X ( X N X + Q Z Y N Y ) (43)Therefore, these ferromagnetic states correspond to thetwo unique NESSs of the original model. One can easilysee that the superposition ( |⇑i + ( − N |⇓i ) / √ |⇑i − ( − N |⇓i ) / √
2) lives in sector (i) (resp. sector(ii)).
V. THE LIOUVILLIAN GAP ON THESELF-DUAL LINE
In this section, we focus on our model on the self-dualline [31] h = J, ∆ = ∆ = ∆ , (44) on which H is invariant under the duality transforma-tion [46] σ zi = e σ xi e σ xi +1 , τ zi = e τ xi e τ xi +1 , σ xi = i Y k =1 e σ zk , τ xi = i Y k =1 e τ zk (45)up to boundary terms. We investigate how the Liou-villian gap, i.e., the inverse relaxation time, and the firstdecay mode behave on the self-dual line. Since the overallscale is not important for the analysis, we set h = J = 1in the following. Let eigenvalues of the Liouvillian L beΛ i ( L ). A Liouvillian gap g is defined as g := − max i Re[Λ i ( L )] =0 Re[Λ i ( L )] , (46)hence, the inverse of the relaxation time. It is clear fromEq. (10) that the Liouvillian gap corresponds to the gapbetween the first and second largest imaginary parts ofeigenvalues of H .As we have seen in Sec. III B and IV, the Hilbert spaceis divided into four sectors (i), (ii), (iii), and (iv), and thetwo ferromagnetic states which correspond to NESSs livein sector (i) and (ii). We now define the local Liouvilliangap g (i) , g (ii) , g (iii) , and g (iv) as follows: g (i) and g (ii) aredefined as the gap between the first and second largestimaginary part of the eigenvalues of e H (+1 , ± N ∆ as we have proved inSec. IV. The other two gaps, g (iii) and g (iv) , are definedas the difference between 2 N ∆ and the largest imaginarypart of the eigenvalues of e H ( − , ± global Liou-villian gap g is then obtained as g = min (cid:0) g (i) , g (ii) , g (iii) , g (iv) (cid:1) . (47)The main result of this section is that we the exact for-mula for g as a function of dissipation strength ∆ in thethermodynamic limit is obtained as g (∆) = (cid:18) < ∆ ≤ √ (cid:19) p ∆ + 1 (cid:18) √ ≤ ∆ (cid:19) , (48)(see also Fig. 3). The result implies a kind of phase tran-sition for the first decay mode. In fact, the first decaymode lives in sector (i) and (ii) in the weak dissipationregion 0 < ∆ ≤ / √
3, while it lives in sector (iii) and(iv) in the strong dissipation region ∆ > / √
3, which wediscuss in the following subsections.
A. In sector (i) and (ii)
First, let us consider the lower bound of g (i) . In Sec.IV, we have shown that the largest imaginary part of H , and also of e H (+1 , +1), is 2 N ∆ and corresponding FIG. 3. (Color online) Liouvillian gap g as a function of ∆. Ablue solid line denotes the exact result Eq. (48) in the thermo-dynamic limit possessing a cusp at ∆ = 1 / √
3, while a orangedashed line is the numerical one with 2 N = 10. eigenstates are ferromagnetic states |⇑i and |⇓i . Then,the state which has the second largest imaginary part of e H (+1 , +1) must have at least one up-spin and at leastone down-spin. Therefore, there are at least two kinks,which together with Lemma IV.1 implies that the secondlargest imaginary part is less than or equal to (2 N − g (i) ≥
4∆ holds. Next, let usexplicitly construct an eigenstate whose imaginary partof the eigenvalue is (2 N − | i, j i thenormalized state which has two down-spins at site i and j , and (2 N −
2) up-spins at the rest of chain. Then, onecan verify that (cid:12)(cid:12)(cid:12) χ (i) E := 1 √ N N X i =1 ( − i − | i, i + 1 i (mod 2 N ) (49)is an eigenstate of e H (+1 , +1) with eigenvalue (2 N − (cid:12)(cid:12) χ (i) (cid:11) is known as a singular state in the con- text of the Heisenberg chain [47–50]. Therefore, we haveproved that g (i) = 4∆ holds and the first decay mode insector (i) is (1 + ( − N Q X ) √ (cid:12)(cid:12)(cid:12) χ (i) E . (50)The Liouvillian gap in sector (ii) g (ii) and the first de-cay mode can also be obtained in a similar way. One cansee that (cid:12)(cid:12)(cid:12) χ (ii) E := 1 √ N N − X i =1 ( − i − | i, i + 1 i + | N, i (51)is an eigenstate of e H (+1 , −
1) with eigenvalue (2 N − g (ii) ≥
4∆ using Lemma IV.1,and therefore, we have proved rigorously that g (ii) = 4∆ , (52)and the corresponding first decay mode is(1 − ( − N Q X ) √ (cid:12)(cid:12)(cid:12) χ (ii) E . (53) B. In sector (iii) and (iv)
It is known that e H ( − , +1) and e H ( − , −
1) give thesame spectrum [51], which leads to g (iii) = g (iv) . Thus,it suffices to consider only e H ( − , +1). In Ref. [36], theenergy of the ferromagnetic XXZ model under this anti-diagonal twisted boundary condition has been studiedby the ODBA method. As a result, the authors obtainedinhomogeneous Bethe ansatz equations and the formulafor the energy E ase i u j N Y l =1 sin( u j − u l + i η )sin( u j + i η/
2) = e − i u j N Y l =1 sin( u j − u l − i η )sin( u j − i η/
2) + 2ie − Nη sin u j − N X l =1 u l ! ( j = 1 , . . . , N ) , (54) E = −
2i sinh η N X j =1 (cid:20) cot (cid:18) u j + i η (cid:19) − cot (cid:18) u j − i η (cid:19)(cid:21) + 2 N sinh η + 2 sinh η, (55)where { u j } are the Bethe roots and cosh η correspondsto the anisotropy parameter. It is important to notethat the ODBA method is also applicable to complexanisotropy parameter, and hence, complex η . The specialcase cosh η = i∆ is particularly relevant to the analysis ofthe gap, as it corresponds to e H ( − , +1). We confirmednumerically the following two facts: 1. The Bethe roots for the eigenvalue of e H ( − , +1)with the largest imaginary part is obtained by ana-lytic continuation of those for the eigenstate of theHermitian XXZ model with the largest eigenvaluein absolute value (Fig. 4 (b)).2. The string hypothesis (3.1) of Ref. [36] may also bevalid for complex anisotropy (Fig. 4 (a)). In partic- ReIm Im Re (a) (b)
FIG. 4. (Color online) Numerical results of (a) Bethe roots and (b) corresponding eigenvalues of P Ni =1 ( X i X i +1 + Y i Y i +1 +∆e iθ Z i Z i +1 ) under the anti-diagonal twisted boundary condition A N +1 = X A X ( A = X, Y, Z ), where 2 N = 8 , ∆ =cosh(1 . iθ varies from θ = 0 (corresponding to the Hermitian XXZ) to θ = π/ e H ( − , +1)). ular, the corresponding Bethe roots for e H ( − , +1) are obtained by analytic continuation as u j = − π (cid:18) N + 12 − j (cid:19) i (cid:20) ln (cid:16) ∆ + p ∆ + 1 (cid:17) + i π (cid:21) + o ( N ) ( j = 1 , . . . , N ) . (56)Then, by using this string hypothesis for the solution, g (iii) in the thermodynamic limit can be obtained as g (iii) = 2 p ∆ + 1 . (57)As a result, we obtain the explicit formula for the globalLiouvillian gap Eq. (48). We have confirmed that thisanalytical result agrees extremely well with the numericalresult obtained by exact diagonalization for 2 N = 10 (seeFig. 3).One can check the validity of Eq. (57) by consideringthe Ising limit ∆ → ∞ . For notational convenience, werescale the Hamiltonian e H ( − , +1) → e H ( − , +1) = H + V (58)with H = i N X i =1 Z i Z i +1 (59)and V = 2∆ N X i =1 ( S + i S − i +1 + S − i S + i +1 ) , (60) where S ± i = ( X i ± Y i ) / Z N +1 = − Z , S ± N +1 = S ∓ . In the Ising limit, one can first ana-lyze H , and then treat the remaining term V as a per-turbation. It is seen by inspection that (2 N − H with the largest imaginary part. Thecorresponding eigenstates with Q Z = − N states |↑↓ . . . ↓i , |↑↑↑↓ . . . ↓i , . . . , |↑ . . . ↑↓i , |↓↑ . . . ↑i , |↓↓↓↑ . . . ↑i , . . . , |↓ . . . ↓↑i . (61)Each state has two kinks (domain walls), one of whichis between sites 2 N and 1, and the other of which isbetween sites 2 j − j ( j = 1 , . . . , N ). Let P bea projection operator onto the subspace spanned by thestates in Eq. (61). We also define R := X n P| ψ n i =0 | ψ n i h ψ n | (2 N − − E n , (62)where | ψ n i are other eigenstates of the non-perturbedHamiltonian with eigenvalue E n . Then the effectiveHamiltonian from the second order perturbation reads H eff = P ( V + V RV ) P . (63)One can see that the first order perturbation term van-ishes, and the second order one leads to the constant shiftof the eigenvalue by4∆ N − − (2 N − − i∆ (64)without lifting the degeneracy. Putting all this together,we have g (iii) = 2 N ∆ − (cid:20) (2 N − − (cid:21) + O (cid:18) (cid:19) (65)= 2∆ + 1∆ + O (cid:18) (cid:19) , (66)which is consistent with Eq. (57). VI. SUMMARY
We have studied a quantum Ising chain with bulk dis-sipation. By vectorizing the density matrix, we showedthat the Liouvillian of the model can be thought ofas a non-Hermitian Ashkin-Teller model. Then usingthe Kohmoto-den Nijs-Kadanoff transformation, we fur-ther mapped it to a staggered XXZ model with pure-imaginary anisotropy parameters. This mapping has en-abled us to prove that the two NESSs are unique. Fur-thermore, we obtained the exact results for the Liouvil-lian gap on the self-dual line corresponding to the uni-form XXZ model, which shows an excellent agreementwith the numerical results even for small system sizes.Though we mostly focused on the self-dual line, it wouldbe interesting to explore the first decay mode in the wholeparameter region of h, J, ∆ , and ∆ . Whether it also un-dergoes the phase transition remains an interesting ques-tion for future research. ACKNOWLEDGMENTS
The authors thank Nobuyuki Yoshioka and MarkoˇZnidariˇc for fruitful discussions. H.K. was supported inpart by JSPS KAKENHI Grant No. JP18H04478 andNo. JP18K03445. N.S. acknowledges support of the Ma-terials Education program for the future leaders in Re-search, Industry, and Technology (MERIT).
Appendix A: Interacting ParafermionRepresentation
Our model has another representation using Z parafermions [52–55]. Neglecting the boundary terms,we start from Eq. (11). First let us define the followingtwo operators: α j = 1 + i2 σ xj + 1 − i2 ( − j − τ xj , (A1) β j = 12 ( σ zj − τ zj ) + i2 ( − j ( σ xj τ yj + σ yj τ xj ) , (A2) which satisfy α j = β j = 1 , α i β j = i δ ij β j α i . (A3)One can verify that α † i α i +1 + α i α † i +1 = σ xi σ xi +1 − τ xi τ xi +1 , (A4) β i + β † i = σ zi − τ zi . (A5)Then, Eq. (11) is rewritten in terms of α i and β i as H bulk = − h N X i =1 ( β i + β † i ) − J N − X i =1 ( α † i α i +1 + α i α † i +1 ) − i∆ N X i =1 β i + β † i − i∆ N − X i =1 ( α † i α i +1 ) + ( α i α † i +1 ) . (A6)Next, we define Z parafermion operators as follows: γ i − = i − Y j =1 β j α i , (A7) γ i = e i π/ i − Y j =1 β j α i β i = e i3 π/ i Y j =1 β j α i . (A8)They satisfy γ i = 1 , γ † i = γ i , γ i γ j = i γ j γ i ( i < j ) . (A9)By noting that β i = e − i π/ γ † i − γ i , (A10) α † i α i +1 = e i3 π/ γ † i γ i +1 , (A11)we obtain H bulk = − h N X i =1 (cid:16) e − i π/ γ † i − γ i + e i π/ γ † i γ i − (cid:17) − J N − X i =1 (cid:16) e i3 π/ γ † i γ i +1 + e − i3 π/ γ † i +1 γ i (cid:17) − i∆ N X i =1 γ i − γ i − i∆ N − X i =1 γ i γ i +1 . (A12)The first and second terms are quadratic, while the thirdand fourth terms are quartic in parafermions. Therefore,our model can be thought of as a Z parafermion chainwith non-Hermitian quartic interactions (up to boundaryterms). [1] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. B¨uchler,and P. Zoller, Nat. Phys. , 878 (2008).[2] T. Prosen and I. Piˇzorn, Phys. Rev. Lett. , 105701(2008).[3] S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller,Phys. Rev. Lett. , 015702 (2010).[4] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Nat.Phys. , 971 (2011).[5] C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico,A. ˙Imamo˘glu, P. Zoller, and S. Diehl, New J. Phys. ,085001 (2013).[6] B. Kraus, H. P. B¨uchler, S. Diehl, A. Kantian, A. Micheli,and P. Zoller, Phys. Rev. A , 042307 (2008).[7] M. J. Kastoryano, F. Reiter, and A. S. Sørensen, Phys.Rev. Lett. , 090502 (2011).[8] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nat.Phys. , 633 (2009).[9] T. Prosen, E. Ilievski, and V. Popkov, New J. Phys. ,073051 (2013).[10] H.-P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems (Oxford University Press, Oxford,UK, 2002).[11] A. C. Y. Li, F. Petruccione, and J. Koch, Sci. Rep. ,4887 (2014).[12] M. ˇZnidariˇc, Phys. Rev. E , 042143 (2015).[13] T. Prosen and M. ˇZnidariˇc, J. Stat. Mech. , P02035(2009).[14] M. ˇZnidariˇc, Phys. Rev. Lett. , 220601 (2011).[15] M. J. Hartmann and G. Carleo, (2019),arXiv:1902.05131.[16] A. Nagy and V. Savona, (2019), arXiv:1902.09483.[17] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti,(2019), arXiv:1902.10104.[18] N. Yoshioka and R. Hamazaki, (2019), arXiv:1902.07006.[19] T. Prosen, Phys. Rev. Lett. , 137201 (2011).[20] D. Karevski, V. Popkov, and G. M. Sch¨utz, Phys. Rev.Lett. , 047201 (2013).[21] T. Prosen, Phys. Rev. Lett. , 030603 (2014).[22] M. V. Medvedyeva, F. H. L. Essler, and T. Prosen, Phys.Rev. Lett. , 137202 (2016).[23] N. Shibata and H. Katsura, Phys. Rev. B , 174303(2019).[24] S. Sachdev, Quantum Phase transitions (Cambridge Uni-versity Press, Cambridge, UK, 1999).[25] A. Dutta, G. Aeppli, B. K. Chakrabarti, U. Divakaran,T. F. Rosenbaum, and D. Sen,
Quantum phase tran-sitions in transverse field spin models: from statisticalphysics to quantum information (Cambridge UniversityPress, Cambridge, UK, 2015).[26] C. Ates, B. Olmos, J. P. Garrahan, and I. Lesanovsky,Phys. Rev. A , 043620 (2012).[27] A. De Luca, J. Viti, D. Bernard, and B. Doyon, Phys.Rev. B , 134301 (2013).[28] B. Horstmann, J. I. Cirac, and G. Giedke, Phys. Rev. A , 012108 (2013).[29] L. M. Vasiloiu, F. Carollo, and J. P. Garrahan, Phys.Rev. B , 094308 (2018).[30] T. Prosen, New J. Phys. , 043026 (2008). [31] M. Kohmoto, M. den Nijs, and L. P. Kadanoff, Phys.Rev. B , 5229 (1981).[32] F. C. Alcaraz, M. N. Barber, and M. T. Batchelor, An-nals of Phys. , 280 (1988).[33] M. Yamanaka, Y. Hatsugai, and M. Kohmoto, Phys.Rev. B , 9555 (1993).[34] F. C. Alcaraz and J. D. de Fel´ıcio, J. Phys. A: Math.Gen. , L651 (1984).[35] M. Yamanaka, Y. Hatsugai, and M. Kohmoto, Phys.Rev. B , 559 (1994).[36] Y. Qiao, Z. Xin, K. Hao, J. Cao, W.-L. Yang, K. Shi,and Y. Wang, New J. Phys. , 073046 (2018).[37] M. ˇZnidariˇc, Phys. Rev. E , 042140 (2014).[38] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Phys.Rev. A , 042118 (2018).[39] In the case where N is odd, our model and the quan-tum Ashkin-Teller model has difference in the boundaryterms.[40] The same model can be constructed from the one-dimensional quantum compass model H = − N/ X i =1 J x σ x i − σ x i − N/ − X i =1 J y σ y i σ y i +1 (A13)with the Lindblad operators L i − = √ ∆ σ x i − σ x i , L i = √ ∆ σ y i σ y i +1 , (A14)although there exist many conserved charges in this case,and hence degenerate NESSs.[41] M. T. Batchelor, R. J. Baxter, M. J. O’Rourke, andC. M. Yung, J. Phys. A: Math. Gen. , 2759 (1995).[42] S. Niekamp, T. Wirth, and H. Frahm, J. Phys. A: Math.Theor. , 195008 (2009).[43] J. Cao, W.-L. Yang, K. Shi, and Y. Wang, Phys. Rev.Lett. , 137201 (2013).[44] Y. Wang, W.-L. Yang, J. Cao, and K. Shi, Off-diagonalBethe ansatz for exactly solvable models (Springer,Berlin-Heidelberg, DE, 2015).[45] ´Angel Rivas and S. F. Huelga,
Open Quantum Systems:An Introduction (Springer, Heidelberg, DE, 2011).[46] P. Di Francesco, P. Mathieu, and D. S´en´echal,
Conformalfield theory (Springer-Verlag, NewYork, US, 1997).[47] L. V. Avdeev and A. A. Vladimirov, Theor. Math. Phys. , 1071 (1986).[48] F. H. L. Essler, V. E. Korepin, and K. Schoutens, J.Phys. A: Math. Gen. , 4115 (1992).[49] J. D. Noh, D.-S. Lee, and D. Kim, Physica A , 167(2000).[50] R. I. Nepomechie and C. Wang, J. Phys. A: Math. Theor. , 325002 (2013).[51] F. C. Alcaraz, M. Baake, U. Grimmn, and V. Rittenberg,J. Phys. A: Math. Gen. , L117 (1988).[52] P. Fendley, J. Stat. Mech. , P11020 (2012).[53] N. Moran, D. Pellegrino, J. K. Slingerland, and G. Kells,Phys. Rev. B , 235127 (2017).[54] A. Calzona, T. Meng, M. Sassetti, and T. L. Schmidt,Phys. Rev. B , 201110(R) (2018).[55] A. Chew, D. F. Mross, and J. Alicea, Phys. Rev. B98