J Noirot
email: jean.noirot@cea.fr

I Aubrun

L Desgranges

K Hanifi

J Lamontagne

B Pasquet

C Valot

P Blanpain

H Cognon

HIGH BURNUP CHANGES IN UO 2 FUELS IRRADIATED UP TO 83 GWd/t IN M5 ® CLADDINGS

Since the 90's, EDF and AREVA-NP have irradiated, up to very high burnups, lead assemblies housing M5 ® cladded fuels. Post-irradiation examination of high burnup UO 2 pellets show an increase in the fission-gas release rate, an increase in fuel swelling, and formation of fission-gas bubbles throughout the pellets. Xenon abundances were quantified, and phenomena leading to this bubble formation were identified. All examinations provided valuable data on the complex state of the fuel during irradiation. They show the good behavior of these fuels, exhibiting various microstructures at very high burnups, none of which is likely to lead to problems during irradiation.

I. INTRODUCTION

Increases in discharge burnup in light water reactor (LWR) fuels, suitable for economic benefits, imply the irradiation of lead fuel rods at very high burnups in order to gather good knowledge of the fuel performance and behavior under these conditions. Since the 90's, AREVA-NP, with the support of Electricité De France (EDF), has launched in-pile research programs on M5 ® cladded fuels used at very high burnups. Detailed post-irradiation examinations of UO 2 fuels having experienced increasing burnups during their irradiations in EDF pressurized water reactors (PWR) have been performed at Commissariat à l'énergie atomique (CEA) laboratories. Even though most of this program focused on cladding behavior, the behavior of the fuel itself was also studied, especially at high burnups. This paper presents an overview of the main observations made in the CEA LECA-STAR hot cells facility within the framework of these programs. Optical microscopy, hydrostatic density measurements, rod diameter measurements, scanning electron microscopy on polished samples as well as on fractures, Electron Probe Micro Analyzes, and Secondary Ion Mass Spectrometry, were performed.

II. FUEL DATA

The fuels described in this paper were fabricated in two batches only. This forced their initial properties to be similar.  Initial 235 U enrichment: 4.5%.  Initial porosities: ≈5%.

 Grains sizes: ≈11 µm.  Resintering tests: 0.4% of densification after one day at 1700 °C, 0.61% after two days, and 0.7% after four days.  The cladding was M5 ® .

All fuels were irradiated in the EDF Gravelines 5 reactor. The highest burnup used here was reached after 7 one-year cycles of irradiation. The highest mean rod burnup was 78.4 GWd/tU, and the highest sample burnup was 83 GWd/tU.

The linear powers during these irradiations differed from rod to rod. The range of linear powers found during these irradiations are displayed in Fig. 1. Meteor fuel behavior code calculations shows, for example, that the fuel centerline temperatures were in the range 800-900 °C during the last 3 cycles for the rods examined after 7 cycles of irradiation.

III. POST IRRADIATION EXAMINATION TECHNIQUES

Examinations used in this study were conducted in the hot cells of the CEA LECA-STAR facility in Cadarache, France and are outlined below.

 A specially designed rod-puncturing device [START_REF] Desgranges | A new puncturing apparatus optimised for free volume determination[END_REF] .

 A LWR dedicated non-destructive testing device.

 Hydrostatic density measurements employing a shielded Mettler precision balance and bromobenzene.

 An Olympus PME3 optical microscope with hot cell adaptations. Image acquisition equipment consisted of an FView II camera integrated into the SIS acquisition system.

 A shielded electronic microprobe (EPMA) Camebax microbeam from CAMECA tailored for acquisition and measurement exploitation using the SAMx system.

 A Philips XL30 scanning electron microscope (SEM) on which the WDX Microspec spectrometer used inhouse shielding. The ADA "SIS" system was used for large field acquisitions, with the AnalySIS software for image analysis.

 A CAMECA IMS 6f secondary ion mass spectrometer (SIMS) capable of analyzing the same samples as the SEM and EPMA.

For observational purposes, samples were prepared either by fracturing the fuel for SEM surface examination or by grinding and polishing the cross sections.

IV. PIE RESULTS

IV. A. Fission-gas release

Fission-gas release, obtained by rod-puncturing, as a function of the rod mean burnup, is displayed in Fig. 2. The results show an increase in fission-gas release rate with mean rod burnup, especially at high burnup.

Figure 2. Fission-gas release

The scattering in fission-gas release around 70 GWd/t is influenced more by the power step between the last two cycles than the actual power level of the last cycles. The specimen referred as the 73 GWd/tU sample in the following was cut in the rod having experienced the largest increase of power between cycles 5 and 6 and exhibiting therefore the highest release for this range of burnup.

IV. B. Fuel swelling

The fuel hydrostatic density measurements show an increase in fuel swelling rate at high burnup.

The swelling rate exceeds 0.09% per GWd/tU at the highest burnups, which is more than 60% higher than its value below 30 GWd/t (Figure 3). This increase is due to an increase of fission-gas precipitation into bubbles at high burnup. Indeed, as will be shown in the following, an increasing proportion of fission-gases form bubbles within the fuel. However, the bubble formation rate is even higher, but the swelling is limited by the low atomic volume of the gas in these bubbles. We study this fissiongas precipitation by analyzing the resulting microstructure via optical microscopy and SEM and through microanalysis by EPMA and SIMS.

Figure 3. Fuel swelling at high burnup

IV. C. Microstructure and microanalysis

In these high burnup fuels, the main microstructural changes are high burnup structure (HBS) formation and spreading, and thermal fission-gas precipitation in the central part of the pellets. Nonetheless, some other interesting features exist all over the pellets.

A macrograph of a radial cut (at 83 GWd/tU local burnup) from a 77 GWd/tU rod is displayed in Fig. 5. This macrograph was taken after chemical etching revealed both grain boundaries and small fission-gas bubbles. Concentric areas are evident at this low magnification. At the periphery (very dark ring) is the HBS. Inside this is a much thicker inhomogeneous (less) dark area that extends down to 0.57R (where R is the radius, and 1R is defined as the periphery). A brighter ring spanning 0.57-0.43R clearly delineates the external limit of the central fission-gas precipitation region. This central area itself is inhomogeneous, being brighter at its centre than at its periphery.

These radial differences in the grey levels are mainly due to factors such as the bubble size and density, some grain defects density, and, especially for the rim area, the size of the grains themselves. Consequences of these variations on the EPMA xenon measurements can be seen in Figure 6. Here, EPMA xenon concentration measurements are compared to a creation profile deduced from the EPMA neodymium measurements along the same radius of the sample. The positions where the measured values are much lower than the creation may partially correspond to xenon release from the pellet, but they are mainly places where much of the fission-gases have precipitated into bubbles. Opening of surface bubbles during polishing significantly influences the measurements because the analysis process involves less than one micrometer at the surface of the fuel 2,3 . EPMA measurements are very low at the periphery and increase rapidly and approach the creation values as we move towards the center-they do not really reach the creation values. The measurements can be very different from one point to the other even in the same area. In the These measurements correspond to the brightest ring in Figure 5. The central area (0.43-0R) also exhibits variation, the peripheral darker area corresponds to a higher EPMA measurement than the brighter centre, with a limit around 0.2R. SIMS measurements [START_REF] Desgranges | A method for the quantification of total xenon concentration in irradiated nuclear fuel with SIMS and EPMA[END_REF] complement the EPMA measurements that give the proportion of xenon precipitated into bubbles and not released from the pellet. In Figure 7 and 8, the black triangles profiles are SIMS results corresponding to Figure 6 EPMA results. All areas where a large part of the xenon is missing in the EPMA measurements are also areas where a large part and a large amount of xenon is measured inside bubbles by the SIMS. Figure 7 also indicates the presence of bubbles throughout the pellet. Shown on these same graphs are the SIMS measurements for both 61 and 73 GWd/tU samples as well as the 83 GWd/tU sample. These measurements show, for all radial positions, an increase (with burnup) in the proportion of xenon inside bubbles and an increase of the amount of xenon inside these bubbles in spite of fissiongas release. These increases are particularly clear in the central area and in the intermediate range of 1-0.57R. We now go into detailed findings working from the periphery towards the centre. Figure 9 first shows the increase (with burnup) of the highly developed HBS area. At 61 GWd/tU, in spite of an excess of 95 GWd/tU at periphery and of 80 GWd/tU 120 µm from the edge of the pellet, the extent of the highly developed HBS is still very limited (≈20 µm). It is wider for the higher burnup samples, but it is still not very wide (≈82 and 92 µm). At these limits, the local burnup is higher than 100 GWd/tU. This has to be understood with the coming high burnup post irradiation examination results. The relatively low HBS area thicknesses at high burnup might be associated with the steep plutonium increase at the periphery-the plutonium has an effect on HBS formation [START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF] .

Figure 10 shows a SEM image of the 83 GWd/tU sample periphery. The interpenetration of the internal zirconia and of the HBS is clear. The strong bonds between the pellet and the clad through this interpenetration induce cracks in the HBS. This is a likely source of fission-gas release, but has not been borne out by SIMS measurements. The atomic volume of the gas inside the HBS bubbles slightly increases with burnup but remains lower than 170 Å 3 /atom-i.e. less than 4.2 times the trivacancy volume.

Figure 9 and 10 show the HBS internal limit to be fuzzy. SEM fractographs (Figure 11) show the transition between the highly restructured HBS and the almost intact microstructure a few micrometers further. They highlight the bubbles measured by the SIMS in the range1-0.57R.

In fact, in this area, two phenomena are evident:  The first is xenon depletion around the grain boundaries. Clearly visible in Figure 9, at 61 GWd/tU, it is much wider in Figure 13 at 83 GWd/tU.  The other phenomenon is partially intragranular bubble formation.

These two phenomena are already active at 61 GWd/t, though the amount of gas involved in the precipitations at this burnup remains limited (≈10% at 0.9R and 5% at 0.7R). They already cover the same areas, reaching the central precipitation area for the first one, and up to 0.57R for the second. Altogether, they imply ≈40% of the local gas in the 83 GWd/tU fuel. Figure 12 displays a xenon X-ray mapping of the 83 GWd/tU cross section near 0.73R. It exhibits both of the above phenomena: the first one delineates the grains and the second leads to larger dark areas. This mapping explains the scattering seen in Figure 6.

In the range 0.57-0.43R, the second phenomenon is apparently, not present, and the amount of xenon inside bubbles greatly decreases. Some intergranular bubbles are detected by SEM, but we suspect an actual partial release from this area as the amount of xenon measured in the bubbles by SIMS does not compensate the EPMA lacks on the xenon profile in this area.

In the central area, the SEM and optical microscopy details (Figure 14-20) explain the differences both on the macrograph after chemical etching (Figure 5) and on the EPMA xenon measurement profile (Figure 6).

