Yang-Baxter integrable Lindblad equations

We consider Lindblad equations for one dimensional fermionic models and quantum spin chains. By employing a (graded) super-operator formalism we identify a number of Lindblad equations than can be mapped onto non-Hermitian interacting Yang-Baxter integrable models. Employing Bethe Ansatz techniques we show that the late-time dynamics of some of these models is diffusive.


Introduction
Weak couplings to an environment can have very interesting effects on the dynamics of manyparticle quantum systems. In particular they can result in desirable non-equilibrium steady states [1][2][3][4][5]. In order to arrive at a tractable theoretical description it is customary to employ a Markovian approximation that assumes that the characteristic times scales associated with the environment are much shorter than those of the many-particle system of interest. The absence of a back action of the system onto its environment then facilitates a well defined mathematical description of open many-particle systems. In the quantum case this a priori results in a Markovian quantum stochastic many-particle system [6][7][8][9], which is however difficult to analyze. The customary approach is therefore to focus on the dynamics averaged over the environment, which leads to a description by the Lindblad master equation [10] for the time-dependent reduced density matrix ρ(t) Here H is the system Hamiltonian, L a are jump operators that encode the coupling to the environment and γ a > 0. While much progress has been made in analyzing Lindblad equations for many-particle systems by employing e.g. perturbative [11,12] and matrix product states methods [13][14][15][16] it clearly is highly desirable to have exact solutions in specific, and hopefully representative, cases. In the context of master equations for classical stochastic many particle systems an example of such a solvable paradigm is the asymmetric simple exclusion process [17][18][19][20][21][22].
In the quantum case it has been known for some time that certain Lindblad equations describing many-particle systems can be represented by Liouvillians that are quadratic in fermionic or bosonic creation and annihilation operators, which makes it possible to solve them exactly by elementary means [23][24][25][26]. Very recently examples of Lindblad equations with Liouvillians related to interacting Yang-Baxter integrable models have been found [27][28][29]. This opens the door for bringing quantum integrability methods to bear on obtaining exact results for the dynamics of open many-particle quantum systems. An obvious question is whether the known cases are exceptional, or whether there are other examples of Yang-Baxter integrable Lindblad equations. In this work we report on the results of a search for integrable cases among a particular class of Lindblad equations for translationally invariant many-particle quantum systems.

Lindblad equations for lattice models
We now turn to the precise definition of the class of quantum master equations we will be interested in. We consider one dimensional lattice models with local Hilbert spaces that can include bosonic as well as fermionic degrees of freedom. A basis of the local Hilbert space is formed by N bosonic and M fermionic quantum states |α j , α = 1, . . . , N + M.
We denote the fermion parity of the state |α j by α An orthonormal basis of the full Hilbert space H L on an L-site chain is then given by the states We define the fermion parity of the states (4) by A basis of the space of linear operators acting on site j is then provided by These are often referred to as Hubbard operators. Their fermion parity is α + β mod 2, i.e. they are fermionic if either the state |α or the state |β is fermionic. The operators E αβ n act on the states |α as Minus signs are acquired when moving fermionic operators past fermionic states. The operators defined in this way either commute or anticommute on different sites For later convenience we define a graded permutation operator on sites j and j + 1 It acts on states as i.e. it permutes the states and generates a minus sign if both states are fermionic.

A useful decomposition for
The general local Hilbert space V N +M introduced above has N bosonic and M fermionic basis states. If N and M are such that they can be expressed as N = n 2 B + n 2 F and M = 2n B n F for integer n B , n F it is possible to express V N +M as a graded tensor product of two n = n B + n Fdimensional spaces V N +M = V n ⊗ V n . Here n B and n F are the numbers of bosonic and fermionic basis states of V n . Denoting the basis of V n by {|1 , . . . , |n } we can express the N + M basis states of V N +M as |α = | α ⊗ |ᾱ , α = 1, . . . , N + M , where 1 ≤ᾱ, α ≤ n are related to α bȳ α = α mod n + nδ α mod n,0 , α = α n + 1 + 1.
We note that α = n( α − 1) +ᾱ and that the fermion parities are related by α = α + ᾱ . Defining operators e α β j = | α j j β| , eᾱβ j = |ᾱ j j β |, we may express E αβ j in the form We will use this decomposition in several models considered below. In the purely bosonic case M = 0 such decompositions are possible for N = n 2 with integer n.

