Alexandre Lazarescu

Vincent Pasquier

Bethe Ansatz and Q-operator for the open ASEP

Keywords: numbers: 05.40.-a, 05.60.-k, 02.50.Ga, 02.30.Ik, ASEP, open boundaries, current fluctuations, integrable systems, Bethe Ansatz

In this paper, we look at the asymmetric simple exclusion process with open boundaries with a current-counting deformation. We construct a two-parameter family of transfer matrices which commute with the deformed Markov matrix of the system. We show that these transfer matrices can be factorised into two commuting matrices with one parameter each, which can be identified with Baxter's Q-operator, and that for certain values of the product of those parameters, they decompose into a sum of two commuting matrices, one of which is the usual one-parameter transfer matrix for a given dimension of the auxiliary space. Using this, we find the T-Q equation for the open ASEP, and, through functional Bethe Ansatz techniques, we obtain an exact expression for the dominant eigenvalue of the deformed Markov matrix.

The asymmetric simple exclusion process (ASEP), a one-dimensional discrete lattice gas model with hard core exclusion and biased diffusion, is one of the most extensively studied models of non-equilibrium statistical physics [START_REF] Derrida | An exactly soluble non-equilibrium system: The asymmetric simple exclusion process[END_REF][START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF][START_REF] Chou | Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport[END_REF][START_REF] Schütz | Exactly solvable models for many-body systems far from equilibrium[END_REF]. There are several good reasons for this. First of all, it is a simple, physically reasonable toy model, yet its behaviour is complex enough to be interesting. Furthermore, it has the mathematical property of being integrable, which makes it a good candidate for exact calculations. It is also related to other models, like growing interfaces [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF], the XXZ spin chain [START_REF] Sandow | Partially asymmetric exclusion process with open boundaries[END_REF], or certain random matrices [START_REF] Ferrari | From interacting particle systems to random matrices[END_REF], which makes its study relevant not only for itself, but for a wide range of physical systems.

One of the fundamental quantities to describe the behaviour of the ASEP, as a particle gas driven out of equilibrium, is the current of particles that flows through it in its steady state. That current is closely related to the entropy production in the system [START_REF] Lebowitz | A Gallavotti-Cohen Type Symmetry in the Large Deviation Functional for Stochastic Dynamics[END_REF], which is the defining characteristic of non-equilibrium systems, and is therefore of particular importance. The statistics of that current can be described through its large deviation function [START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF], which is the rescaled logarithm of its probability distribution, or, equivalently, through the generating function of its cumulants. This generating function can be expressed as the maximal eigenvalue of an integrable deformation of the Markov matrix of the system. That method was used in [START_REF] Derrida | Exact large deviation function in the asymmetric exclusion process[END_REF], in conjunction with the coordinate Bethe Ansatz [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF], to find the generating function of the cumulants of the current in the periodic totally asymmetric simple exclusion process (TASEP), a simpler variant where the particles move only in one direction, and later in the periodic ASEP [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF]. For the open ASEP, unfortunately, the presence of reservoirs makes it impossible to number the particles, which makes the coordinate Bethe Ansatz inapplicable. The generating function of the cumulants of the current was first obtained in [START_REF] De Gier | Large Deviation Function for the Current in the Open Asymmetric Simple Exclusion Process[END_REF] in the large size limit and in certain phases of the system. An exact general expression was then conjectured in [START_REF] Lazarescu | An Exact Formula for the Statistics of the Current in the TASEP with Open Boundaries[END_REF] for the TASEP and [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF] for the ASEP (through calculations [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF] relying on the famous matrix Ansatz [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF] that describes the steady state of the system), and their structure was found to be extremely similar to the periodic case, but no rigorous proof was found, and no precise explanation for that similarity. This is what we propose to do in this paper.

There is a version of the Bethe Ansatz more versatile than the coordinate one, called the algebraic Bethe Ansatz [START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF], which can in principle be used for the open ASEP. Its formulation comes from the fact that the Markov matrix of the ASEP (or the Hamiltonian of the XXZ spin chain), can be expressed in relation to the row-by-row transfer matrix of the six-vertex model [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF]. That transfer matrix commutes with the Markov matrix, and depends on a spectral parameter. It has the form of a product of local tensors, (Lax matrices), traced over an auxiliary space which is usually of dimension 2, like the physical space of a single site, although any dimension can be chosen. That product of Lax matrices, if the trace on the auxiliary space is not taken, is called the monodromy matrix of the system, and creation and annihilation operators for the particles can be extracted from it, that depend on their own spectral parameters. The algebraic Bethe Ansatz then consists in finding a trivial eigenvector of the transfer matrix (a vacuum state), on which the creation operator is then applied to obtain the other eigenvectors. Commutation relations between the transfer matrix and the creation operators, depending on their respective spectral parameters, produce the same Bethe equations as for the coordinate Bethe Ansatz, where the role of the Bethe roots is played by those spectral parameters.

In the case of a periodic system, with a fixed number of particles (or a fixed magnetisation sector), that vacuum state can be chosen as either completely empty, or completely full, which are two trivial eigenvectors of the system. In the case of an open system, there are two extra difficulties. Firstly, the transfer matrix is more complicated, and involves two rows of the vertex model instead of one, with certain reflection operators at the boundaries [START_REF] Kulish | Algebraic structures related to reflection equations[END_REF][START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]. This makes things harder, but not intractable. The second difficulty, however, does: for an open system, with no occupation sectors and no trivial eigenvectors, we do not know, in general, how to find a suitable vacuum state. Such states have been found for certain special boundary conditions, such as triangular boundary matrices [START_REF] Belliard | Algebraic Bethe Ansatz for Open XXX Model with Triangular Boundary Matrices[END_REF][START_REF] Ragoucy | Generalized Coordinate Bethe Ansatz for open spin chains with non-diagonal boundaries[END_REF], or full matrices with constraints on their coefficients [START_REF] Crampé | Eigenvectors of open XXZ and ASEP models for a class of nondiagonal boundary conditions[END_REF][START_REF] Crampé | Matrix coordinate Bethe Ansatz: applications to XXZ and ASEP models[END_REF][START_REF] De Gier | Bethe ansatz solution of the asymmetric exclusion process with open boundaries[END_REF][START_REF] Simon | Construction of a coordinate Bethe ansatz for the asymmetric simple exclusion process with open boundaries[END_REF], for which some pseudo-particles are conserved and the full construction can be performed. Recent progress has also been made for a semi-infinite chain [START_REF] Baseilhac | The half-infinite XXZ chain in Onsagers approach[END_REF][START_REF] Baseilhac | New results in the XXZ open spin chain[END_REF], which has only one boundary, through a method alternative to the Bethe Ansatz. The vacuum state for completely general boundary conditions, however, remains elusive.

There is however a way to obtain the eigenvalues of the Markov matrix (or Hamiltonian, or transfer matrix) we are interested in without having to deal with the eigenvectors at all. This can be done through Baxter's so-called 'Q-operator' method. It was first used as an alternative method to solve the six-vertex model [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF], but was later discovered to be, in fact, a limit of the transfer matrix with an infinite-dimensional auxiliary space [START_REF] Frappat | A complete Bethe ansatz solution for the open spin-s XXZ chain with general integrable boundary terms[END_REF][START_REF] Nepomechie | Functional relations and Bethe Ansatz for the XXZ chain[END_REF]. Certain algebraic relations between the Q-operator and the transfer matrix, called 'T-Q' relations, allow to obtain the functional Bethe Ansatz equations for the eigenvalues directly, without need of the eigenvectors [START_REF] Yung | Exact solution for the spin-s XXZ quantum chain with non-diagonal twists[END_REF][START_REF] Bazhanov | A shortcut to the Q-operator[END_REF][START_REF] Korff | Auxiliary matrices on both sides of the equator[END_REF][START_REF] Korff | A Q-operator for the quantum transfer matrix[END_REF][START_REF] Yang | Q-operator and relation from the fusion hierarchy[END_REF]. However, even with that method, the open case was solved only for certain constraints on the boundary parameters [START_REF] Yang | Q-operator and relation from the fusion hierarchy[END_REF][START_REF] Nepomechie | Bethe ansatz solution of the open XXZ chain with nondiagonal boundary terms[END_REF][START_REF] Murgan | Bethe ansatz derived from the functional relations of the open XXZ chain for new special cases[END_REF] (in the case of the ASEP deformed to count the current, those constraints involve all four boundary parameters, the currentcounting fugacity, and the size of the system). More recently, a variant of the T-Q equation was devised, for the XXZ spin chain [START_REF] Cao | Off-diagonal Bethe ansatz solutions of the anisotropic spin-1/2 chains with arbitrary boundary fields[END_REF][START_REF] Cao | Spin-1/2 XYZ model revisit: general solutions via off-diagonal Bethe ansatz[END_REF][START_REF] Kitanine | Open spin chains with generic integrable boundaries: Baxter equation and Bethe ansatz completeness from SOV[END_REF] as well as in the XXX case [START_REF] Nepomechie | An inhomogeneous T -Q equation for the open XXX chain with general boundary terms: completeness and arbitrary spin[END_REF][START_REF] Cao | Off-diagonal Bethe ansatz solution of the XXX spin chain with arbitrary boundary conditions[END_REF], where an extra inhomogeneous term is added to the equation. This extra term allows to define a polynomial equivalent to the Q-operator which verifies the same relation with the eigenvalue of the corresponding Hamiltonian (the usual Q-operator cannot be polynomial itself for general values of the boundary parameters, unless the aforementioned constraints are verified). While that method applies to the situation we are considering, it is for the moment unclear how it compares with the one presented in this paper. We will come back to this issue in the conclusion.

In this paper, we show that it is in fact possible to treat the most general case without introducing an inhomogeneity, by constructing explicitly the Q-operator (as well as the P-operator for the 'other side of the equator' [START_REF] Pronko | Bethe equations 'on the wrong side of the equator[END_REF]) for the open ASEP with any boundary parameters and current-counting deformation, and obtaining the functional Bethe equations for the eigenvalues of the deformed Markov matrix. We construct a transfer matrix with an auxiliary space of infinite dimension and two spectral parameters instead of one (the second of which is what is usually used as a representation parameter of the U q [SU(2)] algebra [START_REF] Chaichian | Introduction to quantum groups[END_REF] and fixed to a specific value, but we will see that it is essential to us to treat it as a free parameter). This is a natural generalisation of the transfer matrix presented in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF], with the boundary vectors playing the part of the reflection operators we mentioned. We identify that transfer matrix as the product of P and Q, and we find that, taking special values of these two spectral parameters, we can recover the usual one-parameter transfer matrices with any auxiliary space dimension, and the T-Q relations for all of those matrices, as well as the corresponding fusion rules. This is the reverse of the usual construction, where Q is obtained through the fusion of an infinite number of 2-dimensional transfer matrices.

In the first section, we consider the periodic case, where everything is known from the coordinate Bethe Ansatz, as a benchmark for our construction, and we connecting this approach to that of the functional Bethe Ansatz presented in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF][START_REF] Prolhac | Current fluctuations in the exclusion process and Bethe ansatz[END_REF], noticing that the polynomials P and Q that are constructed there are the eigenvalues of the operators P and Q that we consider here. In the second section, we apply the same method to the open ASEP, and see thet the addition of the boundaries does not modify the structure of the Q-operator. This allows us to use the same method as in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF] to obtain the expression of the generating function of the cumulants of the current in the open ASEP that was conjectured in [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. We settle, in passing, the question of how the matrix Ansatz for the steady state of the ASEP [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF] relates to the algebraic Bethe Ansatz, and we show how our construction can also be applied to the spin-1/2 XXZ chain.

I. PERIODIC ASEP

In this first section, we treat the periodic case, for which the coordinate Bethe Ansatz solution is known [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF]. By generalising the tensors X that were defined in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF], as well as the algebraic relations their elements satisfy, we construct a transfer matrix with two free parameters, which commutes with the current-counting deformed Markov matrix of the periodic ASEP for any values of those parameters. We then show that, for certain special values, the transfer matrix decomposes into two independent blocks, one of which is the one-parameter transfer matrix for some dimension of the auxiliary space. We also show that our transfer matrix is in fact the product of two one-parameter operators P and Q. Putting these results together, we are able to recover the functional Bethe equations for the periodic ASEP [START_REF] Prolhac | Current fluctuations in the exclusion process and Bethe ansatz[END_REF].

The matrix that we want to diagonalise here is the Markov matrix of the periodic ASEP, with a current-counting deformation, which is given by:

M = L i=1 M (i) with M (i) =    0 0 0 0 0 -q e µi 0 0 qe -µi -1 0 0 0 0 0   
acting on sites i and i + 1 in basis {00, 01, 10, 11}, and where M (L) connects site L with site 1. If all the µ i are taken to be 0, this gives the Markov matrix of the open ASEP (which is stochastic : the columns sum to 0). Deforming it with fugacities allows to keep track of the number of times a given jumping rate is used in a given realisation of the system, and hence to get the statistics of the particle currents. The largest eigenvalue of that deformed matrix can then be found to be the generating function for the cumulants of those currents, which is the quantity that we want to obtain. It can be shown [START_REF] Lebowitz | A Gallavotti-Cohen Type Symmetry in the Large Deviation Functional for Stochastic Dynamics[END_REF] that the spectrum of that deformed matrix depends only on the sum of all the µ i , so that we can for instance choose µ L = µ and set all the others to 0. We will denote that deformed matrix by M µ , and the corresponding highest eigenvalue by E(µ).

A. Bulk algebra and commutation relations

The starting point for the results of this paper is to realise that the matrices d and e that are defined, for instance, in [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF][START_REF] Prolhac | The matrix product solution of the multispecies partially asymmetric exclusion process[END_REF], with the algebraic relation that they satisfy, deq ed = (1q), correspond to a special representation of the U q [SU(2)] algebra (up to a simple gauge transformation that we present in section II F). In light of this, it seems natural to wonder whether a more general representation might be used, and produce different, and perhaps better, results.

Let us therefore define:

X = n 0 e d n 1 = 1 + xA e d 1 + yA , X = n0 ê d n1 = (1 -q) 2 1 -xA e -d -1 + yA .
in basis {0, 1}, corresponding to the occupancy of one of the sites, and where matrices A, d and e satisfy:

de -q ed = (1 -q)(1 -xyA 2 ), Ae = q eA, dA = q Ad. (1) (2) (3)
For practical purposes, we will be using a specific solution to these equations, given by:

A = ∞ n=0 q n ||n n|| (4) d = ∞ n=1 (1 -q n )||n -1 n|| = S -(1 -A) (5)
and

