Integrability of one-dimensional Lindbladians from operator-space fragmentation

We introduce families of one-dimensional Lindblad equations describing open many-particle quantum systems that are exactly solvable in the following sense: (i) The space of operators splits into exponentially many (in system size) subspaces that are left invariant under the dissipative evolution; (ii) the time evolution of the density matrix on each invariant subspace is described by an integrable Hamiltonian. The prototypical example is the quantum version of the asymmetric simple exclusion process (ASEP) which we analyze in some detail. We show that in each invariant subspace the dynamics is described in terms of an integrable spin-1 / 2 XXZ Heisenberg chain with either open or twisted boundary conditions. We further demonstrate that Lindbladians featuring integrable operator-space fragmentation can be found in spin chains with arbitrary local physical dimensions.


I. INTRODUCTION
In many physical situations, a quantum system weakly coupled to an environment can be described by a Lindblad master equation [1,2]. This description is valid within a Markovian approximation, where the internal bath dynamics is assumed to be much faster than that of the system itself [3]. The study of many-body Lindblad equations is generally quite difficult and even in one spatial dimension the available tools are largely restricted to perturbative [4,5] and sophisticated numerical approaches [6][7][8][9][10]. A natural question is then whether one can find exactly solvable Lindblad equations that can give detailed, nonperturbative insights into representative cases, and provide reliable benchmarks for the development of approximate methods. This question is also highly relevant in light of the important role that exactly solvable models have played in the recent efforts aimed at understand nonequilibrium dynamics in isolated quantum systems [11][12][13][14][15][16].
It has been known for some time that certain Lindblad equations can be cast in the form of imaginary-time Schrödinger equations with non-Hermitian "Hamiltonians" that are quadratic in fermionic or bosonic field operators [17]. This allows one to obtain exact results on the dynamics of correlation functions [18][19][20][21][22][23] as well as entanglement-related quantities [24,25]. Furthermore, it has also been realized that in some models, not necessarily integrable, local correlation functions satisfy closed hierarchies of equations of motion, provided that the Lindbladian satisfies certain conditions [ [26][27][28][29][30][31][32] (although in these cases the full spectrum typically remains out of reach). Very recently, further progress has been made in this direction, with the discovery of Lindbladians that can be mapped onto known Yang-Baxter integrable interacting models [33][34][35][36][37][38][39][40][41][42]. In particular, the approach put forward in Refs. [33,37] is based on a map between Lindblad equations and solvable "two-leg ladder" quantum spin chains.
Here, we highlight a different mechanism that allows us to relate certain Lindbladians to known integrable models. The systems that we study are characterized by the fact that the space of operators "fragments" into exponentially (in system size) many subspaces that are left invariant under the dissipative evolution. This is reminiscent of the fragmentation induced by dipole-conserving Hamiltonians [43,44] or confinement [45], although in our case this mechanism pertains to the space of operators, and not that of states. For the Lindblad equations studied in this paper, the projection of the Lindbladian onto each invariant subspace is mapped onto an integrable Hamiltonian, thus allowing us to obtain the full spectrum analytically. Our construction is carried out for purely dissipative dynamics, although it is possible to add certain coherent terms without spoiling integrability.
We illustrate the above mechanism in a prototypical example: the quantum version of the celebrated asymmetric simple exclusion process (ASEP) [46][47][48][49]. The classical ASEP is known to be integrable [50,51] and this has made it possible to obtain a host of exact results (see, e.g., Refs. [52][53][54][55][56][57][58][59][60][61][62][63][64][65]). This model can be obtained as the averaged dynamics of the so-called quantum ASEP [66], a stochastic quantum model of particles hopping with random amplitudes first introduced in its symmetric form in Ref. [67] (see also Refs. [68][69][70][71]). In this case, we show that the Lindbladian restricted to each invariant subspace can be mapped onto a X X Z Heisenberg Hamiltonian with integrable boundary conditions. We also show that a similar "integrable fragmentation" can be found in spin chains with arbitrary local physical dimensions.
The rest of this paper is organized as follows. We begin in Sec. II, where we introduce the general idea behind the integrable operator-space fragmentation. In Sec. III we work out in detail the interesting case of the ASEP Lindbladian, discussing both its exact spectrum (Sec. III A), and the dissipative dynamics of local observables (Sec. III B). In Sec. IV we present higher-spin generalizations, while we report our conclusions in Sec. V. Finally, we consign the most technical aspects of our paper to a few Appendices.

II. OPERATOR-SPACE FRAGMENTATION
We consider one-dimensional qudit systems defined on the Hilbert space H = H 1 ⊗ · · · ⊗ H L , with H j C d . We denote a basis of the local space H j by |α d α=1 , and set |α = |α 1 1 · · · |α L L . We study the Lindbladian dynamics of the density matrix ρ(t ) of the system defined by with J a > 0, K ∈ N, and where periodic boundary conditions are assumed. In the following, we focus on jump operators L (a) j acting on up to two neighboring sites. We employ a superoperator formalism [72], according to which the density matrix can be viewed as a state in a d 2L -dimensional Hilbert space H ⊗ H, with basis elements |α |β = |α 1 1 · · · |α L L |β 1 1 · · · |β L L . Given any operator O = α,β O α,β |α β| we have the natural mapping O → |O = α,β O α,β |α |β . In this formalism, the Lindblad equation can be cast in the form d|ρ /dt = L|ρ where we introduced the Lindbladian L = −iH + iH + L j=1 D j , and Here, we have employed the notation O = O ⊗ 1, while we definedŌ = 1 ⊗ O T , with α|O T |β = β|O|α . In the following we assume H ≡ 0 (later, we will comment on the possibility of adding simple Hamiltonian terms that do not spoil integrability). For each site j, we take r < d 2 rank-1 projectors P α j , with α = 1, . . . , r, satisfying P α j P β j = δ α,β P α j , and define P 0 can be decomposed as a direct sum of simultaneous eigenspaces K α of the projection operators P α j . These subspaces are labeled by vectors α = {α 1 , . . . , α L }, specifying that K α is associated with the projector P α 1 1 · · · P α L L . We note that by construction K α are left invariant under the action of the Lindbladian, and their number grows exponentially with the size L.
Consider a subspace such that α j = 0 for j ∈ { 1 , . . . q }, and α j = 0 otherwise. Physically, the sites labeled by j can be thought of as "defects," separating different regions in the operator-space fragmentation. Indeed, the Lindbladian restricted to K α can be expressed as with 0 = q and q+1 = 1 . Here, we have defined and , and D j is given in Eq. (2). Our observation is that there exist simple families of jump operators L (a) j such that the action of L α,β [m,n] coincides with that of an integrable Hamiltonian with appropriate boundary terms. When this happens, the full spectrum of the Lindbladian can be obtained analytically, by analyzing each invariant sector individually. Importantly, we stress that L α,β [m,n] acts on a space where the effective local dimension is not necessarily the square of an integer number, allowing for a mapping to different integrable systems with respect to those found in Refs. [33,37].

III. THE ASEP LINDBLADIAN
As a concrete example, we consider the fully dissipative Lindbladian defined by Eq. (2), with d = K = 2, and Here, σ ± = (σ x ± iσ y )/2, and σ α are the Pauli matrices. This Lindbladian admits an equivalent representation in terms of spinless fermionic operators via a standard Jordan-Wigner mapping [73], and in this formulation it can be obtained after averaging over the bath degrees of freedom in the "quantum ASEP" studied in Refs. [66][67][68][69][70]. It was already pointed out in these works that the above Lindblad equation reduces to the classical ASEP in the subspace of density matrices that are diagonal in the position basis introduced above [74]. When combined with the known mapping between the classical ASEP and the quantum X X Z Heisenberg chain [50,51], this gives us the integrability of the Lindbladian restricted to this subspace. In the following, we establish a stronger statement, namely, we show that its full spectrum can be obtained analytically, as a result of the operator-space fragmentation.
Following the general discussion given before, we focus on the Lindbladian restricted to different subspaces K α . Let us first consider the case where α j = 0, for all j. Its restriction to K α acts on the tensor product of two-dimensional local spaces, which are spanned by the states |e 1 j , |e 4 j . Choosing this basis for each site, we show in Appendix A that the restricted Lindbladian can be rewritten as 062210-2 This is precisely the evolution operator of the classical ASEP with periodic boundary conditions, which is known to be integrable [50,51]. Applying a local similarity transformation S = j S j , with S = diag(1, √ J 1 /J 2 ), L ASEP can be mapped to a spin-1/2 Heisenberg Hamiltonian with twisted boundary conditions (cf. Appendix A). The Lindbladian (6) describes the time evolution of quantum states that are by construction unentangled.
When α j = 0 for some j ∈ { 1 , . . . , q }, instead, the Lindbladian can be decomposed as in Eq. (3). If | j − k | 3, the projected Lindbladian is simply a constant. Otherwise, we find where and n] are required in order to describe time evolution of quantum states with nonvanishing entanglement. The operators (6) and (8) are integrable, corresponding to the X X Z Heisenberg chain with diagonal twisted [75,76] and open [77] boundary conditions, respectively.
Before leaving this section, we sketch how a coherent term can be added to the Lindbladian without spoiling integrability. It is easy to see that, for the quantum ASEP, any Hamiltonian whose action preserves the invariant subspaces (namely, that is written entirely in terms of Pauli matrices σ z j ) also preserves integrability. The most general nearest-neighbor Hamiltonian of this form is In Appendix A, we show that the projection of the corresponding term i(H −H ) onto the different invariant subspaces results in diagonal boundary terms, thus preserving integrability.