Super-operator formalism for Lindblad equations
We now consider a Lindblad equation (1) with a Hamiltonian H and jump operators L a acting on H L defined above. We are ultimately interested in cases where the Hamiltonian density and L a have local expansions in terms of the E αβ j . To start with we will assume for simplicity that all jump operators are bosonic. The cases where some of the jump operators are fermionic will be discussed later. The reduced density matrix can be expressed in terms of the basis states defined above as ρ = α,β ρ α,β |α β| .
The matrix elements are related to particular Green's functions of the operators E αβ j ρ α,β = (−1) In terms of components the Lindblad equation reads where we have introduced the following notations for the matrix elements of an operator O We can view the density matrix as a state in a (N + M ) 2L dimensional Hilbert space H S = H L ⊗ H L with basis states In these notations we have and the "wave-functions" ρ α,β correspond to Green's functions in the original problem. The Lindblad equation (17) can be cast in the form where the Liouvillian L for bosonic jump operators L a is given by Here we employ notations such that O = O ⊗ 1 and have defined related One can easily check that taking the scalar product of (21) with the state β| α| precisely reproduces (17). A convenient basis for expanding operators O is constructed in terms of operators E αβ n defined as E αβ n = 1 ⊗ |α n n β| .

Fermionic jump operators
If some of the jump operators are fermionic the super-operator formalism needs to be modified. Let us denote the fermion parity of the jump operator L a by La ∈ {0, 1}. When written in components the Lindblad equation still takes the form (17). However, the Liouvillian (22) is now replaced by The state representing the density matrix is also modified and now takes the form where P ± are projection operators onto states with even and odd fermion parity respectively We have (−1) F |α |β = (−1) α+ β |α |β .
It is straightforward to check that inserting (26) and (27) into the equation and expanding it in a basis of states precisely recovers (17). We stress that in our construction both bosonic and fermionic jump operators can be accommodated as long as any given jump operator has a definite fermion parity.
3 Lindblad equations as non-Hermitian two-leg ladders As we are interested in Liouvillians with local densities we focus on jump operators where the index a runs either over the sites or the nearest-neighbour bonds of a one dimensional ring. In this setting −iH − a γa 2 L † a L a and iH − a γa 2 L † a L a describe interactions along the two legs of the ladder, while a γ a L aL † a play the role of interactions between the two legs.

Single-site jump operators
In translationally invariant situations the most general bosonic single-site jump operator can be written in the form where λ αβ = 0 unless ( α + β ) mod 2 = 0. This generates "interaction terms" between the two legs of the form The other jump operator terms in the Liouvillian generate "generalized magnetic field terms" acting on the two legs † where Λ βγ = α λ * αβ λ αγ .

Single-bond jump operators
The most general bosonic jump operator acting on a bond takes the form This gives rise to quartic, cubic and quadratic "interaction terms" in the Liouvillian. The resulting explicit expression is presented in Appendix A. The extension to fermionic jump operators is straightforward.

General form of the Liouvillian
In the following we will consider Liouvillians of the form where L (a) j are jump operators that act either on site j or the bond (j, j + 1) and γ a > 0. Our aim is to identify cases which are Yang-Baxter integrable. In practice this means that we need to check whether any of the large number of integrable Hamiltonians that can be interpreted as two-leg ladder models can be cast in the particular form (35). An added complication is that we should allow for general similarity transformations, i.e. consider The spatial locality of the Hamiltonian density of the various integrable models imposes strong restrictions on the possible form of S. Transformations of the form where S j acts non-trivially only on site j are always compatible with the aforementioned local structure.

Generalized Hubbard models
The first example of a Lindblad equation that is related to an interacting Yang-Baxter integrable model was presented in Ref. [27], where it was shown that the Lindblad equation for a tightbinding chain with dephasing noise can be mapped onto a fermionic Hubbard model with purely imaginary interactions. We now briefly review some results obtained in that work. We then show that the mathematical structure that underlies the integrability of the Hubbard model quite naturally leads to a connection with a Lindblad equation.