e = ∞ n=0 (1 -xyq n )||n + 1 n|| = S + (1 -xyA) (6)
where S + and S -are simply operators increasing or decreasing n by 1 (not to be confused with spin operators). We recover the simpler versions of these matrices simply by taking x = y = 0. A few remarks need to be made here. First of all, the matrix A that we have just defined plays an important role in building the matrix Ansatz for the multispecies periodic ASEP [START_REF] Prolhac | The matrix product solution of the multispecies partially asymmetric exclusion process[END_REF]. Secondly, we could have chosen for d and e their contragredient representation e = t d and d = t e, which is equivalent to a gauge transformation on d and e. We will be using this fact abundantly in the rest of the chapter. Finally, we can actually define A, S + and S -over Z rather than N, so that S + and S -are the inverse of one another: S + S -= 1 (which wouldn't work on N because of the cut at -1). Because of the term (1q n ) in d, which is 0 between states ||0 and || -1 , we are assured that, if starting from a state ||n with n ≥ 0, we can never go to one with n < 0 through any combination of d, e and A. We just need to make sure that those expressions are always applied to vectors that have non-zero coefficients only for n ≥ 0, which is enforced by the matrix A µ that we will define momentarily, on N alone. This fact will make many future calculations much easier.

All our matrices are now combinations of only A and S, which satisfy a simple algebra:

AS = q SA, (7)
with S + = S and S -= (S) -1 .

We also need to define another diagonal matrix A µ , given by

A µ = ∞ n=0 e -nµ ||n n|| (8)
such that

A µ S = e -µ SA µ . (9)
Note that this doesn't have a well-defined trace for µ = 0. If we want to consider that limit, we will have to multiply A µ by (1e -µ ) first.

Finally, we define the transfer matrix

T per µ (x, y) = Tr[A µ L i=1 X (i) ] (10)
where the product symbol refers to a matrix product in the auxiliary space (i.e. the internal space of matrices A and S + ) and a tensor product in configuration space, and the trace is taken only on the auxiliary space. The superscript (i) only indicates to which physical site each matrix X corresponds to (and their place in the product). In other terms, the weight of that matrix between configurations C = {τ i } and C = {τ i }, where τ i is the number of particles on site i in C, is given by

T per µ (x, y) C ,C = Tr[A µ L i=1 X τ i ,τi ]
Through a calculation which can be found in appendix A 1, we show that each of the local matrices in

M µ satisfy [M (i) , X (i) X (i+1) ] = X(i) X (i+1) -X (i) X(i+1) (11)
and, for M (L) , which contains the deformation,

[M (L) , X (L) A µ X (1) ] = X(i) A µ X (i+1) -X (i) A µ X(i+1) . ( 12
)
Note that these relations are in fact the infinitesimal equivalent of the so called 'RLL equation' for the commutation of matrices X with different parameters, where X is the derivative of X for a well chosen variable.

We can now recover M µ in its entirety by summing over i in [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] and [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF]. The hat matrices cancel out from one term to the next, and we are left with 0, so that:

[M µ , T per µ (x, y)] = 0. ( 13
)
T per µ (x, y) has therefore the same eigenvectors as M µ .

Note that for a general set of fugacities {µ i }, the corresponding transfer matrix is the same as the one we defined here, with matrices A µi inserted at their appropriate place in the matrix product:

T per {µi} (x, y) = Tr[A µ0 L i=1 X (i) A µi ]. (14)

B. Decomposition of the transfer matrix

Considering the representation we chose for matrix e in [START_REF] Sandow | Partially asymmetric exclusion process with open boundaries[END_REF], namely S(1 -xyA), an interesting case to consider is xy = q -k+1 with k ∈ N (which sets one coefficient to 0 in e).

Let us therefore impose y = 1/q k-1 x. The four matrices in X become:

d = ∞ n=1 (1 -q n )||n -1 n|| , e = ∞ n=0 (1 -q n-k+1 )||n + 1 n|| (15)
n 0 = ∞ n=0 (1 + q n x)||n n|| , n 1 = ∞ n=0 (1 + q n-k+1 /x)||n n|| (16)
and the coefficient of ||k k -1|| in e vanishes. This makes all these matrices lower block-triangular (n 0 and n 1 obviously are, since they are diagonal, and d already was, but not e in general) with a block of size k (for n from 0 to k -1 in the four series above) and one of infinite size (for n from k to ∞).

The coefficients of that second block happen to be the same as the coefficients of the whole matrix for x → q k x and y → q/x and in the contragredient representation of d and e (i.e. exchanging and transposing them):

t e = ∞ n=0 (1 -q n+k )||n -1 n|| , t d = ∞ n=1 (1 -q n+1 )||n + 1 n||, (17)
n 0 = ∞ n=0 (1 + q n+k x)||n n|| , n 1 = ∞ n=0 (1 + q n+1 /x)||n n||. (18)
Indeed, taking n → n + k in ( 15) and ( 16), and removing the k negative indices, we get exactly [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF] and [START_REF] Faddeev | How Algebraic Bethe Ansatz works for integrable model[END_REF]. Since the trace of a product of block-diagonal matrices is the sum of the traces of the products of the blocks, this gives us an equation for T per µ , which is one of the two results essential to our derivation of the functional Bethe Ansatz:

T per µ (x, 1/q k-1 x) = t [k] (x) + e -kµ T per µ (q k x, q/x) (19)
where the factor e -kµ is the first coefficient of A µ on the second block. The transfer matrix t [k] is the contribution coming from the first block, which can be written as:

t [k] (x) = Tr[A µ L i=1 X (i) k (x)] ( 20
)
where X k (x) contains the same entries as X(x, 1/q k-1 x), but truncated at n = k-1 (so that the auxiliary space is k-dimensional).

We saw that T per µ commutes with M µ for any values of x and y, so it also commutes with another T per µ at different values of the parameters (this is in fact not certain, because any one of those matrices could have a degenerate eigenspace, but we will assume that it is true, for now; there is a better way to show that two matrices T per µ at different values of x and y commute, and we will come back to it in the next section). This tells us that those matrices also commute with t [k] (x) for any k, and that the t [k] (x)'s with different k's commute together. This matrix equation therefore implies the same relation for the eigenvalues Λ i (x, y) of T per µ (x, y) and the eigenvalues Λ

[k] i (x) of t [k] (x).
To go further, we need to examine the first two cases in this last equation. Note that, to be consistent with our notation for X, all the matrices will be written with the physical space as their outer space and the auxiliary space as their inner space (i.e. separated in 2 × 2 blocks of matrices acting on the auxiliary space), which is opposite to the standard convention in some cases.

For k = 1, the first block of X is given by:

X 1 (x) = 1 + x 0 0 1 + 1 x .
The matrix t [1] (x), which is scalar inside of a given occupancy sector, is then given by:

t [1] (x) = (1 + x) L-N (1 + x -1 ) N = h(x). (21)
This is the same as the function h(x) that is defined in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF], up to a sign.

For k = 2, the 2 × 2 blocks from n 0 , e, n 1 and d are:

X 2 (x) =      1 + x 0 0 0 0 1 + qx 1 -1 q 0 0 1 -q 1 + 1 qx 0 0 0 0 1 + 1 x     
and t [2] (x) = Tr[A [2] µ L i=1

X (i) 2 (x)] (22)
with A

[2] µ = 1 0 0 e -µ . This matrix is, in fact, the standard transfer matrix for the periodic XXZ spin chain, with a twodimensional auxiliary space. To write it in its usual form, we need to make a few transformations and change variables. To that effect, let us consider:

1 1 + x 1 0 0 x X 2 (x) =      1 0 0 0 0 1+qx 1+x q-1 q(1+x) 0 0 x(1-q) 1+x 1+qx q(1+x) 0 0 0 0 1      =      1 0 0 0 0 qλ 1 -λ 0 0 1 -qλ λ 0 0 0 0 1      (23)
with λ = 1+qx q(1+x) , i.e. x = -1-qλ q(1-λ) . We use the symbol to signify a product in configuration space, so that it not be confused with a product in the auxiliary space (for which we use the usual product notation). This is the common form of the Lax matrix for the ASEP:

L (i) (λ) =      1 0 0 0 0 qλ 1 -λ 0 0 1 -qλ λ 0 0 0 0 1      = P (i) (1 + λM (i) )
where P (i) is a permutation matrix which has the effect of exchanging the physical space at site i with the auxiliary space (fig. -1).

M i i a a i i a a L (λ) = +λ (i) FIG. 1. Schematic representation of L (i) (λ).
The first box represents P (i) , which exchanges the occupancies of the auxiliary space a and the physical space i (the matrix is seen as acting from SE to NW). The second box is the same exchange operator, with the local matrix M acting during the exchange.

The matrix we applied to X 2 from the left in [START_REF] Ragoucy | Generalized Coordinate Bethe Ansatz for open spin chains with non-diagonal boundaries[END_REF] multiplies every entry by x for each occupied site (and since this number is conserved between the left and right entries of the transfer matrix, this operation actually commutes with t [2] , so that we haven't modified its eigenvectors). All in all, this operation multiplies t [2] by x N /(1 + x) L = 1/h(x). We therefore define:

t[2] (λ) = t [2] (x) h(x) = Tr[A [2] µ L i=1 L (i) (λ)] (24)
which is the usual transfer matrix for the periodic ASEP with one marked bond. Since L (i) (0) = P (i) is a permutation matrix, and its derivative with respect to λ at 0 is d dλ L (i) (0) = P (i) M (i) , we find that t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] (0) is the matrix that transposes the whole system back one step, and that, for the whole transfer matrix t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] (fig. -2):

M µ = t[2] (λ) -1 d dλ t[2] (λ) λ=0 = d dλ log t[2] (λ) λ=0 (25)
(we have not considered how A µ comes into play, but we can easily check that it gives the correct term in M µ ).

L -1 L 1 2 3 4 M M L -1 L 1 2 3 4 L -1 L 1 2 3 4 L -1 L 1 2 3 4 1 2 1 2 t[2] (0) = = t[2] (0) = d dλ t[2]-1 FIG.
2. Schematic representation of the value and first logarithmic derivative of t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] at 0. The first is a translation matrix which takes each site i to i -1. The second (of which only one part of the sum over sites is represented) gives the local jump matrix M (i) .

Written in terms of t [2] (x), this becomes:

M µ = (1 -1/q) d dx log t [2] (x) h(x) x=-1/q . ( 26
)
Doing the same calculations at λ = ∞ instead of 0, after a few modifications (such as taking the contragredient representation for X and multiplying everything by h(x)/h(qx)) would have instead given us:

M µ = (1 -q) d dx log t [2] (x) h(qx) x=-1 (27)
which is also what we would have obtained if we had considered a system with L-N particles, exchanging Q with P , and kept y = 1/qx as a variable. This identity is usually called the 'crossing symmetry' in the language of quantum spin chains.

C. R matrix

The next step is to show that the eigenvalues Λ i (x, y) of T per µ are in fact a product of a function of x and a function of y, which we will note as P i (x) and Q i (y). This factorisation property is a well known fact for the periodic XXX [START_REF] Derkachov | Factorization of the transfer matrices for the quantum sl(2) spin chains and Baxter equation[END_REF][START_REF] Derkachov | Factorization of the R-matrix and Baxter's Q-operator[END_REF] and XXZ [START_REF] Korff | Auxiliary matrices on both sides of the equator[END_REF][START_REF] Pronko | Bethe equations 'on the wrong side of the equator[END_REF] spin chains. This result will actually be derived in the next section, as there are a few preliminary calculations which need to be done first, mainly in order to find the R matrix of our system.

One way to go about this is through a method used in [START_REF] Bernard | Exchange Algebra and Exotic Supersymmetry in the Chiral Potts Model[END_REF][START_REF] Bazhanov | Chiral Potts model as a descendant of the six-vertex model[END_REF] which consists in introducing two new Lax matrices, defined by:

L 1 = L(a 1 , b 1 , c 1 , d 1 ) = a 1 A 1 b 1 S 1 -c 1 S -1 1 A 1 d 1 , L2 = L(a 2 , b 2 , c 2 , d 2 ) = a 2 A 2 -c 2 S 2 A 2 b 2 S -1 2 d 2 ,
where the operators A i and S j obey (9) for i = j (i.e. if they act on the same space), and commute otherwise.

We then take the product of those matrices (which is a matrix product on the physical two-dimensional space and a tensor product on the infinite-dimensional auxiliary spaces of L 1 and L2 ):

L 1 L2 = a 1 a 2 A + b 1 b 2 S 1 /S 2 b 1 d 2 S 1 -a 1 c 2 S 2 A d 1 b 2 S -1 2 -c 1 a 2 S -1 1 A c 1 c 2 S 2 /S 1 A + d 1 d 2 with A = A 1 A 2 .
We will omit the notation for the product between those new Lax matrices, since the indices are there to signify that their elements act on different spaces. We now need to consider two special cases for the coefficients of L 1 and L2 . Let us first set them as follows: a 1 = x, c 2 = y and the rest is 1. We write the corresponding matrices as L + 1 and L-2 :

L + 1 (x) L- 2 (y) = xA + S 1 /S 2 S 1 -xyS 2 A S -1 2 -S -1 1 A yS 2 /S 1 A + 1
and, by projecting each element on S 1 = S 2 = S (i.e. by applying ||i, j i + j|| to the right and its contragredient to the left), we get:

L + 1 (x) L- 2 (y) = xA + 1 S(1 -xyA) S -1 (1 -A) yA + 1 = X(x, y). (28)
Naturally, we check that A and S satisfy the correct relations.

Let us now set a 2 = x, c 1 = y/q, c 2 = q and the rest to 1. We write the corresponding matrices as L -

1
and L+ 2 :

L - 1 (y) L+ 2 (x) = xA + S 1 /S 2 S 1 -qS 2 A S -1 2 -xy/qS -1 1 A yS 2 /S 1 A + 1
and, through the same operation as before, we get:

L - 1 (y) L+ 2 (x) = xA + 1 (1 -A)S (1 -xyA)S -1 yA + 1 = X(x, y). (29)
In both of those special cases, one of the non-diagonal elements has a factor (1 -A) which allows us to truncate the representation at state ||0 and avoid some convergence issues. It would not be the case, however, if we were to construct

L + 1 (x) L+ 2 (y) or L - 1 (x) L- 2 (
y), which we will therefore avoid at any cost.

We will now use this formalism in order to find the so-called 'R matrix' which is such that:

X(x, y) X(x , y ) R(x, y; x , y ) = R(x, y; x , y ) X(x , y ) X(x, y)
where R acts on the two auxiliary spaces of both X matrices. We will comment on the use of such a matrix at the end of this section.

Considering that X(x, y) X(x , y

) = L + 1 (x) L- 2 (y)L + 3 (x ) L- 4 (
y ), we will perform this exchange of parameters in steps, exchanging the parameters of only two L i 's at a time. The first thing we could try is to exchange y and x , but this would transform L-2 (y) into L+ 2 (x ), so we would have L + 1 (x) L+ 2 (x ) on the left, which we want to avoid. The solution is then to first exchange x and y , then y and y, and finally y and x , to obtain X(x, y ) X(x , y), and then do the same once more to exchange x and x .

We first need to find

f 12 (x, y) such that L + 1 (x) L- 2 (y)f 12 (x, y) = f 12 (x, y)L - 1 (y) L+ 2 (x)
, which is to say:

X(x, y)f 12 (x, y) = f 12 (x, y)X(x, y).
That f 12 may depend on A and S. In terms of the elements of X(x, y), this writes:

[1 + xA, f 12 ] = 0, [1 + yA, f 12 ] = 0, S(1 -xyA) f 12 = f 12 (1 -A)S, S -1 (1 -A) f 12 = f 12 (1 -xyA)S -1 .
The first and second equations tell us that f 12 commutes with A (i.e. it is diagonal), so it should be a function of A alone. The third or fourth equations then give us:

S(1 -xyA)f 12 [A] = f 12 [A](1 -A)S = S(1 -qA)f 12 [qA]
(where the second equality is due to the commutation of S with A), which we can rewrite as

f 12 [A] f 12 [qA] = (1 -qA) (1 -xyA) .
Iterating this last equation, we finally find:

f 12 (xy) = (qA) ∞ (xyA) ∞ , (30)
where (x) ∞ is the infinite q-Pochhammer symbol

(x) ∞ = ∞ k=0 (1 -q k x).
To exchange the parameters back, one simply has to apply f -1 12 (xy).

