N Crampe

K Mallick

E Ragoucy
email: eric.ragoucy@lapth.cnrs.fr

M Vanicat

Open two-species exclusion processes with integrable boundaries

We give a complete classification of integrable Markovian boundary conditions for the asymmetric simple exclusion process with two species (or classes) of particles. Some of these boundary conditions lead to non-vanishing particle currents for each species. We explain how the stationary state of all these models can be expressed in a matrix product form, starting from two key components, the Zamolodchikov-Faddeev and Ghoshal-Zamolodchikov relations. This statement is illustrated by studying in detail a specific example, for which the matrix Ansatz (involving 9 generators) is explicitly constructed and physical observables (such as currents, densities) calculated.

Introduction

Particles that diffuse on a discrete lattice, with anisotropic hopping rates and hard-core exclusion, define one of the simplest examples of a driven diffusive system [START_REF] Katz | Nonequilibrium steady states of stochastic lattice gas models of fast ionic conductors[END_REF][START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF][START_REF] Schmittmann | Statistical mechanics of driven diffusive systems[END_REF], known as the asymmetric simple exclusion process (ASEP). Thanks to its simplicity, the ASEP plays a seminal role in classical low-dimensional transport with excluded volume constraints and appears in numerous phenomenological models in hard and soft condensed matter (see e.g. [START_REF] Chou | Non-equilibrium statistical mechanics: From a paradigmatic model to biological transport[END_REF][START_REF] Derrida | An exactly soluble non-equilibrium system: the asymmetric simple exclusion process[END_REF][START_REF] Schütz | Exactly Solvable Models for Many-Body Systems Far From Equilibrium[END_REF] and references therein). The exclusion process is also a paradigm in the field of non-equilibrium statistical mechanics [START_REF] Bertini | Macroscopic Fluctuation Theory for stationary non-equilibrium states[END_REF][START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF][START_REF] Liggett | Stochastic Models of Interacting Systems: Contact, Voter and Exclusion Processes[END_REF][START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF]: it displays an unexpectedly complex behaviour [START_REF] Varadhan | The complex history of simple exclusion[END_REF], that has stimulated many elaborate analytical studies. A distinctive feature of the ASEP (as compared to many other interacting particles models such as the zero range process) is that it is integrable. This remarkable property explains its tremendous success in mathematical physics, combinatorics, representation theory and probability theory [START_REF] Alcaraz | Reaction-diffusion processes, critical dynamics and quantum chains[END_REF][START_REF] Corwin | The Kardar-Parisi-Zhang equation and universality class[END_REF][START_REF] Duchi | A combinatorial approach to jumping particles[END_REF][START_REF] Isaev | Diffusion algebras[END_REF][START_REF] Kriecherbauer | A pedestrian's view on interacting particle systems, KPZ universality and random matrices[END_REF][START_REF] Uchiyama | Asymmetric simple exclusion process with open boundaries and Ashkey-Wilson polynomials[END_REF]. Besides, many intricate questions about non-equilibrium behaviour, when formulated for the ASEP, allow for precise and quantitative answers, often involving elegant mathematical structures. Hence, the ASEP provides us with benchmarks for a general theory of non-equilibrium systems [START_REF] Bertini | Macroscopic Fluctuation Theory[END_REF][START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF][START_REF] Derrida | Microscopic versus macroscopic approaches to non-equilibrium systems[END_REF].

In the exclusion process, identical particles undergo a stochastic evolution on a lattice; empty sites are called holes. A natural generalization of the original dynamics involves multiple species (or classes) of particles; second class particles were initially used as a tool to locate the shock that may appear at the microscopic level in the density profile of the ASEP [START_REF] Derrida | Exact solution of the totally asymmetric exclusion process: shock profiles[END_REF][START_REF] Derrida | Shock profiles for the asymmetric simple exclusion process in one dimension[END_REF][START_REF] Ferrari | Microscopic shocks in one dimensional driven system[END_REF][START_REF] Ferrari | Microscopic structure of travelling waves in the asymmetric simple exclusion process[END_REF]. More generally, introducing multiple species of particles is a way to couple independent ASEP models, a standard probabilistic technique that allows to gain insight on individual processes [START_REF] Liggett | Stochastic Models of Interacting Systems: Contact, Voter and Exclusion Processes[END_REF][START_REF] Thorisson | Coupling, Stationarity, and Regeneration[END_REF]. Moreover, the ASEP with multiple species has some physical applications, e.g. in relation to traffic flow [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF][START_REF] Evans | Spontaneous symmetry-breaking in a one-dimensional driven diffusive system[END_REF][START_REF] Karimipour | Multispecies asymmetric simple exclusion process and its relation to traffic flow[END_REF].

