Dissipation-driven integrable fermionic systems: from graded Yangians to exact nonequilibrium steady states

Using the Lindblad master equation approach, we investigate the structure of steady-state solutions of open integrable quantum lattice models, driven far from equilibrium by incoherent particle reservoirs attached at the boundaries. We identify a class of boundary dissipation processes which permits to derive exact steady-state density matrices in the form of graded matrix-product operators. All the solutions factorize in terms of vacuum analogues of Baxter's Q-operators which are realized in terms of non-unitary representations of certain finite dimensional subalgebras of graded Yangians. We present a unifying framework which allows to solve fermionic models and naturally incorporates higher-rank symmetries. This enables to explain underlying algebraic content behind most of the previously-found solutions.


Introduction
Remarkable progress in experiments with cold atoms [1][2][3][4][5][6][7] has greatly impacted theoretical research in the area of quantum many-body dynamics [8][9][10][11][12][13][14][15][16]. Quantum systems which reside in the proximity of a quantum integrable point have received a great amount of attention. Non-ergodic character of these systems was revealed through anomalous relaxation and absence of conventional thermalization, and paved the way to study new paradigms in quantum statistical mechanics such as pre-thermalization [17][18][19][20][21][22] and equilibration towards generalized Gibbs ensembles [23][24][25][26][27][28][29][30][31][32]. In an idealized scenario which neglects integrability-breaking perturbations, integrable interacting systems were shown to permit a universal classification of local equilibria [33][34][35] based on a complete set of local conservation laws [36]. Equilibrium statistical ensemble however constitute a fairly small set of quantum manybody states and are outside of perturbative regime insufficient to capture physically interesting situations in which systems support particle and energy currents. An important step towards realizing accessing genuine far-from-equilibrium regimes is to devise an efficient computational framework for accessing regimes of strongly-correlated quantum dynamics which often lie beyond the reach of traditional techniques. Switching from the Hamiltonian approach for closed systems to the open system perspective [37][38][39] offers a promising route to achieve this.
A quantum system is regarded as an open system when as a result of interactions with its surroundings experiences incoherent loss of information (quantum decoherence), making a system evolving according to an effective non-unitary evolution law. In a highly controlled environment however, an irretrievable loss of information due to quantum noise may sometimes even act as a resource [40][41][42][43]. Quantum noise is typically modelled either as a stochastic process, or alternatively via deterministic evolution laws in the form of quantum master equations, where the unitary dynamics is supplemented with additional non-Hamiltonian effective terms. The simplest master equations are Markovian [44,45] and thus entirely discard memory effects between a system and its environment.
Despite quantum many-body systems which undergo dissipation typically evolve to either trivial states, or highly entangled states of prohibitive complexity, there remarkably exist certain non-trivial examples of quantum dissipative Markovian dynamics where an intricate interplay between noise and coherent evolution results in stationary state which are analytically tractable. Integrability of the central model is of central importance here, as it makes it possible to identify Markovian particle reservoirs which, as explained in the manuscript, induce certain 'symmetry protected' nonequilibrium states. The main objective in this regard is to isolate scenarios in which the steady states of dissipative many-body dynamics are of low complexity and permit an exact analytic description in terms of matrix-product states. This programme has been initially employed in classical exclusion processes [47,48] and relatively recently applied to quantum chain of non-interacting fermions [70,71]. The same ideas have been shortly after expanded also to a few representative interacting exactly solvable many-body Hamiltonians, such as the Heisenberg spin chain [72][73][74] and the fermionic Hubbard model [75,76], which led to various applications (see e.g. [72,[77][78][79]). For a historical perspective on the subject the reader is referred to the recent topical review [80].
Despite many promising advancements on the subject, it is rather unsatisfactory that the structure of these solutions still remain elusive and poorly understood. Indeed, no common framework which would explain the origin and meaning of integrable dissipative boundaries and offer a systematic way to extend the results to more general scenarios has been proposed to the date. In particular, all previous attempts to understand the internal structure of these exactly solvable instances based on 'first symmetry principles' and algebraic concepts of Yang-Baxter integrability have been mainly unsuccessful, although a few central insights have been made in [81,82] which unveiled the Lax formulation and highlighted the importance of non-unitary representations of quantum groups. However, a comprehensive group-theoretic approach which would enable to construct a larger class of solutions from first principles remained unknown.
The primary goal of this work is study formal aspects of integrable quantum chains driven far from equilibrium with aid of incoherent Markovian reservoirs attached to their ends. By continuing the bottom-up approach initiated in [81], we shall uncover the symmetry content behind some of the solutions obtained in the previous works, extend these results to models based on higher rank algebras and discuss the paradigm of integrable steady states from the standpoint of representation theory of quantum algebras. This work offers a unifying algebraic construction for an entire class of exact nonequilibrium states belonging to the so-called rational integrable quantum spin chains by making use of tools of quantum integrability theory. Specific instances which have been derived in the previous work with different techniques(cf. [73,74,81,82]) are thus naturally incorporated in a common framework. Moreover, by using graded vectors spaces and Lie superalgebras, we present how to accommodate for fermionic degrees of freedom [83,84] and derive a new class of steady-state solutions for interacting integrable fermionic chains with SU (n|m)-symmetric Hamiltonians. Further quantitative analysis of the constructed solutions go beyond the main scope of the present study and will be thus omitted. The task of computing correlation functions can however be carried out by using standard techniques based on matrix-product states, see e.g. [80].
One of the central insights of our approach is rooted in the universal factorization property of quantum Lax operators. This leads to the so-called 'partonic' Lax operators can be regarded as the elementary constituents of Yang-Baxter integrable systems and are intrinsically related to the notion of Baxter's Q-operators [85][86][87][88], a widely used concept in the Bethe Ansatz diagonalization techniques. The observations that partonic Lax operators which realized over non-unitary irreducible modules may be used as local building units of exact steady-state solutions to certain Lindbladian dynamics is however a curious unconventional feature which displays their proper nonequilibrium character.
Outline. The paper is organized as follows. In the preliminary section 2 we give a quick introduction to the Lindblad master equation and briefly review some basic concepts regarding graded vectors spaces. In section 3, we proceed by introducing an out-of-equilibrium protocol by coupling a quantum chain of interacting particles to incoherent particle reservoirs attached at its ends. We subsequently present the main algebraic structures which are afterwards used in the construction of the steady-state solutions. The notion of graded Yangians is defined in section 4, and identify finite dimensional subalgebras which are intimately related to the Baxter's Q-operators. In section 5 we outline a unifying construction for a class of nontrivial current-carrying steady states, and present a few explicit examples of widely studied integrable spin and fermionic chains. In section 6 we provide some technical remarks on the notion of vacuum Q-operators, and conclude in section 7 by summarizing the main results and providing an outlook.

Lindblad master equation
In the approach of open quantum systems [38,39], the time-evolution of a density operator ρ(t) of the reduced density matrix of a central system (which presently represents a onedimensional system of spins or interacting fermions) is governed by a completely positive and trace-preserving map V(t), reading compactly The Liouville propagator obeys the semi-group property V(t 1 + t 2 ) = V(t 2 )V(t 1 ). Notice that, in contrast to the unitary propagator of the Hamiltonian evolution, the generator V(t) is not invertible. The generator L takes the Lindblad form [44,45] where L 0 ρ ≡ −i[H, ρ] is the ordinary Liouville-von Neumann unitary dynamics generated by the Hamiltonian H of the central system, while D is the dissipator which fully encodes an effective description of the environment and admits a canonical resolution in terms of the Lindblad operators A k , Here each Lindblad 'jump operator' A k acts as an independent incoherent process. 1 In our application, A k will be used to model incoming and outcoming particle flows through the boundaries of the quantum chain. A particular advantange of such a nonequilibrium protocol is to have a simple setup for obtaining exact or approximate results for genuine far-fromequilibrium states, reaching beyond the traditional linear response theory and quasi-stationary regimes described with the hydrodynamic approach [13,92,93].
In this work we shall exclusively restrict our considerations to the steady states. The latter correspond, by definition, to fixed points of the Liouville dynamics, ρ ∞ = lim t→∞ ρ(t). This means that a steady state is an operator ρ ∞ from the kernel (null space) of the generator L, We will also encounter situations when dim ker L > 1, which physically corresponds to degenerate steady states and leads to higher dimensional steady-state manifolds.

Graded vector spaces
In order to incorporate fermionic degrees of freedom in our description we shall make use of graded vector spaces. A local Hilbert space attached to a site in a quantum chain is denoted by C n|m , where integers n and m in the superscript signify the number of bosonic and fermionic states, respectively. Below we briefly recall a few basic notions of graded vectors space and refer the reader for a more detailed exposition to appendix A. The two types of states are distinguished by the Z 2 -parity, The mapping p equips C n+m with a Z 2 -grading: if a belongs to a subset of bosonic (fermionic) indices we assign it a parity p(a) = 0 (p(a) = 1). Gradation is naturally lifted to vectors spaces C n+m and furthermore to the Lie algebra of linear operators acting on C n+m . Specifically, by adopting the distinguished grading in which p(a) = 0 for a ∈ {1, 2, . . . n} and p(a) = 1 for a ∈ {n + 1, . . . m}, the space of (n + m)-dimensional matrices on C n+m block-decompose into the bosonic (even) subspace V 0 and fermionic (odd) subspace V 1 . The two subspaces are typically referred to as the homogeneous components. The fundamental gl(n|m) representation, denoted by V n|m , is spanned by a basis of matrix units E ab , (E ab ) ij = δ ai δ jb . The action of the Lie bracket adjusted to the grading is expressed as Since exchanging two fermionic states results in a minus sign, the presence of fermionic states in a graded tensor product space non-trivially affects the multiplication rule. Namely, for a set of homogeneous elements 2 we have Tensor multiplication can be conveniently recast in the standard from by introducing the graded tensor product ⊛, defined in accordance with (A ⊛ B)(C ⊛ D) = AC ⊛ BD. Further clarifications about the notation can be found in appendix A.
1 Lindbladian flows can be alternatively understood in terms of 'quantum trajectories', i.e. an approach which uses a stochastic differential equation for an ensemble of pure quantum states evolving under an effective non-hermitian Hamiltonian [89][90][91]. 2 Homogeneous elements are linear operators on (C n|m ) ⊗N with a well-defined parity.

Exactly solvable nonequilibrium steady states
The algebraic construction of the solutions which is outlined below consists two steps. The general strategy in some sense reminds of solving a Poisson's equation. Namely, the first step is to identify a space of solutions for the bulk part which only accounts for the unitary part of the generator L 0 . Note that the entire space of bulk solutions is determined purely from the kinematic constraints, i.e. it is determined solely from the quantum symmetry algebra of the spin chain, irrespective of the representation labels. The second step is to impose the dissipative boundary conditions which (when a solution exists) uniquely fixes the physical state at hand. More specifically, this step amounts to chose suitable boundary auxiliary states and subsequently solve a non-linear system of boundary constraints in the space of the free representation parameters. Such a separation of bulk and boundary processes is indeed a characteristic feature of all exactly solvable classical and quantum boundary-driven lattice models.
The aim of this section is to break down the entire procedure into elementary steps and systematically discuss all the necessary ingredients to carry out the algebraic construction for the class of steady-state solutions of integrable quantum chains. The more difficult problem of identifying and classifying the relevant class of subalgebras is postponed to the next section, before finally presenting a few explicit results in Section 5.

Graded Yang-Baxter relation
This work is focused on a particular class of integrable lattice models which involve both bosonic and fermionic states. These models can be systematically derived from the so-called rational solutions to the graded Yang-Baxter relation. On a two-particle space C n|m ⊗ C n|m the latter takes the following form where z 1,2 are two arbitrary complex numbers usually referred to as the spectral parameters.
Here and subsequently we shall use the convention in which bold-faced symbols pertain to operators which act non-identically in the auxiliary space(s). Let us first explain the main objects. The graded R-matrix R n|m (z) acts as an intertwiner on the two-fold space C n|m ⊗ C n|m , i.e. expresses the equivalence of two distinct orderings of the tensor product of two L-operators. Matrices R n|m (z) are simply related to the graded permutation matrices P n|m , where matrix units E ab form the standard basis of linear operator in the fundamental module V n|m . The rational Yang-Baxter relation (8) can be formally understood as the defining relation of an infinite-dimensional associative algebra Y ≡ Y (gl(n|m)) known as the Yangian.
The L-operator from Eq. (8) is in this context interpreted as a Y-valued matrix on C n+m which admits the resolution The class of solutions to the nonequilibrium protocol considered in this work turn out to be related to certain degenerate representation of Y which are defined discussed in Section 4.  (8). From the particle scattering perspective, the relation imposes the equivalence of two apriori district ways of three consecutive pairwise elastic scatterings. Time direction for physical particles flows vertically and is shown by gray trajectories. The R-matrix R(z 1 − z 2 ) acts proportionally to a graded permutation on two fundamental physical particles in a su(n|m) symmetric integrable quantum chain. Lax operators L(z i − z 3 ) on the other hand govern scattering between physical particles with rapidities z i , i ∈ {1, 2}, and a fictitious particle carrying rapidity z 3 whose time-direction runs horizontally. Graded Yangians Y (gl(n|m)) are infinite-dimensional associative algebras for the coefficients of the operator components of the L-operator, with the Yang-Baxter equation on C n|m ⊗ C n|m taking the role of its defining relations.
The presence of non-trivial grading can be seen as a diagonal 'metric tensor' θ on twoparticle space C n+m ⊗ C n+m , θ ac,bd = (−1) ab δ ac δ bd , which allows for an alternative interpretation of Eq. (8) as a braided Yang-Baxter equation on non-graded vector space C n+m ⊗ C n+m , The graded permutation can be expressed in terms of the non-graded permutation P n+m on C n+m as P n|m = θP n+m . The central object of the algebraic Bethe Ansatz solution of integrable quantum models is the fundamental transfer matrix T n|m (z) operating on a N -particle physical space (C n|m ) ⊗N and satisfying the involution property (additional information can be found in appendix B) Commutativity of transfer matrices is ensured by the existence of the R-matrices R n|m , (z) which intertwines two fundamental auxiliary representations V n|m , i.e. it solves the corresponding (graded) Yang-Baxter relation. We need to emphasize however that Eq. (8) is written with the opposite identification of physical and auxiliary degrees of freedom with respect to the form which is most commonly used in the (algebraic) Bethe ansatz technique. The upshot is that our construction necessitates generic auxiliary models and not the conventional fundamental auxiliary representations. Indeed, the ordinary set of transfer matrices T n|m (z) and their fused counterparts which correspond to finite-dimensional auxiliary irreducible representations of gl(n|m) are unitary objects and in fact have no natural place in our application.
Differential Yang-Baxter relation. Taking the derivative of Eq. (8) with respect to z = z 1 − z 2 yields the differential Yang-Baxter relation (sometimes also called the Sutherland relation, cf. [74,81,94,95]), using the short-handed notation L ′ (z) ≡ ∂ z L(z). Equation (14) is simply a consequence of the fact that su(n|m)-symmetric Hamiltonian densities h n|m coincide with graded permutations P n|m over C n|m ⊗ C n|m , i.e.
For the so-called rational spin chains, relation (14) is in fact a simple corollary the zerocurvature property of the Lax connection. 3 What is more important is that the differential Yang-Baxter relation (16) is satisfied on a purely algebraic level, i.e. irrespective of a representations of the auxiliary components of the L-operator. A general solution to Eq. (15) is given by an operator L Λ n+m (z) acting on a product space of a local physical space and an arbitrary auxiliary representations, that is V ⊗ V + Λ n+m . Here V + Λ n+m denotes a generic irreducible highest-weight representation of gl(n|m) Lie superalgebra, characterized by a set of Dynkin labels Λ n+m (cf. appendix B). A key property of algebraic relation (14) is that it remains intact under fusion of auxiliary spaces. This readily makes it possible to extend it to composite (many-particle) auxiliary spaces, namely we may quite generally consider a multi-component Lax operators of the following form acting on V ⊗H aux , with H aux representing an arbitrary ℓ-component auxiliary product space characterized by a set of weight vectors Λ ≡ {Λ 1 n+m , . . . , Λ ℓ n+m } and a vector of complex parameters z ≡ {z 1 , . . . , z ℓ }. It is worthwhile emphasizing at this point that the tensor product in Eq. (16) is written with respect to auxiliary spaces V Λ n+m , and thus differs from the tensor product of two Lax operator from Eq. (8) which multiplies two copies of local physical (fundamental) spaces C n|m . The multi-component Lax operator L Λ (z) obeys an analogue of Eq. (14), where the z-derivative acting on L Λ n+m (z) should be replaced by the chain-rule derivation ∂ z ≡ ℓ i=1 ∂ z i on L Λ (z), as illustrated in Figure 4.

Amplitude factorization
We consider an N -site quantum system with the Hamiltonian and impose open boundary conditions. This class of models describes integrable quantum chains symmetric under su(n|m) Lie superalgebra [96,97] whose interactions take a simple form 4 We adopt the convention for summing over repeated indices throughout the text, unless stated otherwise.
In the boundary-driven setting, the Lindblad dissipator D gets naturally split into two independent incoherent processes assigned to the boundaries of the chain, i.e. D = D L + D R , where D L (D R ) operates only on the first (last) site of the chain. It is perhaps not too surprising that fixed-point solutions ρ ∞ to Eq. (4) for some bulk Hamiltonian H n|m with generic Lindblad boundary dissipators typically yields density matrices lacking any obvious structure. Remarkably however, there exist a classes of dissipative boundary conditions for which one may derive an exact algebraic expression for it. Before presenting the precise form of such integrable dissipative boundaries in section 3.3, we first wish to explain why localizing the dissipators to the chain boundaries plays a vital role in our construction and briefly comment on some important consequences. In simple terms, attaching the dissipators only to the boundaries manifestly ensures that the unitary part of the fixed-point condition (4) annihilates ρ ∞ , modulo some residual terms which stick at the boundary sites of the chain. This neat property motivates to use the algebra of (possibly non-local) commuting operators associated to the Hamiltonian H n|m as a trial space of operators for constructing an appropriate Ω-amplitude introduced in Eq. (18). In other words, by assuming that the steady-state solution of our problem has a well-defined local structure which is related to the symmetry algebra of the Hamiltonian. As explained below, the global symmetry gets broken only due to a mismatch in the boundary conditions, which is essentially the reason why the steady state is non-trivial, [H, ρ ∞ ] = 0.
We now proceed by employing the following amplitude factorization of the density operator ρ ∞ , Let us immediately stress that even though such a decomposition can be applied quite generally, it plays no fundamental role without imposing further restrictions on the amplitude operators Ω N . 5 . Indeed, the factorization property has been originally observed already in the seminal paper [73] where a non-pertrubative steady-state solution of the driven anisotropic, where it is referred to as the'reverse many-body Cholesky factorization'. However, in the class of solutions considered here, Ω N need not be a Cholesky factor of a steady state ρ ∞ , namely there is no requirement that Ω N takes a triangular form when expanded in the standard many-body computational basis of unit matrices spanning (C m+n ) ⊗N . Nonetheless, since the entire class of solutions which are presented below extends the simplest su(2) model to higher dimensional quantum spaces, we shall adopt the factorization property as a starting point of our presentation 6 . Following the above reasoning, the local symmetry of model is manifestly realized by introduce the following homogeneous fermionic matrix-product operator acting on an N -site quantum chain, with symbol ⊛ designating the graded tensor product and takes into account the presence of fermionic states. In the pictorial representation, the amplitude represents the lower leg in Figure 3. The key properties of the amplitude operator are: • Each tensor factor in Eq. (19) is assumed to be a gl(n|m)-invariant Lax operator parametrized by a continuous real parameter g (being the coupling strength parameter associated the Lindblad dissipator D). The L-operator acts (by definition) on a local physical space C n|m and an auxiliary Hilbert space which is at the moment left unspecified. Hence, the auxiliary component of the L-operator can be thought of as a generic representation of the underlying quantum algebra.
• We have introduced the boundary state |vac which will be subsequently referred to as the (auxiliary) vacuum. The vacuum state is determined by the choice integrable dissipative boundaries. In all the instances addressed in this work, |vac is simply an 'empty state', i.e. a product of highest-or lowest-weight vectors from the irreducible components which form a representation of the auxiliary algebra of the L-operator. This uniquely fixes the vacuum state once the representation labels (i.e. Dynkin labels and additional labels to specify the types of modules involved) associated to the L-operator from Eq. (19) are being specified. The role of the auxiliary vacua shall be more carefully explained in Section 5 we treat a few explicit instances.
Let us now return to the differential Yang-Baxter relation. Algebraic property (14) can be readily extended to the entire spin chain Hilbert space H ∼ = (C n+m ) ⊗N by simply expanding out the commutator [H n|m , Ω N ] and iteratively applying Eq. (14) at every pair of adjacent lattice sites. This results in a telescoping cancellation mechanism which on globally almost annihilates the unitary part of the evolution generated by L 0 , leaving behind only residual boundary terms which are an artefact of open boundary condition. This can be formally expressed in the form [74,98] which is can be viewed as the global version of the local condition (14) after contracting with the vacuum |vac at the end. We have written Ξ L,R to denote a pair of 'boundary defect operators' acting only in the boundary particle spaces (their explicit form is not of our interest) . Figure 3: Matrix-product representation of the non-equilibrium steady state ρ ∞ = Ω N Ω † N : the amplitude operator Ω N is represented by the degrees of freedom residing in the bottom row (shown in pink), while its conjugate transpose Ω † N corresponds to the upper row. In terms of an auxiliary scattering process, auxiliary particles are depicted by black lines and propagate in the horizontal direction. They can be viewed as fictitious particles composed of canonical bosons, fermions of complex (super)spins, emanating from the auxiliary vacuum on one end and getting absorbed by the same vacuum at the other end. Physical degrees of freedom (shown by gray vertical lines) are on the other hand associated to N fundamental particles of gl(n|m) Lie superalgebra. The off-shell steady-state density operator admits an interpretation as a vacuum contraction of a homogeneous two-row monodromy operator is a Lax operator which acts on a vertical rung.
Factorization property (18) indicates that the auxiliary Hilbert space H aux associated to the matrix-product representation of the steady-state density matrix ρ ∞ is a two-fold product of auxiliary subspaces which belong to mutually conjugate realizations of the underlying symmetry. Therefore, setting ℓ = 2 in Eq. (16) and writing shortly Λ m+n ≡ Λ, we arrive at the 'off-shell' representation 7 for ρ Λ (z), where is a two-row Lax operator of the form which is represented in Figure 3 by a vertical rung. Similarly, the boundary state |vac represents a factorizable state of two auxiliary vacua, |vac = |0 ⊗ |0 . The internal structure of the vacuum state |vac , which depends on the rank of symmetry algebra and the choice of integrable boundaries, will be detailed out in Section 5.

Boundary compatibility condition
Given the Hamiltonian H n|m and a set of boundary dissipators D, the fixed-point condition (4) imposes a certain type of bulk-boundary matching condition. It can be inferred from expression (20) that the fixed point condition L ρ ∞ = 0 admits a solution ρ ∞ if and only if there exist an Ω-amplitude (which amounts to find the L-operator and the vacuum state |vac ) for which the dissipator D exactly cancels out the right hand-side of Eq. (20). By plugging Figure 4: A schematic depiction of h n|m (L Λ (z) ⊛ L Λ (z)), representing the local action of the Hamiltonian H n|m on the density operator ρ ∞ = Ω N Ω † N (the coloring adopted from Figure 3.1, and spectral and representation parameters are suppressed for clarity). The process of brining the interaction h n|m across the horizontal legs generates terms which can be interpreted as a operator divergence condition for two-row Lax operators L Λ (z).
in a trial off-shell density operator ρ Λ (z) and demanding the on-shell condition one obtains a system of boundary algebraic equations for the undetermined representation parameters which depends also on the physical coupling parameters g of the reservoirs. The solution, when it exists, singles out a unique density operator ρ ∞ (g).
Combining Eq. (21) with a general solution of the bulk condition (16) results in two decoupled sets of boundary compatibility conditions, which can be cast in the compact form [82] vac| The boundary conditions of this form generically yields an overdetermined system of equations for the free parameters of the two-row Lax operator L Λ (z). Indeed, it is not difficult to confirm that in spite of integrability of the bulk interactions generic boundary dissipators do not lead to any solutions of Eqs. (23). In other words, for some general choice of boundary dissipators there exist no off-shell operator ρ Λ (z) which would satisfy the fixed-point condition of Eq. (4). Of course this should not be surprising at all since typical dissipation processes result in a 'non-integrable' Liouvillian dynamics in which a naïve separation of bulk and boundary parts cannot be justified. Needless to say that in such a case there exists no obvious explicit representation of the steady states either. It is therefore quite remarkable that integrable lattice models with su(n|m)-symmetric interactions h n|m do allow for certain elementary (so-called integrable) boundary dissipators which lead to non-trivial solutions to boundary equations (23).

Integrable dissipative boundaries
We consider a pair of dissipative boundary processes which involves any (but arbitrary) pair of states from the local Hilbert space C n|m . Denoting them by |α and |β , we posit the jump operators of the form 8 parametrized by a single reservoir coupling parameter g. Since Lindblad dissipators which enter in Eq. (24) operate non-trivially only on the boundary sites of the chain, the jump Figure 5: Graphical interpretation of the boundary compatibility condition as given by equation (23) displayed for the right boundary at lattice site N . The left-hand side shows schematically the action of the dissipator D on the L-operator decomposed into three terms which the define the action of the jump operator A N . The termination point of the horizontal arrow signifies the contraction with the right auxiliary vacuum state. Note that the boundary condition has to be satisfied for all values of physical indices.
operator from Eq. (24) can be interpreted as a source and drain associated to U (1) particle currents.
In models with multiple states per site such as su(n|m) chains considered here, diagonal 'density operators' E aa i obey the following local continuity equations where j ab i denote partial currents between two level |a and |b locally at lattice site i. Total current densities between two adjacent lattice sites are then obtained by summing over all partial currents, that is j a i,i+1 ≡ n b=1 j ab i,i+1 , and fulfil Integrable su(n|m) symmetric Hamiltonians H n|m conserve each total particle number N a = i E aa i independently, i.e. [H n|m , N a ] = 0. The addition of dissipation however destroys the conservation of N a if a ∈ {α, β}.
To better examine this situation, we notice that a particular choice of dissipative boundary condition as given by Eq. (24) allows to decompose Liouvillian dynamics into invariant subspaces, where orthogonal Hilbert subspaces H ν are defined via N γ H ν = ν γ H ν , with eigenvalues ν γ ∈ {0, 1, . . . , N } for γ ∈ {1, . . . , n + m}. Accordingly, we introduce endomorphisms O ν = End(H ν ), i.e. linear spaces of operators operating on H ν . This means that states from H ν have well-defined values of all particle number operators N γ . When rank(g) > 1, there exist at least one number operator N γ such that This is an example of the so-called 'strong Liouvillian symmetry' [99]. In fact, all N γ correspond to strong symmetries, with the exception of the two distinguished indices which belong to a pair of levels affected by the boundary dissipation, that is γ ∈ {α, β}. This immediately implies degeneracy 9 of the steady states (cf. [82]) and vanishing current expectation values j γ ∞ = 0. Thus only current densities j α ∞ and j β ∞ can take non-vanishing steady-state expectation values. When dealing with degenerate null spaces of the generator L, the steady state operator ρ ∞ naturally decomposes in terms of independent fixed-point components ρ with P µ denoting orthogonal projectors onto subspaces O µ . Because each invariant component ρ µ ∞ satisfies Lρ µ ∞ = 0 for all allowed values of particle numbers µ, they may be combined in any convex-linear combination with c µ representing a (n + m − 2)-component vector with non-negative components. The steady-state operator ρ ∞ as defined by Eq. (30) can be thus regarded as a grand canonical nonequilibrium ensemble with coefficients c µ , which in analogy to grand canonical ensembles have the role of particle chemical potentials. Notice that Eq. (30) can be conveniently cast in the form of a matrix-product operator along the lines of ref. [82] for the su(3) chain.
A heuristic analysing of 'integrable boundaries' specified Eq. (24), e.g. by computing exact solutions for quantum chain of small length, reveals that the fixed point is a non-trivial current-carrying steady states of particularly simple structure. In the remainder of the paper we demonstrate that the steady-state solutions exhibit a particular algebraic representation which directly link to fundamental objects of quantum algebras. It is worthwhile stressing nonetheless that, despite the simplicity of our effective reservoirs, the entire spectrum of L -typically referred to as the Liouville decay modes -remains highly complex and lack any obvious structure. With that said, it is therefore only the fixed-point solutions ρ ∞ of Eq. (4) which admit an exact description. It is also instructive to remark here that even the integrable steady state density operators themselves do not enjoy the full quantum group symmetry of the underlying theory. Indeed, as a consequence of the foliation (27) of the Lindbladian flow, the global residual symmetry of ρ ∞ is merely U (1) ⊗n+m−2 . However, as subsequently demonstrated, the local symmetry of the Ω-amplitude is much larger. The symmetry content of the steady state solution will be carefully examined in the next sections. Particularly, the local symmetry of the Ω-amplitudes will become apparent on the basis of previously discussed Lax representation (see also Section 6 for additional remarks).

Graded Yangians
Yangians are certain infinite-dimensional quadratic associative algebras which belong to a class of (quasi-triangular) Hopf algebras, widely referred to as quantum groups. Yangians can be defined in various equivalent ways [101][102][103]. Here we employ the 'FRT realization' [104] (also known as the 'RTT realization'), in which Yang-Baxter equation (8) takes a role of the defining relation.
We specialize the discussion to Yangians Y ≡ Y (g) of Lie superalgebras g = gl(n|m) [105,106]. Recall that the signature n|m indicates that the local Hilbert space consists of n bosonic and m fermionic states. Generators of Y are given as the operator-valued coefficients of the Lax operator L(z) expanded as a formal Laurent series By imposing Yang-Baxter equation (8) as the defining relation, we obtain an infinite set of quadratic algebraic conditions The level-0 generators L ab (0) are scalars belonging to the center of Y. In the scope of our application, we are only be interested in the class of fundamental rational solutions of Eq. (9) which are of degree one in the spectral parameter z, 10 This choice represents, in mathematical terms, an evaluation homomorphism from the Yangian to the universal enveloping algebra of g, Y(g) → U (g). With this restriction, representations of Y are in one-to-one correspondence with representations of the classical Lie (super)algebra g.
Automorphisms. It is instructive to shorty discuss the gauge freedom due to automorphisms of Y, i.e. transformations which preserve the algebra (32) (cf. [86]). These comprise of (i) rescaling L(z) with an arbitrary complex-valued scalar function f (z), (ii) shifting the spectral parameter z → z + z ′ and (iii) applying a (n + m)-dimensional GL(n|m) gauge transformations which acts in C n|m and is given by two arbitrary invertible matrices G L and G R , In addition, there exist anti-automorphisms of Y, i.e. transformations which only preserve the defining relations (8) up to exchanging the order of tensor factors. Examples of these are transposition of the matrix space L(z) → L t (z), and reflection in the spectral plane L(z) → L(−z). Any composition of two anti-automorphisms is again an automorphism.
Rank-degenerate realizations. The list of transformations given above nonetheless do not exhaust all possibilities of realizing Y. As pointed out in [86,87], equation (8) admits a class of 'degenerate solutions' provided one relaxes the requirement L (0) = 1. This is a viable choice because the level-0 generators L ab (0) are central and can therefore take arbitrary (possibly vanishing) values. We may thus quite generally prescribe modulo equivalent choices which correspond to permutations of 0s and 1s. Such a restriction obviously induces another block structure 11 on C n|m under which the generators of Y split as The ranges of ordinary (undotted) and dotted indices are where p (q) denotes the number of bosonic (fermionic) states in the index set I. Similarly, we shall denote byṗ (q) is the number of bosonic (fermionic) states contained in the complementary set I. The defining relations of the resulting 'hybrid algebra' A I n,m are readily obtained by plugging Eq. (35) in Eq. (32), and select the level-0 generators in accordance with Eq. (34). in accordance with the prescription of Eq. (34), i.e. L ab (0) = δ aI δ bI . Since the generators Dȧ˙b are central, it is convenient to pick a gauge by setting Dȧ˙b = δȧ˙b. The remaining non-trivial commutation relations read

Oscillator realizations
Commutation relations (37) have been derived in [86,87], where the authors provide a realization in terms of gl(p|q) 'super spin' generators J ab (for a, b ∈ I, and with p + q = |I|), and additional |I| · |I| canonical bosonic or fermionic oscillators which obey graded canonical commutation relations where a generator ξ aḃ should be understood as a creation operator of a bosonic (fermionic) oscillator if p(ḃ) = 0 (p(ḃ) = 1), for a ∈ I andḃ ∈ I. The oscillator part of the algebra A I n,m , denoted by osc(p+q|p+q), is associated with a multi-component Fock space B ⊗(p+ṗ) ⊗F ⊗(q+q) , where each factor B (F) belongs to an irreducible bosonic (fermionic) Fock space. In terms of these 'super spins' and 'super oscillators', the level-1 generators L ab (1) take the canonical form 11 We follow the notation of [86,87] and employ the two-index labelling of the Yangian generators. where

Partonic Lax operators
For a Lie superalgebra g of rank(g) = n + m, there are in total 2 n+m distinct types of hybridtype subalgebras A I n,m . The latter are in a bijective correspondence with all possible choices of set I. Notice that this counting we are excluding the various possibilities of choosing the grading. By additionally excluding permutation equivalent choices, we eventually deal with finite-dimensional Lie superalgebras of the type gl(p|q) ⊗ osc(p +ṗ|q +q).
The simplest rank-degenerate solutions of the graded Yang-Baxter algebra (8) belongs to the single-indexed sets |I| = 1 and were dubbed in [85] as the partonic solutions. These consist solely from n + m − 1 oscillators arranged in a distinctive cross-shaped form, with Here the integers a ∈ {1, . . . , n + m} in the subscript of L {a} (z) are being used to indicate the location of the single non-vanishing level-0 generator, cf. Eq. (34). As shown in appendix B, all Lax operators associated to A |I|≥2 m,n can be systematically generated from the partonic solutions which carry A

Exact steady states for integrable quantum spin chains
In this section we finally present a few explicit examples for the steady-state solutions of the boundary-driven su(n|m) quantum chains, subjected to integrable dissipative boundaries given by (24). As the first step, we account for kinematic constraints and construct off-shell density operators which take the universal form of equation (21). Subsequently, the goal is to find an appropriate internal structure of the Lax operator L Λ and the auxiliary vacuum state |vac which fulfil the requirements of the boundary equations (23).
Notice first that there exist in total 2× n+m 2 ways of assigning the dissipators of Eq. (24), representing all the possibilities of selecting a pair of target levels |α and |β . The extra factor of 2 reflects a possibility of exchanging |α with |β which results in a state of opposite chirality, i.e. reverses directions of particle currents. It turns out that every choice of |α and |β leads to a solution to Eq. (4) which is uniquely characterized by specifying a representation and the corresponding labels for the auxiliary algebra A I n,m of the Lax operator L Λ (z).

Fundamental integrable spin models
The significance of partonic Lax operators and the structure of the steady-state solutions is perhaps best illustrated by explicitly working out a few simplest examples. To this end, we first consider the non-graded interactions, and initially examine the most studied case of the su(2) spin chain (the isotropic Heisenberg spin-1/2 model), with interaction density h 2|0 = 2 a,b=1 E ab ⊗ E ba . Let us remark that this particular instance has been considered initially in the seminal paper [73] where the solutions was found with a somewhat different approach, and afterwards re-obtained in a more compact and symmetric form in [74]. The derivation from Yang-Baxter algebra has been presented in [81]. Nevertheless, to uncover the connection with partonic Lax operators and embed this solution in a unified theoretic framework, we shall reproduce it below once again.
In the su(2) spin chain, the local building block of the Ω-amplitude is given by a twoparametric Lax operator L − j (z) acting on V ⊗ V − j , whose auxiliary space V − j represents a lowest-weight sl(2) module spanned by an infinite tower of states {|k } ∞ k=0 . We adopt the sl(2) spin generators obeying algebraic relations whose action on V − j is prescribed by State |0 has the lowest weight, J − |0 = 0, and will be referred to as the vacuum. By recalling that all the solutions factorize in accordance with Eq. (18), the off-shell Lax operator L Λ (z) which defines the steady-state solution ρ ∞ is a product of two copies of auxiliary representations, cf. Eq. (22). This factorization makes is possible to express the final solution by only specifying a pair of complex parameters: a sl(2) weight which is interpreted as a complex spin j, and a complex-valued spectral parameter z. Specifically the two-parameteric off-shell Lax operator which we denote by L j (z) is represented in the following compact form The notation used is as follows. Integer indices in superscript brackets are used to assign an operator L j (z) into the corresponding tensor factor in the multi-component auxiliary space.
In addition, in the superscripts we also employed extra parity signatures which are required to correctly specify the type of the sl(2) module. Namely, while the first auxiliary copy is realized in the lowest-weight module as prescribed by Eq. (45), the second factor in Eq. (46) must be associated with the highest-weight realization of sl(2) algebra V + j . The highestweight Lax operator L + j (z) can be readily obtained from L − j (z) by applying the spin-reversal transformation The highest weight state |0 from V + j is distinguished by J + |0 = 0.
Plugging the off-shell form of Eq. (46) into the boundary conditions (23) yields a system of polynomial equations with a unique solution Notice that the auxiliary vacuum state takes the product form, |vac = |0 ⊗ |0 , and is determined by the internal structure of L j (z). In order to reverse the direction of driving we may simply exchange the target states |α ↔ |β . This amounts to interchange the factors in Eq. (46). Before proceeding with other examples, let us stress again that a proper identification of the internal structure of the two-leg Lax operator L j (z) is crucial. Once the convention for labelling the irreducible sl(2) Verma modules is being fixed, there exist only one correct assignment of irreducible spaces (incorrect assignments produce a system of boundary equations which admits no solution).
Asymmetric driving. Let us mention a simple trick which enables to generalize the solutions to the case of unequal reservoir coupling constants. Considering as an example the su(2) spin chain, we may impose an asymmetric pair of Lindblad jump operators of the form This choice yields an extended class of solutions which is connected to the special case of equal couplings by a diagonal tilting transformation -a one-parameter automorphism of Yby applying on every local spin space C 2 . The solution to the boundary compatibility conditions is then given by Models of higher-rank symmetry. The simplest higher-rank model is the su(3) spin chain (with interaction h 3 ), often called in the literature as the Lai-Sutherland model [107,108]. Solutions to the fixed-point condition (4) have been originally identified and parametrized in [82] and now represent a degenerate manifold of steady states. As discussed earlier in Section 3.4, degeneracy of the null space of L is a consequence of the conservation of the number operator N γ associated to a distinguished noise-protected state |γ (which depends on I).
The Lax operator for the Ω-amplitude now operates on a space of three auxiliary particles, a non-compact sl(2) spin and two species of canonical bosons. Bosonic particles obey canonical oscillator algebra, and live in the canonical Fock space spanned by a tower of states {|k } ∞ k=0 . Similarly as in the case of sl(2) spins, one also has to distinguished two distinct realizations of bosonic Fock spaces B ± , related to each other by an algebra automorphism which is interpreted as the particle-hole conjugation. By assigning the dissipation to states I = {1, 2}, the off-shell Ω-amplitude is constructed from the Lax operator L {1,2} (z) which carries a representation of algebra A I Now is suffices to repeating the logic used before on the su(2) case, and define a factorized off-shell Lax operator L ω j (z) for the steady-state solution ρ ∞ in the form This time we equipped each tensor copy with an additional label ω which, as argued earlier, is needed to supply the information about the types of sl(2) and Fock modules. After determining the right ω, the coupling constant g is linked to the free representations parameters z, j through the boundary equations (23) with the solution The delimiter in ω was used to explicitly distinguish the sl(2) module V ± j (on the left) from the signatures belonging to the product of Fock spaces B ± (on the right, in the ascending order). Specifically, the above instance requires a lowest-weight type sl(2) representation and to assigned B − (B + ) to the first (second) bosonic oscillator.
Before heading on to the more involved examples of fermionic models, let us spent a few more words on the non-trivial structure of the vacuum |vac and, in particular, to inequivalent roles of the highest and lowest type of (auxiliary) representations. As said earlier, in order to construct an off-shell Ω-amplitude it is first required to infer an appropriate 'internal structure' for the auxiliary space of L Λ n+m (z). Only then it is possible to proceed by solving the corresponding finite system of polynomial equations (23). The upshot here is that the module type labels ω are essential to assign ρ ∞ the appropriate chiral structure. For instance, an incorrect assignment of the auxiliary bosons which in expression (56), e.g. by imposing two identical representations B + , would violate the boundary compatibility conditions. Finally, one can easily verify that the Lax operator L {1,2} (z) is in fact equivalent to the Lax operator found previously in [82]. 12

Fermionic models
In this section we generalize the above construction to the steady state solutions which pertain to graded su(n|m) chains, representing the simplest class of interacting integrable models with fermionic degrees of freedom. We retain the dissipative boundaries given by Eq. (24). 13 Besides bosonic oscillators, the auxiliary particle spaces will now also involve canonical fermions from two-dimensional spaces F.
The defining gl(1|1) representation is spanned by two basis states |0 and |1 . The 'highest weight' type representation, denoted by F + , is prescribed by where the generators obey canonical anticommutation relations Similarly, the 'lowest weight' representation F − is obtained from F + by virtue of the particlehole mapping Free fermions. Arguably the simplest fermionic integrable system is a tight-binding model of non-interacting spinless fermions (with a homogeneous chemical potential) whose interaction is invariant under gl(1|1) Lie superalgebra and in terms of canonical fermions reads 14 In spite of its simplicity, it is remarkable find that the corresponding integrable steady states involve auxiliary spaces with belong to non-canonical gl(1|1) representations. The fermionized integrable reservoirs provided by Eq. (24) are interpreted as an inflow (outflow) of spinless fermions at the left (right) boundary with rate g, To find the unique solution to the fixed-point condition (4), we follows the procedure from the previous section and first consider the following two partonic Lax elements, By merging them together using the fusion rule (see appendix B for details) we have The outcome is a two-parameteric Lax operator L λ (z) whose auxiliary space is identified with an indecomposable non-unitary representation denoted here by V λ . The latter can be realized in terms of canonical generators (60) as 13 A comment on the Jordan-Wigner transformation: When expressed in terms of the fermionic generators, the dissipator attached to the first lattice site differs from its non-graded counterpart by a (non-local) Jordan-Wigner 'string operator' W , indicating that fermionization of the boundary-driven spin chain maps to a model of non-local dissipation. This discrepancy between the two formulations which is due to the presence of W is however immaterial as far as only the steady states are of our interest, the reason being that W commutes with both the steady state ρ∞ and the total Hamiltonian H n|m .
14 The interaction h 1|1 can also be expressed in terms of fundamental su(2) generators. Up to boundary terms this yields the XX model in a homogeneous external field, that is where the complex representation parameter λ is the central charge. 15 To distinguish between the highest-and lowest-weight types of representations we shall use an extra superscript label, using the convention that V ± λ is associated with the Fock space F ± , i.e. (L ± λ ) 12 |0 = 0. In close analogy to the su(2) case, the local unit of the amplitude operator Ω N is now built from L-operator L − λ (z). The undetermined representation parameters are finally obtained from the boundary conditions (24), yielding a unique solution The factorization property (18) implies that ρ ∞ (g) is constructed from a two-component Lax operator L λ (z) which explicitly reads The driving may be reversed by first applying the particle-hole transformation on the physical fermions (see Eq. (63)), exchanging the order of factors in Eq. (68), and ultimately setting The auxiliary vacuum state |vac remains the product of the lowest-and highest-weight state |0 from the fermionic Fock modules F ± . We find it instructive to remark that the model of free fermions takes a special place among the gl(n|m) quantum chains which plays nicely with the fact that the model is compatible with a larger set of integrable boundary dissipators. It is rather remarkable however that such an extended set of solutions still admits the Lax representation, albeit the latter does no longer exhibit the usual additive form. We are not sure whether this enlarged set of steady-state solutions is still related to representation theory of Yangians. Further details and explicit results are presented in appendix C.
Example: SUSY t-J model. Integrable spin chains whose interactions coincide with graded permutations have been initially considered in [96,97]. A prominent (and generic) example is the su(1|2)-symmetric integrable spin chain which is mappable to the t-J model at the 'supersymmetric point' [109] (2t = J). The spectral problem of the model has been solved with Bethe Ansatz techniques in [110][111][112][113][114].
The grading can be distributed in various ways. 16 We shall adopt |0| = 0, and regard the empty state |0 as the highest-weight state (vacuum) in the physical Hilbert space. Then, we may consider one of the following three options, depicted by the corresponding Kac-Dynkin diagrams. 17 It is important to remark here that the choice of grading is entirely independent from the set I = {α, β} which specifies a pair of states subjected to the dissipators. Let us first set the grading to |1| = 0 and |2| = |3| = 1, which corresponds to diagram (1). Incoherent conversion processes induced by the dissipators can be described by any of the following three possibilities: Options (a) and (b) represent fermionic driving and physically corresponds to Markovian transitions between two states of opposite parities, namely a spin-carrying electron and the unoccupied state |0 . Option (c) is different, and affects the bosonic sector (i.e. the su(2) doublet) by triggering incoherent spin flips. Let us first address option (a), corresponding to assigning the following pair of Lindblad jump operators This instance pertains to I = {1, 2}, I = {3}, with p = q = 1 andq = 1, which defines the auxiliary algebra A where the generators J ab are associated with the non-unitary gl(1|1) representation V λ , whereas the super oscillators are identified with bosonic and fermionic canonical oscillators in accordance with the rule The solution to Eq. (23) is then given in the form where L ω j (z) now implements the auxiliary algebra A {1,2} 1,2 = gl(1|1) ⊗ osc(1|1), and the representation parameters take the values 16 From the algebraic point of view, distinct inequivalent gradings indicate that Lie superalgebras do not admit unique simple roots. All distinct possibilities are however related under certain boson-fermion duality transformations. In the context of Bethe Ansatz these correspond to inequivalent Bethe vacua and different ways of proceeding to higher levels in the nesting scheme (see e.g. [114][115][116]). 17 By convention we draw an open circle if two adjacent states are of the same parity and a crossed circle when their parities differ (assuming |0| = 0), while moving from left to right.
In particular, the signature labels ω (where the bar in Eq. (76) stands for flipping the signs) indicate that the auxiliary algebra A The same procedure can be repeated for the case of bosonic driving (c), where the jump operators acts as The auxiliary algebra A {2,3} 1,2 = sl(2) ⊗ osc(0|2) now consists of a non-compact sl(2) module V j (with the spin generators denoted by J a ) and a pair of fermionic Fock spaces F ⊗ F, In order to fulfil the boundary constraints, the auxiliary algebra of L ω j (z) should consist of the product V − j ⊗ F − ⊗ F + . Finally, ρ ∞ is cast in the universal form Eq. (21), where now

Vacuum Q-operators
Given that all the solutions are directly related to a particular type of solutions of the graded Yang-Baxter equation (8), it is quite remarkable (and perhaps surprising) that a one-parametric family of density matrices ρ ∞ (g) do not commute for different values of g. On the opposite, an explicit construction of ρ ∞ (g) indicates that (reported first in [117]) One may still wish argue that due to amplitude factorization (18) the nonequilibrium density operators ρ ∞ (g) are not the most 'fundamental' physical objects. Indeed, amplitudes operators Ω N (g) themselves do commute for different values of couplings, This shows that Ω N (g) can be regarded as a family of vacuum highest-weight transfer matrices. However, while the steady states ρ ∞ (g) are diagonalizable objects which encode physical properties of the system, their Ω-amplitudes exhibit a non-trivial Jordan structure. 18 Below we examine this unusual behaviour in more detail and relate it to the vacuum Q-operators.
Baxter's Q-operators. Before introducing the notion of vacuum Q-operators, let us first give some comments on the connection between the conventional Q-operators and Lax operators L I (z) introduced in Section 4. The concept of a Q-operator was originally introduced in Baxter's seminal paper on the 8-vertex model, where it was used as a device to diagonalize the transfer matrix of the problem by solving a suitable second-order difference relation [118], and later revived in the context of Potts model [119] and integrable structure of CFTs [120,121]. For clarity we focus the subsequent discussion entirely on the homogeneous su(2) spin chain, providing only a condensed summary of the main ingredients. For a more comprehensive and pedagogical exposition we refer the reader to [85].
Baxter's TQ-relation is a functional relation for the fundamental transfer operator T of the form 19 The pair of Baxter Q-operators Q ± (z) represents two independent operator solutions to the functional equation (83) and enjoy the involution property Eigenvalues of Q ± (z), denoted by Q ± (z), are (up to a twist-dependent phase which is omitted for brevity) polynomials of the form Their zeros z k (z k ) coincide with the Bethe (dual) roots, and are solutions to the celebrated Bethe quantization condition Polynomiality of eigenvalues of T (z) and Q ± (z) ensure that the TQ-relation (83) is equivalent to Bethe equations (86). Operators Q ± can be conveniently cast as auxiliary traces over quantum monodromies obtained by the lattice path integration of partonic Lax operators. Specifically, in the su(2) case we have Here we have made identifications L {1} (z) ≡ L + (z), L {2} (z) ≡ L − (z), and the trace is with respect to the auxiliary Fock space F. An analogous construction applies to integrable theories based on higher-rank algebras [120][121][122][123] where a complete set of Q-operators is associated to Lax operators L I (z) introduced in Section 4. In the language of Bethe ansatz, this means that eigenvalues of all Q-operators belonging to rational solutions of Yang-Baxter algebra (cf. Eq. (8)) are polynomials whose roots coincide with Bethe roots belonging to different nesting levels. An explicit construction of the full hierarchy of Q-operators for gl(n|m) spin chains can be found in [85][86][87][88] (and in [124,125], using a different approach). The outcome of this procedure is a set of 2 n+m distinct Q-operators which can be arranged on vertices of a hypercubic lattice [126]. In this context, partonic Lax operators L {a} (z) are associated to the distingusihed set of n + m elementary Q-operators Q {a} (z) which can be used to solve the spectral problem by explicit integrating of an auxiliary linear problem [122,123]. An important consequence of this is that eigenvalues of fused transfer matrices which obey the T-system functional identities [127,128] decompose in terms of the elementary Q-functions.
Vacuum Q-operators. Since in open quantum spin chains translational symmetry is manifestly absent, taking (super) traces over auxiliary spaces is no longer a priori justified. As originally noticed in [72], one may instead consider as a meaningful replacement projections onto the highest (or lowest) weight states of auxiliary spaces 20 . To this end we now define the following set of 'vacuum Q-operators' The previous analysis of the steady-state solutions for gl(n|m) spin chains with integrable dissipative boundaries given by Eq. (24) shows that all Ω-amplitudes can indeed by identified with vacuum Q-operators. More specifically, Ω-amplitudes which enter in our nonequilibrium setting always correspond to 'mesonic' Lax operators L I (z) with |I| = 2. We wish to elaborate on a subtle (but important) point in regard to the auxiliary algebra of L I (z) and the structure of the vacuum state |vac . The fact that L I (z) carry (besides Dynkin labels) the information about the types of irreducible components which implement the auxiliary algebra becomes crucial here. For instance, already in the simplest case of the su(2) chain, we had to define and operate with two inequivalent types of vacuum Q-operators (denoted by Q vac,± {a} (z), with a = 1, 2) carrying either of inequivalent bosonic Fock spaces B ± . The explicit structure of the fusion relation for partonic operators L {a} (z) (see appendix B) brings us to the conclusion that the vacuum Q-operators with equal auxiliary modules are still in involution and in turn implies that the same property also holds for the corresponding Ω-amplitudes (as given by Eq. (82)). Conversely, the objects which involve inequivalent auxiliary spaces do not commute, Q vac,± {a} (z), Q vac,∓ {a ′ } (z ′ ) = 0, ∀z, z ′ ∈ C, and a, a ′ ∈ {1, 2}.
Since the steady-state density operators ρ ∞ (g) always consist of two fused mesonic vacuum Q-operators of the opposite type (which is a corollary of property (18)), by virtue of Eq. (90) they do not inherit the involution property (82) from their amplitude operators Ω N (g). 21 It remains an interesting open problem to devise a suitable generalization of the Algebraic Bethe Ansatz procedure to diagonalize ρ ∞ [117].

Conclusion and outlook
In this work we introduced a unifying algebraic description for exact nonequilibrium steady states which belong to an important class of integrable quantum lattice models. We presented an explicit construction of density matrices which appear as non-trivial stationary solutions to a non-unitary relaxation process in which a system is coupled to effective Markovian particle reservoirs attached at its boundaries. We employed a simple set of incoherent particle source 20 In a more general setting, when particle source and drain terms are rotated with respect to the z-axis, the highest-or lowest-weight vacua get replaced by spin-coherent states [129]. 21 It is instructive to remark that tensor products of irreps of mixed types do not admit a resolution in terms of a (finite or infinite) discrete sum over extremal-weight irreps, in contrast to ubiquitous decomposition of tensor products of finite dimensional irreps (or products of extremal-weight irreps of the same type). and drain reservoirs which naturally generalized those used previously in refs. [73,74,76,80,81]. We have shown that such reservoirs partially preserve the integrable structure of the bulk Hamiltonian and permit to obtain analytic closed-form steady-state density operators in a systematic way.
The solutions were presented in the universal form of a homogeneous fermionic matrixproduct operators, and shown to decompose in terms of the vacuum analogues of Baxter's Q-operators. Such a factorization property reflects the chiral structure of the states and also allows to reverse directions of particle currents with aid of suitable particle-hole transformations. The basic building blocks of our construction are the so-called partonic Lax operators which stem from certain degenerate representations of graded Yangians, identified recently in [85][86][87][88]. These rather unconventional algebraic structures admit a canonical realization in terms spins and oscillators. In the context of our application, these appeared as the auxiliary degrees of freedom in the matrix-product operator representation for the steady-state solutions.
The absence of translational symmetry in open quantum chains is profound importance and requires to replace the usual auxiliary traces by the projectors onto highest-or lowestweight auxiliary vacua. The internal algebraic structure of the amplitude operators depends crucially on the parities assigned to particles which experience dissipation. In the case of equal parities (bosonic driving), the amplitude operators always involve a single auxiliary noncompact sl(2) spin, whereas the opposite parities (fermionic driving) require a non-unitary irreducible gl(1|1) representation which are two dimensional. The residual auxiliary degrees of freedom pertain to a finite number of canonical (bosonic or fermionic) oscillators which remain intact upon varying the coupling parameters of the reservoirs. The universal structure of the steady-state solutions signifies that it is the non-unitary part of the auxiliary algebra which ultimately controls their qualitative characteristics: on one end, the presence of sl(2) sectors leads to a universal anomalous (i.e. non-diffusive) j ∼ O(N −2 ) decay of longitudinal currents and cosine-shaped density profiles as already found in [73,76,80,82]. Fermionic driving is on the other hand characterized by gl(1|1) subspaces and triggers ballistic transport with non-decaying currents j ∼ O(N 0 ) and flat density profiles [98]. The solutions at hand can therefore be perceived two particular nonequilibrium universality classes.
The distinguished feature of integrable steady states addressed in this work are the nonunitary representations of Lie (super)algebras . This contrast the conventional approaches to quantum integrable systems which are primarily based on unitary representations and directly relate to physical excitations in the spectrum (described by the formalism of Thermodynamic Bethe Ansatz [130][131][132]). Physical significance of non-unitary representations in integrable theories is on the other hand far less understood and has not been much explored in the literature, although a few prominent examples are worth mentioning. Most notably, the logarithmic CFTs are based on (non-unitary) reducible indecomposable representations of Virasoro algebra [133,134] and are known to capture various phenomena in statistical physics ranging from critical dense polymers [135], symplectic fermions [136,137], critical percolation [138][139][140] to Gaussian disordered systems [141,142]. It is perhaps instructive to add that non-compact spin chains are also found in the hadron scatting in QCD, which is in the Regge regime governed by the s = 0 non-compact isotropic Heisenberg magnet [143,144].
The role of non-unitarity in the present nonequilibrium setting is however different as it is not (at least directly) attributed to physical degrees of freedom, but instead enters on the level of fictitious particles assigned to auxiliary spaces in a matrix-product representation of nonequilibrium steady states. Nevertheless, it has been found that non-unitary represen-tations can sometimes be linked to certain hidden conservation laws which turn out to be responsible for anomalous quantum spin transport (singular DC conductivity) in the linear regime [72,77,145,146].
In the conclusion we wish to highlight a few unresolved aspects of the problem which in our opinion deserve to be better explored and understood. In order to further extend the range of applicability of the present approach, it is of paramount importance to obtain better theoretical understanding of the integrability-preserving dissipative boundaries. In particular, whether there exist a connection between quantum integrability and a special type of Lindblad reservoirs employed here remains unanswered at the moment. Another intriguing open question is to find a field-theoretic version of the Lindbladian evolution which would qualitatively reproduce the scaling regime of integrable quantum lattices (cf. [80]). It is moreover difficult to overlook several discernible similarities with the Caldeira-Leggett approach of modelling a dissipative environment with a boundary-localized friction term [147,148], which has been applied to sine-Gordon theory with an integrable boundary perturbation [149]. In particular, (i) the boundary current is given by the vacuum eigenvalues of CFT analogues of Q-operators, (ii) the reservoir parameters are linked to purely imaginary values of highest weights, and (iii) the particle current is expressed directly in terms of the nonequilibrium partition function Z = Tr ̺ ∞ , which is otherwise common to both the asymmetric classical exclusion processes [60] and their quantum counterparts [80] considered here. In our opinion, these curiosities deserve to be further explored in future studies.

Acknowledgements
The author thanks P. Claeys, V. Popkov, E. Quinn, and especially T. Prosen for valuable and stimulating discussions and/or providing comments on the manuscript.
Funding information. The author acknowledges support by VENI grant by the Netherlands Organisation for Scientific Research (NWO).

A Graded vector spaces and Lie superalgebras
A graded vectors space is a complex vector space C n|m , spanned by basis states {|a } n+m a=1 , which is endowed with a Z 2 map, referred to as the grading: We subsequently adopt (with no loss of generality) the distinguished grading, |a| = 0 a ∈ {1, 2, . . . , n} 1 a ∈ {n + 1, . . . , n + m} .
This assignment induces a block decomposition on End(C n|m ), being the space of matrices acting on C n|m . Specifically, End(C n|m ) = V 0 ⊕ V 1 , where components V 0 (dim V 0 = n) and V 1 (dim V 1 = m) represent bosonic (even) and fermionic (odd) parts, respectively. The subspaces V 0 and V 1 are referred to as the homogeneous components of End(C n|m ). Notice that while V 0 is a subalgebra, the odd part V 1 is not. A vector space End(C n|m ) also constitutes gl(n|m) Lie superalgebra. In particular, any element A admits a block form where sub-matrices A 00 , A 11 , A 01 and A 10 are of dimensions n × n, m × m, n × m and m × n, respectively. The bosonic part decomposes in terms of bosonic subalgebras gl(n|m) 0 ∼ = gl(n) ⊕ gl(m) and corresponds to A 01 = A 10 ≡ 0, whereas the fermionic (odd) part gl(n|m) 1 pertains to elements with A 00 = A 11 ≡ 0. Let E ab denote matrix units, i.e. matrices with the only non-zero element being 1 in the a-th row and b-th column. Basis element E ab are assigned a Z 2 -parity according to p(E ab ) ≡ |E ab | → {0, 1}, with the prescription Element A is of even (odd) parity when non-vanishing blocks A ab are of equal (opposite) parity |a| + |b| = 0 (|a| + |b| = 1). For instance, R-matrices R n|m are always even elements.
Matrix units E ab form a basis of the fundamental representation of gl(n|m) algebra denoted by V n|m . The graded Lie bracket is prescribed by Here and below we simplified the notation by writing (−1) a instead of (−1) |a| . Graded vector spaces and Lie superalgebras are a naturally extended over N -fold product spaces. Product spaces inherit the parity according to the prescription A linear operator A on (C n|m ) ⊗N is called a homogeneous element of parity |A| if it satisfies A product of two homogeneous elements A and B has a good parity and is given by |AB| = |A| + |B|. The presence of non-trivial grading also affects the tensor product. The graded tensor product is denoted by ⊛ and defined as where function r designates the row parity, The advantage of definition (100) is that it preserves the standard tensor multiplication rule, The graded tensor product can be extended to product spaces by introducing (homogeneous) elements E ab j , representing the generators associated to the j-th copy of End(C n|m ). Notice that in contrast to the standard (non-graded) basis 1 ⊗ · · · ⊗ E ab ⊗ · · · ⊗ 1 of End(C n+m ) ⊗N , elements E ab j do not commute at different lattice sites, but we find instead On the same lattice site they however still obey the property of projectors, The last two properties combined yield The graded generators acting on the N -particle space (C n|m ) ⊗N read in terms of the graded tensor product whereas expressed in terms of the standard tensor product they assume the expansion This prescription should be interpreted as the higher-rank version of the Jordan-Wigner transformation [150]. Interaction densities h n|m for a class of the so-called 'fundamental graded models' are identified with graded permutations P n|m on C n|m ⊗ C n|m , Permutations P n|m can be alternatively given also as matrices acting on the two-fold fundamental spaces C n+m ⊗ C n+m , reading The defining su(n|m) representations admit realizations in terms of canonical fermions. In the su(1|1) case, the graded projectors E i act non-identically only on the i-th copy of C 1|1 in the chain, and are realized as a 2 × 2 matrix of spinless fermions Here the generators n i and 1 − n i span the even (bosonic) subalgebra V 0 , while c i and c † i are the fermionic generators which span the odd part V 1 and satisfy canonical anticommutation relations Equation (107) is nothing but the well-known Jordan-Wigner transformation from Pauli spins to canonical spinless fermions In systems with multiple fermionic species (e.g. spin-carrying fermions), the super projectors can be constructed with aid of the fusion procedure [150]. For instance, the local physical space of a su(2|2) spin chain is four dimensional, spanned by states At each lattice site i we thus have 4 × 4 = 16 generators, Flattening the indices readily yields the graded permutation on C 2|2 ⊗ C 2|2 , taking the form Furthermore, the fermionic realization of the graded projector P 2|1 can be obtained from P 2|2 by projecting out e.g. the doubly-occupied state |↑↓ . The local space of states thus consists of the triplet Choosing e.g. the grading as |0| = 0 and |1| = |2| = 1, one finds

B Fusion and factorization properties of Lax operators
In this section we revisit the main factorization and fusion formulae for the gl(n|m) integrable spin chains. A comprehensive and more detailed can be found e.g. in references [85][86][87][88].
In the case of unitary sl(n|m) representations, all µ a for a = n must be non-negative integers, while µ n can take arbitrary real values. The fundamental representation of sl(n|m) is given by the weight vector Λ n+m = (1, 0, . . . , 0). Kac-Dynkin labels and finite-dimensional irreducible representations are in one-to-one correspondence, with Young diagrams corresponding to nonnegative non-increasing weights. Rectangular representations {s, a} with s columns and a rows have a single non-vanishing label µ a = s.

B.1 Factorization of Lax operators
Highest-weight gl(n|m)-invariant transfer operators T + Λ n+m (z) acting on a Hilbert space H ∼ = (C n|m ) ⊗N are given by 22 which due to Yang-Baxter relation (8) enjoy the commutativity property, for all z, z ′ ∈ C and representation labels Λ n+m and Λ ′ n+m . Remarkably however, operators T + Λ n+m do not represent the most elementary objects in the theory. In fact, they factorize 23 into an ordered sequence of Q-operators Q {a} 24 , Here parameters λ ′ a are the 'shifted weights', In the gl(n) case, i.e. when m = 0, the shifts arrange in 'complete n-strings', reading ̺ n = Below we exemplify the factorization procedure on a few concrete instances. 22 Here it is implicitly assumed that the super trace exists. Additional regulators in the form of boundary twists may be needed in general. 23 The algebraic origin of the factorization formula has to do with the U(g)-invariant universal R-matrix which decomposes in terms of tensor products of components from the corresponding Borel subalgebras [120,121]. 24 Geometrically, Q-operators can understood as Plücker coordinates on Grassmannian manifolds [122,151].
Basic example: sl(2) case. The factorization property is best illustrated on the sl(2) case. The corresponding highest-weight Lax operators read and are characterized by a single Dynkin label j which parametrized the gl(2) weight vector Λ 2 = (j, −j). The non-compact spin generators J a act on a sl(2) module V + j , and can be conveniently given in the Holstein-Primakoff form where b and b † are the generators of a bosonic oscillator obeying canonical commutation and h = b † b + 1 2 is the mode number operator. We furthermore define two types of Fock space representations, denoted by The two are related to each other under the particle-hole transformation A pair of partonic Lax operators L {1} (z) and L {2} (z) can be straightforwardly obtained from T + Λ 2 (z) by considering two possible way of taking a (correlated) large-j and large-z limits (cf. [85]). This is achieved by keeping either of the combinations z ± = z ± (j + 1 2 ) fixed, resulting in the 'degenerate' Lax operators of the form which represent two distinct well-defined and valid solutions to the Yang-Baxter equation.

B.2 Fusion of partonic Lax operators
Quantum groups are endowed with a coproduct, ensuring that the algebraic structure gets preserved under tensor multiplication Y → Y ⊛ Y. Partonic Lax operators, as defined in Eq. (42), represent the simplest solutions of the Yang-Baxter equation (8). They serve as irreducible components for obtaining other realizations of Y via fusion. Below we outline the main features of such a procedure, while referring the reading for a more complete and detailed presentation to references [86][87][88]. Let I, J ⊆ {1, 2, . . . , n + m} be two index sets. We shall only consider operators L(z) which are linear in spectral parameter z, requiring the sets I and J to be non-intersecting, I ∩ J = ∅. Set I (resp. J) comprises of p (ṗ) bosonic and q (q) fermionic indices. We furthermore introduce K ≡ I ∪ J, involvingp (q) bosonic (fermionic) indices, such that p +ṗ +p = n and q +q +q = m. Fusion is a process of merging two canonical Lax operators L I and L J of respective ranks |I| = p + q and |J| =ṗ +q, which takes the abstract form for some appropriate choice of shifts z 1 and z 2 . The superscript square brackets were needed here to distinguish inequivalent species. The precise prescription for the fusion rule is entailed by the following form where K is a triangular 'disentangling matrix' and S a suitable global similarity transformation. An implication of fusion formula (132) is that Lax operators which are realized in terms of algebras A K m,n are not elementary, but instead factorize according to A K m,n → A I m,n ⊗ A J m,n . Below we take a closer look at this by inspecting a few explicit examples.
We begin by noticing that the above fusion procedure clearly violates the canonical form given by Eq. (40), as it appears to involve an exceeding number of auxiliary spaces. We shall in turn demonstrate that all redundant auxiliary spaces can be eliminated upon appropriately redefining the generators. In particular, there exist a canonical procedure to reduce the number of oscillators from |I| · |I| + |J| · |J| down to |K| · |K| and expressing the gl(|K|) generators in terms of independent generators of gl(|I|) and gl(|J|), dressed with |K| additional oscillators. This is in practice achieved by virtue of the homomorphisms [87] gl(|K|) → gl(p|q) ⊗ gl(ṗ|q) ⊗ osc(p +ṗ|q +q), (133) in which the post-fusion gl(p +ṗ|q +q) generators J ab are given by the following prescription with J ab 1 and Jȧ˙b 2 denoting the gl(p|q) and gl(ṗ|q) super spins, respectively, whereas the oscillators are to be identified as aḃ , a ∈ I ξ [2] aḃ , a ∈ J .
The latter are either bosonic or fermionic, depending on the grading. Finally, the tridiagonal matrix K is of the form where the double-dotted indices represent the summation over K.
sl (2) case. The basic principle of fusion can be explained on the sl(2) theory. Fusion can be understood as the opposite procedure of factorization which is outlined in the previous section. The partonic pieces given by expressions (130) can be fused in two distinct ways.
To this end we introduce square brackets and assign a boson oscillator to each tensor factor, yielding the following operator identity on C 2 ⊗ B ⊗ B, where the input spectral parameters are given by Notice that in the second line of Eq. (138) the oscillators have been rearranged using the similarity transformation S = exp (b † 1 b 2 ) on B ⊗ B, which reads explicitly On the other hand, fusing in the opposite order yields a similar operator identity where again the parameter constraints (139) are imposed. Formula (138) readily implies factorization property for the highest-weight sl(2)-invariant transfer operator T + Λ 2 (z) ≡ T + j (z) in term of a pair of Q-operators, A sequence of transfer matrices T j (z) with 2j ∈ Z, pertaining to finite-dimensional irreducible su(2) representations, can the be obtained from T + j (z) with aid of the Bernstein-Gelfand-Gelfand resolution of finite-dimensional modules, V j = V + j − V + −j−1 , resulting in Bazhanov-Reshetikhin determinant representation [152] A comment in regard to the so-called vacuum Q-operators is in order here. First, recall that vacuum Q-operators represent a family of transfer operators which are constructed from a path-ordered product of (partonic) Lax operators contracted with respect to a suitable 'vacuum state'. Vacuum state are presently identified with the highest (or lowest) weight state in B. Mutual commutativity of vacuum Q-operators can be inferred from fusion formula (138), after observing that (i) when the two Fock space involved are of same type the product vacuum |0 ⊗ |0 remains inert under the action of the similarity transformation S, i.e. S |0 ⊗ |0 = |0 ⊗ |0 , and (ii) the disentangler K has no global effect due to its triangular form. In the opposite scenario, a fusion of two Lax operators which involve two different types of Fock spaces inevitably excites the vacuum to a coherent state which in turn prevents the vacuum Toperator from decomposing into two vacuum Q-operators. For the very same reason vacuum Q-operators pertaining to inequivalent auxiliary modules are not guaranteed to satisfy the involution property. In fact, it can be explicitly confirmed that they do not.

sl(3) case.
The Ω-amplitudes for the integrable steady states constructed in this work are all formed from 'mesonic' Lax operators, namely objects which result from the fusion of two partonic elements. As an explicit example we consider the SU (3)-symmetric Lai-Sutherland chain [107,108], to which we ascribe auxiliary algebra A with the partonic Lax operators reading The oscillators are disentangled with aid of and an additional similarity transformation S = U V, with Explicitly, these transformations act as and V b respectively. Putting everything together, relabelling the oscillators, and representing the generators of sl(2) spins acting on V j (where 2j = j 1 − j 2 + λ) as [1] , yields precisely the anticipated canonical representation of the mesonic Lax operator gl(1|1) case. It may be instructive to also explicitly spell out the fusion step (132) for the gl(1|1) Lie superalgebra. The latter is spanned by four elements, two bosonic generator N and E, and two fermionic ones ψ ± . Commutation relations read where E is the central element. The defining (fundamental) representation is of dimension 2, and is given by 2 × 2 matrices By setting E = 0, we obtain a trivial one-dimensional irreducible representation with ψ ± ≡ 0. There moreover exists a family of two-dimensional irreducible representations, denoted by n, e (with E = 0), which reads and includes the fundamental representation (155) as 1, 1 . Reducible indecomposable representations of gl(1|1) (see e.g. [133,153]) and not of our interest here. Fusion in the fermionic case works as follows. The partonic Lax operators contain a single fermionic specie and are of the form By employing the universal fusion formula, we readily derive the gl(1|1)-invariant Lax operator which takes the form where the gl(1|1) super spin generators J ab should are identified as together with the following constraint on the gl(1) scalars Therefore, L λ (z) belongs to the two-dimensional representation V λ , with λ being the central charge. A comparison with Eq. (156) shows that V λ ≡ 1, 2λ . In terms of fermionic algebra, operator L λ (z) admits an expansion Recall that for λ = 0, module V λ becomes an atypical indecomposable representation (a short multiplet) which is no longer irreducible; it contains a one-dimensional invariant subspace corresponding to the Fock vacuum. However, these exceptional instances do not seem to be relevant in the context of boundary-driven spin chains. We moreover wish to emphasize that prescription (134) provides an explicit oscillator realization of gl(n|m) Lie superalgebras. Let us consider the gl(2|1) case as an example, and fixing the grading to −− . The bosonic subalgebra is a direct sum gl(2) ⊕ u(1), and is spanned by gl (2) and u(1) generator In addition, there are four fermionic charges which are parametrized as The sl(1|2) case can be obtained by restricting J ab 1 to sl(2) spins acting on V j 1 , while j 2 is the remaining Dykin label.

C Non-interacting fermions
In this section we provide the solution to the problem of boundary-driven non-interacting spinless fermions hopping on a one-dimensional lattice. The Hamiltonian of the model can be seen as a Yang-Baxter integrable spin chain invariant under gl(1|1) Lie superalgebra. 25 25 The model can also be mapped to the XX Heisenberg spin chain Hamiltonian. In the spin picture, we deal with a model invariant under the q-deformed quantum symmetry Uq(sl(2)) for the value of the deformation parameter q = i. In Section 5.2 we constructed steady-state solutions for the simplest fermionic boundary reservoirs with equal coupling strengths. Our aim here is to demonstrate that the problem of free fermions represents a special case which even allows for solutions going beyond those discussed in Section 5.2. Quite remarkably, the operator Schmidt rank 26 of ρ ∞ now equals 4 and thus does not depend on the system size. This is in stark contrasted to the solutions pertaining to higher-rank symmetries which all exhibit Schmidt ranks which grow algebraically with system size. A finite Schmidt rank can be understood as a strong indication that the problem of finding the steady states may be tractable by 'brute force', that is first by explicitly computing the null space of the Liouvillian generator L and subsequently analytically parametrizing the solution (e.g. with help of symbolic algebra routines). An obvious advantage of this approach is that it does not require any prior knowledge of underlying algebraic structures. This allows allows to conveniently parametrize the solutions directly in terms of physical couplings attributed to the boundary reservoirs.
We shall provide an extended four-parametric set of solutions for the asymmetric driving involving coupling rate parameters g, ζ ∈ R, supplemented with two additional boundary external fields, of unequal magnitudes h L , h R ∈ R. We notice that the solutions given below appear to lie outside of gl(1|1)-invariant Lax operators of Eqs. (68) and (69). 27 The L-operator for the Ω-amplitude is now formally linked to a two-dimensional auxiliary representation denoted by V u . Here u is a four-component vector label which involves the boundary parameters, u = (g, ζ, h L , h R ). In terms of Pauli matrices, the L-operator admits the expansion using shorthand notationsh i = h i − 1 and δh = h L − h R . It can easily be verified that the L-operator provides a multi-colored family commuting of transfer matrices, T (u) = Tr Vu L Vu ⊗ · · · ⊗ L Vu , [T (u 1 ), T (u 2 )] = 0 ∀u 1 , u 2 ∈ C 4 .
Involution property of T (u) is ensured by the multi-colored 6-vertex R-matrix R Vu 1 Vu 2 which operates on V u 1 ⊗ V u 2 , 26 The operator analogue of the Schmidt rank characterizes a degree of bipartite entanglement of a mixed quantum state. In the language of matrix-product states it coincides with the bond dimension. 27 In practice it turns out that free fermions allows even more general types of non-perturbative integrable boundaries which then those considered here.
intertwining two copies of auxiliary spaces V u associated to a pair of L-operators L Vu i acting on C 2 ⊗ V u i , R Vu 1 Vu 2 L Vu 1 L Vu 2 = L Vu 2 L Vu 1 R Vu 1 Vu 2 .
By expanding it to the entire spin chain, we obtain the action of the unitary part L 0 which in distinction to the canonical solutions discussed in the paper this time (assuming non-vanishing δh) acquires an additional term,