This was for the exchange of parameters between the first and second or third and fourth matrices in

L + 1 (x) L- 2 (y)L + 3 (x ) L- 4 (
y ) (i.e. inside of a same X matrix). We now need to exchange parameters between the second and third matrices in that product. Let us consider L2 L 3 with c 2 = y, c 3 = y /q, and the rest set to 1:

L- 2 (y)L - 3 (y ) = (1 + yy /qS 2 /S 3 )A 2 A 3 (S 3 -yS 2 )A 2 (S -1 2 -y /qS -1 3 )A 3 S 3 /S 2 + 1 .
We need to find g - 23 (y, y ) such that L-

2 (y)L - 3 (y )g - 23 (y, y ) = g - 23 (y, y ) L- 2 (y )L - 3 (y).
As before, the commutation for each of the four elements of L-

2 (y)L - 3 (y ) is: [(1 + yy /qS 2 /S 3 )A 2 A 3 , g - 23 ] = 0, [1 + S 3 /S 2 , g - 23 ] = 0, (S 3 -yS 2 )A 2 g - 23 = g - 23 (S 3 -y S 2 )A 2 , (S -1 2 -y /qS -1 3 )A 3 g - 23 = g - 23 (S -1 2 -y/qS -1 3 )A 3 .
The first and second equations tell us that g - 23 commutes with S 2 /S 3 and with A 2 A 3 . the third and fourth suggest that g - 23 keeps A 2 and A 3 separated, so it should only depend on S 2 /S 3 . The third equation then gives:

(S 3 -yS 2 )A 2 g - 23 [S 2 /S 3 ] = g - 23 [qS 2 /S 3 ](S 3 -yS 2 )A 2 = g - 23 [S 2 /S 3 ](S 3 -y S 2 )
A 2 which can be rewritten as:

g - 23 [S 2 /S 3 ] g - 23 [qS 2 /S 3 ] = (1 -yS 2 /S 3 ) (1 -y S 2 /S 3 )
and produces, through iteration:

g - 23 (y, y ) = (yS 2 /S 3 ) ∞ (y S 2 /S 3 ) ∞ . (31)
Finally, we look for g

+ 23 (x, x ) such that L+ 2 (x)L + 3 (x )g + 23 (x, x ) = g + 23 (x, x ) L+ 2 (x )L + 3 (x)
. Setting a 2 = x, c 2 = q, a 3 = x , and the rest to 1, in L2 L 3 , we find:

L+ 2 (x)L + 3 (x ) = (xx + qS 2 /S 3 )A 2 A 3 (xS 3 -qS 2 )A 2 (x S -1 2 -S -1 3 )A 3 S 3 /S 2 + 1
and the exact same operations as before produce:

g + 23 (x, x ) = (x S 3 /S 2 ) ∞ (xS 3 /S 2 ) ∞ . ( 32
)
Let us note that g + and g -commute with the projection which we perform on L 1 L2 and L 3 L4 to get the X matrices.

We can finally forget everything about that alternative construction of X(x, y), and simply define those operators f and g ± as what we found them to be. After re-indexing them so that 1 refers to the auxiliary space of the first X matrix, and 2 to the second, we can write: with

X 1 (x, y) X 2 (x , y ) R y (x, y; x , y ) = R y (x, y; x , y ) X 1 (x, y ) X 2 (x , y), X 1 (x, y) X 2 (x , y ) R x (x, y; x , y ) = R x (x, y; x , y ) X 1 (x , y) X 2 (x, y ), (33) (34)
R y (x, y; x , y ) = f 2 (x y )g - 12 (y, y )f -1 2 (x y), R x (x, y; x , y ) = f 1 (xy)g + 12 (x, x )f -1 1 (x y). (35) (36)
where we relabelled f 12 as f 1 , f 34 as f 2 , and g ± 23 as g ± 12 , consistently with the indexes of the X matrices. Applying those one after the other, we find the full R matrix:

X 1 (x, y) X 2 (x , y ) R(x, y; x , y ) = R(x, y; x , y ) X 1 (x , y ) X 2 (x, y) (37)
with

R(x, y; x , y ) = R y (x, y; x , y )R x (x, y ; x , y) = (qA 2 ) ∞ (x y A 2 ) ∞ (yS 1 /S 2 ) ∞ (y S 1 /S 2 ) ∞ (x yA 2 ) ∞ (qA 2 ) ∞ (qA 1 ) ∞ (xy A 1 ) ∞ (x S 2 /S 1 ) ∞ (xS 2 /S 1 ) ∞ (x y A 1 ) ∞ (qA 1 ) ∞ (38)
(cf. figs. -3,4,5), or an equivalent expression if we apply R x to the left of R y instead. There are many things to be said about that R matrix. First of all, it is the rigorous way to prove that T per µ (x, y) commutes with T per µ (x , y ) for any value of those parameters: to do that, we insert R(x, y; x , y )R -1 (x, y; x , y ) at any point in the matrix product expression of T per µ (x, y)T per µ (x , y ), and make R(x, y; x , y ) commute to the left all the way around the trace (fig. -6), exchanging parameters between the two rows along the way. When crossing the marked bond, we just need to note that R commutes with

f g - f f -1 g + f -1
A µ A µ . R -1 R X X X X X X X X R R -1 X X X X X X X X = FIG.
6. Schematic representation of the commutation between R matrices and X matrices with different parameters (represented by different colours).

Equations ( 33) and ( 34) are also of interest by themselves: they tell us that T per µ (x, y) and T per µ (x , y ) can actually exchange only one of their parameters and keep the other, instead of commuting altogether (by applying the procedure we just described, but with R x or R y instead of R). This decomposition of the R matrix is a well known property [START_REF] Derkachov | Factorization of the transfer matrices for the quantum sl(2) spin chains and Baxter equation[END_REF][START_REF] Derkachov | Factorization of the R-matrix and Baxter's Q-operator[END_REF], and will be extremely useful to us in the next section.

Moreover, considering the decomposition ( 19) which we obtained in the previous section for special values of the spectral parameters in T per µ (x, y), by taking y = 1/q k-1 x and y = 1/q l-1 x in R(x, y; x , y ), we should be able to recover the R matrix between auxiliary dimensions k and l as an independent part of the whole matrix. In other words, this R matrix of twice infinite dimension should contain all the smaller R matrices for the XXZ chain. Taking y = 1/qx, in particular, should yield the Lax matrix X (possibly up to a permutation). Because of the complexity of eq.( 38), those facts are yet to be verified.

D. Q-operator and T-Q equations

From the previous section, we now know that T per satisfies:

T per (x, y)T per (x , y ) = T per (x, y )T per (x , y). (39)
Fixing x = y = 0 in that equation (although any other constants would do), we can therefore write:

T per (x, y) = P (x)Q(y) (40)
with:

P (x) = T per µ (0, 0) -1
T per µ (x, 0) , Q(y) = T per µ (0, y). Since the matrices P and Q are combinations of the matrix T per µ taken at various values of its parameters, and considering eq.( 39), we immediately see that P (x) and Q(y) commute for any values of x and y.

This relation is crucial to our reasoning, and, put together with [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF], will allow us to reach our conclusion in just a few more lines.

Using this, we can now rewrite [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF] as:

P (x)Q(1/q k-1 x) = t [k] (x) + e -kµ P (q k x)Q(q/x). (41)
The first and second orders of this equation give:

P (x)Q(1/x) = h(x) + e -µ P (qx)Q(q/x), (42)
P (qx)Q(1/qx) = h(qx) + e -µ P (q 2 x)Q(1/x), (43)
P (x)Q(1/qx) = t [2] (x) + e -2µ P (q 2 x)Q(q/x),

where we have written the first one twice (once at x, once at qx).

Combining those equations as 44) yields the T-Q relation:

Q(1/qx)×(42)+e -µ Q(q/x)×(43)-Q(1/x)×(
t [2] (x)Q(1/x) = h(x)Q(1/qx) + e -µ h(qx)Q(q/x). (45)
Using this equation in conjunction with [START_REF] Simon | Construction of a coordinate Bethe ansatz for the asymmetric simple exclusion process with open boundaries[END_REF] then gives M µ in terms of Q:

M µ = (1 -q) d dx log Q(q/x) Q(1/x) x=-1 . ( 46
)
Let us note that using equs. [START_REF] Kitanine | Open spin chains with generic integrable boundaries: Baxter equation and Bethe ansatz completeness from SOV[END_REF], we can obtain all the equations from the so-called 'fusion hierarchy' [START_REF] Korff | A Q-operator for the quantum transfer matrix[END_REF][START_REF] Yang | Q-operator and relation from the fusion hierarchy[END_REF], which gives equations on the decomposition of products of matrices t [k] , as well as the T-Q equations for any t [k] (which involve products of k -1 matrices Q). We find for instance that t [2] (q 2 x)t [3] (x) = h(q 2 x)t [4] (x) + e -µ h(q 3 x)t [2] (x), and t [3]

(x)Q(1/x)Q(1/qx) =h(x)Q(1/qx)Q(1/q 2 x) + e -µ h(qx)Q(q/x)Q(1/q 2 x) + e -2µ h(q 2 x)Q(q/x)Q(x).
See appendix C for more details on these relations.

At this point, in order to get an explicit result for the largest eigenvalue E(µ) of M µ , we need some additional information, which we get by taking µ to 0.

E. Non-deformed case

We want to recover the expression for E(µ) which was found in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF]. Everything that we have shown so far is valid for all eigenspaces of M µ , but we now need to look for any information that we might have specifically on the dominant eigenvalue or eigenvectors of M µ , so that we can solve equation [START_REF] Prolhac | Current fluctuations in the exclusion process and Bethe ansatz[END_REF] explicitly in that eigenspace. Starting from the coordinate Bethe Ansatz, this was done, in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF], using prior knowledge on the behaviour of Bethe roots for µ → 0. We cannot use that same argument here, for two reasons: we did not start from the coordinate Bethe Ansatz, so that our Q-operator is not defined in terms of Bethe roots, and even though we can, in principle, argue that the two definitions of Q have to be consistent with one another, we do not expect to be able to do the same in the open case. Fortunately, we can get to the same result by purely algebraic calculations.

Let us go back to the definition of T per µ as given in [START_REF] Derrida | Exact large deviation function in the asymmetric exclusion process[END_REF], as well as that of X in terms of S and A:

T per µ (x, y) = Tr[A µ L i=1 X (i) ] with X = 1 + xA S(1 -xyA) S -1 (1 -A) 1 + yA .
Expanding this expression in powers of A, we see that each weight in T per µ is a finite sum of traces of the form

x i y j q k Tr[A µ S l A m ],
where k, l, and m are integers, and the term q k comes from ordering the powers of A and S. Because of the trace, those terms are 0 if l = 0 (which accounts for the conservation of the total number of particles: the number of d's and e's in the trace have to be the same). Otherwise, they are simply equal to x i y j q k /(1q m e -µ ). In particular, there is exactly one term with m = 0, which is equal to 1/(1e -µ ). In the non-deformed limit µ → 0, this term dominates the trace in each entry of the transfer matrix, as it is the only term that diverges, so that:

T per µ (x, y) ∼ 1 1 -e -µ |1 N 1 N |
where |1 N is a vector with all entries equal to 1 if the total number of particles is N , and 0 otherwise. |1 N 1 N | is the projector onto the dominant eigenspace (i.e. the steady state) of M in the sector with N particles.

From this, we draw two important conclusions. First, the prefactor in Q(y) diverges as 1/µ. Secondly, we see that this limit does not depend on x and y, which means that the roots of P (x) and Q(y), in this eigenspace, all go to infinity when µ goes to 0. This allows us to use a contour integral in [START_REF] Nepomechie | An inhomogeneous T -Q equation for the open XXX chain with general boundary terms: completeness and arbitrary spin[END_REF] to separate the contribution from P (x) and P (qx) (whose roots all go to infinity, and are therefore all out of the unit circle for µ small enough) and those from Q(1/x) and Q(q/x) (whose roots, in terms of x, all go to 0, and are inside of the unit circle for µ small enough). We also find that the eigenvalue of t [2] in that eigenspace is

1 N |t [2] |1 N = h(x) + h(qx)
which gives E(0) = 0. Using this, we can solve [START_REF] Nepomechie | An inhomogeneous T -Q equation for the open XXX chain with general boundary terms: completeness and arbitrary spin[END_REF] perturbatively in µ and find Q in terms of h. Putting the result in [START_REF] Prolhac | Current fluctuations in the exclusion process and Bethe ansatz[END_REF], we finally find E(µ). Since we will be doing the same for the open ASEP in the next section with only a few differences, we will give those final steps in detail only in section II E.

II. OPEN ASEP

We will now apply the same procedure to the open ASEP. Considering the same X matrix as before, we first generalise the boundary vectors found in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF] accordingly, and see how our construction relates to the original matrix Ansatz [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF]. We then show the PQ factorisation of the transfer matrix, which is much more straightforward than in the periodic case. Finally, we see how it decomposes into blocks, one of which is a one-parameter transfer matrix with finite auxiliary dimension, according to the same equation as in the periodic case but with a different quantum determinant. We also show, as an appendix, what this all becomes in the language of the XXZ chain with spin 1 2 . The current-counting Markov matrix for the open ASEP is given by:

M µ = m (0) + L-1 i=1 M (i) + m (L) with m (0) = -α γe -µ αe µ -γ , M (i) =      0 0 0 0 0 -q 1 0 0 q -1 0 0 0 0 0      , m (L) = -δ β δ -β
where the matrices m (0) and m (L) represent two particle reservoirs to which the system is coupled, through the first and the last site, respectively. m (0) acts on site 1 and m (L) on site L, both in the canonical occupancy basis {0, 1}, and, as previously, M (i) acts on sites i and i+1 in basis {00, 01, 10, 11}.

As in the periodic case, the bond over which we count the current can be chosen arbitrarily, so for the sake of simplicity we choose the one between the left reservoir and the first site.

A. Boundary algebra and commutation relations

Let us first find out how the presence of boundaries make this case different from the previous one. Most of the results and calculations in this section are generalisations of those found in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF]. We define:

U µ (x) = W ||A µ L i=1 X (i) (x, x)||V , T µ (y) = W ||A µ L i=1 X (i) (y, y)|| Ṽ , (47) (48)
which has the same structure as the transfer matrix from [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF], but with the X matrices replaced by their generalisation. Note that, since we have two transfer matrices, we have, in principle, four free parameters (two in each matrix), but we must in fact put the same parameter twice in each matrix (so that U µ depends only on x, and T µ on y) if we want to be able to find suitable boundary vectors. To simplify notations, we will therefore rewrite X and X as:

X(x) = n x e d n x , X = (1 -q) 2 ñx e -d -ñ x with n x = 1 + xA and ñx = 1 -xA.
Note that for a general set of fugacities, the generalisation ( 14) holds, with the matrices A µi being inserted in both U and T .

The conditions that these boundary vectors must satisfy is:

[β(d + n x ) -δ(e + n x ) -(1 -q)] ||V = 0, W || [α(e + n x ) -γ(d + n x ) -(1 -q)] = 0, [β(d -n y ) -δ(e -n y ) + (1 -q)yA] || Ṽ = 0, W || [α(e -n y ) -γ(d -n y ) + (1 -q)yA] = 0, (49) (50) (51) (52)
where we notice that the first two are the same as in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF][START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF] with n x replacing 1, and the next two also have an extra term (1q)yA. For x = y = 0, we naturally recover the conditions given in [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF]. Note that those conditions were found by trial and error, so that they may not be unique.