SU(2) Hubbard model
The Hubbard Hamiltonian is given by where n j,σ = c † j,σ c j,σ . The model is integrable for any complex value of U/t [41]. In terms of the notations of section 2.1 we can choose a basis such that and concomitantly

Associated Lindblad equation
Let us consider a tight-binding model coupled to an environment by jump operators In the super-operator formalism the corresponding Liouvillian (22) is This is related to the Hubbard Hamiltonian by [27]

Integrable structure of generalized Hubbard models and associated Lindblad equations
The Hubbard model was embedded into the general framework of the Quantum Inverse Scattering Method [40] in seminal work by Shastry [31,32]. This construction was subsequently generalized to other classes of integrable models [35][36][37][38][39]. The construction is based on an R-matrix r 12 (λ) acting on the tensor product of two graded linear vector spaces V ⊗ V and a conjugation matrix C acting on V that fulfil the Yang-Baxter relation as well as the "decorated" Yang-Baxter relation In the cases considered below the r 12 (λ) is given by where Π 12 is a graded permutation operator (9) acting on V ⊗ V and whereπ is a projection operator onto a subspace of V . The R-matrix of an integrable generalized Hubbard model is then obtained by gluing together two copies [37,39,41] Here the function α(λ, µ) is given by where h(µ) is a solution of the equation The local Hamiltonian density of the integrable "fundamental spin model" [40] corresponding to this R-matrix is Here we have generalized the construction of [37] by taking the logarithmic derivative of the transfer matrix at a shifted point u 0 following Ref. [43,44]. Importantly the structure of the Hamiltonians (52) is such that they all can be related to Liouvillians of Lindblad equations. In the following we discuss a number of examples.

USW model
As a first application we consider eqn (52) for the case of the Hubbard model R-matrix [44]. The Hamiltonian of these models was first derived by Umeno, Shiroishi and Wadati in [43] and is of the form where andB j,j+1 is obtained from B j,j+1 by replacing e αβ n →ẽ αβ n . Here u 0 is a free (complex) parameter and the function h(u) is fixed by the requirement We note that the operators e αβ j are related to spinful fermion creation and annihilation operators by (39). The Hamiltonian (53) is SO(4) symmetric [43] and in particular commutes with the total particle numberN

Associated Lindblad equation
The USW model is related to a Lindblad equation with a tight-binding Hamiltonian and jump operators where the parameter u 0 is taken to be purely imaginary. In the super-operator formalism the corresponding Liouvillian (22) is This is related to the USW Hamiltonian by where the unitary transformation U is given by (44) and the parameter u is purely imaginary and related to γ by

Differential equations for correlation functions
As the jump operators are Hermitian the Lindblad equation implies the following time evolution for expectation values of (time independent) operators It is straightforward to verify that the jump operators (58) fulfil This shows that n-particle Green's functions fulfil simple, closed evolution equations. This is analogous to the case of the imaginary-U Hubbard model [27]. For example, the single-particle Green's function has the following equation of motion Here we have defined

Maassarani models
In [33,34] Maassarani introduced a class of integrable 2n-state models that generalize the Hubbard model along the lines set out in section 4.2 above. We now discuss these models in more detail. A basis of the local Hilbert space is given by the tensor product |a ⊗ |ã , a,ã = 1, . . . , n, where all states are bosonic, i.e. a = 0 = ã . While these models a priori are generalized spin models they can be related to interacting fermion models by Jordan-Wigner transformations as is done for a simple case below. A basis of operators acting on these states is then given by e ab jẽãb j . In terms of these (bosonic) operators Maassarani's Hamiltonian reads where x ab e ba j e ab j+1 + x −1 ab e ab j e ba j+1 , Here the two sets A and B form an arbitrary partition of {1, . . . , n} and x ab are arbitrary complex parameters. In the following we will simply set them equal to 1. The operators P (n) j,j+1 and C j are of the same forms as P (n) j,j+1 and C j respectively but with the replacement e ab j → e ab j . Maassarani's models are related to Lindblad equations with Hamiltonians and jump operators where c ∈ R is a free parameter. In the superoperator formalism the corresponding Liouvillian is where H but with e ab j replaced by e ab j . This is related to Maassarani's Hamiltonian by