A. The exact spectrum
We now study the full spectrum of the Lindbladian defined by the jump operators (5). Let us first consider the case where α j = 0 for all j, corresponding to Eq. (6). It is easy to see that the ground-state energy is E GS = 0, with degeneracy L + 1. The associated eigenspace is spanned by |gs M , M = 0, 1, . . . , L, where |gs M is the equal-weight superposition of all spin configurations in the sector of M "down spins." We note that |gs M corresponds to the projection of the nonequilibrium steady state ρ NESS = 1/2 L onto a given magnetization sector. It is worth mentioning that, although ρ NESS is formally an infinite-temperature state, it features a nonvanishing expectation value of the spin current, whose density J z j is defined by tr[ρ(t )σ z j ] = −tr[ρ(t )(J z j+1 − J z j )]. One can study the full spectrum of L ASEP (and the associated eigenstates) via a standard application of the Bethe ansatz formalism. In fact, a detailed analysis has been already reported in Ref. [51], in the context of the classical ASEP, where it was shown that the Hamiltonian features a spectral gap vanishing as L −3/2 . In the following, we sketch how an exact analysis can be carried out also for the spectrum of all other projected Lindbladians.
Let us then consider the case where α j = 0 for some j. We immediately see that, since H X X Z [m,n] 0, the eigenvalues λ k of L [m,n] in Eq. (7) are bounded by where the bound is strict, since the smallest eigenvalue of H X X Z [m,n] is E 0 GS = 0 (corresponding to the two saturated ferromagnetic states). This means that the long-time limit of the Lindbladian evolution is dominated by the subspace of locally diagonal density matrices, containing all the "soft modes" of the dynamics. Once again, the full spectrum of H X X Z [m,n] can be obtained using the Bethe ansatz formalism [77]. According to the latter, each eigenstate is associated with a set of quasimomenta {λ j } M j=1 , satisfying the set of quantization conditions where = n − m + 1 is the length of the open chain, and with η = arccosh( ). The associated energy eigenvalue is then One can explicitly identify the eigenstate associated with the lowest-energy eigenvalue E M GS in each magnetization sector, by following the derivation in Ref. [78], where the analogous problem for periodic boundary conditions was considered. We report this analysis for M = /2 "down spins" (assuming to be even). By numerically solving the Bethe ansatz equations for small systems we find that the smallest eigenvalue in this sector corresponds to a maximal string with real part π/2, namely, to a solution of the Bethe equations of the form λ j = π 2 + i(2 j − 1) γ 2 + δ j , j = 1, . . . , M, where δ j are deviations from a "perfect string." The ground state is doubly degenerate in the thermodynamic limit, and at finite size the exponentially close ground states correspond to different values of the deviations δ j . Assuming δ j to vanish as → ∞, and repeating the steps reported in Ref. [78], one finds where a > 0 is a constant (which depends on ). Note that, as in the case of periodic boundary conditions, low-lying excitations are separated by a finite gap from E M GS . A careful analysis of the deviations δ j for large systems, and the identification of solutions corresponding to higher "excited states" are of interest in particular in light of the complications that are known to arise for = 1 [79,80], but are beyond the scope of this work.
Due to Eq. (10), we see that the invariant subspaces are organized in hierarchies of decreasing energy. Indeed, given α, we can identify the subset {α j } such that α j = 0, and as long as | j − j+1 | > 1, the Lindbladian restricted to the corresponding invariant subspace lowers the ground-state energy by an amount −(J 1 + J 2 )/2. If | j − j+1 | = 1 for all j, then we have a completely fragmented space (α j = 0 for all j), and a dedicated analysis is needed. As we show in Appendix A, in this case, any product state is fixed under the Lindbladian evolution, and thus corresponds to an eigenstate of the Lindbladian. The associated eigenvalue is always bounded by −(J 1 + J 2 )/2, except for the two special states |e 2 = |e 2 ⊗L and |e 3 = |e 3 ⊗L that are annihilated by L, making its ground-state degeneracy L + 3.

