Exact steady state manifold of a boundary driven spin-1 Lai-Sutherland chain

We present an explicit construction of a family of steady state density matrices for an open integrable spin-1 chain with bilinear and biquadratic interactions, also known as the Lai-Sutherland model, driven far from equilibrium by means of two oppositely polarizing Markovian dissipation channels localized at the boundary. The steady state solution exhibits n+1 fold degeneracy, for a chain of length n, due to existence of (strong) Liouvillian U(1) symmetry. The latter can be exploited to introduce a chemical potential and define a grand canonical nonequilibrium steady state ensemble. The matrix product form of the solution entails an infinitely-dimensional representation of a non-trivial Lie algebra (semidirect product of sl_2 and a non-nilpotent radical) and hints to a novel Yang-Baxter integrability structure.


Introduction
Nonequilibrium transport problem in extended low-dimensional (say one-dimensional, 1D) quantum systems is an important current topic in statistical mechanics with possible links to experiments in condensed matter systems [1,2]. Among the most important open issues are (i) classification or identification of possible transport behaviors, ranging from ballistic, via diffusive (normal or anomalous), to insulating, and understanding their microscopic mechanisms [3,4], and (ii) developing nonequilibrium quantum thermodynamics [5] and a theory of nonequilibrium quantum phase transitions (see e.g. [6,7]).
A convenient setup for studying far from equilibrium relaxation dynamics or steady state situations which support macroscopic currents of charge/particles, magnetization, or energy/heat, is to couple a 1D system of strongly interacting quantum particles to incoherent forcing governed by two reservoirs attached at each end of the particle chain and assign them different effective thermodynamic potentials. This can be achieved e.g. by choosing simple Markovian dissipation channels which operate as quantum jumps, i.e. pump-in or absorb-out the elementary excitations at the surface (boundary of the chain). The rest of the system is chosen to be unaffected by the dissipation and hence evolves according to fullycoherent unitary evolution. For a derivation and physical justification of such an approach, see Refs. [5,8].
With the hope of being able to take advantage of their rich and elegant mathematical content, one addresses integrable systems with local interactions first. In this light one obtains a toy model to study dissipative integrable theory with surface 'non-unitary sources'. This model could be in some sense also regarded as a quantum analogue of classical stochastic exclusion processes [9,10,11]. Focusing initially on the steady states alone, the aim is to be able to isolate regimes where complexity of the steady state density operator is drastically reduced, opening a possibility of finding an efficient and exact representation in terms of the matrix product state. Ever since the first solutions in this direction have been presented, addressing quasi-free theory [12] and the paradigmatic strongly interacting case of the (anisotropic) Heisenberg (XXZ) spin-1/2 chain [13,14], the quest for new integrable out-of-equilibrium scenarios continues [15,16], with some recent attempts [17,18,19] of putting these searches under the common roof of theory of integrable quantum systems [20,21].
It has been argued [15,18] that explicit steady state solutions of boundary driven Liouvillian (Lindbladian) flows pertaining to certain integrable models arise as a consequence of the underlying quantum group symmetry of the model. The latter provides a prerequisite condition for the solution in the bulk which needs to be fine-tuned with the form of the quantum noise process applied to the system's boundaries. Two principal insights have been made, namely (i) to re-write the matrix product representation of solution in terms of monodromy matrices with Lax operators arising from solutions of the universal Yang-Baxter equation associated to a symmetry algebra of an interaction, and (ii) to allow for non-unitary irreducible representations over infinitely-dimensional vector spaces [22]. In the prototype case of the Heisenberg spin-1/2 model it actually appears that fundamental (local) building blocks that generate the solution inherit symmetry from the interaction, despite that the latter is finally broken at the level of Liouvillian flow and density operator. Below we demonstrate however, how central objects of our construction could also admit a (non-trivial) continuous symmetry which does not respect that of an integrable bulk interactions.
To this end, we consider an integrable SU(3)-invariant spin-1 chain, commonly referred to as the Lai-Sutherland model [23,24] 1 , and employ a pair of Lindblad jump operators which couple only two extreme levels at the chains end. The intermediate level, which can be viewed as a hole particle, is thus protected from the environment and its number is preserved throughout the (dissipative) evolution. Henceforth, such Lindbladian flow is reducible to a (thermodynamically) infinite number of sectors corresponding to subspaces with fixed 'hole doping'. This allows for a possibility of constructing a grand-canonical steady state ensemble, with chemical potential being an additional parameter which controls the average number of holes (or the filling factor). The paper is organized as follows. In section 2 we introduce the open threestate Lai-Sutherland model and specify suitable 'integrable' boundary dissipative processes in the framework of Markovian (Lindblad) master equations. In section 3 we rigorously construct the solution of the steady state in terms of an infinite rank matrix product ansatz. In section 4 we discuss several important physical properties of the solution: we introduce the grand-canonical nonequilibrium steady state ensemble (subsection 4.1), describe a formal computation of local physical observables (subsection 4.2), discuss graph-theoretic interpretation of the solution in terms of sums over walks (subsection 4.3), characterize the symmetries (subsections 4.4, 4.5, 4.6) and discuss possible connection to quantum inverse scattering method (4.7). Finally, we conclude in section 5.