Commutation relations

As in the periodic case, we will now see how the transfer matrices we have just defined commute with each of the individual jump matrices.

For the bulk matrices M (i) , we can simply use eq. [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF]. Summing all of them, we get

L-1 i=1 M (i) , U µ = W ||A µ X(1) L i=2 X (i) ||V -W ||A µ L-1 i=1 X (i) X(L) ||V , (53)
L-1 i=1

M (i) , T µ = W ||A µ X(1) L i=2 X (i) || Ṽ -W ||A µ L-1 i=1 X (i) X(L) || Ṽ . ( 54
)
Unlike in the periodic case, where the trace over the product of matrices ensures that all the X matrices cancel one another, we here need to consider the action of the boundary matrices m (0) and m (L) as well. We find that we cannot cancel the two boundary terms in [START_REF] Gasper | Basic Hypergeometric Series[END_REF], and those in [START_REF] Prosen | Exterior integrability: YangBaxter form of non-equilibrium steady-state density operator[END_REF], independently. Instead, we have to consider the commutation with the product U µ (x)T µ (y). We show, in appendix A 2, that

[m (0) (µ), U µ T µ ] = -W ||A µ X(1) L i=2 X (i) ||V T µ -U µ W ||A µ X(1) L i=2 X (i) || Ṽ , (55)
[

m (L) , U µ T µ ] = W ||A µ L-1 i=1 X (i) X(L) ||V T µ + U µ W ||A µ L-1 i=1 X (i) X(L) || Ṽ . ( 56
)
which is exactly what is needed to compensate the terms coming from the bulk.

Putting those relations together, we can conclude that, for any values of x and y, we have:

[M µ , U µ (x)T µ (y)] = 0. ( 57
)

Explicit expressions for the boundary vectors

We will need an expression for the boundary vectors in a few calculations later, so we determine them explicitly here.

For that, we need to define four new parameters in terms of the boundary rates:

a = 1 2α (1 -q -α + γ) + (1 -q -α + γ) 2 + 4αγ , ã = 1 2α (1 -q -α + γ) -(1 -q -α + γ) 2 + 4αγ , b = 1 2β (1 -q -β + δ) + (1 -q -β + δ) 2 + 4βδ , b = 1 2β (1 -q -β + δ) -(1 -q -β + δ) 2 + 4βδ ,
and reversely:

α = (1 -q) (1 + a)(1 + ã) , γ = - aã(1 -q) (1 + a)(1 + ã) , δ = - b b(1 -q) (1 + b)(1 + b) , β = (1 -q) (1 + b)(1 + b)
.

In terms of these, equations ( 49)-( 52) become:

[d + b be + (1 + b b)xA -(b + b)] ||V = 0, ( 58
) W || [e + aãd + (1 + aã)xA -(a + ã)] = 0, (59)
[d + b be + (b + b)yA -(1 + b b)] || Ṽ = 0, ( 60
) W || [e + aãd + (a + ã)yA -(1 + aã)] = 0. ( 61
)
We first focus on the right boundary. We saw that a possible representation for d and e is e = S 1 (1 -

x 2 A 1 ) in U µ or e = S 2 (1 -y 2 A 2 ) in T µ , and d = S -1 i (1 -A i ) (
where indices 1 and 2 refer to the auxiliary spaces of U µ and T µ , respectively). Equations ( 58) and ( 59) become:

[S -1 1 (1 -A 1 ) + b bS 1 (1 -x 2 A 1 ) + (1 + b b)xA 1 -(b + b)] ||V = 0, [S -1 2 (1 -A 2 ) + b bS 2 (1 -y 2 A 2 ) + (b + b)yA 2 -(1 + b b)] || Ṽ = 0,
which is to say, multiplying by S i to the left:

[(1 -bS 1 )(1 -bS 1 v) -(1 -xS 1 )(1 -b bxS 1 )A 1 ] ||V = 0, [(1 -S 2 )(1 -b bS 2 ) -(1 -byS 2 )(1 -byS 2 )A 2 ] || Ṽ = 0.
We can write those vectors as generating functions, which will make it easier for us to manipulate them. Let us then write

||V = V (S 1 )||0 = ∞ k=0 V k S k 1 ||0
, which is such that A 1 ||V = V (qS 1 )||0 . Those two last equations now become:

V (S 1 ) V (qS 1 ) ||0 = (1 -xS 1 )(1 -b bxS 1 ) (1 -bS 1 )(1 -bS 1 ) ||0 , Ṽ (S 2 ) Ṽ (qS 2 ) ||0 = (1 -byS 2 )(1 -byS 2 ) (1 -S 2 )(1 -b bS 2 ) ||0 ,
which we can iterate to get:

||V = (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ ||0 , || Ṽ = (byS 2 ) ∞ ( byS 2 ) ∞ (S 2 ) ∞ (b bS 2 ) ∞ ||0 (62) (63)
where we recall that (x) ∞ is the infinite q-Pochhammer symbol:

(x) n = n-1 k=0 (1 -q k x), (x) ∞ = ∞ k=0
(1q k x).

We will sometimes write products of those in a more compact form, with only one symbol, such as

(x, y) ∞ = (x) ∞ (y) ∞ .
As for the left boundary, it is simpler to treat it using the contragredient representation X of X on both vectors (which are then the exact symmetric of ||V and || Ṽ , but with a and ã replacing b and b), and then apply the operator f i that we found in section I C in order to go back to X. This gives us:

W || = 0|| (x/S 1 ) ∞ (aãx/S 1 ) ∞ (a/S 1 ) ∞ (ã/S 1 ) ∞ (x 2 A 1 ) ∞ (qA 1 ) ∞ , W || = 0|| (ay/S 2 ) ∞ (ãy/S 2 ) ∞ (1/S 2 ) ∞ (aã/S 2 ) ∞ (y 2 A 2 ) ∞ (qA 2 ) ∞ , (64) (65)
to which we could add factors (q)∞ (x 2 )∞ and (q)∞ (y 2 )∞ for normalisation.

In fact, in most future calculations, we will use X for U µ and X for T µ , so that the factor

(y 2 A2)∞ (qA2)∞ goes to || Ṽ instead of W ||.
Note finally that scalar products involving those vectors (like, for instance, all the weights of U µ and T µ ) are well defined only for |a|, |b|, |ã| and | b| all smaller than 1. This corresponds to the ASEP being in its so-called 'maximal current' phase. We can assume this to be the case for now, but a simple analytic continuation can be performed in order to access more generic values of the parameters [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF].

We now need to show two things: that the transfer matrix is a product of two commuting matrices, one depending on x and one on y, and that for special values of the parameters it decomposes into two independent blocks.

B. R matrix and PQ factorisation

Unlike what we did for the periodic case, we will first show that U µ (x)T µ (y) factorises into two commuting matrices P (x) and Q(y). Note that this decomposition is not trivial: U µ and T µ do not commute, so P is not simply equal to U µ and Q to T µ .

First, we need to show that U µ (x)T µ (y) commutes with U µ (x )T µ (y ) for any values of x, y, x and y . If we want to be rigorous, the fact that both commute with M µ is not enough (M µ might have a degenerate eigenspace which is not degenerate for U µ (x)T µ (y)), so, as in the periodic case, we need to use R matrices to exchange spectral parameters. We cannot do as easily as in the periodic case, because of the boundaries: instead of taking R around a trace to exchange the top and bottom matrices, we now need to show that when applied to a boundary, R transfers the spectral parameters from a vector to another. We find the two following results (which are proven in appendix B 1 and illustrated on fig. -7):

R -1 (x, x; y, y) ||V (x) || Ṽ (y) = ||V (y) || Ṽ (x) , R -1 (y, y; x, x) || Ṽ (y) ||V (x) = || Ṽ (x) ||V (y) ,
where we use the same R matrix as defined in [START_REF] Murgan | Bethe ansatz derived from the functional relations of the open XXZ chain for new special cases[END_REF].

= V Ṽ x x y y f -1 f g - f -1 g + f V Ṽ x x y y FIG. 7.
Action of the R matrix on the right boundary.

These relations imply, at the level of the transfer matrices, that

U µ (x)T µ (y) = U µ (y)T µ (x), T µ (y)U µ (x) = T µ (x)U µ (y), (66) (67)
which immediately gives that U µ (x)T µ (y) commutes with U µ (x )T µ (y ), but also that

U µ (x)T µ (y) U µ (0)T µ (0) = U µ (x)T µ (0) U µ (0)T µ (y)
so that we can write, dividing by U µ (0)T µ (0),

U µ (x)T µ (y) = P (x)Q(y) (68)
with

P (x) = U µ (x) U µ (0) -1 , Q(y) = U µ (0)T µ (y),
where P and Q commute with each other. Note that, in terms of P and Q and at y = 0, eq.(66) becomes P (x)Q(0) = P (0)Q(x), so that P and Q are the same up to some constant (since P (0) is set to 1 but not Q(0)). We will not need to use this identity in our calculations, so we will keep using distinct names for those two objects, to have notations that are consistent with the periodic case.

C. Decomposition of the transfer matrix

In this section, we investigate whether anything happens for those special values of x and y that we considered in section I B. The calculations involved are slightly more complicated than they were then, not only because of the boundaries, but also because x and y are now separated, and the factors (1-xyA) in e have been replaced by (1x 2 A 1 ) and (1y 2 A 2 ), which do not give anything useful for xy = q 1-k .

The simple solution to this problem is to put x and y back together, using g - 12 (x, y) (as defined in section I C) to go from X(x, x) X(y, y) to X(x, y) X(y, x). We choose to keep the second row of matrices in the contragredient representation to simplify future calculations. We recall that it corresponds to having f (xy) commute through X(x, y), so that

X(x, y) = f (xy)X(x, y)f -1 (xy) = n x t d t e n y
This operation makes the boundary vectors more complicated: they are not a product of two vectors acting each on one row any more (see fig. -8). We will write those vectors as K + for the right boundary and K -for the left one:

K -= 0 1 , 0 2 || (x/S 1 ) ∞ (aãx/S 1 ) ∞ (a/S 1 ) ∞ (ã/S 1 ) ∞ (x 2 A 1 ) ∞ (qA 1 ) ∞ (ay/S 2 ) ∞ (ãy/S 2 ) ∞ (1/S 2 ) ∞ (aã/S 2 ) ∞ (xS 1 /S 2 ) ∞ (yS 1 /S 2 ) ∞ (q) ∞ (x 2 ) ∞ , K + = (q) ∞ (y 2 ) ∞ (yS 1 /S 2 ) ∞ (xS 1 /S 2 ) ∞ (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ (y 2 A 2 ) ∞ (qA 2 ) ∞ (byS 2 ) ∞ ( byS 2 ) ∞ (S 2 ) ∞ (b bS 2 ) ∞ ||0 1 , 0 2 . (69)
(70) 8. Two equivalent representations for Uµ(x)Tµ(y). In the first one, the two rows of X matrices are independent, and the two spectral parameters x (orange) and y (blue) are separated. In the second, x and y are back together in X, but the two rows are now intertwined.

f -1 V Ṽ f -1 g - g - X X W W f -1 V Ṽ f -1 g - g - W W X X X X X X X X X X X X X X FIG.
We will denote by K + i,j the coefficient of ||i, j in K + (i.e. the coefficient of S i 1 S j 2 in the expansion of the ratios of q-Pochhammer symbols in K + ), and by K - j,i the coefficient of i, j|| in K -. Note that, in this section, the normalisation of the boundary vectors is important. The factors (q)∞ (x 2 )∞ and (q)∞ (y 2 )∞ in K + and K -are there so that K - 0,0 = K + 0,0 = 1. Also note that the transpose of K -can be obtained from K + through the transformation {x ↔ y, S 1 → aS 2 , S 2 → aS 1 , b → 1/a, b → ã}, so that we only need to do calculations on K + . Now that x and y have been reunited in the bulk of the system, the same decomposition (into block triangular matrices) happens as did in the periodic case. We have to see whether the same happens to the boundary matrices K ± .

One-way boundaries

We first focus on the simpler case where ã = b = 0, which is to say γ = δ = 0. For this simpler case, we have:

K + = (xS 1 ) ∞ (bS 1 ) ∞ (yS 1 /S 2 ) ∞ (xS 1 /S 2 ) ∞ (y 2 A 2 ) ∞ (qA 2 ) ∞ (byS 2 ) ∞ (S 2 ) ∞ (q) ∞ (y 2 ) ∞
where the factor ||0 1 , 0 2 on which this whole function acts is implicit. This can be expanded into:

K + = (xS 1 ) ∞ (bS 1 ) ∞ ∞ k=0 (y/x) k (q) k (xS 1 /S 2 ) k ∞ n=0 (by) n (y 2 ) n (S 2 ) n
where the first sum comes from (yS1/S2)∞ (xS1/S2)∞ , and the second from everything that is to the right of that.

We will shortly be in need of one of Heine's transformation formulae for basic hypergeometric series, which we give now. For a function 2 φ 1 (a, b; c; z) defined as:

2 φ 1 (a, b; c; z) = ∞ n=0 (a) n (b) n (q) n (c) n z n (71)
we have [START_REF] Gasper | Basic Hypergeometric Series[END_REF]:

2 φ 1 (a, b; c; z) = (zab/c) ∞ (z) ∞ 2 φ 1 (c/a, c/b; c; zab/c). (72)
We will also need this simple identity on q-Pochhammer symbols:

(x) j+k = (x) j (xq j ) k . (73)
Coming back to K + , we have:

K + = (xS 1 ) ∞ (bS 1 ) ∞ ∞ j=0 (by) j (y 2 ) j (S 2 ) j ∞ k=0 (y/x) k (q) k (xS 1 /S 2 ) k (byq j ) k (y 2 q j ) k (S 2 ) k = ∞ j=0 (by) j (y 2 ) j (S 2 ) j (xS 1 ) ∞ (bS 1 ) ∞ 2 φ 1 (y/x, byq j ; y 2 q j ; xS 1 ) = ∞ j=0
(by) j (y 2 ) j (S 2 ) j 2 φ 1 (xyq j , y/b; y 2 q j , bS 1 )

where we used (73) between the first line and the second, and (72) between the second and the third. This finally gives us:

K + i,j (x, y) = (xyq j ) i (y/b) i b i (by) j (q) i (y 2 ) i+j (74)
and, at the left boundary:

K - j,i (x, y) = (xyq i ) j (x/a) i a i (ax) j (q) j (x 2 ) i+j . (75)
We need to compare those values for (x, y) equal to (1/q k-1 y, y) or (q/y, q k y) (which are the same values as before, but keeping y as a variable instead of x). We get:

K + i,j (1/q k-1 y, y) = (q j-k+1 ) i (y/b) i b i (by) j (q) i (y 2 ) i+j , K + i,j (q/y, q k y) = (q j+k+1 ) i (yq k /b) i b i (byq k ) j (q) i (y 2 q 2k ) i+j .
Since (q j-k+1 ) i = 0 for jk + 1 ≤ 0 and i + jk ≥ 0, i.e. for j ≤ k -1 and i ≥ kj, the first matrix is, as we expected, block triangular.

We now consider, for j ≥ k -1, the ratio:

K + i+k,j+k (1/q k-1 y, y) K + i,j (q/y, q k y) = (y/b) k b k (by) k (y 2 ) 2k (q j+1 ) k (q i+1 ) k .
The term (q j+1 ) k (q i+1 ) k accounts for going to the contragredient representation on both lines (consistently with what happens to the X matrices), and can be transferred to the left boundary, where they compensate similar terms, emerging from the same calculation:

K - j+k,i+k (x, 1/q k-1 x) K - j,i (q k x, q/x) = (x/a) k a k (ax) k (x 2 ) 2k (q i+1 ) k (q j+1 ) k .
All the other terms depend only on k, and factor out of the matrix product. We will now put them together. First, we need to make a few transformations on those factors. For y = 1/q k-1 x, we have:

(y/b) k b k = (-x) -k q -k(k-1)/2 (bx) k ,
(by) k = (bq -k+1 /x) k , (x/a) k a k = (-x) k q k(k-1)/2 (aq -k+1 /x) k , which gives us:

(x/a) k a k (ax) k (y/b) k b k (by) k (x 2 ) 2k (y 2 ) 2k = (aq -k+1 /x) k (bq -k+1 /x) k (ax) k (bx) k (x 2 ) 2k (x -2 q -2k+2 ) 2k = h b (q k x)h b (q/x) h b (x)h b (1/q k-1 x) with h b (x) = (x 2 ) ∞ (ax, bx) ∞ . ( 76
)
We see that the boundary matrices, just as the bulk matrices, are upper block triangular, with a first block of auxiliary dimension k, and a second of infinite auxiliary dimension, which is the same as the full matrix at different values of the spectral parameters, up to a global factor which we have just computed. Put into equations, this becomes:

P (x)Q(1/q k-1 x) = t [2] (x) + e -2kµ h b (q k x)h b (q/x) h b (x)h b (1/q k-1 x) P (q k x)Q(q/x)
where the factor e -2kµ comes, as before, from the first coefficient of the matrices A µ (of which there are now 2) in the second block.

To make things simpler, we can redefine all our matrices as:

P (x) = h b (x)P (x) , Q(x) = h b (x)Q(x) , t[k] (x) = h b (x)h b (1/q k-1 x)t [k] (x) (77)
which is equivalent to choosing another normalisation for U µ and T µ . Written in terms of those new transfer matrices, this last equation takes its definitive form:

P (x) Q(1/q k-1 x) = t[k] (x) + e -2kµ P (q k x) Q(q/x). ( 78
)
This has the exact same form as eq.( 41), the only difference being a factor 2µ instead of µ in the right hand side. Note that, since P and Q are the same function up to a constant (as we saw at the end of section II B), we find that t[k] is symmetric under x ↔ 1/q k-1 x.

The rest of the reasoning is the same as in the periodic case. We first consider the case where k = 1. This gives us quite simply: t [1] (x) = (1 + x) L (1 + 1/x) L = h(x). From this, we get:

t[1] (x) = h(x)h b (x)h b (x -1 ) = (1 + x) L (1 + 1/x) L (x 2 , x -2 ) ∞ (ax, bx, a/x, b/x) ∞ = F (x).
This explains how the function F (x) replaces h(x) (from the periodic case) in all the expressions for the cumulants of the current that were found in [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. For k = 2, through the same calculations as in the periodic case, we find the T-Q equation

t [2] (x)Q(1/x) = h(x)Q(1/qx) + e -2kµ h b (qx)h b (q/x) h b (x)h b (1/x) h(qx)Q(q/x) which translates into t[2] (x) Q(1/x) = F (x) Q(1/qx) + e -2kµ F (qx) Q(q/x). (79)
Note that although the equation for t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] is more compact that that for t [2] , the former is not polynomial due to the presence of q-Pochhammer symbols in F and in the normalisation of t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] , whereas the latter is.

We now have to verify that t[2] (x) is related to M µ through some derivative. We will do this directly in the general case (with all four boundary rates), in a few pages. For now, we just give the two-dimensional boundary matrices that we find from K + and K -:

K + 2 = 1 qx(qx-b) (q 2 x 2 -1) x(1-qxb) (q 2 x 2 -1) 0 , K - 2 = 1 (a-x) (1-x 2 ) (ax-1) q(1-x 2 )
0 .

Two-way boundaries

In the general case, where both boundaries have two non-zero rates, the calculations are much more involved, and can be found in appendix B 2. Starting from expression (70) for K + (x, y), and taking xy = q 1-k , we find that

K + i+k,j+k (x, y) K + i,j (q k x, q k y) = (x/b) k b k (xb) k (y/ b) k ( b) k (y b) k (y 2 ) 2k (-y) k q k(k-1)/2 (q j+1 ) k (q i+1 ) k (80)
and the corresponding result for the left boundary:

K - j+p,i+p (x, y) K - j,i (q p x, q p y) = (y/a) p a p (ya) p (x/ã) p (ã) p (xã) p (x 2 ) 2p (-x) p q p(p-1)/2 (q i+1 ) p (q j+1 ) p . ( 81
)
Using this, we do the exact same operations as in the previous case, replacing h b with

h b (x) = (x 2 ) ∞ (ax, ãx, bx, bx) ∞ , (82)
and we get the completely general version of the function F (x) = t[1] (x):

F (x) = (1 + x) L (1 + 1/x) L (x 2 , x -2 ) ∞ (ax, ãx, bx, bx, a/x, ã/x, b/x, b/x) ∞ . ( 83
)
We finally look at the transfer matrix t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] (x) and try to relate it to M µ , as we did in section I B for the periodic case. By analogy with eq.( 24), we will try to rewrite t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] (h)/F (x) as the standard two-dimensional transfer matrix.