B. The dissipative dynamics
Any initial density matrix can be decomposed as |ρ(0) = α ρ αα |α |α + β =α ρ αβ |α |β . Due to Eq. (10), we have that, up to terms that are exponentially small in time, and the late-time behavior is therefore fully determined by the classical ASEP. Here, we have neglected the subspace associated with the fixed points |e 2 , |e 3 , which is expected to yield contributions that are exponentially small in the system size L. On the other hand, for "quantum" operators that have no analogs in the classical ASEP the dynamics takes place entirely in subspaces that are orthogonal to that spanned by |α |α . In order to illustrate this, we consider the transverse spin-spin correlation function In the superoperator formalism, this can be rewritten as ⊗L . This shows that this two-point function is only sensitive to the invariant subspace corresponding to α 1 = 1, α = 2, and α j = 0 otherwise. Assuming for simplicity 3, this leads to where we have defined φ [m,n] | = ⊗ n j=m ( e 1 | j + e 4 | j ), |a = |e a 1 |e a 2 · · · |e a −2 , |b = |e b 1 · · · |e b L− , for a j , b j = 1, 4 and ρ 0 (a, b) = a | b |ρ(0) , where |a = |e 2 |a , |b = |e 3 |b . Equation (16) reduces the problem of calculating the transverse spin-spin correlator to computing φ [m,n] |eL [m,n] t |a . This is still nontrivial because φ [m,n] | is not an eigenstate ofL [m,n] , due to the boundary terms. A simplification occurs in the isotropic limit J 1 = J 2 , where φ [m,n] | becomes a left eigenstate ofL [m,n] with eigenvalue −J. Hence we have We note that this result could be obtained in a direct way from the equations of motion for σ + 1 σ − , which for J 1 = J 2 become linear. For J 1 = J 2 the calculation of φ [m,n] |eL [m,n] t |a is an open question, but we note that similar quantities have been recently determined in problems of real-and imaginary-time evolution from simple product states [81][82][83].