The SEM fractograph of the 73 GWd/tU sample at 0.34R in Fig. 14 presents an intragranular view, whereas Figure 15 is an intergranular view in the same area. Optical microscopy (Figure 19) also shows this periphery of the central area after chemical etching of the 83 GWd/tU cross section. Many small intragranular bubbles and metallic precipitates are present. Larger metallic fission-product precipitates are also evident in intergranular position associated with large bubbles. This microstructure is close to the 61 GWd/tU observations near 0R (Figure 16). The fractographs in Figure 17 18 and the optical micrographs in Fig. 20 show the centers of the 73 GWd/tU and 83 GWd/tU pellets, respectively. Only a few small intragranular bubbles are present, which explains the brighter grey level at the centre on the macrograph. Most of the bubbles are intergranular and large. They are associated with large metallic precipitates. Moreover, these large bubbles are neither spherical nor lenticular, and they do not seem to be widely interconnected. It is not so easy to see the difference between interand intragranular surfaces in these fractographs. Most of the gas retained inside the fuel at the center is in fact located inside these bubbles. Nevertheless, the bubbles are no longer highly over-pressurized, as they can be at 61 GWd/tU 2 . Indeed, our evaluation of the atomic volume of the gas inside the bubbles is about 195 Å 3 /atom at 0.24R and 321 Å 3 /atom at the center (4.7 and 7.9 times the trivacancy volume), whereas our evaluations 2 were only slightly superior to the trivacancy volume.

V. CONCLUSION

With the increasing fuel burnup, a number of observations were made:  An increase in fission-gas release rate and, after 6 cycles, a power-step level effect between cycle 5 and cycle 6 higher than the last power level effect.  An increase in swelling of the fuel borne out by hydrostatic density measurements and through the observed mean increase of the rod diameters after gap closure.  Formation of fission-gas bubbles all over the surface of the pellet.

When compared with other high burnup fuel results [START_REF] Manzel | EPMA and SEM of fuel samples from PWR rods with an average burn-up of around 100 MWd/kgHM[END_REF][START_REF] Goll | UO 2 Fuel Behaviour at Rod Burnups up to 105 MWd/kgHM A Review of 10 Years of High Burn-up Examinations Commissioned by AREVA NP[END_REF][START_REF] Romano | Evolution of porosity in the high-burnup fuel structure[END_REF] , the state and extent of the HBS is far less advanced here. The thickness of the (completely developed) HBS, for instance, is about 92 µm at 83 GWd/tU, whereas more than 350 µm has been reported elsewhere [START_REF] Goll | UO 2 Fuel Behaviour at Rod Burnups up to 105 MWd/kgHM A Review of 10 Years of High Burn-up Examinations Commissioned by AREVA NP[END_REF] . This point has to be understood with the coming high burnup post irradiation examination results.

An increase in porosity in the pellets between 1R and 0.43R was confirmed. It remains very limited and is mainly intergranular in the range 0.43-0.57R. It therefore concerns mainly the area between 0.57R and the HBS area. This increase is attributed to two different phenomena:  The first involves a depletion of fission-gases around the grain boundaries, which brings gas to the grain boundaries. This phenomenon is observed between the central precipitation area and the HBS structure. It is the only phenomenon observed in the grange 0.43-0.57R.

 The second is the formation of additional bubbles between 0.57R and the highly developed HBS. It can be intragranular and seems in continuity with the HBS bubbles, though the initial grains structure remains mostly intact.

At 83 GWd/tU, about 40% of the locally created xenon is found in bubbles between 0.57R and the highly developed HBS.

In the range 0.43-0.57R, the smaller amount of xenon measured in the intergranular bubbles does not seem to compensate for the missing gas at the periphery of the grains, so that a small release can be expected from this area.

At the very center of the pellets, large, irregularly shaped, intergranular bubbles form and retain the main part of the unreleased gases. They are no more as highly overpressurized as they were around 60 GWd/tU for this type of irradiation.

All these examinations have provided valuable data on the complex state of fuel during irradiation. They are critical for a good understanding of following tests results from test reactors or annealing devices.

They especially show the good behavior of these fuels, and exhibit various microstructure at very high burnup, none of which is likely to lead to particular problems during irradiation.

Figure 1 .

 1 Figure 1. Linear power histories

Figure 4 .

 4 Figure 4. Zirconia formation corrected rod mean diameter in span 5. Fuel swelling was also detected in the non-destructive rod diameter measurements. Mean span 5 rod diameters as measured are shown in Fig. 4. Subtracted from those data are slight increases in diameter due to zirconia formation on the external surface of the clad. Surface increases associated with diameter increases indicate pellet surface increases of 0.058% and 0.086% per GWd/tU in the range 60-70 GWd/tU and 70-83 GWd/tU, respectively.

Figure 5 .

 5 Figure 5. Macrograph of a radial cut of a 83 GWd/tU sample. Chemical etching revealed small bubbles and grain boundaries.

Figure 6 .

 6 Figure 6. Xenon concentration as a function of radius for the 83 GWd/tU sample as measured by EPMA.

  57-0.43R, some measurements are close to the creation values, and others are not-the scatter is large.

Figure 7 .

 7 Figure 7. Proportion of the retained xenon precipitated into bubbles as a function of radius for the 61, 73 and 83 GWd/tU samples as measured by SIMS.

Figure 8 .

 8 Figure 8. Local xenon weight% measured in the bubbles as a function of radius for the 61, 73 and 83 GWd/tU samples as determined by SIMS.

Figure 9 .

 9 Figure 9. Increase of the highly developed HBS with the increasing burnup, EPMA xenon X-ray mappings at the periphery of the 61, 73 and 83 GWd/tU fuels.

Figure 10 .

 10 Figure 10. Backscattered electrons SEM image of the 83 GWd/tU fuel periphery.

Figure 11 .

 11 Figure 11. SEM fractograph of the 73 GWd/tU sample periphery.

Figure 12 .

 12 Figure 12. EPMA xenon X-ray mapping at 0.73R of the 83 GWd/tU sample.

Figure 13 .

 13 Figure 13. EPMA xenon X-ray mapping at 0.51R of the 83 GWd/tU sample.

Figure 14 .

 14 Figure 14. SEM fractograph of the 73 GWd/tU sample at 0.34R, intragranular view.

Figure 15 .

 15 Figure 15. SEM fractograph of the 73 GWd/tU sample at 0.34R, intergranular view.

Figure 16 .

 16 Figure 16. Optical microscopy after chemical etching of the 61 GWd/tU sample at 0R.

Figure 17 .

 17 Figure 17. SEM fractograph of the 73 GWd/tU sample at center of the pellet, intragranular view.

Figure 18 .

 18 Figure 18. SEM fractograph of the 73 GWd/tU sample at the center of the pellet.

Figure 19 .

 19 Figure 19. Optical microscopy after chemical etching of the 83 GWd/tU sample at 0.25R

ACKNOWLEDGMENTS

The results presented here have all been obtained within the framework of collaboration between AREVA, Electricité De France (EDF) and the CEA.


Jean-Baptiste Courbot
email: jean-baptiste.courbot@uha.fr

Vincent Mazet

Emmanuel Monfrini

Christophe Collet

Pairwise Markov fields for segmentation in astronomical hyperspectral images

Keywords: Markov Random Fields, Pairwise Markov Fields, Bayesian Unsupervized Segmentation, Blurred Hyperspectral Image Segmentation

come

Introduction