Open Lai-Sutherland model far from equilibrium
Consider a finite chain of n sites and let H 1 C 3 be a local quantum ('physical') space associated to each spin-site x ∈ {1, 2, . . . , n}. The entire 3 n dimensional many-body quantum space H s is constructed as n-fold tensor product of local spaces, H s = H ⊗n 1 . Using the Weyl matrix basis {e i j = |i j| ; i, j = 1, 2, 3} of End (H 1 ) = gl 3 , we define a full set of local generators of the matrix algebra F = End (H s ) as e i j ½ d defining a d−dimensional unit matrix, satisfying the Lie algebra relations The spin-1 Lai-Sutherland model [23,24] for a chain of n sites is given by the Hamiltonian H ∈ F, where form independent spin-1 variables (local s = 1 representations of su 2 ) satisfying Straightforward inspection shows that the local Hamiltonian h x,x+1 -the interaction -is in fact just the permutation operator between neighboring sites The local Hilbert state basis is therefore given by a triple of states |1 ≡ |↑ , |2 ≡ |0 , |3 ≡ |↓ , which can be interpreted as three different particle species; respectively, as spin-up particles, holes, and spin-down particles. The model then becomes equivalent to the so-called supersymmetric t-J model [26]. Lai-Sutherland chain is a multi-component quantum model and we may associate with it a skew-symmetric tensor of particle currents, with two-site density which, by construction, satisfies the following continuity equation J i j can be considered as a partial current of the particle of species i into particles of species j. The total current of particles of species i, then fulfills the continuity equation where e i i x can be considered as the operator of particle density of species i. We shall now open the Lai-Sutherland chain and couple it to the environment via Markovian processes which act only on local quantum spin spaces at the boundary, i.e., at x = 1 and x = n. The many-body density operator ρ t , t ∈ R + , considered as an element of F which may be here considered as a Liouville vector space of operators, then evolves according to Liouvillian semigroup with time-independent generator -the LiouvillianL ∈ End (F) being split into non-dissipative partL 0 (ρ) ≡ −i[H, ρ] governing unitary Liouville-von Neumann evolution, and a dissipatorD ∈ End (F) describing the incoherent, dissipative (non-unitary) processes of overall strength ε. The latter is given in terms of a set of jump operators {A α ∈ F} and takes a general canonical Lindblad [27,28] form In particular, we install a single local jump operator at each end of the chain: Two dissipation channels, interpreted as the left and right magnetization bath, perform the processes |↑ → |↓ and |↓ → |↑ , respectively, with the rates ε. Both processes keep the hole state |0 unaffected. Since also the bulk dynamics generated byL 0 conserves the number of particles of each species, it follows that the whole Liouvillian dynamics (master equation) preserves the number of holes. More precisely, defining the hole-number operator N 0 ∈ F as we have that the set of all, Hamiltonian and jump operators, commute with N 0 which implies that N 0 generates a strong [29] U(1) symmetry of the Liouvillian flow (11). N 0 foliates the physical space into n + 1 orthogonal eigenspaces, The theorem A.1 of Ref. [29] then guarantees that the full Lindblad dynamics (11) is closed on F (ν) = End (H (ν) s ),L (ν) =L| F (ν) , and that a fixed point ρ (ν) ∞ = lim t→∞ exp(tL (ν) )ρ (ν) 0 -nonequilibrium steady state (NESS)exists for each symmetry subspace flow 2 , The theorem by Evans [31] can then be used to show uniqueness of NESS ρ (ν) ∞ for each fixed ν. In the next section we shall outline a simple algebraic procedure for actual explicit construction of density operators ρ (ν) ∞ .