IV. HIGHER-SPIN GENERALIZATIONS
Lindbladians featuring similar types of fragmentation can be constructed also for higher local dimensions. We consider a Hilbert space H = H 1 ⊗ · · · ⊗ H L , with H j C N and a purely dissipative Lindblad equation of the formρ = − j D j (ρ), with D j as in Eq. (2). We choose K = N 2 , J a = 1, and jump operators [84] L (α,β ) Here, E α,β is the operator with matrix elements (E α,β ) α ,β = δ α,α δ β,β . In the superoperator formalism, it is readily seen projects onto the subspace spanned by the "diagonal" states |e α j = |α j ⊗ |α j , α = 1, . . . , N. This implies that the Hilbert space H ⊗ H splits into exponentially many invariant subspaces K α . We start by analyzing the one corresponding to α j = 0 for all j. In this case the local physical dimension is N and choosing the local diagonal basis introduced above, the projection of L can be rewritten as where j, j+1 is the operator swapping neighboring sites (cf. Appendix B). We see that L | α=0 is indeed integrable, as its action coincides with that of the SU(N )-invariant Sutherland model, first solved in Ref. [85]. Next, we consider the generic case where α j = 0 for some j ∈ { 1 , . . . , q }. If | j − k | 3, the projection of L vanishes; otherwise, we have the factorization (7) This is the SU(N )-invariant chain with diagonal boundary conditions, which is once again integrable [86][87][88] and, as a result, the full spectrum can be analyzed analytically via the Bethe ansatz. We mention that, although here we have focused on J a = 1 for all a, we expect that one could also choose J a such that each projected Lindbladian is mapped onto an integrable deformation of the isotropic Sutherland model.