This paper considers the problem of astronomical Hyperspectral Image (HSI) segmentation, which consists in estimating the classes X from an observation Y = y. This kind of images can be provided, e.g., by the Multi-Unit Spectroscopic Explorer (MUSE) instrument [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey: I. Survey description, data reduction and source detection[END_REF]. Among others, the MUSE instrument targets the distant Universe to study very faint objects. Considering that these objects are blurred by a spread function and buried in the noise, both originating from the observational setting (instrument and atmosphere), their segmentation is particularly difficult.

The problem of faint source detection in astronomical HSI was recently addressed by dedicated methods within an hypothesis testing framework [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF][START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF], in which the absence and presence of signal are two competing hypothesis. In [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF], we accounted for spatial and spectral priors to perform a local contextual detection, and errors were handled locally with the expected Probability of False Alarm (PFA). In [START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF], the authors use a matched-filter preprocessing to tackle the spatial structure, and the errors are globally handled through robust False Discovery Rate (FDR) control. Both methods perform locally, and the handling of the Field Spread Function (FSF) relies on contextual information. An alternative can be found in [START_REF] Meillier | Nonparametric Bayesian extraction of object configurations in massive data[END_REF], where the authors use a point marked process modeling within a Monte Carlo Markov chain procedure to detect sources in astronomical HSI. This model however relies on the expected shapes of objects (ellipses), which does not matches the general cases of sources detection.

In a more general fashion, the detection can bee seen as a segmentation task. For this purpose, a Markov field modeling is relevant in order to bolster the regularity of the segmentation results. Numerous works were published for Bayesian HSI segmentation, using hidden Markov fields (HMF) models [START_REF] Schweizer | Hyperspectral imagery: Clutter adaptation in anomaly detection[END_REF][START_REF] Rellier | Texture feature analysis using a Gauss-Markov model in hyperspectral image classification[END_REF][START_REF] Li | Spectral-spatial hyperspectral image segmentation using subspace multinomial logistic regression and Markov random fields[END_REF][START_REF] Eches | Adaptive Markov random fields for joint unmixing and segmentation of hyperspectral images[END_REF][START_REF] Li | Hyperspectral image classification using Gaussian mixture models and Markov random fields[END_REF][START_REF] Xia | Spectral-spatial classification for hyperspectral data using rotation forests with local feature extraction and Markov random fields[END_REF]. Other works specifically handle astronomical multi-spectral images [START_REF] Salzenstein | Fuzzy Markov random fields versus chains for multispectral image segmentation[END_REF][START_REF] Vollmer | Simultaneous multi-band detection of low surface brightness galaxies with Markovian modeling[END_REF]. However, as far as we know, there is no existing work on the Bayesian segmentation of faint objects in astronomical HSI handling the FSF. The FSF cannot be modeled within an HMF framework. This is why we turn, in this paper, to the pairwise Markov field (PMF) model [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF] which is a generalization of HMF models.

The contributions of this paper are the following. We introduce a PMF model which allows the modeling and the unsupervized Bayesian segmentation of blurred images. This model is applied for the detection of signals in extremely noisy HSI (Section 2 and Section 3). The proposed model is evaluated on synthetic images, and it is shown that it outperforms the state-of-the-art methods from [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF][START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF] shows that these latter are outperformed by the proposed method (Section 4). Finally, the segmentation of real MUSE images provide satisfying results with regards to the target application (Section 5).

Modeling

Along this paper, a random variable (resp. vector) is noted A (resp. A), its realization a (resp. a), and its distribution p(a) (resp. p(a)).

Let Y = (Y s ) s∈S and X = (X s ) s∈S be respectively the observation process and the class process, which has to be retrieved by the segmentation task. S is the lattice of the sites s in the image, which is ruled by a 8-neighborhood system denoted (N s ) s∈S and associated to two-elements cliques. For all s ∈ S, X s takes values in the finite set Ω and Y s ∈ R Λ , Λ being the number of spectral channels of Y. 1: Dependency graph corresponding to the CoPMF models introduced in this paper. The graphs represent the distribution of (X, Y) given the center node (X s , Y s ). Roughly speaking, eliminating a graph edge is equivalent to not account for a variable in one of the conditional distributions ruling the process (X, Y). The graph correspond to a local neighborhood, restricted to 4 sites for clarity purpose.

Within the pairwise Markov field (PMF) modeling [START_REF] Pieczynski | Pairwise Markov random fields and segmentation of textured images[END_REF], the pair (X, Y) is a Markov field :

p(x, y) ∝ s∈S p(x s , y s |x N s , y N s ). (1)
We assume that (X, Y) is stationary and that Y s and Y N s are independent given (X s , X N s ):

p(x s , y s |x N s , y N s ) = p(y s |x s , x N s )p(x s |x N s y N s ) (2)
The two resulting distributions are detailed in the following subsections.

Observation model

We assume the noise to be normally distributed:

p(y s |x s , x N s ) ∼ N(γ s , Σ); (3)
where γ s ∈ R Λ and Σ ∈ R Λ×Λ are the mean and covariance matrix respectively. The parameter γ s models the convolution by the FSF f1 :

γ s = f 0 µ s + r∈N s f r µ r ; (4)
where µ s ∈ R Λ is the mean associated to the class of x s . We assume in the remaining of this paper that f is known prior to the segmentation, e.g. from instrumental calibration. In this paper, the covariance matrix is designed to model the MUSE observations (cf. Section 5): it is pentadiagonal, with the same parameters for each spectral component:

Σ = σ 2 ρ 1 ρ 2 0 0 ρ 1 ρ 2 0 0 ρ 2 ρ 1 0 0 ρ 2 ρ 1 σ 2                                                                         (5)
This configuration of Σ is set with respect to the application on MUSE observations (cf. Section 5), and can be adapted to other applications.

Convolutional pairwise Markov fields

In this paper, we investigate two possible choices for the second term of (2) :

• in the first case, we assume that X s and Y N s are independent given X N s , and that X has a Markov field distribution. We choose to use an Ising potential as it is a widespread prior for image regularization [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]:

p(x s |x N s , y N s ) = p(x s |x N s ) ∝ exp         - s ∈N s α 1 -2δ x s (x s )         (6)
where δ x s is the Kronecker function of x s and α is a model parameter controlling the "granularity" of the realization x. This model will be hereafter referred as "Convolutional Pairwise Markov Field" (CoPMF). • In the second case, we choose to penalize the gaps between the spectral mean at the site s, and the observed spectra in the neighboring sites within N s . To do so, we add a Gaussian term to (6), yielding:

p(x s |x N s , y N s ) ∝ exp         - s ∈N s α 1 -2δ x s (x s ) - 1 2         γ s - 1 |N s | s ∈N s y s         Σ -1         γ s - 1 |N s | s ∈N s y s                 (7)
where | • | is the cardinal function. This model will hereafter be referred as "Multivariate Gaussian CoPMF" or MGCoPMF. Fig. 1 illustrates the dependency graphs yielded by the two models, in comparison with the classical HMF model. The two models will be evaluated in Section 4.

Unsupervized Bayesian segmentation

Bayesian segmentation

We compare two estimators for segmentation, namely the Maximum A Posteriori (MAP) [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] :

xMAP = arg max ω∈Ω |S| p(X = ω|Y = y); (8)
The exact distributions p(x|y) is known only up to a constant. However, the segmentation remains possible, using e.g. adaptations of Iterated Conditional Modes (ICM) [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] for the MAP and Marroquin's algorithm [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF] for the MPM. These methods rely on simulation algorithms, such as the Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF].

Parameter estimation

In an unsupervised context, the unknown parameters θ = {µ, σ, ρ 1 , ρ 2 , α} are estimated prior to segmentation with dedicated methods. It is possible to use, e.g., methods inspired from the Expectation-Maximization (EM) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF], the Stochastic EM (SEM) [START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF] or the Iterative Conditional Estimator (ICE) [START_REF] Delignon | Estimation of generalized mixtures and its application in image segmentation[END_REF] algorithms. For robustness reasons [START_REF] Monfrini | Estimation de mélanges généralisés dans les arbres de Markov cachés, application à la segmentation des images de cartons d'orgue de barbarie[END_REF][START_REF] Dias | An empirical comparison of EM, SEM and MCMC performance for problematic Gaussian mixture likelihoods[END_REF], we use an adaptation of the SEM method with, at a given step q:

• simulate x (q) along p θ(q-1) (x|y)

• estimate θ (q) from (x (q) , y) with the estimators described below. Let w = (w s ) s∈S = x * f be the convolution of x by the FSF f. For µ and Σ, we use the following maximum likelihood estimators:

μ = s∈S w s y s s∈S w 2 s (10) Σ = 1 |S| -1 s∈S (y s -w s μ)(y s -w s μ) (11)
The proof for ( 10) is given at the end of the paper, and the derivation of ( 11) is straightforward. Thanks to the function invariance property of the maximum likelihood estimator [START_REF] Casella | Statistical inference[END_REF], the estimators of σ, ρ 1 , ρ 2 are obtained by :

σ = 1 Λ Λ i=1 Σi,i ; ρ1 = 1 Λ-1 Λ-1 i=1 Σi,i+1 ; ρ2 = 1 Λ-2 Λ-2 i=1
Σi,i+2 .

(12) Finally, α is estimated with the least-squares estimator of [START_REF] Derin | Modeling and segmentation of noisy and textured images using Gibbs random fields[END_REF]. To sum up, the Bayesian unsupervized detection within blurred and noisy HSI is made possible by both CoPMF and MGCoPMF models, using either the MPM or the MAP es()timators.

Numerical results

Settings and alternatives

In this section, we evaluate the method for detection of faint sources in hyperspectral images, with Ω = {0, 1}. The numerical experiments are conducted on HSI containing 60 × 60 spectra and Λ = 10 spectral channels. The noiseless image is a 2D Gaussian where the faintest intensities are set to zero (cf. Fig. 2a). The intensities are then convolved by a 3 × 3 Moffat FSF (cf. figures 2c, 2b and also [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF]). Finally, µ contains 3 non-zero coefficients, as shown in Fig. 2e.

Six methods are evaluated: 5. the detection method introduced in [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF], based on an hypothesis-testing framework. We use the same experimental settings proposed as in [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF] concerning the choice of the catalog and the P FA , and use the same FSF as in the synthetic images ;

6. the detection method proposed in [START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF] and extended in [START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF]. We perform the preprocessing with the same FSF as in the images, and set the FDR value to 0.1.

All methods are evaluated under a varying signal-to-noise ratio (SNR), defined as

SNR = 10log 10       µ 2 2 Λσ 2       ; ( 13
)
where µ the product of µ by the average of non-zero intensities in the intensity map. An insight of the results for the six evaluated method is provided in Fig. 3.

MGCoPMF and CoPMF

In Fig. 4, we report the performances of the CoPMF-based methods in terms of false classification, false positive and false negative rates. These results yield the following findings: • the four methods yield very close averaged error rates, providing for all of them less than 10% error rate for SNR > -14dB. In most cases, the MGCoPMF with the MAP yields slightly better results than the other methods, and for both MGCoPMF and MPM, the MAP estimator yields smaller error rates than the MPM estimator; • these two estimators have however distinct behaviors in terms of false positive rates. One one hand, the MAP provides less false negatives that any other methods (less than 10% when RSB > -14 dB), at the cost of having small, but non-zero false positive rates. On the other hand, the MPM provides almost no false alarm (less than 3% at all SNR), but has a non-zero false negative rate. We interpret this observation as an insight of the stability of the MPM criterion with respect to the MAP in practical implementation.

• in average, the MGCoPMF yields better results with either the smaller error rates for the MGCoPMF with the MAP, or the smallest error rate with zero-valued false alarms for the MGCoPMF with the MPM.

4.3.

Comparison with state-of-the-art alternatives Fig. 5 reports the performance comparison of the MGCoPMF model with [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF] and [START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF]. The investigation of these results reveals the following points:

• the CoPMF segmentations yield, at all SNR, the best error rate with respect to the alternatives, with similar error rates for the MAP and MPM estimators; False negative rate Figure 5: Results for the methods from [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF], from [START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF] and with the MAP [START_REF] Delignon | Estimation of generalized mixtures and its application in image segmentation[END_REF] and the MPM (9) within the MGCoPMF introduced in this paper (the red and black curves are the same as in Fig. 4). The legend is the same as in Fig. 4.

• the method from [START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF] yields fairly high error rate when SNR < -8 dB, due to an excess of non-detection required for the FDR control; • the method from [START_REF] Courbot | Extended faint source detection in astronomical hyperspectral images[END_REF] provides no clear advantage over the other detection methods. This is partly caused by the use of the 3 × 3 pixels FSF from the simulations.

Real MUSE images