3-state Maassarani model
The simplest Maassarani model is obtained by considering a local Hilbert space of three bosonic states. Choosing a decomposition In order to fermionize this model we embed it into an enlarged Hilbert space with four states per site, and then employ the results of section 2.1. This gives e 12 j = e 12 j e 11 j , e 13 j = e 11 j e 12 j .
Finally we carry out a Jordan-Wigner transformation After these transformations the Hamiltonian H 0 can be written in the form where is a projection operator that ensures that all sites are at most singly occupied. The Hamiltonian (77) can be viewed as the U → ∞ limit of the Hubbard model and is sometimes referred to as the t − 0 model. In terms of the fermionic operators the jump operator takes the form Choosing c = 1 we have which shows that the bath acts on the charge degrees of freedom. The Hamiltonian part H has a free fermionic spectrum [45,46], but the creation operators of the non-interacting fermion degrees of freedom are related to the c † j,σ in a non-local way [47,48]. As a result the single-particle Green's function does not obey a simple evolution equation. The time evolution is again given by the general expression (62), where the relevant commutators are

4-state Maassarani model
In the 4-state case we can express the e ab j in terms of two species of Pauli operators, cf.
The jump operators become (setting again c = 1 in (71))

Bethe Ansatz solution
The Maassarani models have been solved by Bethe Ansatz in Ref. [49]. Without loss of generality we restrict our discussion to the case where the sets A and B in (69) are given by commute with H Ma,n (U ) and with one another. Hence their eigenvalues N a ,Ñ a can be used as good quantum numbers. Following Ref. [49] we introduce integers and N ≥ N A + N B +Ñ A +Ñ B . We then define sets and finally introduce two non-intersecting ordered sets of integers 1 ≤ a j ≤ N ≤ L By ordered we mean that a j < a j+1 if a j , a j+1 ∈ A A and similarly forÃ where we require arg(b ) < arg(b +1 ) and the phases Φ and Ψ are given by The corresponding eigenvalues of L Ma,n (γ) are

String solutions and vanishing of the Liouvillian gap in the thermodynamic limit
The first two sets (89) of the Bethe Ansatz equations are the same as for the Hubbard model with imaginary interactions strength and twisted boundary conditions. This ensures that the "k-Λ string solutions" constructed in [27] are valid solutions for the n-state Maassarani models as well. A k-Λ string of length m corresponds to the following pattern of rapidities Here the string centres Λ In the framework of the string hypothesis the equation that fixes the allowed positions of the string centres Λ (m) α is obtained my "multiplying out the string" [50], which gives Taking logarithms this can be cast in the form where we have defined ϕ = m(2Φ + Ψ) mod 1.
For even lattice lengths L the J (m) α are integers with range We now focus on the particular sequence of string states characterized by integers In the limit of large system sizes L 1 the corresponding string centres follow from (96) Substituting this into our expression (94) for the eigenvalue of the Liouvillian gives This shows that in the large-L limit we have a band of Liouvillian eigenstates with eigenvalues that scale as L −2 . This establishes that the Liouvillian gap vanishes in the thermodynamic limit. Moreover, the scaling with system size suggests that the corresponding eigenmodes are diffusive. Our calculation does not rule out the existence of eigenstates with gaps that approach zero faster than L −2 .

GL(N, M ) Maassarani models
As we already mentioned above in section 4.2 the Shastry-Maassarani construction can be generalized to graded magnets based on GL(N, M ). Following Ref. [37] we consider the class of Hamiltonians where We can relate this to a Lindblad equation with Hamiltonian and jump operators