Matrix product solution
LetP (ν) ∈ End (F) be an orthogonal projector to F (ν) . We define a universal density matrix of NESS as a direct sum of non-trivial solutions of (16) for all ν, being solution of the fixed point equation (16) as well. The state ρ ∞ shall be sought for in terms of Cholesky factorization (in analogy to previous solutions of XXZ [14] and Hubbard [16] models) where S n (ε) ∈ End (H s ) is some yet unknown operator which is represented by an upper triangular matrix in the computational basis |i 1 , . . . , i n . Introducing an auxiliary Hilbert space H a -separable, but of infinite dimensionality as will become clear later -we define the monodromy operator M(ε) ∈ End (H s ⊗ H a ) as a spatially-ordered product of some local Lax operators 3 L x (ε) ∈ End (H s ⊗ H a ), Throughout the paper, the upright-boldface notation designates object which are not scalars in auxiliary space H a . Index free Lax operator can be defined as L(ε) ∈ We further assume existence of a special state |vac ∈ H a , such that Cholesky factor writes as the auxiliary expectation value of monodromy operator, or equivalently, as a matrix product operator (MPO) Fixing an arbitrary, fixed orthonormal basis we denote the second copy of auxiliary space carrying conjugate representation of L i j as H a . One can then write MPO formulation of NESS density operator , and Note the transposition in the quantum space of the conjugated factor of (22). Here a two-leg monodromy operator and a product of a pair of vacua vac| = vac| ⊗ vac|, |vac = |vac ⊗ |vac have been introduced, so that (23) is merely a formal rewriting of (18). These definitions become particularly handy when we consider evaluation of expectation values of local observables with respect to NESS ρ ∞ (ε). Let η := iε be a complex-rotated coupling parameter and let us (for convenience) relabel the quantum space matrix elements of the L-operator as The key results of this paper are the following: Theorem 1. Suppose that 9 matrix elements {L i j } generate the Lie algebra g defined by commutation relations, with a representation over the Hilbert space H a satisfying the following conditions Then, the universal solution (17) to NESS fixed point condition (16) is given via Cholesky factorization (18) with explicit MPO expression (21) for S n (ε) with η = iε.

Theorem 2.
A possible irreducible explicit representation of Lie algebra g (26) satisfying (27) is given as in terms of three auxiliary degrees of freedom with a three dimensional lattice and a complex spin (Verma module of sl 2 ) with |vac = |0, 0, 0 being the highest-weight-state. The complex spin parameter p should be linked to dissipation parameter via Proof. The proof of the theorems is based on verifying that the Lie algebra g, given by (26), can be equivalently defined by means of an identity over End (H s ⊗ H a ) in the form of local operator divergence (LOD) condition (customary referred to as the Sutherland equation which is equivalent to zero curvature/Lax condition), with the-so-called boundary operator B x (ε) ∈ End (H s ⊗ H a ) -operating nontrivially only in the local quantum space Identification of (26) Multiplying LOD by a string L 1 · · · L x−1 from the left and a string L x+2 · · · L n from the right, summing over x and taking vacuum expectation value yields the global almost conservation condition for the Cholesky factor (the so-called defining relation, analogous to similar relations in other integrable nonequilibrium models [14,16]), Consequently, by expanding the unitary part of LiouvillianL 0 , in conjunction with (35), and employing the definition (22), the steady state condition (16) yields a decoupled system of boundary equations where two-leg boundary operators have been defined. Note that, due to (33), b x = iεs 3 x = −b x for ε ∈ R. The last two lines of (26) indicate that pairs of auxiliary operators (t + , t − ) and (u + , u − ) span the Weyl-Heisenberg algebra. In conjunction with the highest weight conditions (27) this fixes (uniquely, up to unitary transformations) the representation of (t + , t − ) and (u + , u − ) to be that of a Fock space of two canonical bosonic (oscillator) modes, specified by creation/annihilation operators, suggesting that the auxiliary space H a is perhaps just a two-mode boson Fock space. While realization for all the other generators consistent with the bulk algebra g is not difficult to construct (e.g. v ± , l ↑ + l ↓ can be just the Schwinger boson representation of su 2 -see 5th line of (26)), it turns out not to be consistent with the boundary conditions (27) 4 . Therefore the auxiliary space H a has to contain (at least) one additional degree of freedom.
Ultimately, in order to fulfill (37), a straightforward calculation shows that it is enough to add a Verma module S of complex spin representation (30) of sl 2 and consider a triple-product space H a B ⊗ B ⊗ S = lsp{| j, k, l ; j, k, l ∈ Z + }, and find a representation of the algebra (26) which is compliant with conditions with vacuum being given by the ground state |vac ≡ |0, 0, 0 . These requirements are all satisfied by choosing representation (28,29,30) with p being fixed (31) as required by the conditions in the first two lines of (27). The last two lines of (27) hold due to highest-weight-property of |vac . As such a representation is clearly irreducible, this concludes the proof of theorems 1 and 2.
Remark. All MPO (21) amplitues, i.e., matrix elements of the Cholesky factor of the density operator are polynomials (of order not more than n) in η = iε with integer coefficients. This is a simple consequence of Wick theorem, or representation of Theorem 2.