Detection in astronomical hyperspectral images may reveal spatially extended, spectrally-located light emission called Lyman-alpha, tracing the circum-galactic medium (CGM). The CGM may be occurring in the vicinity of known galaxies, making possible to extract from a MUSE observation (see [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey: I. Survey description, data reduction and source detection[END_REF]) a smaller HSI centered spatially and spectrally on the emission of interest. Such HSI may be contaminated by brightest objects presenting a "continuum" emission. These contamination are efficiently removed by a spectral median filtering prior to the segmentation. Fig. 6 illustrates the unsupervized segmentation of 4 MUSE HSI. In the absence of ground truth, the qualification of the results must be conducted thoroughly by experts. Nevertheless, two points are worth mentioning :

• in most cases, the MAP-based segmentation yields a larger detected region than its MPM counterpart. This is consistent with the results obtained on synthetic images; • the inspection of the average spectra reveals that their shape is consistent with a possible Lyman-alpha emission. Besides, the spatial shapes make likely the presence of an extended emission wider than a FSF blur on a point source, bolstering the hypothesis of the detection of an extended Lyman-alpha emission.

Conclusion

In this paper, we introduced the CoPMF and the MGCoPMF model, allowing to handle the convolution in the Markovian modeling of images. This model was specified for the application to detection in astronomical hyperspectral images, and the unsupervized segmentation in this context was presented. Then, comparative results with state-of-the art method indicate that the CoPMF methods outperform their state-of-the-art alternatives, with in addition complimentary results provided by the MAP and MPM estimators. We finally presented results on real MUSE images, showing that the both segmentation method are robust to a real-world use.

Acknowledgement

J.-B. Courbot acknowledges support from the ERC advanced grant 339659-MUSICOS. This work was funded in part by the DSIM project under grant ANR-14-CE27-0005. The authors would like to thank R. Bacon (CNRS-CRAL) for providing the MUSE data, and R. Bacher (CNRS-Gipsa-lab) for providing the source code for [START_REF] Bacher | Robust Control of Varying Weak Hyperspectral Target Detection With Sparse Nonnegative Representation[END_REF][START_REF] Bacher | Global error control procedure for spatially structured targets[END_REF] and for the fruitful discussions on the detection problem.

Proof : maximum likelihood estimator of µ

Let us derive the expression of the maximum likelihood estimator (MLE) of µ, stated in [START_REF] Derin | Modeling and segmentation of noisy and textured images using Gibbs random fields[END_REF]. The log-likelihood is written: log (x, y; µ, Σ) = c -1 2 s∈S (y sw s µ) Σ -1 (y sw s µ). [START_REF] Li | Spectral-spatial hyperspectral image segmentation using subspace multinomial logistic regression and Markov random fields[END_REF] where c = -|S| 2 log (2π) Λ det(Σ) . The MLE of µ, noted μ, verifies:

∂ log (x, y, µ, Σ) ∂µ µ= μ = 0 [START_REF] Li | Hyperspectral image classification using Gaussian mixture models and Markov random fields[END_REF] which can be rewritten using [START_REF] Li | Spectral-spatial hyperspectral image segmentation using subspace multinomial logistic regression and Markov random fields[END_REF] as:

0 = - 1 2 ∂ ∂µ         s∈S y s Σ -1 y s -2w s µ Σ -1 y s + w 2 s µ Σ -1 µ         µ= μ 0 = - 1 2 s∈S -2w s Σ -1 y s + 2w 2 s Σ -1 μ. ( 16
)
Assuming that

Figure

  Figure1: Dependency graph corresponding to the CoPMF models introduced in this paper. The graphs represent the distribution of (X, Y) given the center node (X s , Y s ). Roughly speaking, eliminating a graph edge is equivalent to not account for a variable in one of the conditional distributions ruling the process (X, Y). The graph correspond to a local neighborhood, restricted to 4 sites for clarity purpose.

  (a) Intensity map. (b) Intensity map after convolution. (c) 3 × 3 Moffat FSF. (d) Spectral mean of a realization Y = y. Noiseless (red) and noisy (gray) realizations of y s .

Figure 2 :

 2 Figure 2: Synthetic image formation.

1 .

 1 the proposed CoPMF model with the MAP (

Figure 3 :

 3 Figure 3: Results insights under low SNR for synthetic images. Red circles represent the ground truth (cf. Fig. 2a).

Figure 4 :

 4 Figure 4: Average performances of the MGCoPMF and CoPMF models using the MAP (8) and the MPM (9). Each point is an average results, and the shaded regions are between the first and fourth quartiles. The results were obtained with 100 realizations Y = y.

  MGCoPMF + MAP Courbot et al. [2] Bacher et al. [4] -20 -18 -16 -14 -12 -10 -8 -6 -4

Figure 6 :

 6 Figure 6: Results on 4 real MUSE images. Left : spectral average of Y = y (inverse video) after the spectral median subtraction has removed the brightest objects. Center : detection results with the MGCoPMF model, using the MPM (gray) and MAP (red) criteria. Right : average spectra in the two detected regions, with intensities in 10 -20 erg s -1 cm -2 Å -1 as a function of the wavelength (in Ångström).

Note that since since (X, Y) is stationary, f does not depend on s.


N Crampe

M R Evans

K Mallick
email: kirone.mallick@cea.fr

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

M Vanicat

Matrix product solution to a 2-species TASEP with open integrable boundaries

We present an explicit representation for the matrix product ansatz for some two-species TASEP with open boundary conditions. The construction relies on the integrability of the models, a property that constrains the possible rates at the boundaries. The realisation is built on a tensor product of copies of the DEHP algebras. Using this explicit construction, we are able to calculate the partition function of the models. The densities and currents in the stationary state are also computed. It leads to the phase diagram of the models. Depending on the values of the boundary rates, we obtain for each species shock waves, maximal current, or low/high densities phases.

Introduction

The totally asymmetric exclusion process (TASEP) has proved itself to be a paradigmatic model for non-equilibrium, current-carrying states [START_REF] Chou | Non-equilibrium statistical mechanics: From a paradigmatic model to biological transport[END_REF][START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF]. The model comprises particles hopping stochastically along a one-dimensional lattice with hard-core exclusion interactions. The version with open boundary conditions, where particles enter and leave at the boundaries, has revealed the existence of boundary induced phase transitions [START_REF] Krug | Boundary-induced phase transitions in driven diffusive systems[END_REF][START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF][START_REF] Schütz | Phase transitions in an exactly soluble one-dimensional exclusion process[END_REF]. An exact solution of the stationary state is provided by a matrix product Ansatz [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF][START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF] which allows the stationary density profiles and all equal-time correlation functions to be computed. The combinatorial nature of this solution has been explored in various works [START_REF] Angel | The stationary measure of a 2-type totally asymmetric exclusion process[END_REF][START_REF] Brak | A combinatorial derivation of the PASEP stationary state[END_REF][START_REF] Duchi | A combinatorial approach to jumping particles[END_REF][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]. Moreover, recent works have explored the integrability of the model and allowed the spectrum of the Markov matrix governing the dynamics to be established via the Algebraic Bethe Ansatz [START_REF] De Gier | Bethe Ansatz Solution of the Asymmetric Exclusion Process with Open Boundaries[END_REF][START_REF] Crampe | Algebraic Bethe Ansatz for the totally asymmetric simple exclusion process with boundaries[END_REF].

The generalisation to several species of particles has been under intense study, initially motivated by the case of second-class particles. Second-class particles hop forward when there is a vacancy immediately in front of them but may be overtaken by a first-class particle immediately behind them, in which case they hop back one lattice site. Second-class particles have proven to be a useful tool in tracking the positions of microscopic shocks in the system [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]. In the case of periodic boundary conditions the exact stationary state of a system of first and second class particles has been established in a matrix product formulation [START_REF] Derrida | Exact solution of the totally asymmetric exclusion process: shock profiles[END_REF]. Inspired by work in the probabilistic literature which invoked a queueing interpretation [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], the solution was later generalised to a multi-species hierarchy [START_REF] Evans | Matrix Representation of the Stationary Measure for the multi-species TASEP[END_REF][START_REF] Mallick | Exact solution of an exclusion process with three classes of particles and vacancies[END_REF][START_REF] Prolhac | Matrix product solution of the multi-species partially asymmetric exclusion process[END_REF][START_REF] Arita | Recursive structures in the multi-species TASEP[END_REF]. The resulting solution is constructed from tensor products of the fundamental matrices used in the two-species solution. The Bethe integrability of the multi-species problem has been established through co-ordinate Bethe Ansatz for the two-species case [START_REF] Alcaraz | Exact solution of the asymmetric exclusion model with particles of arbitrary size[END_REF][START_REF] Derrida | Bethe Ansatz Solution for a Defect Particle in the Asymmetric Exclusion Process[END_REF][START_REF] Cantini | Algebraic Bethe Ansatz for the two species ASEP with different hopping rates[END_REF] and by Algebraic Bethe Ansatz for the multi-species case [START_REF] Arita | Spectrum of a multi-species asymmetric simple exclusion process on a ring[END_REF][START_REF] Arita | Generalized matrix Ansatz in the multi-species exclusion process -the partially asymmetric case[END_REF].

Various examples of open boundary conditions with two or more species of particles have been considered initially motivated by boundary-induced symmetry breaking transitions [START_REF] Evans | Spontaneous Symmetry Breaking in a One Dimensional Driven Diffusive System[END_REF]. Matrix product solutions have been found in the two-species case where the stationary current of secondclass particles vanishes, which occurs under certain symmetric boundary conditions [START_REF] Evans | Asymmetric exclusion model with two species: spontaneous symmetry-breaking[END_REF][START_REF] Derrida | Chapter 14 in Nonequilibrium Statistical Mechanics in One Dimension[END_REF] or when the second-class particles are constrained to remain in the system through semi-permeable boundaries [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] Uchiyama | Two-Species Asymmetric Simple 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].

However, the case of general open boundary conditions for a system of first and second-class particles has remained an elusive problem. Recently progress was made in identifying classes of boundary conditions for which the system is integrable, in the sense of the Algebraic Bethe Ansatz [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF]. In that work the algebraic structure of the stationary state was worked out in a particular case of integrable boundary conditions.

In the present work we consider the manifold of integrable boundary conditions established in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF] and show that these boundary conditions lead to a non-trivial phase diagram with phases that manifest non-zero currents of the second-class particles. We find an explicit solution for the stationary state in terms of tensor products of the fundamental matrices which appear in the two-species solution for periodic boundary conditions. Interestingly, the stationary state has various factorisation properties and these allow us to compute exactly the partition function for the system. The proof of the stationary state is simplified by the factorisation properties and allows us to make the connection with techniques used in the study of integrable systems such as Zamolodchikov-Faddeev and Ghoshal-Zamolodchikov equations.

The paper is organised as follows. In section 2 we review the integrable boundary conditions 1 established in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF] and derive the corresponding phase diagrams. We then present in section 3 the tensor product solution for the stationary state and present various factorisation properties, and use this solution to compute the partition function for various cases. In section 4 we prove that the tensor product state is stationary under the Markov matrix that generates the dynamics.

In section 5, we elucidate the connection to techniques of Algebraic Bethe Ansatz. We conclude in section 6 with an overview of open problems. Finally, supplementary technical details are presented in appendices.

Definition of the models and Phase Diagrams

The two-species totally asymmetric exclusion process (2-TASEP) is a stochastic dynamical system, defined on a one-dimensional lattice with L sites in contact with two boundary reservoirs, where each site i = 1, . . . , L can be in one of three states σ i = 0, 1 or 2. State 0 may be considered as an empty site or hole. State 1 corresponds to a first class particle, and state 2 corresponds to a second class particle. In the bulk, at each pair of nearest neighbor sites, the rates of exchange are 1 0

1 -→ 0 1 , 2 0 1 -→ 0 2 , 1 2
1 -→ 2 1 . (2.1)
(We remark that various other conventions for the labelling of the three particle states have been employed in the literature e.g. [START_REF] Evans | Spontaneous Symmetry Breaking in a One Dimensional Driven Diffusive System[END_REF][START_REF] Evans | Asymmetric exclusion model with two species: spontaneous symmetry-breaking[END_REF].) The sites 1 and L are in contact with boundary reservoirs and particles are exchanged at different rates at the boundaries. For generic values of these boundary rates, the system is not integrable (in contrast with the 1-species TASEP, which is integrable for arbitrary boundary rates); finding an exact solution looks hopeless. However, in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF], using a systematic procedure, all possible boundary rates for the 2-species TASEP that preserve integrability were classified. Amongst such models, some had been studied earlier: the first open two-species matrix product solutions were derived in [START_REF] Evans | Asymmetric exclusion model with two species: spontaneous symmetry-breaking[END_REF]; in [START_REF] Derrida | Chapter 14 in Nonequilibrium Statistical Mechanics in One Dimension[END_REF] the boundary conditions for which the stationary state may be expressed using the matrices D, E, A of [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF][START_REF] Derrida | Exact solution of the totally asymmetric exclusion process: shock profiles[END_REF] were deduced; in [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] the restricted class of semi-permeable boundaries, in which second class particles can neither enter nor leave the system was studied. In all of these cases a matrix product representation of the stationary state was found involving the quadratic algebra used by Derrida, Evans, Hakim and Pasquier [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF] in their exact solution of the 1-species exclusion process with open boundaries.