3-state GL(1, 2) model
The simplest example is the 3-state model based in GL(1, 2). Like in the case of the 3-state Maassarani model considered above we may represent the Hamiltonian in terms of canonical spinful fermion creation and annihilation operators by identifying the three states per site as Then H 0 can be represented as where P is the projection operator on singly occupied sites (78) and S + j = c † i,↑ c j,↓ . This describes correlated hopping of the up fermions, whereas the down fermions can only move through spin-flip processes. The jump operator is L j = 2n j,↑ − 1 .

Other integrable two-leg ladder models
The generalized Hubbard models considered above are all related to Lindblad equations with a single jump operator on each bond by virtue of their integrability structure. There are many other integrable models that can be represented as two-leg ladders and a question we have investigated at some length is whether some of them can be associated with Lindblad equations as well.
The Hamiltonian H is GL(N 2 ) symmetric and hence

Representation as a 2-leg ladder
The permutation models can be viewed as 2-leg ladders by employing the decomposition of section 2.1 for M = N . This provides a representation of the permutation operator as a tensor product Noting that L † j L j = 1 we conclude that the corresponding Liouvillian is where By construction the first term in (117) commutes with the second, which is γH GL(N 2 ) . As H GL(N 2 ) is invariant under all global GL(N 2 ) rotations U we conclude that (117) is integrable for choices of γ αβ and λ αβ such that

Twisting the boundary conditions
As we have mentioned above, in general we need to consider similarity transformations when trying to ascertain whether a Lindblad equation is related to an integrable Hamiltonian. A simple example is provided by considering a Lindblad equation with vanishing Hamiltonian and jump operators The corresponding Liouvillian is where we have used the decomposition 2.1 and fixed the phases φ α by where α, β,ᾱ,β, α, β are related by (12). To relate this to the GL(N 2 ) Hamiltonian we consider the canonical transformation under which the Liouvillian transforms as where we have imposed twisted boundary conditions We conclude that the Liouvillian is related to the integrable GL(N 2 ) Hamiltonian with twisted boundary conditions The integrability of twisted boundary conditions in the GL(N 2 ) models is well known [53][54][55].
The related Lindblad equation has no Hamiltonian and two sets of jump operators The corresponding Liouvillian is After a local basis rotation around the y-axis this maps onto (128) (up to a constant contribution) if we identify γ = J/4 and h = −γ .