The two-dimensional blocks from K + and K -for y = 1/qx are:

K + 2 = 1 qx(qx+qxb b-b-b) (q 2 x 2 -1) x(1+b b-qxb-qx b) (q 2 x 2 -1) -b bx , K - 2 = 1 (a+ã-x-aãx) (1-x 2 ) (ax+ãx-1-aã) q(1-x 2 )
aã qx and the corresponding blocks from the bulk are:

X 2 (x) =      1 + x 0 0 0 0 1 + qx 1 -1 q 0 0 1 -q 1 + 1 qx 0 0 0 0 1 + 1 x      , X 2 (x) =      1 + 1 qx 0 0 0 0 1 + 1 x 1 -q 0 0 1 -1 q 1 + x 0 0 0 0 1 + qx     
.

Consider these transformations, with λ = 1+qx q(1+x) , i.e. x = -1-qλ q(1-λ) :

L (i) (λ) = 1 1 + x 1 0 0 -q X 2 (x) 1 0 0 -1/q 1 0 0 x =      1 0 0 0 0 qλ 1 -qλ 0 0 1 -λ λ 0 0 0 0 1      and
L (i) (λ) = 1 1 + x x 0 0 1 0 1 1 0 X 2 (x) 0 1 1 0 =      1 0 0 1 -qλ 0 λ 0 0 0 0 qλ 0 1 -λ 0 0 1     
where the matrix products inside the parentheses are done in the auxiliary space, on each element of X 2 or X 2 , and the third product is done on the physical space at each site. Notice that the inner products cancel out between one site and the next, and that the outer products are done between X 2 and X 2 and amount to a global factor x (1+x) 2 on each site. Taking a product of L matrices X 2 X 2 , this transformation gives, apart from the inner products at each end of the chain, a global factor x L (1+x) 2L , which accounts for the bulk part h(x) of F (x). Also note that the matrices from L (i) need to be transposed if multiplied from right to left. Considering that t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] has a factor h b (x)h b (1/qx), and that

h(x)h b (x)h b (1/qx) = F (x) (1 -1/q 2 x 2 )(1 -1/qx 2 ) (1 -a/qx)(1 -ã/qx)(1 -b/qx)(1 -b/qx)
the transformations we need to do on the boundary matrices are:

K+ 2 (λ) = (1 -1/q 2 x 2 ) (1 -b/qx)(1 -b/qx) 1 0 0 -q K + 2 0 1 1 0 , K- 2 (λ) = (1 -1/qx 2 ) (1 -a/qx)(1 -ã/qx) 0 1 1 0 K - 2 1 0 0 -1/q ,
of which we will not write the full expression in terms of λ. Their values and first derivatives at λ = 0, which is all we will need, are, in terms of the original boundary parameters:

K+ 2 = 1 0 0 1 , d dλ K+ 2 = -2δ 2β 2δ 1 -q -2β , ( 84
) K- 2 = 1-α+γ 1+q γ q α α-γ+q 1+q , d dλ K- 2 = -α + γ + A -B γ(2q-α-γ) q α(1 + q -α -γ) -2γ -q -A + B , ( 85
) with A = 2+(α-2)α-γ 2 1+q and B = 2(1-α+γ) (1+q) 2 .
We now put all these matrices together, obtaining t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] (h)/F (x) as desired. The Lax matrices L (i) are not exactly the same as those we had for the periodic case. Their values at λ = 0 are L (i) (0) = P (i) and

L (i) (0) = P (i)
, where P (i) is the permutation matrix exchanging the auxiliary space with the physical space when applied to the right in the matrix product, and P v is the same but when applied to the left. We see that the two boundaries do not play the same role. The right boundary matrix is simply the identity at λ = 0, and serves to connect the Lax matrices on the last site. The left boundary matrix is traced, at λ = 0, because of the Lax matrices on the first site, and we see that its trace is 1. All in all, at λ = 0, the whole transfer matrix is simply the identity (see fig.

-9-a). As for its first derivative with respect to λ, we find that each pair of Lax matrices L (i) and L (i) have a total contribution of 2M (i) + (1q)(τ iτ i+1 ), where τ i is the number of particles on site i. The right boundary matrix, as can be seen in ( 84), gives a term 2m (L) + (1q)τ L . At the left boundary, we have three terms to consider, one involving the derivative of K -and the first Lax matrices taken at 0, and two more where we take the derivative of the first Lax matrices, and the boundary matrix at 0. The sum of those terms gives a total contribution of 2m (0) -(1q)τ 1 . If we now sum all the terms that we have found, the parts proportional to (1q) all cancel out, and we are simply left with 2M µ .

At the end of the day, we find that:

M µ = 1 2 1 - 1 q d dx log t[2] (x) F (x) x=-1/q (86)
which is the same as eq.( 26), with h(x) replaced by F (x), and an extra factor 1 2 . We haven't mentioned the matrices A µ in those last calculations, because their behaviour is trivial, and it is left to the reader, as an exercise, to check that adding one between two sites gives the correct deformation for M µ .

As for the periodic case, we could have done all those calculations around λ = ∞ instead of 0, which would have given an equivalent result:

M µ = 1 2 (1 -q) d dx log t[2] (x) F (qx) x=-1 . ( 87
)
In order to obtain the results found in [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF], we need one last piece of information: how P and Q behave as µ goes to 0.

D. Non-deformed case: matrix Ansatz

As in the periodic case, the final step is to consider the µ → 0 limit in our transfer matrices, in order to get some information specific to the dominant eigenspace of M µ and the corresponding eigenvalue E(µ). Once more, this is made much more difficult by the presence of the boundaries, but the principle remains the same.

We first need to consider the behaviour of T µ (y) alone. As we recall, it is defined by:

T µ (y) = W ||A µ L i=1 X (i) || Ṽ where X = 1 + yA S(1 -y 2 A) S -1 (1 -A) 1 + yA and W || = 0|| (ay/S) ∞ (ãy/S) ∞ (1/S) ∞ (aã/S) ∞ (y 2 A) ∞ (qA) ∞ , || Ṽ = (byS) ∞ ( byS) ∞ (S) ∞ (b bS) ∞ ||0 .
As in the periodic case, we can expand each entry in T µ as a finite sum of terms of the form y i q k W ||A µ S l A m || Ṽ (but this time, there is no constraint on l, as the number of particles is no longer conserved). In each entry, there is exactly one term with m = 0, equal to W ||A µ S l || Ṽ , where l is the difference between the number of particles in the initial and final configurations. We show, in appendix D, that for µ → 0, this term behaves as 1/(1e -µ ) times a prefactor which does not depend on l, and all the others remain finite, so that:

T µ (y) ∼ 1 1 -e -µ
(ay, ãy, by, by

) ∞ (y 2 ) ∞ (aã, b b, q) ∞ |1 1|.
Notice that the prefactor is equal to 1/h b (y) as defined in (82), up to a constant.

We can do the same calculations for U µ , and see that there is no divergence there, but it is not necessary. Instead, we just apply U µ (x) on this last result. As we recall, U µ (x) is defined by

U µ (x) = W ||A µ L i=1 X (i) ||V with X = n x e d n x
and where W || and ||V verify:

[β(d + n x ) -δ(e + n x ) -(1 -q)] ||V = 0, W || [α(e + n x ) -γ(d + n x ) -(1 -q)] = 0.
Since the vector |1 can be written as the tensor product of the constant vector |1 i on each site i, applying each X to |1 gives the vector

X|1 i = e + n x d + n x which we will note E D .
Since, as we recall, d, e, and A satisfy

de -q ed = (1 -q)(1 -x 2 A 2 ), Ae = q eA, dA = q Ad,
we find that D and E satisfy

DE -q ED = (1 -q)(D + E) (88)
and the conditions on the boundary vectors simply write

[βD -δE -(1 -q)] ||V = 0, W || [αE -γD -(1 -q)] = 0. (89) (90)
The matrix U 0 (x)T 0 (y) becomes, up to a prefactor:

U 0 (x)T 0 (y) ∼ |P 1|
where, for any configuration C = {τ i },

P (C) = W || L i=1 [(1 -τ i )E + τ i D]||V (91)
(where we used the fact that A µ = 1 for µ = 0), which does not depend on x. Moreover, since we know that [M, U 0 (x)T 0 (y)] = 0 (by considering eq.( 57) at µ = 0), and that 1|M = 0 (because M is a stochastic matrix), we find that

M |P = 0 (92)
which is to say that |P is the steady state of the open ASEP. We therefore recover the original 'matrix Ansatz', as presented in [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF]. Note that a matrix Ansatz also exists for the steady state of the periodic multi-species ASEP [START_REF] Prolhac | The matrix product solution of the multispecies partially asymmetric exclusion process[END_REF], which one should be able to recover from an appropriate transfer matrix with spectral parameters and current-counting deformations, and which would be related to the U q [SU(k)] algebra (for k different types of particles, counting holes as one of those types).

The normalisation of U 0 (x)T 0 (y) still needs to be determined. Since we know that the only dependence in y is a prefactor 1/h b (y), and that U µ (x)T µ (y) = U µ (y)T µ (x) (as shown in eq.( 66)), we conclude that, up to a constant term which we include in |P :

U 0 (x)T 0 (y) = 1 1 -e -µ 1 h b (x)h b (y) |P 1|
which is to say that, similarly to the periodic case, the prefactor in Q(y) (as defined in (77)) diverges as 1/µ, and that h b (x)h b (y)U 0 (x)T 0 (y) is independent of x and y, so that the roots and poles of P (x) and Q(y) all go to infinity when µ goes to 0.

Since, in that limit, Q(y) is a constant, we also find, from eq.( 79), that the eigenvalue of t [START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF] in that eigenspace is

1| t[2] |P = F (x) + F (qx).
It is straightforward to check that eq.( 86) then gives E(0) = 0, which we knew from the fact that M 0 is a stochastic matrix.

Starting from µ = 0, where, as we just saw, we know Q explicitly, we can expand everything in series in µ, and find an explicit expression for E(µ) perturbatively in µ. This is what we do in the next section, where we put everything together and give a summary of the whole procedure.

E. Summary -Functional Bethe Ansatz for the open ASEP