Discussion
The formulae (18,21,25,28,29,30,31) are the main result of this paper: They generate explicit construction of a many-body density matrix of a family of degenerate NESSes ρ (ν) ∞ =P (ν) ρ ∞ for any number of holes ν ∈ {0, 1 . . . n}. The computational complexity of obtaining any locality-based information about the state ρ ∞ , say to compute its matrix elements of the type i 1 , . . . , i n | ρ ∞ | j 1 , . . . , j n or local observables, is at most polynomial in n. Since the eigenspaces H (ν) of number-ofholes operator N 0 or orthogonal, one can also split decompose the Cholesky factors Projected Cholesky factor satisfies a projected defining relation (35) and can be expressed in terms of a constrained or microcanonical MPO Note that since [S (ν) , N 0 ] = 0, the Kronecker-δ constraint can just as well be replaced by δ ( x δ jx,2 ),ν as only operators e i 1 j 1 ⊗ · · · ⊗ e i n j n for which x δ i x ,2 = x δ j x ,2 appear in MPO expansion (21). We note two limiting cases of our new solution. For zero hole sector ν = 0 one obtains exactly the fully polarized boundary driven isotropic (XXX) Heisenberg spin-1/2 chain and reproduces the solution of Ref. [14] as formulated in [17]. The other extreme case (ν = n) is the so-called dark state, i.e. a pure state ρ (ν=n) ∞ = (e 22 ) ⊗n = |2, 2 . . . 2 2, 2, . . . 2| which is unaffected by the dissipation, i.e. it simultaneously annihilated byL 0 andD,L 0 ρ (n) ∞ =Dρ (n) ∞ = 0.

Grand-canonical nonequilibrium steady state ensemble
Any convex mixture of states ρ ∞ = ν c ν ρ (ν) ∞ , c ν ∈ R + , is a valid NESS density operator as well, which factorizes (18) with a Cholesky factor S n = ν √ c ν S (ν) n . Microcanonical constraint in (44) seems cumbersome as it prevents facilitating transfer matrices for computation of local observables. There seems to be a particularly attractive option which overcomes this problem. Namely, one may define a grand canonical nonequilbrium steady state (gcNESS) ensemble by taking a hole chemical potential µ with c ν = exp(µν): Clearly, the addition theorem for exponential function erases the constraint in MPO expansions: where the chemical potential only modifies the components of the Lax operators as Moreover, introducing a transfer vertex operator we define the nonequilibrium partition function and express it via the transfer matrix method Z n (ε, µ) = tr (ρ ∞ (ε, µ)) = vac| ( (ε, µ)) n |vac .
As usual, we expect the fluctuations ν 2 /n 2 − r 2 to be thermodynamically small. We can make a simple assertion about the thermodynamic behavior of Z n . In the regime, n → ∞, one can write an asymptotic expansion where f j (n) are all possible -perhaps non-analytic -super-linear dependencies satisfying lim n→∞ n f j (n) = 0 (as we shall argue later the most typical being f (n) = n log n), and o(n) is the standard 'little-o' notation. Here we have assumed that the chemical potential µ is an intensive quantity, i.e., independent of n. According to the definition (51), the doping should be confined to the unit interval, 0 ≤ r ≤ 1, ∀n, so the following identities follow in the thermodynamic limit i.e., coefficients in front of all super-linear dependencies can not depend on chemical potential.