GL(n
We now turn to particular graded magnets, where we have n 2 B +n 2 F bosonic and 2n B n F fermionic states at a given site of the lattice, where n B,F ∈ N 0 . A much studied family of integrable models is given by [52, 56- where Π j,j+1 is a graded permutation operator (9) and λ α are generalized chemical potentials. The case n B = n F = 1 gives the EKS model (a.k.a. supersymmetric extended Hubbard model). We now employ the decomposition 2.1 and choose a tensor product basis for the local Hilbert space as where α = (n B + n F )( α − 1) +ᾱ. The E αβ j 's can then be expressed as which in turn leads to the following decomposition of the graded permutation operator

Consider now a Lindblad equation with no Hamiltonian and Hermitian jump operators
Noting that we conclude that the corresponding Liouvillian is We can slightly generalize this by following the construction for the GL(N 2 ) case, e.g. we can add a Hamiltonian

Associated Lindblad equation
The Hamiltonian (140) is related to a Lindblad equation with no Hamiltonian part and a set of Hermitian jump operators After a local basis rotation τ a j → τ y j τ a j τ y j , a = x, y, z and setting the γ parameters to be equal for all jump operator terms we arrive at a Liouvillian

A comment on scaling limits
A standard way of generating integrable QFTs is by taking appropriate scaling limits of integrable lattice models. A paradigmatic example is the scaling limit of the Hubbard model, which gives rise to the integrable Yang-Gaudin model. An interesting question is whether we can carry out an analogous construction for our integrable Liouvillians and arrive at non-unitary integrable QFTs.
The answer seems to be negative. Let us consider a lattice model with Hamiltonian where c j and c † j are annihilation and creation operators of spinless fermions, and jump operators These give rise to a Liouvillian whereH 0 is of the same form as H 0 but written in terms of fermion annihilation and creation operatorsc j andc † j . The sign difference betweenH 0 and H 0 can be removed by a canonical transformationc In analogy of what we do in order to obtain the Yang-Gaudin model from the Hubbard model we now consider the scaling limit t → ∞ , a 0 → 0, ta 2 0 fixed.
In this limit lattice fermion operators are replaced by continuum fields The Liouvillian becomes A problem now occurs in the first term. If γ were purely imaginary, as is the case for the Hubbard model, we could tune the chemical potential in such a way to ensure that the prefactor remains finite in the scaling limit. This would leave us with a Yang-Gaudin model at a finite density, but with imaginary interaction strength. However, given that γ is real and positive we cannot take γ → ∞, but must keep it finite in order to describe states with finite real parts of their "energies". This means that the only scaling limit is trivial as the interaction term disappears. This would appear to be a more general feature, independent of integrability.

Some unsuccessful maps
Most of the integrable ladder models we have considered cannot be associated in a straightforward way with Lindblad equations. In the following we present some representative examples.

Perk-Schultz models
As an example we consider the N = 4 Perk-Schultz model [64,65] This can be viewed as a q-deformation of the GL(4) Hamiltonian considered above. Using the decomposition 2.1 we can rewrite H PS as As the spectra of σ z j+1 − σ z j and 1 + τ z j τ z j+1 are different the term in the second line cannot be related to a jump operator structure in this representation.

Higher conservation laws
A well-known way of obtaining integrable spin-ladder models is by considering higher conservation laws [40,66]. In case of the spin-1/2 Heisenberg XXX chain higher conservation laws H (k+1) can be obtained from the transfer matrix by taking logarithmic derivatives at the "shift point". By construction we have [H (k) , H (l) ] = 0. The Hamiltonian we want to consider here is H(b) = H (2) + bH (4) + const [66][67][68], which takes the form This can be viewed as a zig-zag ladder model by associating all even (odd) sites with the first (second) leg, which gives This is asymmetric under leg exchange in a way that precludes a direct relation with a Lindblad equation.

Alcaraz-Bariev model
The Alcaraz-Bariev two-parameter families of integrable models [69] come in two classes denoted by A ± and B ± respectively. The B ± family contains the Hubbard model as a special limit and this is the only case in which we succeeded in obtaining an interpretation in terms of a Lindblad equation. We now discuss why such a relation does not seem to exist in general for the A ± family of models. The Hamiltonian of the A ± family can be cast in the form where g = (1 + )(1 − sin θ) and T j,j+1 = −e 21 j e 12 j+1 + e 21 j e 12 j+1 + h.c., T Here we have carried out a unitary transformation on the Hamiltonian given in [69] in anticipation of relating it to a Liouvillian on a Lindblad equation. We start by noting that we require g = 0 for such an interpretation to be possible. The reason is that the only way to generate T j,j+1 is as a "cross-term" in j † j with j = ae 21 j e 12 j+1 + be 12 j e 21 j+1 + ce 22 However, such jump operators would also generate an unwanted contribution As this cannot be cancelled by introducing additional jump operators and does not feature in H ( ) A we conclude that we must have g = 0. Next we turn to the cubic terms T j,j+1 . These must arise from jump operators of the form L j = ae 21 j e 12 j+1 + be 12 j e 21 j+1 + ce 22 j .
In order to produce the cubic terms in H ( ) Combining these with the requirement that g = 0 leads to = sin θ = 1.
In this case the A ± model reduces to free fermions. We have also investigated whether carrying out a similarity transformation SH may facilitate a Lindblad interpretation. The answer appears to be negative.

Discussion
In this work we have reported our findings for a search for Yang-Baxter integrable Lindblad equations. We have focused on translationally invariant situations where jump operators act on bonds or sites of a one dimensional chain. We have derived a superoperator representation for lattice models with both fermionic and bosonic degrees of freedom, and jump operators which can be bosonic or fermionic. In this representation the Lindblad equation takes the form of a imaginary time Schrödinger equation with a non-Hermitian "Hamiltonian" with local density, which can be thought of in terms of a two-leg ladder model of interacting spins or fermions. We have then investigated which Yang-Baxter integrable two-leg ladder models can be related to such Lindblad equations in a "direct" way. Our main result is that a wide class of generalized Hubbard models can be interpreted as Liouvillians of Lindblad equations. We traced this back to their integrability structure, which is based on gluing together certain solutions of the Yang-Baxter equation in a particular way. Some of the corresponding dissipative models are physically meaningful, an example being the infinite-U Hubbard model subject to on-site dephasing noise. As the jump operators in this class of models are Hermitian, the completely mixed state is a steady state in all cases. Using the Bethe Ansatz solution we have shown for a subclass of generalized Hubbard models that the Liouvillian gap vanishes like L −2 as the thermodynamic limit is approached. The corresponding eigenstates correspond to particle-like "excitations" with quadratic dispersions, which suggests that the late-time behaviour in these models is likely to be diffusive. We have identified a few Yang-Baxter integrable Lindblad equations that are not generalized Hubbard models by showing that certain known integrable Hamiltonians can be cast in the form of Liouvillians associated with a Lindblad equation. However, in most cases we have considered such mappings are not possible. As this is often difficult to see we have presented a non-trivial case of such a failure in the Alcaraz-Bariev two-parameter family of integrable models.
We stress that in this work we have focused on a particular "direct" relation between Liouvillians of Lindblad equations and Hamiltonians of Yang-Baxter integrable models. There are known cases where it is possible to establish such relationships by means of more complicated (non-local) maps [29]. Moreover, as we pointed out in section 3.3, one ought to allow for similarity transformations that maintain locality of the Hamiltonian density in integrable models when trying to establish relations with Lindblad equations. A systematic way of doing this is by considering invariances of the Yang-Baxter equation, cf. Chapter 12.2.5 of Ref. [41]. For example, given a solution R(λ, µ) ∈ End(C ⊗ C) of the Yang-Baxter equation other solutions can be obtained as where V (λ) is an invertible n × n matrix. This allows one to introduce additional free parameters in the resulting Hamiltonian. The latter will generally be non-Hermitian, but this is not a problem in the present context of Lindblad equations. It would be interesting to pursue this line of enquiry further and a good starting point will be the models successfully related to Lindblad equations in this work.
In this we work we focused on identifying integrable Lindblad equations and only briefly explored using methods of quantum integrability to obtain physical properties. A good starting point for this is to determine the spectrum of the Liouvillian, which is given in terms of the solutions of the relevant Bethe Ansatz equations. It is well understood that the nature of solutions to Bethe Ansatz equations changes quite substantially when a parameter is made complex, as this results in the "scattering phases" acquiring magnitudes different from unity. In practice this means that the structure of solutions to the Bethe Ansatz equations, which is usually encoded in appropriate string hypotheses, must be revisited and typically becomes more involved. Even in the simplest case of the Hubbard model the structure of Bethe Ansatz roots for Liouvillian eigenstates with eigenvalues that have large real parts and non-zero imaginary parts appears to be non-trivial. We plan to report on this issue in a future publication. Ultimately one would like to determine the dynamics of general Green's functions Tr ρ(t)E α 1 β 1 j 1 . . . E αnβn jn (167) for evolution from a given initial density matrix ρ(0). In some of the cases discussed above this is relatively simple because the equations of motion for these Green's functions decouple and for two-point functions can thus either be integrated numerically or determined from the exact Liouvillian eigenstates in the two-particle sector [70]. In cases like the 3-state Maassarani model a more involved analysis is required and it would be interesting to investigate this case in more detail.
A Structure of the Liouvillian for the most general jump operator acting on a bond The most general two site bosonic jump operator with nearest-neighbour interactions is This gives rise to interaction terms between the two legs of the ladder L j L † j = I (2) j + I where I (n) j involves n Hubbard operators E αβ j , E αβ j . The interaction along a single rung of the ladder is while the three and four point interactions on a given plaquette are given by There are also interaction terms along the two legs of the ladder where f αβγδ = λ * βα λ γδ + η λ * ηα µ ηβγδ + λ * ηγ µ αβηδ + 1 2 ην (−1) ( α+ β )( η + γ ) µ νβηδ µ * ναηγ . (173)