In this section, we collect all the results we found for the open ASEP, and show how they lead to the expressions for the cumulants of the current obtained in [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. The first step is to construct two transfer matrices U µ (x) and T µ (y):

U µ (x) = h b (x) W ||A µ L i=1 X (i) (x, x)||V , (93)
T µ (y) = h b (y) W ||A µ L i=1 X (i) (y, y)|| Ṽ (94)
(involving the function h b defined in eq.( 82)), such that, for any x and y, we have:

[M µ , U µ (x)T µ (y)] = 0. ( 95
)
Using these, we can construct two commuting matrices P (x) and Q(y) as:

P (x) = U µ (x) U µ (0) -1 , Q(y) = U µ (0)T µ (y) (96)
such that:

U µ (x)T µ (y) = P (x)Q(y). (97)
We can show that those matrices verify, for any positive integer k:

P (x)Q(1/q k-1 x) = t [k] (x) + e -2kµ P (q k x)Q(q/x) (98)
(where we omit the tildes, since we have correctly normalised P and Q from the start). The matrix t [k] (x) is the one-parameter transfer matrix with a k-dimensional auxiliary space. The first two of those relations write:

P (x)Q(1/x) = F (x) + e -2µ P (qx)Q(q/x), (99)
P (x)Q(1/qx) = t [2] (x) + e -4µ P (q 2 x)Q(q/x), (100)
where F (x) is a scalar function given by

F (x) = (1 + x) L (1 + 1/x) L (x 2 , x -2 ) ∞ (ax, ãx, bx, bx, a/x, ã/x, b/x, b/x) ∞ (101)
and t [2] is such that:

M µ = 1 2 (1 -q) d dx log t [2] (x) F (qx) x=-1 . (102)
Using eq.( 99) at x and qx, and eq.( 100), we can find the T-Q equation:

t [2] (x)Q(1/x) = F (x)Q(1/qx) + e -2µ F (qx)Q(q/x) (103)
which allows us to express M µ in terms of Q instead:

M µ = 1 2 (1 -q) d dx log Q(q/x) Q(1/x) x=-1 . (104)
We now consider:

B = -e 2µ Q(0) -1 = -e 2µ U µ (0)T µ (0) -1 . (105)
As we saw in section II D, the first eigenvalue of B goes to 0 for µ → 0, which is not a priori the case for the others. Moreover, all the roots and poles of P (x) and Q(y) go to infinity in that eigenspace.

From here on, we restrict ourselves to that specific eigenspace, so that P , Q and B refer to functions rather than matrices. The next step is to rewrite eq.( 99) in a different way, and see that all the information we have about P and Q makes it solvable. This was done in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF] for the periodic case, and we reproduce it here, under a slightly different form, for the open case. Let us therefore define a function W as:

W (x) = - 1 2 log P (x)Q(1/x) e -2µ P (qx)Q(q/x) , (106)
and a convolution kernel K, as:

K(x, x) = 2 ∞ k=1 q k 1 -q k (x/x) k + (x/x) -k (107)
along with the associated convolution operator X:

X[f ](x) = c1 dx ı2πx f (x)K(x, x). ( 108
)
Using those, as we show in appendix E, one can rewrite eq.(99) in terms of only one unknown function W :

W (x) = - 1 2 ln 1 -BF (x)e X[W ](x) (109)
which is the same as eq.( 59) in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF]. The last step is to take eq.( 104), and eq.( 106) at x = 0, to find:

E(µ) = 1 2 (1 -q) d dx log Q(q/x) Q(1/x) x=-1 , µ = -W (0) (110)
Considering what we said before about the roots and poles of P (x) being outside of the unit circle, and those of Q(1/x) being inside, we can replace 1 2 log Q(q/x) Q(1/x) by -W (x) when expressing E(µ) as a contour integral over the unit circle (since P will not contribute), and obtain:

µ = - c1 dz ı2πz W (z) (111)
and