Computation of local observables
Expectation values of (local) observables can be extracted by facilitating auxil- be a generic local observable supported on a sublattice between sites x and y. Then, a formal expression can be calculated from the MPO representation of ρ ∞ (ε, µ) by tracing out the physical space H s and associating to each observable X [x,y] a corresponding vertex operator via a mapping Λ ℓ : H ⊗ℓ 1 → H a ⊗ H a , where ℓ = y − x + 1, using the prescription For a complementary part of a lattice, i.e. where X [x,y] operates trivially, one has the transfer vertex operator = Λ 1 (½ 3 ), eq. (49), so the final expectation value reads For example, for on-site observables we have auxiliary vertex operators Λ 1 (e i j ) = j i , e.g. for magnetization density Λ 1 (s 3 ) = 11 − 33 . As for two point observables, we consider an interesting example of the current density tensor Stationarity (time-independence) of NESS and continuity equation (8) imply spatialindependence of current expectation values. In auxiliary transfer matrix formulation (49) this implies commutation of transfer vertex operator with current vertex operators when when projected onto subspace of states created upon action of on the vacua, namely Additionally, using representation given in Theorem 2 and highest weight nature of the vacuum, one can with some effort express the expectation values of total current operators (9) in terms of the nonequilibrium partition function (50) Using parametrization of thermodynamic scaling (52) we can express large n asymptotics of the spin current J s = J 1 − J 3 as For example, in the limiting case r → 0 of XXX spin 1/2 chain, we have [14] a single term in the sum of (52) with f 1 (n) = n log n and β 1 = 2, implying a subdiffusive scaling J s x ∝ n −2 . We claim that such scaling may be quite generic, yielding a power-law scaling of the current, β 1 being the power-law exponent.
In order to obtain more precise, or explicit results on the thermodynamics of observables in our nonequilibrium model one would need to have a better understanding of the algebra of auxiliary vertex operators generated by i j and of analytic properties of the partition function Z n , such as in the case of XXX model [14,18,32]. A very attractive question would be two investigate ε − µ phase diagram of the open Lai-Sutherland chain and to analyze possibilities of nonequilibrium phase transitions.
For example, one may define the minimal and maximal doping, accessible by an intensive (n−independent) chemical potential in the non-equilibrium grandcanonical state, as r ± := lim µ→±∞ lim n→∞ n −1 ∂ µ log Z n (ε, µ) (note the importance of the order of the limits!). Depending on the tails of the ν−dependence of tr ρ (ν) ∞ , one may have r − = 0, or r − > 0 (and r + = 1, or r + < 1). In the latter of the case(s) one may hence expect a phase transition at r = r ± , whereas the rest of the doping range, [0, r − ] (or [r + , 1]), is only accessible by considering a carefully chosen n−dependent chemical potential µ(n). Eqs. (53) imply that the current scaling exponent β 1 (r) is constant on the entire range [r − , r + ], nevertheless it may be different than the XXX exponent β 1 | r=0 = 2, which can be obtained from our solution of the Lai-Sutherland chain via different order of the limits lim n→∞ lim µ→−∞ . It is thus in principle possible to find even a normal diffusive exponent β 1 = 1 and/or transitions to other, say super-diffusive or ballistic behaviors with changing the doping r. Investigating these exciting questions will be a subject of intense future work.

The solution as a walking graph state
In Ref. [16] a universal interpretation of NESS density operators of integrable boundary driven chains have been given in terms of walking graph states (WGS). WGS can be considered as an appealing and compact formulation of matrix product state with infinite dimensional matrices having a simple local structure.

Characterization of the Lie algebra.
The Lie algebra g, eq. (26), has a non-trivial structure. It can be decomposed however (according to Levi theorem) as a semi-direct product of a solvable ideal (radical) and semi-simple part, In our case a is given by lsp{v + , v − , l + }, writing l ± := l ↑ ± l ↓ , i.e. a sl 2 is isomorphic to spin algebra, whereas r = lsp{t ± , u ± , l − , l 0 }, generates a (non-nilpotent) radical. The element l 0 lies in the center of g. It is worth noticing also, that parameter η can be fully removed from the algebra (26) by diving all generators by η, except t + , u − .