In the present work, we construct a matrix Ansatz for integrable 2-TASEP with open boundaries that allow all species of particles to enter and leave the system. The algebraic structures required will be much more involved than the fundamental quadratic algebra of [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF].

We shall study two classes of 2-species TASEP models with the following boundary rates left boundary right boundary

2 α ---→ 1 2 β ---→ 0 (P 1 ) : 0 α ---→ 1 1 β ---→ 0 0 1-α --→ 2 1 1-β --→ 2 (2.2) or left boundary right boundary 2 α ---→ 1 2 β ---→ 0 (P 2 ) : 0 α ---→ 1 1 β ---→ 0 0 1-α --→ 2 (2.3)
Hereafter, the two different models will be denoted by (P 1 ) and (P 2 ). Note that in the classification of [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF] the left boundary conditions were referred to as L 2 and the right hand boundary conditions for (P 1 ) or (P 2 ) were referred to as R 2 and R 3 respectively. It is a simple matter to translate our results for (P 2 ) to the case of right boundary R 2 and left boundary L 3 . The final case of right boundary R 3 and left boundary L 3 leaves the stationary state devoid of holes and thus reduces to a one-species TASEP.

The physical interpretation of the boundary conditions is as follows. In both models (P 1 ), (P 2 ) the left-hand boundary conditions correspond to a boundary reservoir containing only first and second class particle with densities α and 1 -α respectively, with no holes. In model (P 1 ) the right-hand boundary conditions correspond to a reservoir containing second-class particles and holes with densities 1 -β and β respectively, with no first-class particles. In model (P 2 ) the right-hand boundary conditions correspond to a reservoir containing first-class particles and holes with densities 1 -β and β respectively, with no second-class particles.

The 2-TASEP is a finite Markov process that reaches a unique steady-state in the long time limit, in which each configuration has the stationary probability (or weight) P (σ 1 , σ 2 , . . . , σ L ). The column-vector P of length 3 L , whose components are the probabilities P (σ 1 , σ 2 , . . . , σ L ), satisfies the stationary master equation

M (3) P = 0, (2.4)
where M (3) is the 3 L × 3 L Markov matrix for the 2-TASEP system. Finding the steady state thus amounts to solving a linear system that grows exponentially with the size of the system. Basically, the matrix product representation of the stationary weights, based on integrability, will allow us to reduce this exponential complexity to a polynomial computation.

Phase diagrams

The stationary state of the exclusion process can exhibit different qualitative features and different analytical expressions for macroscopic quantities in the infinite size limit, L → ∞. The system is said to exhibit various phases, that depend on the values of the boundary exchange rates. These different phases can be discriminated by the values of the currents and by the shapes of the density profiles. More refined features, such as correlations length or even dynamical behaviour, can even lead us to define subphases (see [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF] for details and references). The phase diagram of the one-species TASEP has been well-known for a long time; first determined using a mean-field approximation [START_REF] Krug | Boundary-induced phase transitions in driven diffusive systems[END_REF][START_REF] Derrida | An exact solution of a one-dimensional asymmetric exclusion model with open boundaries[END_REF], it was rigorously established and precisely investigated after the finding of the exact solution [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF][START_REF] Schütz | Phase transitions in an exactly soluble one-dimensional exclusion process[END_REF][START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF]. We recall that the dynamical rules of the one-species TASEP are given by left boundary bulk right boundary

0 α ---→ 1 1 0 1 --→ 0 1 1 β ---→ 0 (2.5)
The phase diagram is determined by the behavior of the stationary current J and bulk density of particles in the limit L → ∞ [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF]. The different phases are detailed in table (2.6).

Region Phase Current J Bulk density

α < β, α < 1 2 Low-density (LD) α(1 -α) α β < α, β < 1 2
High-density (HD)

β(1 -β) 1 -β α > 1 2 , β > 1 2 Maximal current (MC) 1 4 1 2
(2.6)

The phase diagrams of the (P 1 ) and (P 2 ) models can be determined rigorously without having to compute exactly the steady-state probabilities. Indeed, the various phases of these two species models can extracted from the knowledge of the one-species TASEP phase diagram, by using an identification procedure formalised in [START_REF] Ayyer | On Some Classes of Open Two-Species Exclusion Processes[END_REF]. Then, 2-TASEP models (P 1 ) and (P 2 ) can both be mapped to the one-species TASEP model using two possible identifications:

1. One can identify holes and species 2 to get a one-species TASEP model for which the phase diagram is given in table (2.6). The boundary conditions read left boundary right boundary

(P 1 ) : (0, 2) α ---→ 1 1 1 --→ (0, 2) (P 2 ) : (0, 2) α ---→ 1 1 β ---→ (0 , 2) (2.7)
2. One can identify species 1 and 2 to get another version of the one-species TASEP model. In that case the two models (P 1 ) and (P 2 ) produce the same boundary conditions: left boundary right boundary

(P 1 ) & (P 2 ) : 0 1 --→ (1, 2) (1, 2) β ---→ 0 (2.8)
For the two-species TASEP, we denote by J 1 , J 2 the particle currents in the stationary state for the particles of species 1 and 2 respectively. The currents are counted positively when particles flow from the left to the right. In the same way, ρ 1 and ρ 2 denote the densities of particles of species 1 and 2 respectively. Identification 1. allows us to compute the current J 1 and the density ρ 1 , while identification 2. yields the current J 1 + J 2 ≡ -J 0 and the density ρ 1 + ρ 2 . Gathering these results, we obtain the phase diagrams depicted in figure 1.

Phase diagram of the (P 1 ) model

The phase diagram of the (P 1 ) model is displayed in Figure 1(a). It comprises four phases. Using the identification procedure, we observe that ρ 1 behaves as the density for the one-species TASEP with boundary rates (α, 1), while ρ 1 + ρ 2 behaves as the density for the one-species TASEP with boundary rates (1, β). The values of the currents in each phase (see fig 1(a)) are readily found by this identification. The behaviour of the density profile (ρ(j)) in each phase can be found in [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF].

α β 1 2 1 2 IV I III II J 1 = α(1 -α) J 2 = (α -1 2 ) 2 J 1 = 1 4 J 2 = 0 J 1 = α(1 -α) J 2 = J (III) 2 J 1 = 1 4 J 2 = -(β -1 2 ) 2 (a) Boundary conditions (P1) α β 1 2 1 2 IV I III II J 1 = α(1 -α) J 2 = (α -1 2 ) 2 J 1 = 1 4 J 2 = 0 J 1 = α(1 -α) J 2 = J (III) 2 J 1 = β(1 -β) J 2 = 0 (b) Boundary conditions (P2)
J (III) 2 = (β -α)(1 -α -β)).
Phase I: For α > 1 2 and β > 1 2 , first-class particles exhibit a maximal current, whereas the current of the second-class particles vanishes. The bulk density of first-class particles and of holes is equal to 1/2, while the number of second class particles in the bulk is vanishingly small. The density profiles of first and second class particles are characterised by power law decays to the bulk values:

ρ 1 = 1 2 + 1 2 √ πj + O 1 j 3/2 and ρ 2 = O 1 j 3/2 , (2.9)
where j is the site position on the lattice. The system is similar to the one-species TASEP in its maximal current phase.

Phase II: For α > 1 2 and β < 1 2 , none of the currents J 1 , J 2 and J 0 vanishes. This is a genuine 2-TASEP phase with boundaries permeable to all the species. The two species and the holes coexist in the bulk with non-zero bulk densities and density profiles characterised by power-law decays:

ρ 1 = 1 2 + 1 2 √ πj + O 1 j 3/2 and ρ 2 = 1 2 -β - 1 2 √ πj + O 1 j 3/2 .
(2.10)

First-class particles are in their maximal current phase. Boundary effects are long-range for species 1 and 2.

Phase III: For α < 1 2 and β < 1 2 , we obtain a 'massive' phase in which boundary effects are localized: after a finite correlation length, the system reaches its bulk behaviour,

ρ 1 = α + O 1 j 3/2 exp - j ξ and ρ 2 = 1 -α -β + O 1 j 3/2 exp - j ξ . (2.11)
The current of second-class particles J 2 vanishes along the line α = β < 1 2 and changes its sign across this line.

Phase IV: This phase, obtained for α < 1 2 and β > 1 2 , is massive for first-class particles but 'massless' (exhibiting long-range correlations characterised by power laws) for second-class particles and holes. Here again, the two species and the holes coexist in the bulk:

ρ 1 = α + O 1 j 3/2 exp - j ξ and ρ 2 = 1 2 -α - 1 2 √ πj + O 1 j 3/2 .
(2.12)

Holes are in their maximal current phase J 0 = -1/4.

Phase diagram of the (P 2 ) model

The phase diagram of the (P 2 ) model also comprises four phases, displayed in Figure 1(b). The diagram is qualitatively different from that of the (P 1 ) model. Here, ρ 1 behaves as the density for the one-species TASEP with boundary rates (α, β), while ρ 1 + ρ 2 behaves as the density for the one-species TASEP with boundary rates (1, β).

Phase I: For α > 1 2 and β > 1 2 , first-class particles exhibit a maximal current. This phase is similar to Phase I of model (P 1 ).

Phase II: This phase is obtained for β < α < 1 2 . First-class particles are in their high density phase. The bulk density of second-class particles in the bulk vanishes; moreover, the probability to find a second class particle at a distance larger than the correlation length ξ away from the boundaries, is exponentially small:

ρ 1 = 1 -β + O 1 j 3/2 exp - j ξ and ρ 2 = O 1 j 3/2 exp - j ξ .
(2.13)

Phase III: For α < β < 1 2 , the two species and the holes are simultaneously present with nonvanishing currents. The current of second-class particles J 2 is strictly positive. This phase is massive for the two classes of particles and the holes:

ρ 1 = α + O exp - j ξ and ρ 2 = 1 -α -β + O exp - j ξ . ( 2

.14)

Shock Line: This line corresponds to α = β < 1 2 . The density profiles ρ 1 and ρ 2 display a linear behaviour that reflect a coexistence between a low density and a high density regions:

ρ 1 = α + j L (1 -2α) and ρ 2 = (1 -α)(1 - j L ). (2.15)
The density profile of first-class particles takes the values α and 1 -α with a discontinuous shock between the two regions. The second-class particles have a plateau density of 1 -α to the left of the shock and zero density to the right shock. This means effectively that in the stationary state in the infinite system limit only the left reservoir is active as far as second-class particles are concerned.