V. CONCLUSIONS
We have studied families of Lindbladians that can be mapped onto known Yang-Baxter integrable models. These systems are characterized by a fragmentation of the operator space into exponentially many invariant subspaces, where the projected Lindbladian acts as an integrable "Hamiltonian". We have worked out in detail the case of the quantum ASEP Lindbladian [66], and exhibited explicit further examples with arbitrary local dimensions. Our study raises several questions. Most prominently, one could wonder what are the consequences of the operator-space fragmentation and integrability beyond correlation functions. For instance, it would be particularly interesting to investigate the dynamics of entanglement-related quantities [25]. Here, the block-diagonal structure of the evolved density matrix gives us a promising starting point for the study of the entanglement negativity [89]. We note that, while using Eq. (10) the latter can be seen to be proportional to e −(J 1 +J 2 )t/2 at late times, its short-time dynamics is expected to be nontrivial. We plan to come back to these questions in future research [84].

ACKNOWLEDGMENTS
We thank Denis Bernard and Aleksandra Ziolkowska for helpful discussions. This work was supported in part by the EPSRC under Grant No. EP/S020527/1 (F.H.L.E.) and the Alexander von Humboldt foundation (L.P.).

APPENDIX A: THE ASEP LINDBLADIAN
We consider the Lindbladian defined by Eq. (2). As in the main text, we define the operators E α,β by their action E α,β |e γ = δ β,γ |e α . Explicitly, we have With these notations, the Lindbladian superoperator reads We see that L commutes with P 1 j = E 2,2 j , P 2 j = E 3,3 j , so that P 0 j = 1 j − (P 1 j + P 2 j ) is the projector onto the subspace associated with density matrices that are diagonal at site j. We now proceed to analyze the Lindbladian restricted to the different subspaces K α , defined in the main text.
Let us first consider the case where α j = 0, for all j, corresponding to the subspace of locally diagonal density matrices. The Lindbladian restricted to K α acts on the tensor product of two-dimensional local spaces, which are spanned by the states |e 1 j , |e 4 j . Choosing this basis for all sites, we can rewrite the projected Lindbladian as In the general case, instead, α j = 0 for j ∈ { 1 , . . . , q }, and the Lindbladian can be decomposed as in Eq. (3). Suppose first that | k − k+1 | 3. In this case, it is straightforward to write Note thatL [m,n] does not depend on α and β. As stated in the main text, we can map these Lindbladians onto X X Z Heisenberg Hamiltonians, by considering the similarity transformation where Q = √ J 1 /J 2 . Applying this to both Eqs. (A5) and (A7), we finally obtain and where and We are now left to consider the two special cases L α,β [m,m+1] and L α,β [m,m+2] . A simple calculation shows that in both cases the restricted Lindbladian is proportional to the identity. Specifically, we find and Finally, let us discuss the coherent term generated by the Hamiltonian (9). It is straightforward to rewrite From this expression, it follows immediately that the projection of this operator onto each invariant subspace results in diagonal boundary terms, thus preserving integrability.
This implies that the Hilbert space H ⊗ H splits into exponentially many invariant subspaces K α . We start by analyzing the subspace corresponding to α j = 0 for all j. First, it is straightforward to show that the local terms in the Lindbladian L act as SWAP operators, when applied to the states (B6). Explicitly, we have N α,β=1 where j, j+1 is the operator permuting sites j and j + 1. Thus, we simply obtain which is Eq. (19) in the main text. Next, suppose α j = 0 for j ∈ { 1 , . . . , q }. Then, the Lindbladian L can be decomposed as in Eq. (3). First, we consider | k − k+1 | 3. In this case, we see that the factorization condition (7)  Choosing the diagonal basis (B6), and making use of Eq. (B7), we obtain Eq. (20). Finally, when | j − j+1 | 3, it is easy to show that the projection of L is vanishing.