The steady state distribution of the two-species exclusion process on a ring was obtained in [START_REF] Derrida | Exact solution of the totally asymmetric exclusion process: shock profiles[END_REF] by using the matrix product representation introduced by B. Derrida et al. in [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF]. This exact solution of the two-species exclusion process had multiple outcomes: shock profiles were calculated analytically [START_REF] Derrida | Exact solution of the totally asymmetric exclusion process: shock profiles[END_REF]; the invariant measure was proved to be non-Gibbsean [START_REF] Speer | The two-species asymmetric simple exclusion process[END_REF]; quadratic algebras involved in the 'matrix Ansatz' were studied in [START_REF] Aneva | Deformed coherent and squeezed states of multiparticle processes[END_REF][START_REF] Essler | Representations of the quadratic algebra and partially asymmetric diffusion with open boundaries[END_REF][START_REF] Isaev | Diffusion algebras[END_REF]; combinatorial interpretations of the weights generated by the matrices were found [START_REF] Angel | The stationary measure of a 2-type totally asymmetric exclusion process[END_REF][START_REF] Duchi | A combinatorial approach to jumping particles[END_REF][START_REF] Ferrari | Invariant measure for a two-species asymmetric process[END_REF] etc. These results were extended to models with an impurity [START_REF] Mallick | Shocks in the asymmetry exclusion model with an impurity[END_REF] and to systems defined on the infinite line [START_REF] Derrida | Shock profiles for the asymmetric simple exclusion process in one dimension[END_REF] (many references can be found in the standard review [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF]). The solution of the exclusion process with three or more species on a periodic or on the infinite lattice turned out to be more involved; it relies on an astute combinatorial construction [START_REF] Ferrari | Multiclass processes, dual points and M/M/1 queues, Markov Processes[END_REF][START_REF] Ferrari | Stationary distributions of multi-type totally asymmetric exclusion processes[END_REF] (see also [START_REF] Ayyer | Correlations in the multispecies TASEP and a conjecture by Lam[END_REF][START_REF] Aas | Continuous multi-line queues and the TASEP[END_REF]) that can be encoded by a tensor matrix Ansatz [START_REF] Evans | Matrix Representation of the Stationary Measure for the Multispecies TASEP[END_REF][START_REF] Mallick | Exact solution of an exclusion process with three classes of particles and vacancies[END_REF][START_REF] Prolhac | The matrix product solution of the multispecies partially asymmetric exclusion process[END_REF]. These multi-species exclusion processes remain integrable and systems defined on a periodic lattice can, in principle, be solved by a nested Bethe Ansatz [START_REF] Alcaraz | Exact solution of asymmetric diffusion with second-class particles of arbitrary size[END_REF][START_REF] Arita | Spectrum of a multi-species asymmetric simple exclusion process on a ring[END_REF][START_REF] Arita | Remarks on the multi-species exclusion process with reflective boundaries[END_REF][START_REF] Cantini | Algebraic Bethe Ansatz for the two species ASEP with different hopping rates[END_REF][START_REF] Derrida | Bethe Ansatz Solution for a Defect Particle in the Asymmetric Exclusion Process[END_REF].

Models defined on a finite lattice with open boundaries connected to reservoirs, that allow for the exchange of particles, represent a major challenge. From the physicist's point of view, these models are more realistic and therefore the most interesting. Indeed, the fundamental picture often used to illustrate a non-equilibrium system is a pipe connecting two (or more) reservoirs at different temperatures, or chemical potentials. The discovery of boundary-induced phase transitions in the open exclusion process [START_REF] Krug | Boundary-induced phase transitions in driven diffusive systems[END_REF], followed by the matrix Ansatz solution [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF] were the seminal contributions that did trigger the whole field. However, for the case of the multispecies ASEP with open boundaries many question remain to be clarified [START_REF] Arita | Exact analysis of two-species asymmetric exclusion process with open boundary conditions[END_REF][START_REF] Arita | Phase transitions in the two-species totally asymmetric exclusion process with open boundaries[END_REF][START_REF] Ayyer | On the Two Species Asymmetric Exclusion Process with Semi-Permeable Boundaries[END_REF][START_REF] Ayyer | On Some Classes of Open Two-Species Exclusion Processes[END_REF][START_REF] Derrida | Nonequilibrium Statistical Mechanics in One Dimension[END_REF][START_REF] Uchiyama | Two-Species Asymmetric Simple Exclusion Process with Open Boundaries[END_REF] (see section 2.3 for the state of the art).

In the present work, we wish to investigate the asymmetric exclusion process with different species of particles. This model, which will be called in brief the '2-ASEP', is defined on a finite lattice of L sites where boundary sites are connected to reservoirs. The precise dynamical rules of the 2-ASEP will be given in Section 2. At the first and the last sites, particles can either enter or exit: thus, at each boundary, there are 6 possible choices for the corresponding rates. The aim of this work is to give a complete classification of the boundary rates that lead to integrable models. In other words, we find the conditions on the boundary rates that preserve the integrability of the 2-ASEP. This is done by solving simultaneously the Yang-Baxter equation in the bulk and the Sklyanin reflection equations at the boundaries [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF][START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF][START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]. The models thus obtained, are the only ones solvable by Bethe Ansatz and we are convinced that they are also analytically tractable by Matrix Ansatz [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF].

The outline of this work is as follows. In section 2, we present all the dynamical rules for the 2-ASEP with open boundaries that preserve integrability; we also discuss the relation between these integrable models and previous studies. In section 3, we apply the theory of integrable systems to the open 2-ASEP, starting from the (bulk) R-matrix and two boundary K matrices. The operators R and K have to satisfy non-linear algebraic relations (the Sklyanin reflection equations) for which all solutions are found and classified. Each solution allows us to define a one-parameter family of commuting transfer matrices, which yields, finally, to an integrable Markov operator of the 2-ASEP with open boundaries. In the subsection 3.4, we explain, using the results of [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF], how a matrix Ansatz can be derived from this formalism. Section 4 is devoted to a specific example, perhaps the simplest 2-TASEP model with open boundaries in which all types of particles can enter and leave the system. We construct a Matrix Ansatz for this model: the underlying algebra contains nine generators and differs from the quadratic algebra that solves the single species ASEP with open boundaries (hereafter, we shall refer to this quadratic algebra, introduced in [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF], as the DEHP algebra). The last section is devoted to some concluding remarks.

Integrable dynamics for the open 2-ASEP

The ASEP case

We start by recalling the dynamical rules of the 2-ASEP. Each site on the lattice is occupied by a single particle, and there are three different types (or classes) of particles, denoted by 3, 2 and 1. The occupation variable at a site i of the system, with 1 ≤ i ≤ L, is denoted by τ i : we have τ i = 1, 2 or 3. The bulk dynamics is defined as follows: a bond (i, i + 1), with 1 ≤ i ≤ L -1, between two neighboring lattice sites, is updated between time t and t + dt by swapping the particles at i and i + 1; if τ i > τ i+1 , the exchange occurs with rate 1, whereas if τ i < τ i+1 the exchange occurs with rate q < 1.