Symmetries of the Lax and transfer operators
In contrast to situation with XXX or XXZ spin 1/2 chain [14,15], it might seem surprising here that the fundamental local unit, the Lax operator L, does not exhibit the full SL(3) symmetry of the interaction. However, the dissipative driving breaks the SL(3) symmetry, resulting in only remaining U(1) global symmetry of the Liouvillian flow generated by It remains to be investigated whether other gauges exist in which Lax operators exhibit non-Abelian symmetry. Furthermore, transfer vertex operator exhibits U(1) × U(1) symmetry, i.e. there exist two conserved auxiliary-space operators ± ∈ End (H a ⊗ H a ), Conserved operators look particularly useful in the auxiliary spin-boson representation of Theorem 2 (note that η = −η for physical (real) dissipation) where we have now four independent bosonic modes b ↑↓ , b ↑↓ and two complex spins s α and s α with representation parameters p = 1/2 − 1/η and p = 1/2 + 1/η Computation of nonequilibrium partition function (50) should hence be performed on a 4D sub-lattice of a 6D lattice (a basis of H a ⊗ H a , {| j, k, l, j, k, l ; j, k, l, j, k, l ∈ Z + }) where, say j, k can be eliminated using constraints: This is analogous to 'diagonal reduction' of the transfer matrix for the open XXZ chain proposed in Refs. [14,19].

Symmetries of the Liouvillian flow and its fixed point
Besides the strong U(1) symmetry generated by N 0 , the full Liouvillian flow has another U(1) symmetry generated by magnetization operator M (65), as noted in subsect. 4.5. This is a weak symmetry in the sense of Ref. [29] and can be formally written as The Liouvillian flow and NESS display additional Z 2 -parity symmetry which is a composition of lattice reversalR ∈ End (F), R(e i 1 j 1 ⊗ e i 2 j 2 ⊗ · · · ⊗ e i n j n ) = e i n j n ⊗ e i n−1 j n−1 ⊗ · · · ⊗ e i 1 j 1 , and a product of local mirror symmetriesŜ ∈ End (F) which exchange spin-up and spin-down particles,Ŝ =Ŝ ⊗n 1 ,Ŝ 1 (e i j ) = e 3−i+1,3− j+1 , namely [RŜ,L] = 0 andRŜ ρ ∞ = ρ ∞ .
Cholesky factor S n (ǫ) however acquires another Z 2 -parity symmetry. By means of transposition mapT ∈ End (F), one findsRŜ S n =TŜ S n = S n .
Notice thatTŜ pertains to the symmetry with respect to exchange of the bosonic modes in auxiliary space.

Transfer matrix property of Cholesky factors
Similarly to XXX [17] and Hubbard [16] chains the Cholesky factor is found, empirically by checking explicitly systems of small size n, to exhibit a transfer matrix properly, namely [S n (ε), S n (ε ′ )] = 0, ∀ε, ε ′ ∈ C.
This property justifies calling the L-operator a quantum Lax matrix with M-operator being a corresponding monodromy matrix and S n (ε) = vac| M(ε) |vac the corresponding transfer matrix where the trace in infinitely dimensional auxiliary space is replaced by a ground state expectation value [17]. It remains an open issue though to prove that L(ε) belongs to solutions of the quantum Yang-Baxter equation. Clearly, the property (79) can be extended to grand-canonical objects due to orthogonality of subspaces H (ν) , namely [S n (ε, µ), S n (ε ′ , µ ′ )] = 0.

Conclusions
We have presented an explicit infinite rank matrix-product state construction of an exact solution for a grand-canonical nonequilibrium steady state of boundarydriven integrable SU(3)-symmetric spin-1 chain. Beside the external chemical potential, controlling an average filling factor of conserved "hole particles", the NESS (continuously) depends on also on the bath coupling parameter, describing strength of incoherent processes at the boundaries. Quite remarkably, the elements of the main building block (the L-operator) generate a Lie algebra of non-trivial structure whose simple part is given by classical sl 2 algebra. Despite the fact that L-operator does not exhibit invariance with respect to any continuous non-Abelian symmetry, empiric evidence clearly suggests that it generates a quantum transfer matrix of an (abstract) integrable system, indicating that a Yang-Baxter structure is sitting underneath. Another central aspect to the problem is that the auxiliary space can be factored into three-fold product of infinite-dimensional quantum spaces -a Fock space of two independent bosonic modes and one generic representation of sl 2 (Verma module) -depending on one complex continuous representation (spin) parameter p. In order to fulfill the boundary system of equations which guarantee solutions to our problem, the value of spin parameter must be chosen according with the dissipation rate. The solution contains, as a special extreme case, previously known NESS for symmetrical driving of the spin-1/2 (isotropic) Heisenberg model. It remains an open issue whether presented solution admits any integrable continuous deformation (q−deformation), enabling generalization to anisotropic versions of the Lai-Sutherland model, say for the Perk-Schultz model [33]. Another interesting open issue is generalizing our solution to more than three (N > 3) component Sutherland models.