E(µ) = -(1 -q) c1 dz ı2π(1 + z) 2 W (z), (112)
in which we recognise ( 13) and ( 14) from [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. All this is done for a < 1 and b < 1, but can then be generalised to any a and b through the same reasoning as in [START_REF] Blythe | Nonequilibrium steady states of matrix-product form: a solver's guide[END_REF] for the mean current, replacing the unit circle c 1 by small contours around Γ = {0, q k a, q k ã, q k b, q kb }.

In this last result, E(µ) is expressed as an implicit function of µ, through the variable B. In order to get an explicit expression, one has to invert (111) to get B in terms of µ, and inject the result into (112). Since B is of order µ around µ = 0, this can be done perturbatively in B and µ, to yield the coefficients of E(µ) expanded in powers of µ. Those are, up to a factorial, the cumulants of the current in our system. In the general case, they can be expressed as combinations of multiple contour integrals around powers of F convolved through K [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF]. In the simpler case of the TASEP, where q = 0, that convolution kernel vanishes, and (111) and (112) simplify to

µ = - 1 2 ∞ k=1 Γ dz ı2πz F (z) k B k k and E(µ) = - (1 -q) 2 ∞ k=1 Γ dz ı2π(1 + z) 2 F (z) k B k k ,
where F (z) is also much simpler, and has at most 9 poles. This makes calculating the cumulants of the current much easier, especially in the case where a = b = 0 [START_REF] Lazarescu | An Exact Formula for the Statistics of the Current in the TASEP with Open Boundaries[END_REF].

F. XXZ spin chain with general boundary conditions

In this section, we explain how our construction for the open ASEP can be translated for the spin-1 2 XXZ chain with non-diagonal boundary conditions [START_REF] Sandow | Partially asymmetric exclusion process with open boundaries[END_REF].

Let us first define the bulk Hamiltonian of the XXZ spin chain of length L:

H b = 1 2 L-1 k=1 h (i)
with h (i) acting as:

h (i) =      ∆ 0 0 0 0 -∆ 1 0 0 1 -∆ 0 0 0 0 ∆     
on sites i and i + 1 (in basis {00, 01, 10, 11}, as usual), and as the identity on the rest of the chain. We define ∆ as 1 2 (q -1/2q 1/2 ), which is not the usual definition for the XXZ chain (that can be obtained simply by replacing q by q 2 ).

Let us also write the deformed Markov matrix M {µi} for the special choice of weights defined by:

{µ 0 = 1 2 log γ α + ν 0 , µ i = 1 2 log (q), µ L = 1 2 log δ β + ν L }
which is on the line µ = 1 2 log γδ αβ q L-1 + ıR if ν 0 and ν L are imaginary numbers (in which case M {µi} is Hermitian). The deformed local matrices become:

m (0) (µ 0 ) = -α √ αγ e -ν0 √ αγ e ν0 -γ , M (i) (µ i ) =      0 0 0 0 0 -q √ q 0 0 √ q -1 0 0 0 0 0      , m (L) (µ l ) = -δ √ βδe ν L √ βδe -ν L -β .
It is straightforward to check that in this case, we have M {µi} = √ qH + , where is a constant, with the boundary matrices being equal to:

h (0) = 1 2 √ q (1 -q -α + γ) 2 √ αγ e -ν0 2
√ αγ e ν0 (-

1 + q + α -γ) = q -1/2 -q 1/2 2 1 (1 + a)(1 + ã) a + ã 2 √ -aãe ν0 2 √ -aãe -ν0 -a -ã , h (L) = 1 2 √ q (-1 + q + β -δ) 2 √ βδe ν L 2 √ βδe -ν L (1 -q -β + δ) = q -1/2 -q 1/2 2 1 (1 + b)(1 + b) -b -b 2 -b be -ν L 2 -b be ν L b + b .
Since we have three nontrivial parameters in each of those matrices, they are completely general: we can write (without restricting ourselves to hermitian matrices)

h (0) = a z σ z + a + σ + + a -σ -= a z a - a + -a z , h (L) = b z σ z + b + σ + + b -σ -= b z b - b + -b z ,
with a z = (1qα + γ)/2 √ q , a + = αγ/q e ν0 , a -= αγ/q e -ν0 , b z = (-1 + q + βδ)/2 √ q , b + = βδ/q e -ν L , b -= βδ/q e ν L , which is to say

ν 0 = -2 log(a + /a -), α = ( √ qa z -(1 -q)/2) 2 + a + a -- √ qa z + (1 -q)/2, γ = ( √ qa z -(1 -q)/2) 2 + a + a -+ √ qa z -(1 -q)/2,
and

ν L = -2 log(b -/b + ), β = ( √ qb z + (1 -q)/2) 2 + b + b -+ √ qb z + (1 -q)/2, δ = ( √ qa z + (1 -q)/2) 2 + b + b -- √ qb z -(1 -q)/2.
Those are well defined for any values of a z , a + , a -, b z , b + and b -.

Considering expression [START_REF] Lazarescu | An Exact Formula for the Statistics of the Current in the TASEP with Open Boundaries[END_REF], or its equivalent for an open chain, and noting that A µ i 2 = A -1/4 , we can rewrite U (x) and T (y) in a way better suited to this situation:

U (x) = φ|| L i=1 Y (i) (x)||ψ , T (y) = φ|| L i=1 Y (i) (y)|| ψ , with Y (x) = N x Σ + Σ -N x
where

N x = A -1/2 -xA 1/2 Σ + = A -1/4 eA -1/4 = q -1/4 S + (A -1/2 -x 2 A 1/2 ), Σ -= A -1/4 dA -1/4 = q 1/4 S -(A -1/2 -A 1/2 ).
The boundary vectors become:

φ|| = W ||A µ0 A 1/4 , ||ψ = A µ L A 1/4 ||V , φ|| = W ||A µ0 A 1/4 , || ψ = A µ L A 1/4 || Ṽ .
Matrices N x , Σ + and Σ -satisfy the U q [SU (2)] algebra [START_REF] Chaichian | Introduction to quantum groups[END_REF]:

[Σ + , Σ -] = (q -1/2 -q 1/2 )(A -1 -x 2 A), Σ -A = q AΣ -, AΣ + = q Σ + A,
and the conditions on the boundary vectors become:

φ|| [a + Σ + -a -Σ --2a z N x -(q -1/2 -q 1/2 )xA 1/2 ] = 0, [b -Σ --b + Σ + + 2b z N x -(q -1/2 -q 1/2 )xA 1/2 ] ||ψ = 0, φ|| [a + Σ + -a -Σ --2a z N y + (q -1/2 -q 1/2 )A -1/2 ] = 0, [b -Σ --b + Σ + + 2b z N y + (q -1/2 -q 1/2 )A -1/2 ] || ψ = 0.
We may note that the structure of this solution bears a strong resemblance to that of the Lindblad master equation found in [START_REF] Prosen | Exterior integrability: YangBaxter form of non-equilibrium steady-state density operator[END_REF][START_REF] Karevski | Exact Matrix Product Solution for the Boundary-Driven Lindblad XXZ Chain[END_REF]. In that case, the algebraic relations satisfied by the boundary vectors are different from ours, as there is only one vector per boundary but two equations per vector, which may constrain the values of the boundary parameters. The connection between integrable spin chains and certain boundary-driven Lindblad models is also investigated in [START_REF] Prosen | Eigenvalue Statistics as an Indicator of Integrability of Nonequilibrium Density Operators[END_REF][START_REF] Ilievski | Exact steady state manifold of a boundary driven spin-1 LaiSutherland chain[END_REF], and a matrix product structure of the density matrix was observed in [START_REF] Prosen | Exact Nonequilibrium Steady State of an Open Hubbard Chain[END_REF] for an open Hubbard chain, but the precise conditions for such systems to be integrable have yet to be understood.

III. CONCLUSION

In the present work, we treat the case of the asymmetric simple exclusion process with generic open boundaries and a current-counting deformation. We construct a transfer matrix with an infinitedimensional auxiliary space and two free parameters, which commutes with the deformed Markov matrix of our system. We then derive the two essential properties of that transfer matrix: it is a factor of two matrices P and Q which commute and carry one of the free variables each, and for certain values of the product of those parameters it breaks into two independent matrices, one of which is the usual oneparameter transfer matrix for a given dimension of the auxiliary space. From these results, we derive the T-Q equations for the transfer matrices, which allows us to identify our Q matrix with Baxter's Q-operator [START_REF] Baxter | Exactly Solved Models in Statistical Mechanics[END_REF], although it is obtained through an entirely different method. Writing the algebraic relations between the P and Q matrices in a given eigenstate of our system then yields the functional Bethe Ansatz equations, which we can use, through the same method as in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF], to find an explicit expression for the leading eigenvalue of our deformed Markov matrix, and thus give a rigorous proof of the results conjectured in [START_REF] Gorissen | Exact Current Statistics of the Asymmetric Simple Exclusion Process with Open Boundaries[END_REF].

We have mainly focused on one specific eigenvalue, which is of particular importance as to the physics of the ASEP, and which we were able to obtain explicitly, for any system size and boundary parameters, and perturbatively in the deformation parameter µ, using specific knowledge of the corresponding eigenvector for the non-perturbed case. It is however to be noted that only a small part of our derivation depends on that specific knowledge, and that, in principle, with similar information on other eigenstates (i.e. the positions of the roots and poles of Q for µ = 0, which can in principle be obtained through the usual Bethe equation [START_REF] De Gier | Bethe ansatz solution of the asymmetric exclusion process with open boundaries[END_REF][START_REF] De Gier | Exact spectral gaps of the asymmetric exclusion process with open boundaries[END_REF]), we would be able to obtain similar results. The density of the spectrum was obtained, in the thermodynamic limit, in the context of the periodic TASEP in [START_REF] Prolhac | Spectrum of the totally asymmetric simple exclusion process on a periodic latticebulk eigenvalues[END_REF][START_REF] Prolhac | Spectrum of the totally asymmetric simple exclusion process on a periodic lattice -first excited states[END_REF], using the coordinate Bethe Ansatz, and it would be interesting to know whether that method can be extended to the the open case or to the ASEP, and be written in terms of the functional Bethe Ansatz (for the open ASEP, only the first excited state without deformation was analysed, in [START_REF] De Gier | Slowest relaxation mode of the partially asymmetric exclusion process with open boundaries[END_REF]). Our construction can also be applied in more detail to the XXX and the XXZ spin chain, which will be the subject of a future paper.

It should be noted that, when solving the T-Q equation for µ = 0, Q is usually taken to have the so-called 'crossing symmetry' x ↔ 1/qx, which is a symmetry of the equation [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF][START_REF] Sklyanin | Boundary conditions for integrable quantum systems[END_REF]. In our case, Q(y) is explicitly constructed as a series in powers of y, and the fact that it is entire and without zeros inside of the unit circle (which is impossible with the aforementioned symmetry) plays an important part in finding the explicit form of E(µ). It would be interesting to analyse further the precise relation between our construction and the usual approach to the T-Q equation.

There remains the problem of finding an expression for the eigenvectors, possibly in a form similar to that which exists for the periodic case (as a generalised determinant, or equivalently as a product of creation operators on a vacuum state). This was done for special cases where the boundary parameters satisfy certain constraints [START_REF] Crampé | Eigenvectors of open XXZ and ASEP models for a class of nondiagonal boundary conditions[END_REF][START_REF] Crampé | Matrix coordinate Bethe Ansatz: applications to XXZ and ASEP models[END_REF][START_REF] De Gier | Bethe ansatz solution of the asymmetric exclusion process with open boundaries[END_REF][START_REF] Simon | Construction of a coordinate Bethe ansatz for the asymmetric simple exclusion process with open boundaries[END_REF], and also more recently in the general symmetric case (or XXX spin chain) [START_REF] Belliard | Heisenberg XXX Model with General Boundaries: Eigenvectors from Algebraic Bethe Ansatz[END_REF], but no complete solution for the ASEP has yet been found.

Let us also comment on the alternative method of constructing the functional Bethe equations devised recently [START_REF] Cao | Off-diagonal Bethe ansatz solutions of the anisotropic spin-1/2 chains with arbitrary boundary fields[END_REF][START_REF] Cao | Spin-1/2 XYZ model revisit: general solutions via off-diagonal Bethe ansatz[END_REF][START_REF] Kitanine | Open spin chains with generic integrable boundaries: Baxter equation and Bethe ansatz completeness from SOV[END_REF][START_REF] Nepomechie | An inhomogeneous T -Q equation for the open XXX chain with general boundary terms: completeness and arbitrary spin[END_REF][START_REF] Cao | Off-diagonal Bethe ansatz solution of the XXX spin chain with arbitrary boundary conditions[END_REF] in the context of the XXX and XXZ spin chains, which consists in adding an extra inhomogeneous term to the T-Q equation, in order to make polynomial solutions possible (it is easy to see, from the T-Q equation obtained here, that it is not in general the case without that extra term), and without modifying the relation between Q and the eigenvalues of the Hamiltonian. We do not know at present exactly how that method compares with ours, and whether it solves the same aspects of the problem. As far as we can tell, both methods have their merits. The alternative construction has polynomial solutions and allows for a finite set of Bethe equations, which are easier to obtain explicitly, but only when that is possible (e.g. for small system sizes), while our construction allows to use the PQ equation, which seems to be essential in obtaining a general but perhaps less explicit expression of the desired eigenvalue regardless of the system size, and which we would not be able to obtain in the same way with an inhomogeneous term. Moreover, although our construction requires some heavy calculations involving q-series, it is done explicitly and at the level of the whole matrices rather than for individual eigenvalues only, which spares us from having to make assumptions on Q. It would be interesting to see whether those two methods can be combined in any way.

Finally, we may note that certain steps in our derivation seem somewhat sub-optimal (in particular the proof of the block decomposition of our infinite dimensional K matrices, which can be found in appendix B 2), and we believe that simpler purely algebraic derivations might exist. The commutation relation for the periodic case or for the bulk of the open case involves a calculation almost identical to that which can be found in the appendixes of [START_REF] Lazarescu | Matrix ansatz for the fluctuations of the current in the ASEP with open boundaries[END_REF].

Let us recall:

X = n 0 e d n 1 = 1 + xA e d 1 + yA , X = n0 ê d n1 = (1 -q) 2 1 -xA e -d -1 + yA .
The parts of the transfer matrix and of the Markov matrix involving sites i and i + 1 are:

X (i) X (i+1) =      n 0 n 0 n 0 e en 0 ee n 0 d n 0 n 1 ed en 1 dn 0 de n 1 n 0 n 1 e dd dn 1 n 1 d n 1 n 1      , M i =      0 0 0 0 0 -q 1 0 0 q -1 0 0 0 0 0      .
The commutator of those two gives:

[M (i) , X (i) X (i+1) ] =      0
q(n 0 een 0 ) en 0n 0 e 0 dn 0q n 0 d deq ed

(1q)ed n 1 eq en 1 q n 0 ddn 0 (q -1)de q edde q en 1n 1 e 0 q(dn 1n

1 d) n 1 d -dn 1 0      =      0 n0 e -n 0 ê ên 0 -en 0 0 n0 d -n 0 d n0 n 1 -n 0 n1 êd -e d ên 1 -en 1 dn 0 -dn 0 de -dê n1 n 0 -n 1 n0 n1 e -n 1 ê 0 dn 1 -dn 1 n1 d -n 1 d 0      = X(i) X (i+1) -X (i) X(i+1) . ( A1
)
where we got from the second to the third line using:

deqed = (1q)(1 -xyA 2 ) = n0 n 1n 0 n1 , q(n 0 een 0 ) = (q -1)xAe = n0 en 0 ê, en 0n 0 e = (1q)xeA = ên 0en 0 , n 1 eq en 1 = (1q)e = ên 1en 1 , q en 1n 1 e = (q -1)e = n1 en 1 ê,

q(dn 1 -n 1 d) = (q -1)ydA = dn 1 -dn 1 , n 1 d -dn 1 = (1 -q)yAd = n1 d -n 1 d, dn 0 -q n 0 d = (1 -q)d = n0 d -n 0 d, q n 0 d -dn 0 = (q -1)d = dn 0 -dn 0 .
(A2)

Boundaries

We recall, for the left boundary vectors:

W || [α(e + n x ) -γ(d + n x ) -(1 -q)] = 0, W || [α(e -n y ) -γ(d -n y ) + (1 -q)yA] = 0.
We consider the commutator of the left boundary matrix from M µ with the relevant part of U µ and T µ (omitting the matrix A µ here, because, as in the periodic case, its action is trivial), and use those relations. We get:

m (0) , W ||X (1) = W || γd -αe -(α -γ)e (α -γ)d αe -γd = (α -γ) W || n x -e d -n x + (1 -q) W || -1 0 0 1 m (0) , W ||X (1) = W || γd -αe -(α -γ)e (α -γ)d αe -γd = (α -γ) W || -n y -e d n y + (1 -q) W || yA 0 0 -yA .
The first terms in each of those two equations cancel one another:

n x -e d -n x X(y) + X(x) -n y -e d n y = 0
We note that if we had had two different spectral parameters for the two diagonal terms in X(x) or X(y), this relation wouldn't have been possible (or we would have needed two equations on each boundary vector instead of one).

As for the second terms, a straightforward calculation gives:

(

1 -q) -1 0 0 1 X(y) + X(x) yA 0 0 -yA = -X(x) X(y) -X(x) X(y). so that [m (0) (µ), U µ T µ ] = -W ||A µ X(1) L i=2 X (i) ||V T µ -U µ W ||A µ X(1) L i=2 X (i) || Ṽ .
The exact same calculations can be done at the right boundary, and yields: The integral representation of the q-Beta function:

[m (L) , U µ T µ ] = W ||A µ L-1 i=1 X (i) X(L) ||V T µ + U µ W ||A µ L-1 i=1 X (i) X(L) || Ṽ .
(a) n (b) m (ab) n+m = (a) ∞ (b) ∞ (ab) ∞ (q) ∞ 1 0 d q t t (qt) ∞ (aq n t) ∞ t log q (b)+m (B9)
where q log q (x) = x and

1 0 dqt t f (t) = ∞ k=0 f (q k ).
The q-derivative:

D x f (x) = f (x)-f (qx)
x , and the q-Leibniz formula:

D n x f (x)g(x) = n k=0 (q) n (q) k (q) n-k D k x f (x) x→q n-k x D n-k x g(x) (B10)
Finally, a useful relation:

k (a) k (b) n-k (q) k (q) n-k x k y n-k = k (ax/y) k (by/x) n-k (q) k (q) n-k y k x n-k (B11)
which can be found by expanding (axt)∞(byt)∞ (xt)∞(yt)∞ in two different ways.

Action of the R matrix

We start from

I 0 = R -1 (x, x; y, y) ||V (x) || Ṽ (y) = (xy)∞ (q)∞ (qA2)∞ (xyA2)∞ K + ,
as defined in ...

I 0 = (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ (qA 2 ) ∞ (xyA 2 ) ∞ (yS 1 /S 2 ) ∞ (xS 1 /S 2 ) ∞ (y 2 A 2 ) ∞ (qA 2 ) ∞ (byS 2 ) ∞ ( byS 2 ) ∞ (S 2 ) ∞ (b bS 2 ) ∞ (xy) ∞ (y 2 ) ∞ (B12)
where the term ||0 1 , 0 2 is implicit. We expand all but the leftmost two ratios:

I 0 = (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ n,m,k (xy) k (y/x) n+m-k (q) k (q) n+m-k (xS 1 ) n+m-k S k 2 (q) n+m (y 2 ) n+m (yb) n (y/b) m (q) n (q) m (b b) m .
Using equation (B1), and relabelling nr and ms as n and m, we get:

I 0 = (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ n,m,r,s (xy) r+s (y/x) n+m (y 2 ) n+m+r+s (yb) n+r (y/b) m+s (q) n (q) m (q) r (q) s (q s xS 1 ) n (b bxS 1 ) m S r 2 (b bS 2 ) s = (xS 1 ) ∞ (b bxS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ r,s (yb) r (y/b) s (q) r (q) s (xy) r+s (y 2 ) r+s S r 2
(b bS 2 ) s Φ 1 y/x, q s y/b q r yb ; q r+s y 2 ; b bxS1 q s xS1 )

We then use equation (B7):

I 0 = (xS 1 ) ∞ ( bS 1 ) ∞ r,s (yb) r (y/b) s (q) r (q) s (xy) r+s (y 2 ) r+s S r 2 (b bS 2 ) s (q r+s xy) n (q s y/b) n (y b) n (q r+s y 2 ) n (q) n (bS 1 ) n (q n+s xy bS 1 ) ∞ (q s xS 1 ) ∞ = (xS 1 ) ∞ ( bS 1 ) ∞ n,s (y/b) s+n (q) s (xy) s+n (y 2 ) s+n (b bS 2 ) s (y b) n (q) n (bS 1 ) n (q n+s xy bS 1 ) ∞ (q s xS 1 ) ∞ 2 φ 1 ( yb,xyq s+n y 2 q s+n ; S 2 )
Then we use (B5) on 2 φ 1 :

I 0 = (xS 1 ) ∞ (xbS 2 ) ∞ ( bS 1 ) ∞ (S 2 ) ∞ n,r,s (y/b) r+s+n (y 2 ) r+s+n (y/x) r (q) r (xy) s+n (y b) n (q) n (q) s (q n+s xy bS 1 ) ∞ (q s xS 1 ) ∞ (xbS 2 ) r (b bS 2 ) s (bS 1 ) n = (xS 1 ) ∞ (xbS 2 ) ∞ ( bS 1 ) ∞ (S 2 ) ∞ r,s
(y/b) r+s (y 2 ) r+s (y/x) r (q) r (xy) s (q) s (xbS 2 ) r (b bS 2 ) s (q s xy bS 1 ) ∞ (q s xS 1 ) ∞ 3 φ 2 ( y/bq r+s ,xyq s ,y b y 2 q r+s ,q s xy bS1 ; bS 1 ) and (B8) on 3 φ 2 :

I 0 = (xS 1 ) ∞ (xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ n,m,r,s (y/b) r+s (yb) n+m (y 2 ) n+m+r+s (xy) n+s (y/x) m+r (q) n (q) m (q) r (q) s ( bS 1 ) n (q s xS 1 ) m (xbS 2 ) r (b bS 2 ) s
We then use (B1) in reverse (and exchange k and n + mk for convenience) to find:

I 0 = (xS 1 ) ∞ (xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ n,m,k (y/x) k (xy) n+m-k (q) k (q) n+m-k bn+m-k x k (q) n+m (y 2 ) n+m (y/b) n (yb) m (q) n (q) m (xbS 2 ) n (xS 1 ) m (B13)
In this expression, we can see that, after rescaling S 1 to S 1 /b, we have the symmetry {S 1 ↔ S 2 , b ↔ 1/b}, so that

I 0 = (xbS 2 ) ∞ (x bS 2 ) ∞ (S 2 ) ∞ (b bS 2 ) ∞ (qA 1 ) ∞ (xyA 1 ) ∞ (yS 2 /S 1 ) ∞ (xS 2 /S 1 ) ∞ (y 2 A 1 ) ∞ (qA 1 ) ∞ (yS 1 ) ∞ (b byS 1 ) ∞ (bS 1 ) ∞ ( bS 1 ) ∞ (xy) ∞ (y 2 ) ∞
which is to say, after rearranging a few terms:

R -1 (x, x; y, y) ||V (x) || Ṽ (y) = ||V (y) || Ṽ (x) .
Through the exact same calculations, we can also show that R -1 (y, y; x, x) || Ṽ (y) ||V (x) = || Ṽ (x) ||V (y) .

Truncation at xy=q 1-p

What we want to prove here, taking xy = q 1-p , is twofold. First, we need to show that K + i,j (x, y) = 0 for j ≤ p -1 and i ≥ p. Then, we need to find the relation between K + i+p,j+p (x, y) and K + i,j (q p x, q p y), which we will do by calculating

I 0 = D p S1 D p S2 I 0 .
The first part is rather straightforward. Considering equation (B14), and since

K + = (q)∞ (xy)∞ (xyA2)∞ (
qA2)∞ I 0 , we see that K + has a factor (xyA1)∞ (xyA2)∞ which results in a factor (xy)i (xy)j in K + i,j . Taking xy = q 1-p , this factor is 0 precisely when j ≤ p -1 and i ≥ p, and since the rest of the expression for K + has no pole for xy = q 1-p , we conclude that the corresponding terms in K + are equal to 0. Note that, even though that prefactor is infinite for i ≤ p -1 and j ≥ p, the absence of poles in K + tells us that it is compensated by the rest of the expression.

We now need to calculate I 0 = D p S1 D p S2 I 0 . In order to do this, we first notice that

(q) n+m (q) n+m-k bn+m-k = D k b bn+m
and we use equation (B9) on (y/x) k (xy) n+m-k (y 2 )n+m in (B13) to write:

I 1 = (xS 1 ) ∞ (xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ n,m,k 1 0 d q t t (qt) ∞ (y/xq k t) ∞ t log q (xy)+n+m-k x k (q) k D k b bn+m (y/b) n (yb) m (q) n (q) m (bS 2 ) n (S 1 ) m = k 1 0 d q t t (qt) ∞ (y/xq k t) ∞ t log q (xy) (x/t) k (q) k D k b (xS 1 ) ∞ (yb btS 1 ) ∞ (bS 1 ) ∞ ( btS 1 ) ∞ (xbS 2 ) ∞ (y btS 2 ) ∞ (S 2 ) ∞ (b btS 2 ) ∞
where I 0 = I 1 (xy)∞(y/x)∞ (y 2 )∞(q)∞ . That prefactor will be put back in later when we re-sum the integral. Note that the term (xy) ∞ gives 0 when we take xy = q 1-p , but since that term is part of I 0 and not of K + , it merely serves to simplify the calculation, and will be taken away at the end. For that reason, we will not apply xy = q 1-p on that term, but keep it in its generic form.

We will now differentiate with respect to S 1 and S 2 . Consider the following relation, which is easy to prove by recursion: if ab/cd = q 1-p , then

D p X (a) ∞ (b) ∞ (c) ∞ (d) ∞ = (a/c) p (b/c) p c p (q p a) ∞ (q p b) ∞ (c) ∞ (d) ∞ .
This relation applies to both groups of q-Pochhammer symbols at the right in the integral (remember that xy = q 1-p ), and gives:

D p S1 (xS 1 ) ∞ (yb btS 1 ) ∞ (bS 1 ) ∞ ( btS 1 ) ∞ = (x/b) p (y bt) p b p (q p xS 1 ) ∞ (q p yb btS 1 ) ∞ (bS 1 ) ∞ ( btS 1 ) ∞ D p S2 (xbS 2 ) ∞ (y btS 2 ) ∞ (S 2 ) ∞ (b btS 2 ) ∞ = (xb) p (y bt) p (q p xbS 2 ) ∞ (q p y btS 2 ) ∞ (S 2 ) ∞ (b btS 2 ) ∞
We then have:

I 1 = (x/b) p b p (xb) p (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ k 1 0 d q t t (qt) ∞ (y/xq k t) ∞
t log q (xy) (x/t) k (q) k D k b (q p y bt) p (q p y bt) p (yb btS 1 ) ∞ (y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞

We use the formula for the q-derivative of a product (B10) on the derivative with respect to b, with f ( b) = (y bt) p and the rest in g. Considering that D l b(y bt) p = (q p yt) l (q -p ) l (y btq l ) p-l , we find:

D k b (y bt) p (y bt) p (yb btS 1 ) ∞ (y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞ = k l=0 (q) k (q) l (q) k-l (q p yt) l (q -p ) l (y btq k ) p-l D k-l b (y bt) p (q p yb btS 1 ) ∞ (q p y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞ so that, renaming kl as k, we have:

I 1 = (x/b) p b p (xb) p (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ k,l 1 0 d q t t (qt) ∞ (y/xq k+l t) ∞
t log q (xy) (q p xy) l (q -p ) l (q) l (y btq k+l ) p-l (x/t) k (q) k D k b (y bt) p (q p yb btS 1 ) ∞ (q p y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞

Consider the part that depends on l (writing 1 (y/xq k+l t)∞ = (y/xq k t) l (y/xq k t)∞ ): l (y/xq k t) l (q p xy) l (q -p ) l (q) l (y btq k+l ) p-l = l (y/xq k t) l ((y btq k+p-1 ) -1 ) p-l (q) p (q) l (q) p-l (xy) l (y btq k ) p-l (-1) p q p(p-1)/2 which, through (B11), is equal to: l (y/ b) l (q 1-p /xy) p-l (q) l (q) p-l (y btq k ) l (xy) p-l (-1) p q p(p-1)/2 (q) p = l (y/ b) l (xyq l ) p-l (q) l (q) p-l (y btq k ) l (-1) l q l(l-1)/2 (q) p Using that xy = q 1-p , only the term for l = p survives, namely (y/ b) p (y btq k ) p (-1) p q p(p-1)/2 , so that: I 1 = (x/b) p b p (xb) p (y/ b) p (y b) p (-1) p q p(p-1)/2 (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ k,l 1 0 d q t t (qt) ∞ (y/xq k t) ∞ t log q (xy)+p (q p x/t) k (q) k D k b (y bt) p (q p yb btS 1 ) ∞ (q p y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞ Using again the formula for the derivative of a product and relabelling kl as k again, we find:

I 1 = (x/b
) p b p (xb) p (y/ b) p (y b) p (-1) p q p(p-1)/2 (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ k,l 1 0 d q t t (qt) ∞ (y/xq k+l t) ∞ t log q (xy)+p (q 2p xy) l (q -p ) l (q) l (y btq k+l ) p-l (q p x/t) k (q) k D k b (q p yb btS 1 ) ∞ (q p y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞

We have almost the same sum as before over l, only with an extra term q pl which changes the result to: l (y/xq k t) l (q 2p xy) l (q -p ) l (q) l (y btq k+l ) p-l = l (q p y/ b) l (xyq p+l ) p-l (q) l (q) p-l (y btq k ) l (-1) l q l(l-1)/2 (q) p so that, this time, the whole sum survives to taking xy = q 1-p . We now have to re-sum the integral over t. First, we expand the two ratios of q-Pochhammer symbols on the right, and apply the derivative with respect to b:

D k b (q p yb btS 1 ) ∞ (q p y btS 2 ) ∞ ( btS 1 ) ∞ (b btS 2 ) ∞ = t n+m n,m
(q) n+m (q) n+m-k bn+m-k (q p y/b) n (q p yb) m (q) n (q) m (bS 2 ) n (S 1 ) m so that I 1 = (x/b) p b p (xb) p (y/ b) p (y b) p (-1) p q p(p-1)/2 (q) p (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ n,m,k,l 1 0 d q t t (qt) ∞ (y/xq k t) ∞ t log q (xy)+p+n+m+l-k (q p y/ b) l (xyq p+l ) p-l (q) l (q) p-l (y bq k ) l (-1) l q l(l-1)/2 (q p x) k (q) k (q) n+m (q) n+m-k bn+m-k (q p y/b) n (q p yb) m (q) n (q) m (bS 2 ) n (S 1 ) m

We only have to consider the part that depends on k and l, as we need to keep the terms that depend on n and m intact. We re-sum the integral, going back from I 1 to I 0 : k,l (y/x) k (xy) n+m-k+l+p (y 2 ) n+m+l+p (q p y/ b) l (xyq p+l ) p-l (q) l (q) p-l (y bq k ) l (-1) l q l(l-1)/2 (q p x) k (q) k (q) n+m (q) n+m-k bn+m-k Noticing that (xy) n+m-k+l+p (xyq p+l ) p-l = (xyq p+l ) n+m-k (xy) 2p , we find: bn+m (q) n+m (xy) 2p k,l (y/x) k (xyq p+l ) n+m-k (q) k (q) n+m-k (q p+l x/ b) k (q p y/ b) l (q) l (q) p-l (y b) l (y 2 ) n+m+l+p (-1) l q l(l-1)/2

Using (B11) on the sum over k, we get: bn+m (q) n+m (xy) 2p k,l (y b) k (q p+l y/ b) n+m-k (q) k (q) n+m-k (q p+l x/ b) k (q p y/ b) l (q) l (q) p-l (y b) l (y 2 ) n+m+l+p (-1) l q l(l-1)/2 (B14)

We now write:

(q p+l y/ b) n+m-k (q p y/ b) l (y 2 ) n+m+l+p = (q p y/ b) n+m-k+l (y 2 ) n+m+l+p = (q p y/ b) ∞ (y 2 ) ∞ (y 2 q n+m+l+p ) ∞ (q n+m-k+l+p y/ b) ∞ and (y 2 q n+m+l+p ) ∞ (q n+m-k+l+p y/ b) ∞ = r (y bq k ) r (q) r (q p+n+m-k+l y/ b) r

Putting this into (B14), we get: bn+m (q) n+m (xy) 2p (q) p (q p y/ b) ∞ (y 2 ) ∞ k,l,r (y b) k+r (q) k (q) n+m-k (q) r (q p x/ b) k (q p+n+m-k y/ b) r (q) p (q) l (q) p-l q l(l-1)/2 (-q k+r y b) l Since (a) p = l (q)p (q) l (q) p-l q l(l-1)/2 (-a) l , this simplifies to bn+m (q) n+m (xy) 2p (q) p (q p y/ b) ∞ (y 2 ) ∞ k,r (y b) k+r (q) k (q) n+m-k (q) r (q p x/ b) k (q p+n+m-k y/ b) r (q k+r y b) p = bn+m (q) n+m (xy) 2p (q) p (q p y/ b) ∞ (y 2 ) ∞ k,r (q p y b) k (q p+k y b) r (q) k (q) n+m-k (q) r (q p x/ b) k (q p+n+m-k y/ b) r (y b) p

We can now get rid of the sum over r. We have:

(q p y/ b) ∞ (y 2 ) ∞ r (q p+k y b) r (q) r (q p+n+m-k y/ b) r = (q p y/ b) n+m-k (y 2 ) n+m+2p which gives (y b) p bn+m (xy) 2p (q) p (q) n+m (y 2 ) n+m+2p k (q p y b) k (q p y/ b) n+m-k (q) k (q) n+m-k (q p x/ b) k

Using (B11) one final time on the sum over k gives:

(y b) p bn+m (xy) 2p (q) p (q) n+m (y 2 ) n+m+2p k (y/x) k (q 2p xy) n+m-k (q) k (q) n+m-k (q p x/ b) k

We can at last put this back into K + . We get: (q p xS 1 ) ∞ (q p xbS 2 ) ∞ (bS 1 ) ∞ (S 2 ) ∞ n,m,k (y/x) k (q 2p xy) n+m-k (q) k (q) n+m-k (q p x) k bn+m-k (q) n+m (y 2 q 2p ) n+m (q p y/b) n (q p yb) m (q) n (q) m (bS 2 ) n (S 1 ) m which is to say (-y) p q p(p-1)/2 (xy) 2p I 0 (q p x, q p y).

I 0 = (x
In terms of K + , this gives us, for xy = q 1-p : K + i+p,j+p (x, y) K + i,j (q p x, q p y) = (x/b) p b p (xb) p (y/ b) p ( b) p (y b) p (y 2 ) 2p (-y) p q p(p-1)/2 (q j+1 ) p (q i+1 ) p .

(B15)

The equivalent result for K -is:

K - j+p,i+p (x, y) K - j,i (q p x, q p y) = (y/a) p a p (ya) p (x/ã) p (ã) p (xã) p (x 2 ) 2p (-x) p q p(p-1)/2 (q i+1 ) p (q j+1 ) p . (B16) so that

t [2] 0 t [k] 1 -t [k+1] 0 h 1 = e -(k+1)µ P k+1 Q 1 P 1 Q 1 +e -µ P 0 Q k P 2 Q -1 -e -2µ P 2 Q k P 1 Q -1 -e -kµ P 0 Q 0 P k+1 Q 1 = e -µ (P 0 Q 0 -e -µ P 1 Q -1 )(P 2 Q k -e -(k-1)µ P k+1 Q 1 ) = e -µ t [k-1] 2 h 0 .
We obtain the fusion equations:

t [2] 0 t [k] 1 = h 1 t [k+1] 0 + e -µ h 0 t [k-1] 2 (C5)
which, for k → ∞, allows to recover the T-Q equation [START_REF] Chaichian | Introduction to quantum groups[END_REF].

Through the same reasoning, we also find

t [k] 0 t [2] k-1 = h k-1 t [k+1] 0 + e -µ h k t [k-1] 0 , ( C6
)
and combining the two, we find

t [2] 0 t [k] 1 t [2] k = h 1 h k t [k+2] 0 + e -µ h 0 h k t [k]
2 + e -µ h 1 h k+1 t

[k]

0 + e -2µ h 0 h k+1 t

[k-2] 2

.

(C7)

FIG. 3 .

 3 FIG. 3. Schematic representation of f12 and f -112 . The exponents + and -of the L matrices on one side and the other are represented by, respectively, right and left arrows (the top row being L1, and the second L2), and their arguments are represented by colours (blue for x and orange for y).

FIG. 4 .

 4 FIG. 4. Schematic representation of g + 23 and g - 23 .The top row represents L2, and the second L3, and their arguments are represented by colours (red for x and y and green for x and y ).

FIG. 5 .

 5 FIG. 5. Schematic representation of the R matrix for the periodic ASEP.

FIG. 9 .

 9 FIG. 9. Schematic representation of the value and first logarithmic derivative of t[2] at 0. The first is the identity matrix (a). The second is a sum of terms adding up to Mµ: three terms at the left boundary (b), two for each bond in the bulk (c) and one term at the right boundary (d).

3 φ 2

 2 ( a,b,c d,e ; de/abc) = (e/b) ∞ (e/c) ∞ (e) ∞ (de/abc) ∞ Φ 1 d/a, b d/b ; d; e/b e/c ). (B8)

  /b) p b p (xb) p (y/ b) p ( b) p (y b) p (-y) p q p(p-1)/2 (xy) 2p (y 2 ) 2p

I 0

 0 (x, y) = (x/b) p b p (xb) p (y/ b) p ( b) p (y b) p (y 2 ) 2p

The authors would like to thank K. Mallick, R. Blythe, E. Ragoucy, P. Di Francesco, V. Terras, G. Misguich, T. Prosen and P. Baseilhac for their helpful comments and discussions. A. L. gratefully acknowledges financial support from the Interuniversity Attraction Pole -Phase VII/18 (Dynamics, Geometry and Statistical Physics) at the K. U. Leuven.

Appendix B: Boundary vectors: R matrix and truncation

We start this appendix with a few definitions and formulae that we will need during our calculation. Most of those can be found in [START_REF] Gasper | Basic Hypergeometric Series[END_REF].

First, the q-binomial recursion formula at order k:

(q) n+m (q) n (q) m = r+s=k (q) n+m-k (q) n-r (q) m-s (q) k (q) r (q) s q s(n-r) .

(B1)

Then, some useful basic hypergeometric functions:

and the q-Appell function:

along with special relations that they verify: the q-Euler formulae

and relations that link Φ 1 and 3 φ 2 :

; az 2 /b 1 ), (B7)

Appendix C: T-Q equations and fusion hierarchy

In this appendix, we show how the set of PQ equations

allows to recover all the T-Q equation for each t [k] , as well as fusion equations, where products of t [k] 's are expressed as linear combinations of t [l] 's with different dimensions l. We will consider the periodic case, where t [1] (x) = h(x). For the open case, one simply has to replace h(x) by F (x) and µ by 2µ.

To make our notations more compact, we will make the dependence in x implicit, and write

and h(q k x) = h k . In a given equation, one can shift all the indices by any integer k, which corresponds to taking the equation at q k x instead of x.

T-Q equations

This first calculation is rather straightforward. Consider the PQ of order 1

and multiply the PQ equation of order k

by the product of all Q l 's for l between 0 and k -2:

The right-hand side of the equation can then be rewritten to give

Q l in which we can simply use eq.(C1) shifted by n to obtain the T-Q equation:

As we see, this equation involves products of k-1 matrices Q. By transferring all the Qs to the right-hand side, we find an alternative expression:

Fusion equations

Let us consider t

1 . Using the PQ equations of order 2 and k, we find:

The first two terms of this sum can also be obtained by considering t

[k+1] 0

Appendix D: Non-deformed limit: asymptotic calculations

We want to estimate terms of the form

for µ → 0, assuming that |a|, |ã|, |b| and | b| are all smaller than 1. In particular, we want to show that it diverges if i = 0, and is finite otherwise. Let us assume that j ≥ 0 (the calculations for j < 0 are similar). We have:

Since we are looking for possible divergences when e -µ goes to 1, we should focus on the asymptotic behaviour of the terms of the series, as n + m becomes large. As we will need to approximate q-Pochhammer symbols, let us first remark that, for 0

which can be proven by recursion by considering that (xq

From this, we deduce that

Consider now, for N large, the term n+m=N (ya)n(y/a)m (q)n(q)m (aã) m . It can be approximated, for 0 ≤ n ≤ N/2, by

and, for N/2 < n ≤ N , by

The same can be done with

, which leaves us with

As we can see, if i = 0, I ij diverges as 1/µ, with the prefactor that is given here, and is otherwise finite. We can check that, in all cases, the error in this estimation is finite.

Appendix E: Self-consistent integral representation of the PQ equation

We consider eq.(99):

and define a function W (x) as in (106):

x) e -2µ P (qx)Q(q/x) .

Combining them, along with B = -e 2µ Q(0) -1 , we see that we also have

Since for µ small enough, all the roots and poles of P (x) and Q(y) are outside of the unit circle, we can expand their logarithms on and inside of the unit circle in powers of x and y, respectively. We write:

where we recall that P (0) was set to 1.

In those terms, the two expressions for W (x) become

Notice that in the first equation, each of the coefficients a k and b k is connected to a different power of x, so that W (x) determines those completely. We can therefore, in principle, write the argument of the exponential in the second equation in terms of W (x) alone. In this case, we can even do it explicitly: consider that the coefficients of f (x)f (qx) are of the form a k (1q k ), and those of f (qx) are a k q k . To express the latter in terms of the former, each x k has to be replaced by q k 1-q k x k , which can be done through a simple convolution:

The same can be done for g with negative powers of x, and, putting everything together, we find:

f (qx)g(q/x) = X[W ](x) where

Notice that the constant term 2µ in W disappears in the convolution.

This allows us to finally write:

This same method was used in [START_REF] Prolhac | Tree structures for the current fluctuations in the exclusion process[END_REF] for the periodic ASEP.