These bulk rules are summarized in the following table: These rules show that particles of type 3 have the highest priority, followed by 2's and lastly by 1's (we emphasize that our notation differs from that of some previous works [START_REF] Evans | Matrix Representation of the Stationary Measure for the Multispecies TASEP[END_REF]. Type 3 particles having highest priority were called first class particles, type 2 second class and 1's play the role of holes). These bulk rules define an integrable dynamical system, that can be encoded by a R-matrix, that satisfies the Yang-Baxter equation [START_REF] Alcaraz | Exact solution of asymmetric diffusion with second-class particles of arbitrary size[END_REF][START_REF] Arita | Spectrum of a multi-species asymmetric simple exclusion process on a ring[END_REF][START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF]. The R-matrix will be given explicitly in the next section.

Particles are allowed to enter or to exit from both boundaries and the corresponding entrance/exit rates may depend on the type of the particle that was previously located at the boundary. More precisely, both on the left and on the right boundary, we can have a transition of the type I → J, where I and J represent particles of different classes, i.e. 1 ≤ I, J ≤ 3 and I = J. This leads to 12 independent rates (6 on each side). However, for arbitrary choices of these rates, the models will not be integrable.

Hence, in order to preserve integrability, these 12 rates cannot be chosen in an arbitrary manner. They have to satisfy integrable boundary conditions that are associated to the bulk R-matrix. These boundary conditions are expressed as a set of reflection equations [START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF], which will be solved in Section 3 and their solutions, the K-matrices, which will be fully determined.

From this theory of integrable systems with boundaries, we shall derive four sets of integrable rules at the left boundary: out of the six possible rates, only two are independent, the other ones either vanish or are simple functions of the two independent rates. Calling by α and γ two arbitrary positive real numbers, we find that the four set of left-rules are given by

L 1 L 2 L 3 L 4 1 f (α,γ) -----→ 2 1 γ -----→ 2 1 α --→ 3 1 α --→ 3 2 α -----→ 1 1 α -----→ 3 2 γ --→ 1 3 γ --→ 1 3 α -----→ 1 2 α -----→ 3 2 α --→ 3 3 γ -----→ 2 3 g(α,γ) -----→ 2 3 γ --→ 1 (2.2)
where we have defined f (x, y) = y(x + y + 1 -q)

x + y and g(x, y) = y(x + y + q -1)

x + y .

(2.3)

Remark: Note that in the L 2 case, we must have α + γ > 1 -q to ensure g(α, γ) > 0

Similarly, we obtain four integrable rules at the right boundary that can be encoded in the following table

R 1 R 2 R 3 R 4 1 δ -----→ 2 1 g(β,δ) -----→ 2 1 δ --→ 3 1 δ --→ 3 1 β -----→ 3 2 β -----→ 1 2 β --→ 1 3 β --→ 1 2 β -----→ 3 3 β -----→ 1 2 δ --→ 3 3 f (β,δ) -----→ 2 3 δ -----→ 2 3 β --→ 1 (2.4)
where β and δ are two positive real numbers and f (x, y), g(x, y) were defined in (2.3).

The TASEP case

We emphasize here the TASEP limit (q → 0). For the TASEP, the bulk dynamics is given by:

21 1 --→ 12 ; 31 1 --→ 13 ; 32 1 --→ 23. (2.5)
The integrable boundary transition rates compatible with the 2-TASEP dynamics are obtained by taking the limit q → 0 in (2.2) and (2.4). However, this limit allows for some 'unphysical' transitions, where particles can be injected or extracted in the direction opposite to the natural orientation defined by the bulk. In order to avoid these situations, we have to fix the values of the parameters γ and δ. We take γ = 1 -α and δ = 1 -β for L 2 and R 2 (so that g(α, γ) = 0 and g(β, δ) = 0), and γ = δ = 0 for boundaries of type L 3 , L 4 , or R 3 , R 4 . Besides, in the cases L 1 and R 1 , we need first to rescale α = α ′ q, γ = γ ′ q, β = β ′ q and δ = δ ′ q and then take the limit q → 0.

Altogether, we end up with four sets of integrable transition rates on the left boundary:

L 1 L 2 L 3 L 4 1 µ --→ 2 1 1-α ----→ 2 1 α --→ 3 1 α --→ 3 1 α ---→ 3 2 α --→ 3 µ = γ ′ α ′ +γ ′ 2 α ---→ 3 (2.6)
There are also four sets of integrable transition rates on the right boundary:

R 1 R 2 R 3 R 4 3 ν --→ 2 2 β ---→ 1 2 β --→ 1 3 β --→ 1 3 β ---→ 1 3 β --→ 1 ν = δ ′ β ′ +δ ′ 3 1-β ----→ 2 (2.7)
The transition rates can be chosen among this four classes independently for each boundary so that we have at the end 4 × 4 = 16 different integrable models. The local boundary operators and the corresponding K-matrices will be given in the next section.

Comparison with previous works

We now relate the integrable boundary conditions that we have found to previously known results about the 2-TASEP with open boundaries. In [START_REF] Arita | Exact analysis of two-species asymmetric exclusion process with open boundary conditions[END_REF][START_REF] Arita | Phase transitions in the two-species totally asymmetric exclusion process with open boundaries[END_REF], C. Arita found a matrix Ansatz for the nonequilibrium stationary state (NESS) of a 2-TASEP, using a variant of the quadratic algebra of [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF]. His solution was then generalized to the partially asymmetric case [START_REF] Uchiyama | Two-Species Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. The structure of this NESS, the density profiles, the existence of fat shocks and some mathematical properties of the invariant measures were thoroughly investigated by A. Ayyer at al. [START_REF] Ayyer | On the Two Species Asymmetric Exclusion Process with Semi-Permeable Boundaries[END_REF]. These authors emphasized that all the models hitherto considered had the very special property that the boundaries were semi-permeable: second class particles (2) can neither enter nor leave the system. These semi-permeable boundaries correspond to our L 4 and R 4 . This is a drastically restrictive condition: it implies, in particular, that the current J 2 of second class particles vanishes. In fact, it was shown in [START_REF] Ayyer | On the Two Species Asymmetric Exclusion Process with Semi-Permeable Boundaries[END_REF] that the distribution of these second class particles is given by an equilibrium ensemble in fixed volume. We emphasize that the integrable boundary conditions that we have found, are not restricted solely to semi-permeable boundaries. Some of the models we defined are genuinely out of equilibrium, in all respects.

In a second work [START_REF] Ayyer | On Some Classes of Open Two-Species Exclusion Processes[END_REF], A. Ayyer et al. considered more general classes of open two-species exclusion processes. They showed that a matrix solution based on a quadratic algebra could also be used for some permeable boundaries, if the rates obey some particular relations. The set of solvable 2-TASEP models was thus extended. However, these authors also proved that although the boundaries were permeable, the current J 2 still identically vanishes in these models (in the limit of large system sizes). It seemed that all models, solvable by a quadratic DEHP algebra with 3 generators, similar to the ones introduced in [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF], had J 2 = 0: this corresponds to the choices

L 4 /R 2 , L 2 /R 4 , L 4 /R 1 , L 1 /R 4 , L 2 /R 2 and L 1 /R 1 .
Note that for the choices

L 4 /R 3 , L 3 /R 4 and L 3 /R 3 (resp. L 1 /R 2 , resp. L 2 /R 1 )
, the density of 2 (resp. 3, resp. 1) identically vanishes for all system sizes. The stationary state can be computed using exactly the same algebra as in [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF].

Therefore, we are left with 4 possible integrable boundaries L 2 /R 3 , L 3 /R 2 , L 3 /R 1 and L 1 /R 3 with non-vanishing currents for each species of particles. Building a matrix Ansatz for these models requires a more complicated algebra, using the construction introduced in [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF]. In Section 4, we shall analyze an example, that cannot be solved by a DEHP algebra: it requires at least 9 generators (or a more complex tensor-product construction based on a 3 generators algebra [25]).

Finally, in [START_REF] Ayyer | On Some Classes of Open Two-Species Exclusion Processes[END_REF] other types of models, called colorable, were studied; these models are not necessarily solvable (e.g., a matrix representation for the steady state does not exist in general) but currents and density profiles can be easily extracted from the known results for the single species TASEP. The idea is to identify two different classes of particles (by painting them with the same color) and reduce the 2-TASEP to the original TASEP. There are two possible colorings: either we group the 1's and the 2's together (and treat them as holes) and leave the particles 3 alone; or we group the 2's and the 3's (colored particles) and leave the 1's alone (holes). Both colorings are compatible with bulk dynamics; the question is whether the boundary conditions also permit these identifications. The boundary conditions L 2 , L 3 and R 2 , R 3 are colorable. Therefore, a model with these boundaries conditions is integrable (by construction) and the rates α, β that appear in (2.6) and (2.7) can be chosen so that none of the currents vanishes: there exists 2-ASEP models with open boundaries which are exactly solvable and in which all particles genuinely behave out of equilibrium.

Commuting transfer matrices for the 2-ASEP

In this section, we recall the integrability properties of the 2-ASEP and explain how the boundary dynamical rules described in section 2 are derived.

Integrability of the open 2-ASEP

A Markov process can be defined by a master equation

d|P t dt = M |P t
where |P t is the vector containing the probabilities of the different configurations of the system at time t and M is the Markov matrix. For the exclusion process only local jumps are permitted and M can be written as:

M = B 1 + L-1 ℓ=1 w ℓ,ℓ+1 + BL . (3.1)
The operator M acts in the configuration space C 3 ⊗L , where each component of the tensor space stands for a site on the lattice. The subscripts in (3.1) indicate in which component of the tensor space the operators are acting on non trivially. Each w ℓ,ℓ+1 is a 9 by 9 matrix that encodes the dynamics in the bulk (i.e. the update rules on the bond (ℓ, ℓ + 1)). The matrix B (respectively B) is 3 by 3 and gives the dynamics on the left (respectively right) boundary. This stochastic model is said to be integrable when the Markov matrix M can be embedded in a one-parameter family of commuting operators t(z), called transfer matrices. A systematic method of constructing integrable models with open boundaries was given by Sklyanin [START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]. It rests on three algebraic objects: the (bulk) R-matrix and two boundary matrices K and K.

The R-matrix obeys the Yang-Baxter equation [START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF]:

R 12 x 1 x 2 R 13 x 1 x 3 R 23 x 2 x 3 = R 23 x 2 x 3 R 13 x 1 x 3 R 12 x 1 x 2 . ( 3.2)
The matrix K has to satisfy the so-called reflection equation [START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]:

R 12 x 1 x 2 K 1 (x 1 ) R 21 (x 1 x 2 ) K 2 (x 2 ) = K 2 (x 2 ) R 12 (x 1 x 2 ) K 1 (x 1 ) R 21 x 1 x 2 . (3.3)
Finally, the matrix K obeys a dual-reflection equation:

K 2 (x 2 ) R t 1 21 (x 1 x 2 ) -1 t 1 K 1 (x 1 ) R 21 x 2 x 1 = R 12 x 2 x 1 K 1 (x 1 ) R t 2 12 (x 1 x 2 ) -1 t 2 K 2 (x 2 ) (3.4)
where t 1 and t 2 represent transpositions with respect to the first and the second space, respectively.

Remark: K(x) can be constructed [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF] from an auxiliary operator K(x), using the relation

K 1 (x) = tr 0 K0 1 x ((R 01 (x 2 ) t 1 ) -1 ) t 1 P 01 . (3.5)
This auxiliary K(x) obeys the reflection equation (3.3) associated with the matrix R = R -1 :

R 21 x 2 x 1 K1 (x 1 ) R 12 1 x 1 x 2 K2 (x 2 ) = K2 (x 2 ) R 21 1 x 1 x 2 K1 (x 1 ) R 12 x 2 x 1 . (3.6)
The symmetry between left and right boundaries is now much more manifest (compare (3.3) and (3.6)).

Using R, K and K, a one parameter family of transfer matrices t(z) is defined as follows [START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]:

t(z) = tr 0 K 0 (z) T 0 (z) K 0 (z) T 0 (1/z) -1 with T 0 (z) = R 0,L (z) R 0,L-1 (z) • • • R 01 (z). (3.7)
The parameter z that labels these transfer matrices is known as the spectral parameter. The Yang-Baxter and the reflection equations imply that [t(x), t(y)] = 0, for all x, y [START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF][START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]. The transfer matrices t(z) form a one-parameter family of commuting operators that can be co-diagonalized. Their commutation ensures that there exists enough conserved quantities to imply integrability. The last step is to relate this formal construction to the ASEP Markov matrix M . The connection is given by

M = q -1 2 t ′ (1) = B 1 + L-1 ℓ=1 w ℓ,ℓ+1 + BL (3.8) with B = q-1 2 K ′ (1), B = -q-1 2
K′ (1) and w = (q -1)P R ′ (1), where the prime denote derivative w.r.t. the spectral parameter x and P is the permutation matrix. The Markov matrix M belongs to the integrable commuting family t(z); it is itself integrable.

We observe that the bulk transition rates are given by the choice of the R-matrix, while the boundary matrices B and B are obtained from the K-matrices. Therefore, once R is given, the problem of finding integrable boundaries is equivalent to solving the reflection and the dualreflection equations (3.3) and (3.4).

Classification of the Markovian K-matrices for the 2-ASEP

The R-matrix of the multi-species ASEP, satisfying the Yang-Baxter equation is well known [START_REF] Alcaraz | Exact solution of asymmetric diffusion with second-class particles of arbitrary size[END_REF][START_REF] Arita | Spectrum of a multi-species asymmetric simple exclusion process on a ring[END_REF]. For the 2-ASEP it is given by

R(x) =                  1 0 0 0 0 0 0 0 0 0 (x-1)q qx-1 0 (q-1)x qx-1 0 0 0 0 0 0 0 (x-1)q qx-1 0 0 0 (q-1)x qx-1 0 0 0 q-1 qx-1 0 x-1 qx-1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 (x-1)q qx-1 0 (q-1)x qx-1 0 0 0 q-1 qx-1 0 0 0 x-1 qx-1 0 0 0 0 0 0 0 q-1 qx-1 0 x-1 qx-1 0 0 0 0 0 0 0 0 0 1                  . ( 3

.9)

Remark: The ASEP R-matrix is a twisted version of the U q ( gl 3 ) R-matrix R(x) = F 12 R XXZ (x) F -1 21 . However, the twist F 12 is not factorisable (i.e. one cannot write F 12 = U 1 V 2 ). Therefore, the classification [START_REF] Abad | Nondiagonal solutions to reflection equations in SU(N) spin chains[END_REF] of K-matrices for U q ( gl 3 ) cannot be used here and has to be performed ab initio.

The problem consists in finding the matrices K and K obeying the reflection equations (3.3) and (3.4). We focus on unitary Markovian matrices K (i.e. unitary matrices whose entries of each column sum up to 1) such that the entries of the column of B = q-1 2 K ′ (1) sum up to 0 (and similarly for B).

The different steps to derive this classification are the following: we take a generic K-matrix with 9 unknown functions as entries, we project on the different entries the reflection equation (3.3) to get 81 relations, we solve these functional equations and among the solutions we retain only the Markovian ones. Altogether, we have found four different classes of solutions (each depending on α and γ, two free parameters):

K (1) (x) =       x((α 2 -γ 2 )(x-1)+(γx+α)(q-1)) (αx+γ)((α+γ)(x-1)+q-1) (x 2 -1)(α+γ)α (αx+γ)((α+γ)(x-1)+q-1) (x 2 -1)(α+γ)α (αx+γ)((α+γ)(x-1)+q-1) (x 2 -1)(α+γ+1-q)γ (αx+γ)((α+γ)(x-1)+q-1) -(α 2 -γ 2 )(x-1)+(αx+γ)(1-q) (αx+γ)((α+γ)(x-1)+q-1) (x 2 -1)(α+γ)γ x(αx+γ)((α+γ)(x-1)+q-1) 0 0 -(α+γ)(x-1)+x(1-q)
x((α+γ)(x-1)+q-1)

      (3.10) K (2) (x) =       -x((α+γ)(x-1)+1-q)
(α+γ)(x-1)+x(q-1) 0 0

x(x 2 -1)(α+γ)γ (γx+α)((α+γ)(x-1)+x(q-1)) -x((α 2 -γ 2 )(x-1)+(γx+α)(1-q)) (γx+α)((α+γ)(x-1)+x(q-1))

(x 2 -1)(α+γ+q-1)γ (γx+α)((α+γ)(x-1)+x(q-1)) (x 2 -1)(α+γ)α (γx+α)((α+γ)(x-1)+x(q-1)) (x 2 -1)(α+γ)α (γx+α)((α+γ)(x-1)+x(q-1)) (α 2 -γ 2 )(x-1)+(αx+γ)(q-1) (γx+α)((α+γ)(x-1)+x(q-1))       (3.11) K (3) (x) =       -x((α-γ)(x-1)+1-q) (γx+α)(x-1)+x(q-1) (x 2 -1)γ (γx+α)(x-1)+x(q-1) (x 2 -1)γ (γx+α)(x-1)+x(q-1) 0 -(αx+γ)(x-1)+x(1-q) (γx+α)(x-1)+x(q-1) 0 (x 2 -1)α (γx+α)(x-1)+x(q-1) (x 2 -1)α (γx+α)(x-1)+x(q-1) (α-γ)(x-1)+x(q-1) (γx+α)(x-1)+x(q-1)       (3.12)
and

K (4) (x) =      x((γ-α)(x-1)+q-1)
(γx+α)(x-1)+x(q-1) 0 (x 2 -1)γ (γx+α)(x-1)+x(q-1) 0 1 0

(x 2 -1)α
(γx+α)(x-1)+x(q-1) 0 -(γ-α)(x-1)+x(1-q) (γx+α)(x-1)+x(q-1)

     (3.13)
They correspond respectively to the following left boundary matrices, obtained by taking the derivative at x = 1 as in (3.8):

B (1) =      -γ(α+γ+1-q) α+γ α α γ(α+γ+1-q) α+γ -α γ 0 0 -α -γ      , B (2) =      -α -γ 0 0 γ -α (α+γ+q-1)γ α+γ α α -(α+γ+q-1)γ α+γ      (3.14) B (3) =     -α γ γ 0 -α -γ 0 α α -γ     , B (4) =     -α 0 γ 0 0 0 α 0 -γ     . (3.15)
This corresponds using (3.8) to the transition rates at the left boundary given by L 1 to L 4 in the Table (2.2). The K matrices, that obey (3.6) can be obtained directly using the formula

K(i) (x) = U K (i) (1/x) U -1 with U =   0 0 1 0 1 0 1 0 0   , (3.16)
and replacing the parameters α and γ by β and δ. The associated right boundary matrices B are given by:

B(1) =      -β -δ 0 0 δ -β (β+δ+1-q)δ β+δ β β -(β+δ+1-q)δ β+δ      , B(2) =      -δ(β+δ+q-1) β+δ β β δ(β+δ+q-1) β+δ -β δ 0 0 -β -δ      , (3.17) B(3) =     -δ β β 0 -β -δ 0 δ δ -β     , B(4) =     -δ 0 β 0 0 0 δ 0 -β     (3.18)
where β and δ are free parameters. The transition rates at the right boundary given in Table (2.4) are retrieved. Note the exchange q ↔ 1 between the solutions for B and those for B.

The K-matrices for the 2-TASEP

One can perform the same classification of unitary Markovian integrable boundary conditions for TASEP models. The R-matrix for these models is the limit q → 0 of the ASEP R-matrix (3.9):

R(x) =               1 0 0 0 0 0 0 0 0 0 0 0 x 0 0 0 0 0 0 0 0 0 0 0 x 0 0 0 1 0 1 -x 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 x 0 0 0 1 0 0 0 1 -x 0 0 0 0 0 0 0 1 0 1 -x 0 0 0 0 0 0 0 0 0 1               (3.19)
However, the solutions to the corresponding reflection equation include unphysical solutions, where the corresponding boundary matrix B has unwanted transition rates as discussed previously. This unphysical solutions are exactly the cases for which it is impossible to construct the transfer matrix (the RHS of (3.5) diverges when q → 0). Hence, in the following, we give only the solutions leading to physical boundary matrix B, which correspond to triangular K matrices. The triangular Markovian K-matrices solutions to the reflection equation are given by

K (1) (x) =      x(µ(x-1)+1) (1-µ)x+µ 0 0 -µ(x 2 -1) (1-µ)x+µ 1 0 0 0 1      , K (2) (x) =      x 2 0 0 -x(x 2 -1)(α-1) α(x-1)-x -x(α(x-1)+1) α(x-1)-x 0 α(x 2 -1) α(x-1)-x α(x 2 -1) α(x-1)-x 1      (3.20) K (3) (x) =       -x(α(x-1)+1) α(x-1)-x 0 0 0 -x(α(x-1)+1) α(x-1)-x 0 α(x 2 -1) α(x-1)-x α(x 2 -1) α(x-1)-x 1       , K (4) (x) =      -x(α(x-1)+1) α(x-1)-x 0 0 0 1 0 α(x 2 -1) α(x-1)-x 0 1      (3.21)
They correspond respectively to the following left boundary matrices:

B =     -µ 0 0 µ 0 0 0 0 0     ,     -1 0 0 1 -α -α 0 α α 0     ,     -α 0 0 0 -α 0 α α 0     ,     -α 0 0 0 0 0 α 0 0     (3.22)
where α, µ are free parameters. The transition rates given in Table (2.6) are thus obtained. Using the same symmetry as for the ASEP case, we can construct the K matrices and the B operators:

B =     0 0 0 0 0 ν 0 0 -ν     ,     0 β β 0 -β 1 -β 0 0 -1     ,     0 β β 0 -β 0 0 0 -β     ,     0 0 β 0 0 0 0 0 -β    
(3.23) where β, ν are free parameters. This corresponds to the transition rates given in Table (2.7).

Matrix Ansatz

We now apply the results of [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF] to construct a matrix product representation of the steady state using the fundamental R, K and K matrices found above. We believe that integrability provides a natural framework to obtain a bona fide matrix ansatz. As explained in [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF], this claim is based on the analysis of ZF algebras and the study of their representations as a starting point. In section 4, we shall give an example of a 2-TASEP for which integrability leads to a matrix ansatz.

More precisely, we define a 3-component column vector A(x) depending on a spectral parameter x and with algebraic entries, satisfying the following equations:

R 12 x 1 x 2 A 1 (x 1 )A 2 (x 2 ) = A 2 (x 2 )A 1 (x 1 ). (3.24) W | K(x)A 1 x -A(x) = 0, and K(x)A 1 x -A(x) |V = 0. (3.25)
Equation (3.24), valid in the bulk, is known as the Zamolodchikov-Faddeev (ZF) relation [START_REF] Zamolodchikov | Factorized S-matrices in two dimensions as the exact solutions of certain relativistic quantum field theory models[END_REF]. The boundary equations of (3.25) are known as Ghoshal-Zamolodchikov (GZ) relations [START_REF] Ghoshal | Boundary S-matrix and boundary state in twodimensional integrable quantum field theory[END_REF]. These equations allow us to construct a matrix Ansatz in the sense of [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF]. The link between ZF and the matrix product representation was first found by [START_REF] Sasamoto | Stationary states of integrable systems in matrix product form[END_REF] and the precise relation between ZF and GZ algebras and the matrix Ansatz was given in [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF]. The key point to derive the matrix ansatz from the ZF algebra and the GZ relations is to take the derivative of (3.24) at x 1 = x 2 = 1 and of (3.25) at x = 1. This leads to

w A 1 (1)A 2 (1) = (q -1) A 1 (1)A ′ 2 (1) -A ′ 1 (1)A 2 (1) , (3.26)
W | BA(1) -(q -1)A ′ (1) = 0, and BA(1) + (q -1)A ′ (1) |V = 0. (3.27)

Using these relations, on can show [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF] that the vector defined by

|S = 1 Z L W |A 1 (1) . . . A L (1)|V , with Z L = W |C L |V and C = (1, 1, 1).A (1), (3.28)
is the steady state of the process, i.e. M |S = 0. Writing the components of this vector |S , the weights of any configuration at the stationary state are obtained. The crucial step is to construct the vector A(x). As explained in [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF], this vector is a Laurent series with respect to the spectral parameter x, but this series can be truncated and only a finite number of terms retained. We shall illustrate this procedure in the next section, devoted to the analysis of a specific example.

A specific example

In this section, we shall consider a specific example for which we shall carry out a more detailed analysis. We consider a 2-TASEP and take the set of transition rates L 2 with α = 1/2 and the set of transition rates R 3 with β = 1 (see Figure 1). In brief, the example we consider is given by the following rules:

Left Bulk Right 1 1/2 ---→ 2 21 1 --→ 12 2 1 --→ 1 1 1/2 ---→ 3 31 1 --→ 13 3 1 --→ 1 2 1/2 ---→ 3 32 1 --→ 23 (4.1)
Note that for these rules, the boundaries are permeable for all the species.

1 2 1 2 1 1 1 1 1 1
Figure 1: The 2-TASEP with special, integrable, boundary rates. The empty sites stands for species 1, circles for species 2 and bullets for species 3. On the left boundary the continuous line means injection of bullets whereas the dashed line means injection of circles.

We now construct an algebra which allows us to compute the weights of the stationary state. We will then give several quantities computed with this Ansatz.

The stationary state of the process can be expressed in a matrix product form, i.e., the weight of the configuration C = (τ 1 , . . . , τ L ) in the stationary state can be written as

P (τ 1 , . . . , τ L ) = 1 Z L W |X τ 1 X τ 2 • • • X τ L |V , (4.2)
where Z L = W |(X 1 + X 2 + X 3 ) L |V is the normalization factor. Using (3.28), we find that the operators X 1 , X 2 and X 3 are given by

A(1) =   X 1 X 2 X 3   , (4.3)
where the vector A(x) satisfies the ZF and GZ relations (3.24) and (3.25). The exchange relations among the X i 's are obtained from relations (3.26). However, the generators that appear in A ′ (1) are not necessarily scalars (they will not be scalars in the present case, see below). Therefore, more relations are required to close the algebra generated by the X i 's. A systematic way to deal with this question [START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF] is to assume the following expansion for the vector A(x):

A(x) =   x 2 + G 9 x + G 8 + G 7 /x G 6 x + G 5 + G 4 /x G 3 x + G 2 + G 1 /x + 1/x 2   , (4.4)
where the G i 's belong to a non-commuting algebra. We have observed that 9 generators is the minimal choice that ensures that no word of length 3 built from the X i 's identically vanishes. A general proof of this fact for words of arbitrary length is still missing. A way to proceed would be to construct an explicit representation of the algebra generated by the X i 's. We have obtained some partial results, representing the X i 's in terms of tensor products of DEHP generators [25], but the general scheme remains to be found. Below, we give general formulas for some special words of arbitrary length. From (4.3), we have

X 1 = 1 + G 9 + G 8 + G 7 X 2 = G 6 + G 5 + G 4 X 3 = 1 + G 3 + G 2 + G 1 (4.5)
The principle is the similar as the one introduced in [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF] except that in our case the algebra is more involved since it contains 9 fundamental generators G 1 , G 2 , . . . , G 9 from which the three matrices, corresponding to the three types of particles, X 1 , X 2 and X 3 are expressed. The algebra satisfied by the nine generators G i is found by writing each component of the ZF relation and identifying the coefficients of the polynomials in x 1 and x 2 . The generators G i satisfy a quadratic algebra, given by the following exchange relations:

[G 1 , G 2 ] = 0, [G 1 , G 3 ] = 0, [G 2 , G 3 ] = 0, G 1 G 4 = G 5 , G 2 G 4 = G 6 , G 3 G 4 = 0, [G 1 , G 5 ] = G 6 -G 4 G 2 , G 2 G 5 = G 1 G 6 , G 3 G 5 = 0, [G 4 , G 5 ] = 0, [G 1 , G 6 ] = -G 4 G 3 , [G 2 , G 6 ] = -G 5 G 3 , G 3 G 6 = 0, [G 4 , G 6 ] = 0, G 1 G 7 = G 8 , G 2 G 7 = G 9 , G 3 G 7 = 1, G 4 G 7 = 0, [G 1 , G 8 ] = G 9 -G 7 G 2 , G 2 G 8 = G 1 G 9 , G 3 G 8 = G 1 , [G 4 , G 8 ] = -G 7 G 5 , [G 1 , G 9 ] = 1 -G 7 G 3 , [G 2 , G 9 ] = G 1 -G 8 G 3 , G 3 G 9 = G 2 , [G 4 , G 9 ] = -G 7 G 6 , [G 5 , G 6 ] = 0, (4.6)
G 5 G 7 = 0, G 6 G 7 = 0, G 5 G 8 = G 4 G 9 , G 6 G 8 = G 4 , [G 7 , G 8 ] = 0, [G 5 , G 9 ] = G 4 -G 8 G 6 , G 6 G 9 = G 5 , [G 7 , G 9 ] = 0, [G 8 , G 9 ] = 0.
Remark: From the knowledge of the exchange relations for the G i 's, the algebra generated by the X i 's can be obtained using (4.5) and (4.6). We need six more generators Xi and Xi corresponding to

A ′ (1) =   X1 X2 X3   , A ′′ (1) =   X1 X2 X3   . (4.7)
Note that these generators are not scalar. The algebra generated by X i , Xi and Xi is the same as the one generated by the G i 's, written in a different basis. We chose to present the commutation relations in the G i basis because they are simpler. The action of the G i 's on the boundary vectors is derived from the GZ relations (3.25):

W | (G 4 -1) = 0, G 3 |V = 0, W | G 7 = 0, G 5 |V = 0, W | (G 8 -1) = 0, G 6 |V = 0, W | (G 1 -G 3 -G 5 ) = 0, (G 2 -1) |V = 0, W | (G 6 + G 9 -G 5 -1) = 0, (G 1 + G 4 + G 7 -G 9 ) |V = 0. (4.8)
Note that the relations are not all independent. For instance, G 3 G 5 = 0 and (G 2 -1) |V = 0 can be deduced from the other exchange and boundary relations. Using the algebraic relations (4.6)-(4.8), any expression containing X 1 , X 2 and X 3 and placed between W | and |V can be reduced to a multiple of W |V . Therefore, the weights of any configuration and the partition function of the model can be calculated. A computer program allowed us to obtain all the weights for systems of small sizes. In the rest of this section, we present some results and defer technical details and proofs to a forthcoming publication [25].

• A word containing only the letters X 2 and X 3 is evaluated as

W |Y (k) (X 2 , X 3 )X p 3 |V = k + 2 p + k + 2 2p + k + 1 p W |V , for p, k = 0, 1, 2, . . . , (4.9)
where Y (k) (X 2 , X 3 ) is the empty word for k = 0 and is a word of length k containing only the letters X 2 and X 3 and ending with X 2 .

• A word with the letter X 1 on the left is calculated as

W |X p 1 Y (k) (X 2 , X 3 )|V = 2p + 1 p + 1 2p p W |V , for p, k = 0, 1, 2, . . . , (4.10)
where Y (k) (X 2 , X 3 ) is as above.

• The partition function is given by

Z L = W |C L |V = (2L + 1)A L A L+1 W |V , (4.11)
where

C = X 1 + X 2 + X 3 and A L = 1 L+1 2L L
is the Catalan number.

• The average density of particles of type i at site k is given by n

(k) i = 1 Z L W |C k-1 X i C L-k |V .
We obtain:

n (k) 1 = 1 A L+1 k-1 i=0 A i A L-i , (4.12) n (k) 2 = 1 A L+1 L i=k L -i + 1 L + 2 A i A L-i , (4.13
)

n (k) 3 = 1 A L+1 L i=k i + 1 L + 2 A i A L-i . (4.14)
• The current j i of the particles of type i is defined as

j 1 = - 1 Z L W |C k-1 (X 2 + X 3 )X 1 C L-k-1 |V , (4.15)
j 2 = 1 Z L W |C k-1 (X 2 X 1 -X 3 X 2 )C L-k-1 |V , (4.16
)

j 3 = 1 Z L W |C k-1 X 3 (X 1 + X 2 )C L-k-1 |V . (4.17)
Using algebraic relations (4.6)-(4.8) we can show that the currents are independent of k. They have the following simple expressions

j 1 = - L + 2 2(2L + 1) , j 2 = 1 2(2L + 1)
and

j 3 = L + 1 2(2L + 1) . (4.18)
We note that densities and particle currents can be computed using the identification method of [START_REF] Ayyer | Correlations in the multispecies TASEP and a conjecture by Lam[END_REF], a procedure that can also be well-defined at the level of the algebra [25]. However, the algebraic relations (4.6)-(4.8) allow us to compute all the individual weights and correlations between different type of particles: this cannot be obtained from the identification procedure. Finally, we mention that similar algebras, with at most nine generators, can be defined to compute the weight of the stationary state of the other integrable 2-ASEP model found in this work.

Conclusion

In this work, we give a complete classification of the integrable unitary Markovian boundary conditions for the 2-ASEP with open boundaries. This classification is obtained by solving the Sklyanin reflection equations for the boundary K-matrices, and then demanding the solutions to be unitary and Markovian. We emphasize that for some choices of these boundary conditions, all particle currents are non-vanishing. Hence the system is genuinely out-of-equilibrium and can exhibit phase transitions when the boundary parameters are varied. We claim that the stationary state of all these models can be built in a matrix product form. This statement is illustrated by a detailed study of a particular open 2-TASEP model. For this model, we construct explicitly a matrix Ansatz algebra by using two key ingredients, the Zamolodchikov-Faddeev (ZF) and the Ghoshal-Zamolodchikov (GZ) relations. The algebra has nine generators and two boundary vectors; the ZF relation gives the bulk exchange relations between these generators whereas GZ gives their action on the two boundary vectors. It can be shown that these rules allow us to calculate all stationary weights. Several quantities are computed thanks to this algebra, revealing the beautiful combinatorial structure that underlies this integrable model.

Many problems remain to be solved. We have not performed here an exhaustive analysis of all the types of boundaries obtained in the classification. Some of the models may display a rich phenomenological behaviour and deserve further study [25]. Besides, it would be useful to relate the algebra generated from ZF and GZ with the tensor products of q-deformed oscillator algebras used in [START_REF] Prolhac | The matrix product solution of the multispecies partially asymmetric exclusion process[END_REF] to study the N -ASEP (it is plausible that these tensor products give explicit representations of the former algebra). More generally, integrability should allow us to calculate various physical quantities; for example, the spectral gap of the standard ASEP is known [START_REF] De Gier | Bethe Ansatz Solution of the Asymmetric Exclusion Process with Open Boundaries[END_REF]. Can one determine the gap of the 2-ASEP? Recent progress has been made in computing the eigenvectors and the eigenvalues using the algebraic Bethe Ansatz for integrable models with generic boundaries [START_REF] Belliard | Algebraic Bethe ansatz for open XXX model with triangular boundary matrices[END_REF][START_REF] Belliard | Heisenberg XXX Model with General Boundaries: Eigenvectors from Algebraic Bethe Ansatz[END_REF][START_REF] Belliard | Modified algebraic Bethe ansatz for XXZ chain on the segment-I-triangular cases[END_REF][START_REF] Crampe | Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries[END_REF]. These techniques should be generalized to the models investigated here.

Finally, a systematic investigation of integrable boundaries for various generalizations of the exclusion processes remains an open problem. Some particularly interesting cases are the following: the 2-TASEP model with unequal hopping rates [START_REF] Ayyer | On Some Classes of Open Two-Species Exclusion Processes[END_REF] (the bulk dynamics is known to be integrable [START_REF] Cantini | Algebraic Bethe Ansatz for the two species ASEP with different hopping rates[END_REF]), the homogeneous N-species ASEP (for generic N a brute force resolution of the reflection equation seems difficult but partial results were obtained using a baxterisation of a boundary Hecke algebra [20]) and the bridge model which is variant of the 2-TASEP. The bridge model, which for some choices of boundary rates displays a fascinating spontaneous symmetrybreaking [START_REF] Evans | Spontaneous symmetry-breaking in a one-dimensional driven diffusive system[END_REF][START_REF] Jona-Lasinio | Spontaneous symmetry breaking: variations on a theme[END_REF], has eluded an exact solution. We hope that the formalism developed here will lead to some progress in solving this problem.