Integrability of $1D$ 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 dimension.

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, non-perturbative 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 non-equilibrium 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 interacting integrable 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], although in our case this mechanism pertains the space of operators, and not that of states. For the Lindblad equations studied in this work, the projection of the Liouvillian 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 Liouvillians, 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 simple asymmetric exclusion process (ASEP) [45][46][47][48]. The classical ASEP is known to be integrable [49,50] and this has made it possible to obtain a host of exact results, see e.g. [51][52][53][54][55][56][57][58][59][60][61][62][63][64]. This model can be obtained as the averaged dynamics of the so-called quantum ASEP [65], a stochastic quantum model of particles hopping with random amplitudes first introduced in its symmetric form in Ref. [66], see also Refs. [67][68][69][70]. In this case, we show that the Liouvillian restricted to each invariant subspace can be mapped onto a XXZ Heisenberg Hamiltonian with integrable boundary conditions. We also show that a similar "integrable fragmentation" can be found in spin chains with arbitrary local physical dimension.
We consider onedimensional qudit systems defined on the Hilbert space 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 [71], 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 Liouvillian L = −iH + iH + L j=1 D j , and Here we have employed the notation In the following we assume H ≡ 0 (see however [72] for a discussion of simple Hamiltonian terms that can be added without spoiling 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 j = 1 1 j − r α=1 P α j . Now, suppose [L, P α j ] = 0 for all α, j. Since [P α j , P β k ] = 0, the space H ⊗ H 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 Liouvillian, and their number grows exponentially with the size L.
Consider a subspace such that α j = 0 for j ∈ { 1 , . . . q }, and α j = 0 otherwise. Then, the Liouvillian restricted to K α can be expressed as with 0 = q and q+1 = 1 . Here, we have defined and L α,β , 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 Liouvillian 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].
The ASEP Lindbladian. As a concrete example, we consider the fully dissipative Liouvillian defined by Eq. (2), with d = K = 2, and Here σ ± = (σ x ± σ 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. [65][66][67][68][69]. 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 XXZ Heisenberg chain [49,50] this gives us the integrability of the Liouvillian restricted to this subspace. In the following, we establish a stronger statement, namely we show that the full spectrum of the Liouvillian can be obtained analytically, as a result of the operator-space fragmentation.
It is convenient to order the basis elements generating the local Hilbert space H j ⊗ H j as |e 1 = |1 ⊗ |1 , |e 2 = |2 ⊗ |1 , |e 3 = |1 ⊗ |2 , |e 4 = |2 ⊗ |2 and define, accordingly, the operators E α,β by their action E α,β |e γ = δ β,γ |e α . With these notations, it is straightforward to see that L commutes with P 1 j = E 2,2 j and P 2 j = E 3,3 j , so that P 0 j is the projector onto the subspace spanned by |α j |α j Following the general discussion given before, we focus on the Liouvillian restricted to different subspaces K α . Let us first consider the case where α j = 0, for all j. The Liouvillian restricted to K α acts on the tensor product of 2-dimensional local spaces, which are spanned by the states |e 1 j , |e 4 j . Choosing this basis for each site the restricted Liouvillian can be rewritten as [72] This is precisely the evolution operator of the classical ASEP with periodic boundary conditions, which is known to be integrable [49,50]. 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 [72]. The Liouvillian (6) describes time evolution of quantum states that are by construction unentangled.
The exact spectrum. 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 non-equilibrium steady state ρ NESS = 1 1/2 L onto a given magnetization sector. It is worth mentioning that, although ρ NESS is formally an infinite-temperature state, it features a non-vanishing 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. [50], 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 Liouvillians.
Let us then consider the case where α j = 0 for some j. We immediately see that, since H XXZ [m,n] ≥ 0, the eigenvalues λ k ofL [m,n] in Eq. (7) are bounded by where the bound is strict, since the smallest eigenvalue of H XXZ [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 XXZ [m,n] can be obtained using the Bethe Ansatz formalism [77]. According to the latter, each eigenstate is associated with a set of quasi-momenta {λ 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(2j − 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. (9), 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 Liouvillian 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 [72]. In this case, any product state is fixed under the Lindbladian evolution, and thus corresponds to an eigenstate of the Liouvillian. 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.
Dissipative dynamics. Any initial density matrix can be decomposed as |ρ(0) = α ρ αα |α |α + β =α ρ αβ |α |β . Due to Eq. (9), we have that, up to terms that are exponentially small in time, |ρ(t) = α ρ αα e LASEPt |α |α + O(e −(J1+J2)t/2 ) , (13) 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 analogues 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 a1 |e a2 · · · |e a −2 , |b = |e b1 · · · |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 . Eq. (15) reduces the problem of calculating the transverse spin-spin correlator to computing φ [m,n] |eL [m,n] t |a . This is still non-trivial 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].
Here, E α,β is the operator with matrix elements (E α,β ) α ,β = δ α,α δ β,β . In the superoperator formalism, it is readily seen [72] that L = L − L1 1, where L commutes with the d(d − 1) rank-1 projectors P (α,β) j = E α,α j ⊗ E β,β j for α = β, α, β = 1, . . . , d, so that P 0 j = 1 1 j − α =β P (α,β) j 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 [72] where Π j,j+1 is the operator swapping neighboring sites. 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 [72]; otherwise, we have the factorization (7) with 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 Liouvillian is mapped onto an integrable deformation of the isotropic Sutherland model.

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 Liouvillian acts as an integrable Hamiltonian. We have worked out in detail the case of the quantum ASEP Lindbladian [65], and exhibited explicit further examples in arbitrary local dimension. 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. (9) the latter can be seen to be proportional to e −(J1+J2)t/2 at late times, its short-time dynamics is expected to be non-trivial. We plan to come back to these questions in future research.

SUPPLEMENTARY MATERIAL
Here we present all the calculations leading to the results reported in the main text.

The ASEP Lindbladian
We consider the Liouvillian 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 Liouvillian 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 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 Liouvillian 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 Liouvillian restricted to K α acts on the tensor product of 2-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 Liouvillian as In the general case, instead, α j = 0 for j ∈ { 1 , . . . , q }, and the Liouvillian 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 Liouvillians onto XXZ Heisenberg Hamiltonians, by considering the similarity transformation where Q = J 1 /J 2 . Applying this to both Eqs. (24) and (26), 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 Liouvillian is proportional to the identity. Specifically, we find and Adding a coherent Hamiltonian term As we have mentioned in the main text, we can add a coherent term 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 yielding the following additional term to the Liouvillian In turn, it is straightforward to see 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 Liouvillian L act as swap operators, when applied to the states (41). Explicitly, we have N α,β=1 E α,β j ⊗ E α,β j E β,α j+1 ⊗ E β,α j+1 |e γ j |e δ j+1 = |e δ j ⊗ |e γ j = Π j,j+1 |e γ j |e δ j+1 , where Π j,j+1 is the operator permuting sites j and j + 1. Thus, we simply obtain which is Eq. (18) in the main text. Next, suppose α j = 0 for j ∈ { 1 , . . . , q }. Then, the Liouvillian 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) holds, withL Choosing the diagonal basis (41), and making use of Eq. (42), we obtain Eq. (19). Finally, when | j − j+1 | ≤ 3, it is easy to show that the projection of L is vanishing.