Phase IV This phase, obtained for α < 1 2 and β > 1 2 , is similar to Phase IV of the (P 1 ) model.

Summary of matrix product solutions

In this section, the stationary state of the 2-TASEP with open boundaries is represented as a matrix state, for the models (P 1 ) and (P 2 ). We show that the steady-state probability of finding a given configuration (σ 1 , σ 2 , ..., σ L ) can be written as a contraction on two vectors W | and |V over a suitable algebra

P (σ 1 , σ 2 , ..., σ L ) = 1 Z L W |X σ 1 X σ 2 • • • X σ L |V , (3.1)
where Z L is a normalisation constant ensuring that P (σ 1 , σ 2 , ..., σ L ) is a probability.

In the following, we give explicit formulas for the operators X σ that generate the algebra and for the boundary vectors W | and |V . The operators X σ and the left vector W | are the same for both models (P 1 ) and (P 2 ); only the right vectors |V differ. We shall also present an important factorisation property of the Matrix Ansatz that will allow us to derive explicit expressions for the normalisation constant Z L . The proofs of these results will be given in section 4.

Explicit representation of the matrices

The Matrix Ansatz for the 2-TASEP will be constructed in terms of tensor products of the fundamental operators A, δ and ε that appear in the solution of the one-species exclusion process [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF]. These operators A, δ and ε define a quadratic algebra and satisfy

δ ε = 1 , A = 1 -ε δ , δ A = 0 , A ε = 0 . (3.2)
The relation with the operators D and E of [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF] is

δ = D -1, ε = E -1 and A = DE -ED.
We also define the following parameters

a = 1 -α α and b = 1 -β β . (3.3)
Finally, we shall need four commuting copies of the algebra (3.2), (ε n , δ n , A n ), n = 1, 2, 3, 4. A simple way to achieve this is to make four-fold tensor products:

ε 1 = ε ⊗ 1 ⊗ 1 ⊗ 1, ε 2 = 1 ⊗ ε ⊗ 1 ⊗ 1, ε 3 = 1 ⊗ 1 ⊗ ε ⊗ 1, ε 4 = 1 ⊗ 1 ⊗ 1 ⊗ ε, (3.4)
and similarly for δ n and A n .

We are now in a position to present explicit matrices for the 2-TASEP with open boundaries

X 0 = 1 + aA 1 A 2 + ε 2 δ 3 (1 + ε 4 ) + ε 2 + ε 3 + aA 1 A 2 ε 3 + ε 1 A 3 (1 + δ 4 ), (3.5)
X 2 = aδ 1 A 2 (1 + ε 4 ) + aδ 1 A 2 ε 3 + aA 2 A 3 (1 + δ 4 ), (3.6
)

X 1 = δ 2 + δ 3 (1 + ε 4 ) + 1 + δ 2 ε 3 + ε 1 δ 2 A 3 (1 + δ 4 ). ( 3

Expressions of the boundary vectors

To construct the vectors W | and |V , we first define the elementary vectors x| and |x that obey

x| ε = x x| and δ |x = x |x . (3.8)
It is known [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF] that explicit representations of such elementary vectors exit. Here we use a representation where

x|y = 1 1 -xy . (3.9)
The left boundary vector reads:

W | 1234 = 1| 1 0| 2 0| 3 0| 4 , (3.10)
where the indices indicate again which copy of the (A, δ, ε) algebra acts on the vector. To make the notation less cluttered, we shall simply write W | instead of W | 1234 . Note that the left vector is the same for the models (P 1 ) and (P 2 ).

The right boundary vector depends on the choice of the dynamics at the right boundary (i.e. on the choice of the model (P 1 ) or (P 2 )). We have We shall simply write |V for the right vector, without specifying the indices and which model we consider. This should be unambiguous from the context.

|V (P 1 ) 1234 = | b a 1 |0 2 |1

A factorisation property of the matrix ansatz

The expressions (3.5)-(3.7) for the X σ 's can be written in a factorized form which will be useful to compute the normalisation and for the proof of the matrix ansatz in section 4.

Let us consider the following matrices

L (3) =   1 + λA 1 A 2 + ε 2 δ 3 ε 2 + ε 3 + λA 1 A 2 ε 3 + ε 1 A 3 λδ 1 A 2 λδ 1 A 2 ε 3 + λA 2 A 3 δ 2 + δ 3 1 + δ 2 ε 3 + ε 1 δ 2 A 3   , L (2) = 1 + ε 4 1 + δ 4 , (3.13)
where L (3) is a 3 × 2 matrix which contains a parameter λ, while L (2) is a 2 × 1 matrix. Then, for λ = a, the following, important, relation is satisfied

  X 0 X 2 X 1   = L (3) L (2) (3.14)
where X 0 , X 2 and X 1 are the operators that perform the matrix Ansatz for the open 2-TASEP. This identity can readily be checked using equations (3.5)-(3.7). Furthermore, the operator L (3) , defined above, can be factorized into the product of a 3 × 3 matrix by a 3 × 2 matrix, as follows

L (3) = L (3) L (3) (3.15)
with

L (3) =   1 + λA 1 A 2 ε 1 ε 2 λδ 1 A 2 λA 2 0 δ 2 ε 1 δ 2 1   and L (3) =   1 ε 3 0 A 3 δ 3 1   . (3.16)
A similar type of factorisation holds for L (2) :

L (2) = L (2) L (2) (3.17)
with

L (2) = 1 ε 4 δ 4 1 and L (2) = 1 1 . (3.18)
The origin of these factorisations will be clarified in Section 4. Note that the factorisations are valid for arbitrary values of the parameter λ. We present also in Appendix A a relation between L and L.

Factorisation property of the steady-state probabilities

A compact way of writing the steady-state probabilities is to define a vector :

X =   X 0 X 2 X 1   . (3.19)
Then, the matrix Ansatz for the stationary probability vector reads

P = 1 Z L W |X ⊗ X ⊗ • • • ⊗ X |V . (3.20)
The factorisation (3.14) leads to

P = 1 Z L P (3) P (2) , (3.21)
where

P (3) = W | 123 L (3) ⊗ • • • ⊗ L (3) |V 123 and P (2) = W | 4 L (2) ⊗ • • • ⊗ L (2) |V 4 . (3.22)
Here, P (3) is a 3 L × 2 L matrix and P (2) is a 2 L -component vector so that P is a 3 L -component vector as expected. We also remark that P (2) (up to a normalisation) is identical to the steadystate vector of the one species TASEP with open boundaries. Therefore, we have

M (2) P (2) = 0 , (3.23)
where M (2) is the Markov matrix of the one-species TASEP. Finally, thanks to the factorisation properties of (3.15)-(3.18) of L (3) and L (2) , the stationary state can be further decomposed as

P = 1 Z L P (3) P (2) = 1 Z L P (3) P (3) P (2) P (2) , (3.24)
with

P (3) = W | 12 L (3) ⊗ L (3) ⊗ • • • ⊗ L (3) |V 12 (3.25)
P (3) = W | 3 L (3) ⊗ L (3) ⊗ • • • ⊗ L (3) |V 3 (3.26) P (2) = W | 4 L (2) ⊗ L (2) ⊗ • • • ⊗ L (2)
|V 4 (3.27)

P (2) = 1 1 ⊗ 1 1 ⊗ • • • ⊗ 1 1 . (3.28)
Let us note that

P (3) is a 3 L × 3 L matrix, P (3) is a 3 L × 2 L matrix, P (2)
is a 2 L × 2 L matrix and P (2) is a 2 L -component vector with constant components.

Calculation of the normalisation

We may now use the factorisation properties of the previous subsection to calculate the normalisation Z L of the stationary probabilities (3.1). The results we obtain are

Z L = a a -b Z L (α, 1) Z L (1, β)
for (P 1 )

Z L = (1 -ab) Z L (α, β) Z L (1, β) for (P 2 ) (3.29)
where Z L (α, β) is the partition function of the open one-species TASEP with injection rate α and extraction rate β. Its exact expression [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF] is given by

Z L (α, β) = a|(2 + ε + δ) L |b = a|(D + E) L |b = L p=0 p (2L -p -1)! L!(L -p)! 1 α p+1 -1 β p+1 1 α -1 β a|b .
(3.30) From the matrix Ansatz, we know that

Z L = W | (X 0 + X 1 + X 2 ) L |V . (3.31)
Using the factorisations (3.14) and (3.15), we obtain

X 0 + X 1 + X 2 = (1, 1, 1) •   X 0 X 2 X 1   = (1, 1, 1)L (3) L (3) L (2) . (3.32)
We first compute

(1, 1, 1) • L (3) = 1 + a(A 1 + δ 1 )A 2 + δ 2 , ε 1 (1 + δ 2 ) + aA 2 , 1 + ε 2 . (3.33)
Then, from the relations 1|(A + δ) = 1| and 1|ε = 1|, we deduce

1| 1 (1, 1, 1) • L (3) = 1 + aA 2 + δ 2 , 1 + aA 2 + δ 2 , 1 + ε 2 1| 1 (3.34)
This implies that the space 1 drops out (because neither L (3) nor L (2) act on it). Remarking that

1 + aA 2 + δ 2 , 1 + aA 2 + δ 2 , 1 + ε 2 L (3) = 1 + aA 2 + δ 2 , 1 + ε 2 1 1 0 0 0 1 L (3) = 1 + aA 2 + δ 2 , 1 + ε 2 1 A 3 + ε 3 δ 3 1 (3.35)
and using (A + ε)|1 = |1 and δ|1 = |1 , we have

1 A 3 + ε 3 3 1 |1 3 = |1 3 1 1 1 1 = |1 3 1 1 (1, 1) (3.36)
so that space 3 also drops out. Gathering the different results, we obtain

   Z L = 1|b/a 1 0|(2 + aA 2 + δ 2 + ε 2 ) L |0 2 0|1 3 0|(2 + ε 4 + δ 4 ) L |b 4 for (P 1 ) Z L = 1|0 1 0|(2 + aA 2 + δ 2 + ε 2 ) L |b 2 0|1 3 0|(2 + ε 4 + δ 4 ) L |b 4 for (P 2 ).
(3.37)

We conclude the derivation of (3.29) by using (3.30) and by observing that

0|(2 + aA + δ + ε) L |b = 0|b a|b Z L (α, β) (3.38)
because the operators ε = aA + ε and δ obeys the same algebraic rules as ε and δ, but now 0| is a left eigenvector of ε with eigenvalue a.

Proof of the matrix ansatz

In this section, we give an algebraic proof that the matrix Ansatz given in the previous section is indeed a representation of the steady-state probabilities of the models (P 1 ) and (P 2 ). We shall use the method of auxiliary matrices [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF][START_REF] Hinrichsen | On matrix product ground states for reactiondiffusion models[END_REF][START_REF] Rajewsky | The asymmetric exclusion process: Comparison of update procedures[END_REF] that set out a general cancellation scheme and led Krebs and Sandow [START_REF] Krebs | Matrix Product Eigenstates for One-Dimensional Stochastic Models and Quantum Spin Chains[END_REF] to a general proof (albeit not constructive) of the matrix-product form for a general class of stochastic processes.

Local update operators

The evolution rules of the exclusion process are local: a particle moves to one of its neighbouring sites. Hence, the Markov matrix of the process can be written as the sum of local operators. For the one-species TASEP with open boundaries, we have

M (2) = B (2) 1 + L-1 ℓ=1 m (2) ℓ,ℓ+1 + B (2) L , (4.1)
where the local bulk Markov matrix between site ℓ and ℓ + 1 and the boundary matrices are given by6

m (2) =     0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0     ; B (2) = -1 0 1 0 ; B (2) = 0 β 0 -β . (4.2)
Similarly, the dynamics of the 2-species TASEP is governed by the Markov matrix M (3) . It can be decomposed as

M (3) = B 1 + L-1 ℓ=1 m (3) ℓ,ℓ+1 + B L , (4.3)
with the local bulk update operator acting on nearest neighbour sites

m (3) =               . . . .
              (4.4)
where the points in the matrix stand for vanishing entries. The boundary operators read

B =   -1 0 0 1 -α -α 0 α α 0   , B =   0 β β 0 -β 1 -β 0 0 -1   , B =   0 β β 0 -β 0 0 0 -β   . (4.5)
These operators are written in the local state basis (0, 2, 1) which is the natural choice corresponding to increasing order of priority in the update rules. In equation (4.3), the subscripts indicate on which sites of the lattice the local operators act non-trivially, and the right boundary matrix B corresponds to B for the processes (P 1 ) and B for (P 2 ). As a rule, the superscripts in (4.2) and (4.3) indicate the number of possible states at a site, i.e. the number of species plus one. However, to lighten the notation we do not put a superscript (3) in the boundary matrices for the 2-TASEP, defined in (4.5).

Auxiliary matrices

We want to prove that (3.20) is a representation of the stationary vector of the 2-TASEP with open boundaries, i.e. that the master equation (2.4) is satisfied. The matrix Ansatz has a straightforward algebraic proof [START_REF] Hinrichsen | On matrix product ground states for reactiondiffusion models[END_REF][START_REF] Krebs | Matrix Product Eigenstates for One-Dimensional Stochastic Models and Quantum Spin Chains[END_REF][START_REF] Rajewsky | The asymmetric exclusion process: Comparison of update procedures[END_REF]: if one can find auxiliary operators

X ′ =   X ′ 0 X ′ 2 X ′ 1   such that X and X ′ satisfy m (3) X ⊗ X = X ′ ⊗ X -X ⊗ X ′ (4.6) W |B X = -W |X ′ and BX|V = X ′ |V (4.7)
where we recall that B is either B or B, then the stationary master equation (2.4) is satisfied for for the stationary probability vector given by (3.20).

Before giving an explicit realisation of these new operators X ′ 0 , X ′ 1 and X ′ 2 , we want to explain the notations. The auxiliary generators X ′ 0 , X ′ 2 and X ′ 1 are often denoted by a hat or a bar in the literature (and nicknamed 'hat-matrices'). However, we have purposely written them with a prime, because we shall show in section 5 that the X ′ can be constructed by taking the derivative with respect to a spectral parameter of the Zamolodchikov-Faddeev relation. Recalling that the tensor product in equation (4.6) is given by

X⊗X ′ =               X 0 X ′ 0 X 0 X ′ 2 X 0 X ′ 1 X 2 X ′ 0 X 2 X ′ 2 X 2 X ′ 1 X 1 X ′ 0 X 1 X ′ 2 X 1 X ′ 1               , (4.8)
we can spell out the formula (4.6) to obtain the quadratic relations that couple X 0 , X 2 and X 1 with the auxiliary matrices X ′ 0 , X ′ 2 and X ′ 1 .

[X i , X ′ i ] = 0 , i = 0, 1, 2 (4.9)

X 1 X 0 = X ′ 0 X 1 -X 0 X ′ 1 = X 1 X ′ 0 -X ′ 1 X 0 , (4.10)
X 2 X 0 = X ′ 0 X 2 -X 0 X ′ 2 = X 2 X ′ 0 -X ′ 2 X 0 , (4.11)
X 1 X 2 = X ′ 2 X 1 -X 2 X ′ 1 = X 1 X ′ 2 -X ′ 1 X 2 . (4.12)
In Appendix B, we give the connection between the bases presented here and the ones introduced in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF].

Explicit formulas for the auxiliary matrices

As was done for the matrices X i , we wish to express the auxiliary matrices X ′ i in terms of tensor products of the fundamental operators A, δ and ε satisfying the defining relations (3.2). We have seen, in (3.14), that X = L (3) L (2) . A similar representation for X ′ is given by

X ′ = L (3) ′ L (2) + L (3) L (2) ′ (4.13)
where L (3) and L (2) have been defined in (3.13) and

L (3) ′ =   1 ε 1 A 3 + ε 3 0 0 -δ 3 -1   and L (2) ′ = 1 -1 . (4.14)
As in (3.4), the generators (ε n , δ n , A n ), n = 1, 2, 3, 4 are commuting copies of the algebra (3.2). At first sight, the sum expression (4.13) for X ′ may seems a bit arbitrary but this form that is reminiscent of the derivative of a product will appear natural in section 5. Now, using (3.14) and (4.13), one can verify by a direct but lengthy calculation that the relations (4.6) and (4.7) are satisfied and thus prove the matrix Ansatz. However, the calculation can be simplified using some factorisations, as explained in the next subsection.

Synthetic proof of the auxiliary algebra

We have shown in Section 3.3 that the matrices L (3) and L (2) can be factorized, see equations (3.15)- (3.18). We have a corresponding property that holds for L (3) ′ :

L (3) ′ = L (3) ′ L (3) + L (3) L (3) ′ (4.15)
with

L (3) ′ =   1 ε 1 ε 2 0 0 0 0 0 0   and L (3) ′ =   0 0 0 0 -δ 3 -1   (4.16)
and where L (3) and L (3) are given in (3.16).

Using the explicit form of these operators, one can verify the following relations

m (3) L (3) ⊗ L (3) -L (3) ⊗ L (3) m (3) = L (3) ′ ⊗ L (3) -L (3) ⊗ L (3) ′ (4.17) m (3) L (3) ⊗ L (3) -L (3) ⊗ L (3) m (2) = L (3) ′ ⊗ L (3) -L (3) ⊗ L (3) ′ (4.18)
where m (2) is the local operator for the one-species TASEP, given in (4.2). Recalling that L (3) = L (3) L (3) and

L (3) ′ = L (3) ′ L (3) + L (3) L (3) ′ , we conclude that m (3) L (3) ⊗ L (3) -L (3) ⊗ L (3) m (2) = L (3) ′ ⊗ L (3) -L (3) ⊗ L (3) ′ . (4.19)
Similarly, from the factorisation of L (2) , (3.17) and (3.18), we can check that

m (2) L (2) ⊗ L (2) = L (2) ′ ⊗ L (2) -L (2) ⊗ L (2) ′ , (4.20)
where L (2) ′ has been defined in (4.14).

Combining equations (4. [START_REF] De Gier | Bethe Ansatz Solution of the Asymmetric Exclusion Process with Open Boundaries[END_REF]) and (4.20) with expressions (4.13) ends the proof of the bulk relation (4.6).

The following identities hold for L (3) W | 123 BL (3) -

L (3) B (2) = -W | 123 L (3) ′ (4.21) BL (3) -L (3) B (2) |V 123 = L (3) ′ |V 123 . (4.22)
Similar relations also exist for L (2) , we deduce that M (3) P (3) = P (3) M (2) , so that we obtain M (3) P = M (3) P (3) P (2) = P (3) M (2) P (2) = 0 , (

0| 4 B (2) L (2) = -0| 4 L (2) ′ (4.23) B (2) L (2) |b 4 = L (2) ′ |b 4 . ( 4
where we have used (3.23). This concludes the proof that P = 1 Z L P (3) P (2) is the stationary state of the 2-TASEP with open boundaries.

Zamolodchikov-Faddeev algebra and Ghoshal-Zamolodchikov relation

In this section, we show that the relations used in section 4.4 can be obtained from a more general framework. The relations introduced in this section will depend on an additional parameter, called the spectral parameter. The relations of the previous section will be recovered by setting this parameter to a specific value.

The main objects necessary in this section are similar to those used to prove the integrability of the Markov matrix M (3) (see [START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF] for historical paper or [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF] for the case treated here). We need the R-matrix encoding the bulk dynamics and the K-matrices encoding the boundaries rates.

For the 2-species TASEP, the braided R-matrix reads Ř(3) (x) = 1 + (1 -x) m (3) with the property -Ř(3) ′ (1) = m (3) . The K-matrix for the left boundary is

K(x) =      x 2 0 0 - ax(x 2 -1) xa+1 x(a+x) xa+1 0 -x 2 -1 xa+1 -x 2 -1 xa+1 1      (5.1)
and the ones for the two choices of right boundary are

K(x) =    1 x 2 -1 (b+x)x x 2 -1 (b+x)x 0 bx+1 (b+x)x b(x 2 -1) (b+x)x 2 0 0 x -2    and K(x) =    1 x 2 -1 (b+x)x x 2 -1 (b+x)x 0 bx+1 (b+x)x 0 0 0 bx+1 (b+x)x    . (5.2) One obtains - 1 2 K ′ (1) = B , 1 2 K ′ (1) = B and 1 2 K ′ (1) = B. (5.3)
We also define the Lax operators by

L (3) (x) = L (3) (x) L (3) (x)
where

L (3) (x) =   x + λA 1 A 2 xε 1 xε 2 λδ 1 A 2 λA 2 0 δ 2 ε 1 δ 2 1   and L (3) (x) =   1 ε 3 0 A 3 δ 3 /x 1/x   , (5.4)
and

L (2) (x) = L (2) (x) L (2) (x)
where

L (2) (x) = x xε 4 δ 4 1 and L (2) (x) = 1 1/x . (5.5)
Finally, we introduce the 3-component vector

X(x) = L (3) (x)L (2) (x).
(5.6)

The matrices introduced in this section depend on the supplementary parameter x (the spectral parameter). They lead to the matrices used previously to construct the matrix ansatz by remarking that

X = X(1) , X ′ = d dx X(x) x=1 , L (3) = L (3) (1) , L (3) ′ = d dx L (3) (x) x=1 , ... (5.7)
As mentioned previously, the introduction of this spectral parameter makes some of our previous definitions now appear natural. For example, the factorized form of

L (3) (x) = L (3) (x) L (3) (x) implies L (3) ′ = d dx L (3) (x) x=1 = L (3) ′ L (3) + L (3) L (3) ′ (5.8)
and reproduces relation (4.15).

The results obtained in section 4 can be deduced from the two main relations (proved in appendix C):

Ř(3) (x 1 /x 2 )X(x 1 ) ⊗ X(x 2 ) = X(x 2 ) ⊗ X(x 1 ), (5.9)
W |K(x) X(1/x) = W |X(x) and K(x)X(1/x)|V = X(x)|V , (5.10)
where K(x) is K(x) or K(x) depending on the right boundary considered. Taking the derivative of these relations w.r.t. x 1 and setting x 1 = x 2 , we recover the main relations (4.6) and (4.7) used to prove the matrix ansatz. In the context of integrable quantum field theory, an equation such as (5.9) is usually called a Zamolodchikov-Faddeev (ZF) relation and equations (5.10) are called Ghoshal-Zamolodchikov (GZ) relations. For more details about the use of these relations in the context of Markov chains see [START_REF] Sasamoto | Stationary state of integrable systems in matrix product form[END_REF][START_REF] Crampe | Integrable approach to simple exclusion processes with boundaries. Review and progress[END_REF].

Conclusion

In this paper, the stationary state of the open 2-species TASEP in the case of integrable boundary conditions (2.2), (2.3) is computed using a matrix product ansatz. We find that the 'matrices' are in fact four-fold tensor products whose generators are expressed in terms of more fundamental generators introduced to study the 1-species TASEP. However in the stationary state there are various factorisation properties which reduce the complexity of the calculations. We also show the utility of these expressions by computing exactly the normalisation factor of the stationary state.

The two-species models we have considered are distinguished by the the integrability of the boundary conditions [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF]. Here, we have further shown that the proof of the matrix product state involves relations that can be recovered from equations of the Zamolodchikov-Faddeev and Ghoshal-Zamolodchikov types.

In the case of the periodic boundary conditions, a similar tensor-product construction to that employed here has been used for the totally asymmetric and partially asymmetric simple exclusion process with N-species [START_REF] Evans | Matrix Representation of the Stationary Measure for the multi-species TASEP[END_REF][START_REF] Prolhac | Matrix product solution of the multi-species partially asymmetric exclusion process[END_REF][START_REF] Arita | Recursive structures in the multi-species TASEP[END_REF][START_REF] Cantini | Matrix product formula for Macdonald polynomials[END_REF]. The present paper takes a first step to generalizing these results to the case of integrable open boundary conditions (see [START_REF] Crampe | Integrable boundary conditions for multispecies ASEP[END_REF] for integrable boundaries of the N-species ASEP). In particular, we believe that the factorisation scheme proposed in this paper, using L and K matrices of decreasing sizes, see (3.14) and (C.8), will remain valid in the N-species case and allow the matrix ansatz for an N-species model to be constructed from that for an (N-1)-species model. Let us also mention that the queueing interpretation [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][START_REF] Evans | Matrix Representation of the Stationary Measure for the multi-species TASEP[END_REF][START_REF] Prolhac | Matrix product solution of the multi-species partially asymmetric exclusion process[END_REF][START_REF] Arita | Recursive structures in the multi-species TASEP[END_REF] provides a clear understanding of the role of the different spaces used for the fundamental generators in the case of the periodic boundary conditions. It would be of interest to see whether a similar analysis may also be carried out for the open boundary case presented here, which would then clarify the raison d'être of each of the spaces needed in our construction.

Besides, the matrix ansatz for the one-species ASEP was used in [START_REF] Lazarescu | Bethe Ansatz and Q-operator for the open ASEP[END_REF] to obtain the Baxter's Q-operator of the model. We believe that the matrix ansatz found here may be useful to construct the Q-operator for the multi-species TASEP.

Finally, another interpretation of the matrix ansatz was proposed recently for the multi-species TASEP with periodic boundary conditions [START_REF] Kuniba | Multispecies TASEP and the tetrahedron equation[END_REF], leading to new relations between the tetrahedron equation, the stochastic R-matrix and the matrix ansatz. It would be of interest to see whether the construction of [START_REF] Kuniba | Multispecies TASEP and the tetrahedron equation[END_REF] may be generalized to the cases with integrable boundaries. In turn this might reveal a 3D integrability in the matrix ansatz proposed here. choose the trivial representation for the generators in the space 1 (i.e. ε 1 = δ 1 = 1 and A 1 = 0). Changing the name of space 2 to space 3 and putting λ = 1, one establishes a link with the matrix L (3) :

L (3) ε 1 =δ 1 =1,A 1 =0,λ=1 =   1 ε 3 ε 3 0 A 3 A 3 δ 3 1 1   = L (3) 1 0 0 0 1 1 . (A.5)
The procedure to choose the trivial representation to get a simpler matrix has been used previously for the periodic case in [START_REF] Cantini | Matrix product formula for Macdonald polynomials[END_REF].

B Relation with the algebra found in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF] In [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF], the stationary state of the model (P 2 ) with α = 1 2 and β = 1 was constructed using an algebra based on 9 generators G i (i=1,2,. . . ,9). The algebra generated by the G i 's was shown to be well defined. However, no explicit realisation of that algebra was given. Using the matrix Ansatz found here, we give, in this appendix, an explicit representation of the G i 's. Using the following identification [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF]

X(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   , (B.1)
the generators X 0 , X 1 and X 2 are given by

X 0 = 1 + G 7 + G 8 + G 9 (B.
G 7 = λA 1 A 2 ε 3 + ε 2 ; G 8 = λA 1 A 2 ε 4 + ε 2 δ 3 ε 4 + λA 1 A 2 ε 3 δ 4 + ε 2 δ 4 + ε 1 A 3 + ε 3 (B.7) G 9 = λA 1 A 2 + ε 2 δ 3 + ε 4 + ε 1 A 3 δ 4 + ε 3 δ 4 (B.8)
We checked using a symbolic calculation program [START_REF] Vermaseren | New features of FORM[END_REF] that the representation presented here indeed obeys the commutation relations given in [START_REF] Crampe | Open two-species exclusion processes with integrable boundaries[END_REF].

C Algebraic proof of relations (5.9) and (5.10)

We can prove relations (5.9) and (5.10) by direct computations. However, using the factorisation (5.6), we can split the proof of these relations into simpler ones.

C.1 ZF relations

One can show that the following relations hold Ř(3) (x 1 /x 2 )L (3) (x 1 ) ⊗ L (3) (x 2 ) = L (3) (x 2 ) ⊗ L (3) (x 1 ) Ř(3) (x 1 /x 2 ), (C.1) Ř(2) (x 1 /x 2 )L (2) (x 1 ) ⊗ L (2) (x 2 ) = L (2) (x 2 ) ⊗ L (2) (x 1 ) Ř(2) (x 1 /x 2 ), (C.2) Ř(3) (x 1 /x 2 ) L (3) (x 1 ) ⊗ L (3) (x 2 ) = L (3) (x 2 ) ⊗ L (3) (x 1 ) Ř(2) (x 1 /x 2 ), (C.3) Ř(2) (x 1 /x 2 ) L (2) (x 1 ) ⊗ L (2) (x 2 ) = L (2) (x 2 ) ⊗ L (2) (x 1 ), (C.4)

where we used the braided R-matrix for the single-species TASEP, built on the local operator m (2) (see (4.2)): Ř(2) (x) = 1 + (1 -x) m (2) . (C.5)

These identities imply Ř(3) (x 1 /x 2 )L (3) (x 1 ) ⊗ L (3) (x 2 ) = L (3) (x 2 ) ⊗ L (3) (x 1 ) Ř(2) (x 1 /x 2 ), (C.6) Ř(2) (x 1 /x 2 )L (2) (x 1 ) ⊗ L (2) (x 2 ) = L (2) (x 2 ) ⊗ L (2) (x 1 ).

(C.7)

Let us remark that taking the derivative of these relations w.r.t. x 1 and setting x 1 = x 2 , we recover the relations we used previously. For instance, (C.1) implies (4.17), and (C.6) implies (4.19). Finally, using (C.6) and (C.7), we prove (5.9).

C.2 GZ relations

The following relations hold

Figure 1 :

 1 Figure 1: Phase diagrams of the 2-TASEP with open boundaries for the boundary conditions (P 1 ) and (P 2 ) (we have used the notation J

. 7 )

 7 N. B.It is important to realise that the integer suffices on the left and right hand sides of (3.5)-(3.7) are unrelated: on the left σ = 0, 1, 2 corresponds to particle species whereas on the right n = 1, 2, 3, 4 labels the tensor product as in the example (3.4).

2 )X 2 = G 4 + G 5 + G 6 (B. 3 )X 1 = 1 + G 1 + G 2 + G 3 . (B. 4 )

 22463111234 We get the following realisation for the 9 generators GG 1 = δ 3 ε 4 + ε 1 δ 2 A 3 + δ 2 ε 3 + δ 4 ; G 2 = δ 3 + δ 2 ε 4 + δ 2 ε 3 δ 4 + ε 1 δ 2 A 3 δ 4 ; G 3 = δ 2 (B.5) G 4 = λ(δ 1 A 2 ε 3 + A 2 A 3 ) ; G 5 = λ(δ 1 A 2 ε 4 + δ 1 A 2 ε 3 δ 4 + A 2 A 3 δ 4 ) ; G 6 = λδ 1 A 2 (B.6)

W

  | 123 K(x)L(3) (1/x) = W | 123 L(3) (x)K (2) (x) ; W | 4 K (2) (x)L (2) (1/x) = W | 4 L (2) (x) K(x)L (3) (1/x) |V 123 = L (3) (x)K (2) (x) |V 123 ; K (2) (x)L (2) (1/x) |V 4 = L (2) (x) |V 4 ,(C.8) where we have introduced the K-matrices for the single-species open TASEP K (2) (x) = that taking the derivative of relations (C.8) w.r.t. x and setting x = 1, we recover the relations (4.21)-(4.24).Finally, relations (C.8) imply equations (5.10).

  .24) Combining equations (4.21)-(4.24) leads us to the boundary relations (4.7). Finally, from the identities (4.19), (4.21) and (4.22)

We present here only the boundary matrix B[START_REF] Angel | The stationary measure of a 2-type totally asymmetric exclusion process[END_REF] with an injection rate α = 1 that is needed for our purposes, see relations (4.21) and(4.22). Remark that the matrices B[START_REF] Angel | The stationary measure of a 2-type totally asymmetric exclusion process[END_REF] and B[START_REF] Angel | The stationary measure of a 2-type totally asymmetric exclusion process[END_REF] correspond to identification (2.8).

Acknowledgement: M.E. thanks the C.N.R.S. and LAPTh for partial support during the completion of this work.

A Relation between L and L

In this section, we show how L can be obtained from L.

We define the transposition in the space of generators as follows

Let us remark that starting from an L (3) solution to the relation (4.17), the matrix

is also a solution of (4.17). We have used the following property of the local operator m

Starting from the realisation (5.4) for L (3) , one gets

The trivial representation for the ε, δ, A algebra is defined as ε = δ = 1 and A = 0. These values are consistent with the relation (3.2) and the definition of A. In the L (3) matrix, we may