Integrable Floquet dynamics, generalized exclusion processes and"fused"matrix ansatz

We present a general method for constructing integrable stochastic processes, with two-step discrete time Floquet dynamics, from the transfer matrix formalism. The models can be interpreted as a discrete time parallel update. The method can be applied for both periodic and open boundary conditions. We also show how the stationary distribution can be built as a matrix product state. As an illustration we construct a parallel discrete time dynamics associated with the R-matrix of the SSEP and of the ASEP, and provide the associated stationary distributions in a matrix product form. We use this general framework to introduce new integrable generalized exclusion processes, where a fixed number of particles is allowed on each lattice site in opposition to the (single particle) exclusion process models. They are constructed using the fusion procedure of R-matrices (and K-matrices for open boundary conditions) for the SSEP and ASEP. We develop a new method, that we named"fused"matrix ansatz, to build explicitly the stationary distribution in a matrix product form. We use this algebraic structure to compute physical observables such as the correlation functions and the mean particle current.

product form. We illustrate the method with two examples: the Symmetric Simple Exclusion Process (SSEP) and the Asymmetric Simple Exclusion Process (ASEP). In section 3 we use the method to introduce new integrable models with generalized exclusion rules where several particles are allowed on the same site. The building blocks of the models, the R and K matrices, are constructed through the fusion procedure. We present a precise description of the transition probabilities of the Markov matrix. The associated stationary state is exactly computed in a matrix product form, by introducing a new technique that we call fused matrix ansatz. We finally end up in section 4 with some concluding remarks and interesting open questions.
2 Integrable Floquet dynamics from transfer matrix formalism.

Integrability and transfer matrix formalism.
We give in this subsection a very brief review of integrability in the context of exclusion processes. We introduce the main objects that are needed in our construction of discrete time models. The reader is invited to refer to [13,48] for details.
Preliminaries. We are interested in describing a system of interacting particles on a one dimensional lattice with L sites. Each site of the lattice can carry at most s particles. Note that when s = 1, there is at most one particle per site and this class of model is known as exclusion processes. Paradigmatic examples of exclusion processes are the ASEP and the SSEP [14]. For s > 1 the exclusion constraint is partially relaxed and the corresponding class of models is called generalized exclusion process. For each site i of the lattice, we define a local state variable τ i ∈ {0, 1, . . . , s} denoting the number of particles lying on site i. A configuration of particles on the lattice is thus efficiently encapsulated in a L-uplet τ = (τ 1 , τ 2 , . . . , τ L ). We will specify later on a stochastic dynamics on this configuration space, which motivates the introduction of the quantity P t (τ ), which stands for the probability to observe the system in the configuration τ at time t.
The time evolution of P t (τ ) obeys a master equation. It will be particularly efficient to formulate this master equation in a matrix form, introducing a Markov matrix and a probability vector. The first step toward this goal is to associate to each value τ = 0, . . . , s of the local state variable with a local basis vector of V = C s+1 |τ = (0, . . . , 0 τ , 1, 0, . . . , 0 (2.1) It allows us to define a probability vector where we have constructed a basis vector of V ⊗L associated to each configuration τ as tensor product of elementary vectors |τ = |τ 1 ⊗ |τ 2 ⊗ · · · ⊗ |τ L . 3) The time evolution of the probability vector can be encoded in a master equation that takes slightly different forms depending if the stochastic model describes continuous-time or discretetime dynamics. For continuous-time dynamics, the master equation reads where M is a continuous-time Markov matrix (acting on the space V ⊗L ) with non-negative offdiagonal entries and whose sum of the entries of each column vanishes.
In the case of a discrete-time dynamics, the master equation reads where M is a discrete-time Markov matrix with non-negative entries whose sum over each column is one. The present paper will be essentially devoted to the construction of integrable discrete-time Markov matrices M from the transfer matrix approach. We will also be interested in computing exactly in a matrix product form the stationary state 2 3 |S , satisfying Periodic boundary conditions. The building block of integrable models with periodic boundary conditions is a matrix R(z) depending on a spectral parameter z and acting on two components V ⊗ V of the tensor space V (i.e on two sites of the lattice).
The key equation for integrability in the periodic boundary case is the Yang-Baxter equation (2.7) The previous equation holds in V ⊗ V ⊗ V. The subscripts indicate on which components of the tensor space the matrix R is acting non-trivially. For instance R 12 (z) = R(z) ⊗ 1, R 23 (z) = 1 ⊗ R(z),... We also require that the matrix R(z) satisfies further properties. We are in this paper interested in constructing integrable discrete time Markov matrices, i.e matrices with nonnegative entries whose sum on each column is equal to one. We will therefore consider R-matrices satisfying the Markovian property σ| ⊗ σ|R(z) = σ| ⊗ σ| (2.8) where σ| = s v=0 v| is the row vector with all entries equal to one. The latter property is nothing else than imposing that the sum of the entries of each column is equal to one. In our construction the R-matrix will indeed play the role of a local Markovian operator, see for instance (2.29). We will also require that R-matrices satisfy the regularity property R(1) = P (2.9) where P is the permutation operator in V ⊗ V, i.e we have P |τ ⊗ |τ ′ = |τ ′ ⊗ |τ for all τ, τ ′ = 0, . . . , s. In the continuous time setting, it is known that this property ensures that one obtains a continuous time Markov matrix (or a quantum Hamiltonian in the context of spin chains) acting locally on the lattice from the transfer matrix constructed below. In the discrete time setting, we will show that this property allows us to simplify the expression of the transfer matrix to end up with a Markov matrix which has an easy interpretation in the parallel update picture. We also require the R-matrix to satisfy the unitarity property where R 21 (z) = P R 12 (z)P . This property is used to show the commutation of the transfer matrix below 4 and it appears as a consistency relation when constructing the stationary state in a matrix product form, see subsection 2.3.

Remark 2.1
Up to now, we wrote all the definitions and properties about the R-matrix assuming that the spectral parameter is multiplicative. We can easily translate them to an additive spectral parameter: • the Yang-Baxter equation reads • the Markovian property reads σ| ⊗ σ|R(z) = σ| ⊗ σ| • the regularity property reads R(0) = P • the unitarity property reads R 12 (z)R 21 (−z) = 1 We are now equipped to construct the inhomogeneous transfer matrix for a model with periodic boundary conditions The transfer matrix acts on the tensor space V ⊗L (i.e on the whole lattice). The subscripts indicate on which components of the tensor space (i.e on which sites of the lattice) the operators are acting non-trivially. Note that an auxiliary space 0 is introduced and traced out at the end. The main feature of this transfer matrix is that it commutes for different values of the spectral parameter [t(x|z), t(y|z)] = 0. (2.13) It can be shown using the Yang-Baxter equation (2.7) and the unitarity property (2.10). The transfer matrix is commonly used in the context of exclusion processes (respectively of quantum spin chains) without inhomogeneity parameters, i.e with all z i = 1. In this setting the first logarithmic derivative of the transfer matrix w.r.t the spectral parameter generates the continuous time Markov matrix of the exclusion process (or respectively the Hamiltonian of the quantum spin chain). The higher order logarithmic derivatives generates other local charges which commute with the Markov matrix (respectively with the Hamiltonian) because of the commutation relation (2.13).
In the next subsection 2.2, we will use the inhomogeneous transfer matrix with very particular inhomogeneity parameters (different from 1) to define a discrete time Markov matrix.

Remark 2.2
In the case of an additive spectral parameter, the inhomogeneous transfer matrix reads (2.14) Open boundary conditions. In the case of a model with open boundary conditions, the building blocks of integrable models (together with the R-matrix already introduced in the last paragraph) are the reflection matrices K(z) and K(z), both depending on the spectral parameter z and acting on a single copy V of the tensor space (i.e on a single site of the lattice). The matrix K(z) is associated to the left boundary and satisfies the reflection equation The previous equation holds in V ⊗ V. Once again the subscripts indicate on which components of the tensor space the operators are acting non-trivially. For instance K 1 (z) = K(z) ⊗ 1 and K 2 (z) = 1 ⊗ K(z). The matrix K(z) is associated to the right boundary and satisfies the reversed reflection equation We require also that the K-matrices satisfy further properties. We recall that we are interested in constructing discrete time Markov matrices. We therefore impose that the K-matrices fulfill the Markovian property σ|K(z) = σ| and σ|K(z) = σ| (2.17) which is nothing else than demanding that the sum of the entries on each column of the boundary matrices is equal to one. In our construction the K-matrices will indeed play the role of local Markovian operators on the boundaries. We also require the boundary matrices to satisfy the regularity property K(1) = 1 and K(1) = 1 (2.18) and the unitarity property

Remark 2.3
Once again we wrote all the definitions and properties about the K-matrices assuming that the spectral parameter is multiplicative. We can easily translate them to an additive spectral parameter (we give them for the matrix K(z) but similar ones hold also for the matrix K(z)): • the reflection equation reads • the Markovian property reads σ|K(z) = σ| • the regularity property reads K(0) = 1 • the unitarity property reads K(z)K(−z) = 1 We are now in position to define the inhomogeneous transfer matrix for a model with open boundary conditions [46] where the dual boundary matrix K(z) is defined as 5 22) or equivalently The matrix K(z) satisfies the dual reflection equation In the context of continuous time exclusion processes or of quantum spin chains, the transfer matrix is often used with inhomogeneity parameters z i = 1. In this case the derivative with respect to the spectral parameter gives a continuous time Markov matrix (respectively a quantum spin chain Hamiltonian) with two-site local updating rules (respectively with nearest neighbor interactions). We are going to see in the next subsection 2.2 that the inhomogeneity parameters can be efficiently used to define discrete time processes from the transfer matrix. We stress that the construction is quite general in the sense that it does not rely on a particular choice for the R-matrix and the K-matrices (they only have to fulfill the properties mentioned above).

Remark 2.5
In the additive spectral parameter picture, the inhomogeneous transfer matrix reads and the matrices K and K are related for instance through the equation (2.27) 5 We recall that · t i denotes the usual matrix transposition in the i-th tensor space component

Definition of the process.
The goal of this subsection is to use the transfer matrix to define a discrete time Markov matrix. The general idea is to specify particular values of the spectral parameter z and of the inhomogeneity parameters z i in the transfer matrix to obtain Markovian dynamics on the lattice with simple stochastic rules. Several different approaches have already been proposed in the literature. Among the recent developments we can mention the work [34] where a large class of integrable discrete time exclusion processes were introduced using particular values for the inhomogeneity parameters (this was defined for periodic boundary conditions). It would be interesting to investigate what is the physical interpretation, in terms of transition probabilities on the lattice, of all these models. In [12] the transfer matrix was used, without specifying any inhomogeneity parameters, to define a discrete time process in the specific case of the totally asymmetric exclusion process for both periodic and open boundary conditions. A physical interpretation in terms of a sequential update was provided. Note that a similar result has been obtained earlier in the periodic homogeneous case in [23].
In fact the use of the transfer matrix approach in the definition of discrete time exclusion processes follows the pioneering work of Baxter [5]. It was soon realized that the inhomogeneous transfer matrix of the six-vertex model can be used to define an asymmetric exclusion process with discrete time parallel update on the periodic lattice [29,45,43,44]. Our construction is inspired by these works and aims to provide a general framework (independent of the R matrix considered) to define an integrable discrete time Markov matrix from the transfer matrix machinery of integrable systems. Moreover we generalize it also for open boundary conditions using the Sklyanin transfer matrix [46]. The discrete time dynamics that we obtain for both periodic and open boundary conditions appears to be two-step Floquet dynamics. Note that another approach has been developed recently in [24] to construct integrable continuous time Floquet dynamics for quantum systems from the transfer matrix formalism.
Periodic boundary conditions. To define our process on the periodic lattice, we need to consider a lattice with an even number of sites L. We fix z 1 = z 3 = · · · = z L−1 = 1 κ and z 2 = z 4 = · · · = z L = κ. The inhomogeneous transfer matrix then takes the following staggered form Note that a similar staggered construction already appeared in a monodromy matrix in [19]. Let us introduce the local Markovian operator (acting on two sites) whereŘ(z) = P.R(z). Note that U also fulfill the Markovian property σ|U = σ|. A straightforward computation 6 yields We can thus define a Markov matrix M in the following way The property that expression (2.32) together with (2.33) fulfills the requirement of a Markov matrix, i.e non-negativity and sum to one property of the entries, follows from the fact that the entries of the local operator U are non-negative and sum to one on each column (thanks to the Markovian property (2.8) of the R-matrix). Note that the local operators U involved in the definition of U o all commute one with each other and thus the order in the product is not important. The same is true for U e . From the expression (2.32) it is clear that the Markov matrix M can be interpreted as describing a Floquet dynamics composed of two steps: the first one embodied in U e and the other one in U o . We stress that U e and U o do not commute. The operator U e updates every pair of consecutive sites with the first site located at an even position. The local stochastic update is encoded in the matrix U acting on V ⊗ V. We will give precise examples of such matrices in subsections 2.4 and 3.1. The operator U o updates every pair of consecutive sites with the first site located at an odd position. A pictorial representation of the discrete time dynamics is given in figure 1. In the exclusion process literature this kind of dynamics is sometimes referred to as discrete time parallel update, see for instance [43,44]. The same structure also appeared in the propagator of the quantum Hirota model [19] and in a deterministic cellular automaton [35]. The model is said to be integrable because the Markov matrix M commutes with a whole family of operators generated by t(z) which is direct consequence of the commutation relation (2.13) and of the definition (2.32).
Open boundary conditions. To define our process on a lattice with open boundaries, we need to consider a lattice with an odd number of sites L. The idea is essentially the same as for the periodic boundary conditions case. We fix z 1 = z 3 = · · · = z L = κ and z 2 = z 4 = · · · = z L−1 = 1 κ . The inhomogeneous transfer matrix then takes the following staggered form where B = K(κ) and B = K(1/κ). Once again we can define a Markov matrix as The indices in the formula have to be understood modulo L, for instance 0 ≡ L and L + 1 ≡ 1. with operators U e and U o having similar expression as in the periodic case but with additional boundary terms The local operators involved in U o (respectively in U e ) commute one with each other because they are acting on different sites of the lattice. U o realizes a stochastic update on every pair of consecutive sites with the first site located at an odd position, and also a stochastic update on the last site L (describing an interaction with reservoir at the right boundary). U e similarly realizes a stochastic update on every pair of consecutive sites with the first site located at an even position, and also a stochastic update on the first site 1 (describing an interaction with reservoir at the left boundary). Note that the operators U o and U e do not commute and can be interpreted as describing a two-step Floquet dynamics. A pictorial representation of the discrete time dynamics is given in figure 2. This discrete time update is very similar to the sublattice-parallel update [39]: the difference is that in the present model the size of the lattice L has to be odd instead of even in the sublattice-parallel update with open boundaries. Once again the model is said to be integrable because the Markov matrix M commutes with a whole family of operators generated by the transfer matrix The goal is now to use the integrable machinery to handle exact computations for these models.

Stationary distribution in a matrix product form.
It turns out that the integrability framework can be efficiently exploited to construct exactly in a matrix product form the stationary state of the models. The general method has been proposed in [42,13] and we sketch here the main points. We will focus on the open boundaries case but the periodic boundaries case can be treated similarly and we just provide the solution at the end of the subsection. We introduce a (s + 1)-component vector A(z) depending on a spectral parameter and with algebraic-valued (i.e non-commutative) entries. These algebraic elements will be basically the matrices entering the matrix product construction of the stationary state. The construction relies essentially on two key relations. The Zamolodchilov-Faddeev (ZF) relation encodes the commutation relations of the matrix ansatz algebrǎ The associativity of the algebra is ensured by the Yang-Baxter equation (2.7). Another consistency relation is ensured by the unitarity property (2.10). The Ghoshal-Zamolodchilov (GZ) relations encode the action of the matrix ansatz algebra on the boundary vectors W | and |V (which are used to contract the matrix product into a scalar) The consistency between these actions and the commutation relations of the matrix ansatz algebra (i.e the ZF relation) is ensured by the reflection equations (2.15) and (2.16).
We can now introduce the matrix product state where the algebraic element C(z) is defined as the sum of the components of the vector A(z) It has been shown in [13] that if the ZF and GZ relations (2.40) and (2.41) are fulfilled, then where t(z) is the inhomogeneous transfer matrix defined in (2.21) with inhomogeneity parameters z 1 , . . . , z L . On specific models, it is even possible to prove using symmetry and degree arguments that |S(z 1 , z 2 , . . . , z L ) is an eigenvector of t(z) for all z (but with an eigenvalue possibly different from 1), see for instance [13,48]. For our process we recall that we took a specific staggered picture for the inhomogeneity parameters, see subsection 2.2. The steady state of the process is thus given in the staggered form (similar to the transfer matrix (2.35)) where Z L = Z L (κ, 1/κ, . . . , κ). Let us briefly show that it provides the correct stationary state.
The following relations indeed hold (they are a direct consequence of (2.40) and (2.41)) (2.46) Using these properties, it is then easy to check that where It thus implies M |S = |S . Note that this kind of staggering mechanism for the stationary state has already appeared in [37].

Remark 2.6
For an additive spectral parameter, the stationary state is given by We recall also that the construction of the stationary state exposed in this subsection holds only for models with open boundary conditions. For models with periodic boundary conditions, we would need to adapt the construction of the stationary state where tr(·) denotes a trace operator (i.e satisfying the cyclic property) on the auxiliary space of the matrix ansatz algebra. π N is the projector (in the physical space of configurations) on the sector with N particles on the lattice. Indeed with periodic boundary conditions, the number of particles in the system is often conserved and the configuration space can be decomposed into different sectors (labeled by the number of particles) that are stable under the action of the Markov matrix. Z L,N is the appropriate normalization (so that the sum of the components of |S N is 1).
To ease the presentation, and to avoid dealing with these different sectors, we will focus on systems with open boundary conditions in the rest of the paper. We will provide explicit examples of the construction of discrete time models and of the associated stationary state.

Examples: the SSEP and ASEP case.
In this subsection we apply the generic construction introduced above in subsection 2.2 to the R and K matrices associated with the SSEP and the ASEP. We define the discrete time stochastic dynamics and provide explicitly the transition probabilities between the configurations. The stationary states are then computed following the general procedure presented in 2.3. The matrix ansatz algebra is found to be closely related to the one used to solve the continuous time models.
The SSEP case. We are first interested in the SSEP, which is a stochastic model where particles can jump to the left or to the right with equal probability (i.e there is no driving force in the bulk) and with an exclusion constraint. We introduce an R-matrix which reads in the basis |0 ⊗ |0 , |0 ⊗ |1 , |1 ⊗ |0 and |1 ⊗ |1 (ordered this way) It satisfies the Yang-Baxter equation (2.11), the Markovian, regularity and unitarity properties. The R-matrix will be used below to define a local Markovian operator in the bulk of the system. The jumping probabilities to the left and to the right will be equal because of the following property of the R-matrix: P R(z)P = R(z). We also introduce two reflection matrices which read in the basis |0 , |1 (ordered this way) (2.52) They satisfy the reflection equation (2.20) and the reversed reflection equation (2.16) respectively (with additive spectral parameters). They also fulfill the Markovian, regularity and unitarity properties. They will be used to define local Markovian operators acting on the first and last sites, which describe the interaction with particle reservoirs at the boundaries. Following the general procedure introduced in subsection 2.2 we define the local Markovian operators using the R and K matrices (2.55) It may for instance describe the interaction with a particle reservoir at fixed density. On the right boundary the dynamics encoded by the matrix B reads (2.56) Now that the discrete time model is precisely defined, we can look for its stationary distribution in a matrix product form following the general method presented in subsection 2.3. We only need to construct a vector A(z) satisfying the ZF (2.40) and GZ (2.41) relations (with additive spectral parameters).
For such a purpose, we introduce the vector A(z) with algebraic entries The matrices E and D satisfy the following commutation relation We stress that, despite the similitude of the underlying algebraic relations, the stationary distribution of the present discrete time SSEP model is not the same as the one of the continuous time model because of the staggered construction 8 We can use this matrix product solution to compute physical quantities. The first step is to compute the normalization Z L . We recall that it is defined as It is exactly the expression of the normalization for the continuous time SSEP [14,7]. We thus have where the function Gamma satisfies Γ(z + 1) = zΓ(z). We can also compute the mean particle density on site i. We have to be careful because of the two-step Floquet dynamics: in the stationary regime the system is for one out of two time steps described by the probability vector |S and for one out of two time steps described by the probability vector |S ′ (defined in (2.48)).
In the stationary regime the mean particle density (averaged also over the two time steps of the Floquet dynamics) is thus given by where the last equality is obtained using again the results derived for the continuous time SSEP matrix ansatz algebra [14,7]. It appears to be identical to the continuous time SSEP density. Let us mention here that the similarity of the physical observables between the discrete time and continuous time SSEP is specific to this model. We will see in the next paragraph that the results are different between discrete and continuous time in the ASEP case. Note that the result (2.63) can be obtained without using the matrix product structure because the mean densities satisfy a set of closed equations, similarly to the continuous time case. The matrix product solution provides nevertheless an efficient tool to compute multi-point correlation functions. The computation of the mean particle current can also be handled. We have to be careful again because of the Floquet two steps dynamics: particles can only jump from site i to site i + 1, where i is odd (respectively even), during the first time step (respectively the second time step) encoded by U o (respectively by U e ). It implies that when i is odd we have to investigate the action of U o on sites i and i + 1 of the vector |S , whereas when i is even we have to investigate the action of U e on sites i and i + 1 of the vector |S ′ . It turns out that both cases reduce consistently to the same following computation Due to the factor 2κ it differs slightly from the corresponding result for the continuous time SSEP.
The ASEP case. We are now interested in the ASEP, which is a stochastic model where particles can jump to the left or to the right with an asymmetric probability (which can describe a driving force in the bulk) and with an exclusion constraint. We introduce a R-matrix which reads in the basis |0 ⊗ |0 , |0 ⊗ |1 , |1 ⊗ |0 and |1 ⊗ |1 (ordered this way) It satisfies the Yang-Baxter equation (2.7), the Markovian (2.8), regularity (2.9) and unitarity (2.10) properties. The R-matrix will be again used below to define a local Markovian operator in the bulk of the system. Note that the R-matrix of the SSEP (2.51) can be recovered from the R-matrix of the ASEP (2.68) taking the limit R SSEP (z) = lim h→0 R ASEP (e hz )| t 2 =e h . We also introduce two reflection matrices which read in the basis |0 , |1 (ordered this way) (2.69) They satisfy the reflection equation (2.15) and the reversed reflection equation (2.16) respectively. They also fulfill the Markovian (2.17), regularity (2.18) and unitarity (2.19) properties. They will be again used to define local Markovian operators acting on the first and last sites, which describe the interaction with particle reservoirs at the boundaries.
Following one more time the general procedure introduced in subsection 2.2 we define the local Markovian operators using the R and K matrices It allows us to define an integrable discrete time Markov matrix defined by (2.37). We now give an explicit description of the local stochastic rules of the model. In the bulk the dynamics encoded by the matrix U explicitly reads (2.71) Note that the dynamics is asymmetric as announced: the particles jump to the left and to the right with different probabilities. It can for instance describe a driving force in the bulk. Note that the SSEP dynamics is recovered by setting κ → e hκ , t 2 → e h and taking the limit h going to zero. On the left boundary the dynamics encoded by the matrix B reads (2.72) It may for instance describe the interaction with a particle reservoir at fixed density. On the right boundary the dynamics encoded by the matrix B reads (2.73) Now that the discrete time dynamics has been introduced, we can look for its stationary distribution in a matrix product form following the general method presented in subsection 2.3. Again, we only need to construct a vector A(z) satisfying the ZF (2.40) and GZ (2.41) relations. To achieve this goal, we define the algebraic-valued vector A(z) by the following expression The generators e and d satisfy the following commutation relation Similarly to the symmetric case, we stress that, despite the similitude of the underlying algebraic relations, the stationary distribution of the present discrete time ASEP model is not the same as the one of the continuous time model because of the staggered construction 9 We can in principle use this matrix product solution to compute physical quantities. However in the present case of the ASEP, we will not handle the computations up to the final result but rather give an idea on how to use the matrix product formalism in this discrete time dynamics framework. Performing a complete derivation of the physical quantities would require dealing with the explicit representation of the algebraic elements e, d , W | and |V and introducing q-calculus and q-hypergeometric functions [41,7] which is beyond the scope of this paper. The first step is to compute the normalization Z L . We recall that it is defined as because C(z) = σ|A(z) = C(1/z). It depends explicitly on the parameter κ and thus differs from the corresponding expression for the continuous time ASEP. It would be interesting to compute the asymptotic behavior of this normalization for large system size L and see how it is modified in comparison to the continuous time [15,41]. As we will see below the asymptotic behavior of the normalization plays an important role in the determination of the phase diagram of the model. We can also study the mean particle density on site i. Similarly to the symmetric case, taking into account the two steps Floquet dynamics, We are finally interested in the mean particle current. We have first to compute the quantity which is related to the mean number of particles jumping over a particular bond per time step. Hence we deduce that It has a similar structure as the mean particle current of the continuous time ASEP but with an additional factor 1/κ − κ. 9 The stationary state of the continuous time ASEP is given by (with appropriate change of parameters) |S = 1 Z L W |A(1) ⊗ · · · ⊗ A(1)|V

Integrable generalized exclusion processes from the fusion procedure.
This section is devoted to the construction of integrable generalized exclusion processes, i.e models for which the exclusion constraint of the SSEP and ASEP (there is at most one particle per site, s = 1) is a bit relaxed: there is at most s particles per site with s > 1. The general idea to construct such models is the following. The R-matrix associated to the ASEP (respectively to the SSEP) can be viewed as the fundamental representation (also known as spin 1/2 representation) of the universal R-matrix associated to the quantum group U q (ŝl 2 ) (respectively to the Yangian Y(ŝl 2 )). We can construct the R-matrices associated to generalized exclusion processes by taking higher dimensional irreducible representations (i.e higher spin representations) of the universal R-matrix. A well-known technique, called fusion procedure [30,33,32,28], allows us to construct easily the higher dimensional representations of the R-matrix starting from the R-matrix in fundamental (spin 1/2) representation. We summarize briefly the key points of the procedure in the next subsection 3.1 because it will shed some light on the fused matrix ansatz construction of subsection 3.2.

Fusion of R and K matrices and description of the new models.
We wish to give in this small paragraph some insight on the fusion procedure. It is not meant to be fully rigorous. The reader is invited to refer to [30,33,32,28] and references therein or to the specific examples developed hereafter to get into the technical details. Let us first introduce some notations and properties associated with quantum groups. The discussion holds for the quantum group A = U q (ŝl 2 ), corresponding to the asymmetric case, and also for A = Y(ŝl 2 ) corresponding to the symmetric case. The quantum group is endowed with an algebra homomorphism ∆ : A → A ⊗ A called the coproduct, which allows us to construct easily tensor products of representations (see below). There exists an element R ∈ A ⊗ A called universal R-matrix, satisfying the Yang-Baxter equation (3.1) It satisfies also convenient relations with the coproduct, such as (1 ⊗ ∆)(R) = R 12 R 23 . The irreducible finite dimensional representations π (m) z : A → End(V m ) are labeled by a half-integer m called spin, and depend on a spectral parameter z. The dimension of the representation is equal to 2m + 1. Taking finite dimensional representation of the element R allows us to construct a finite dimensional matrix acting on V m ⊗ V l (i.e on two sites with spin m and j respectively) satisfying the Yang-Baxter equation Note that the R-matrix (3.2) depends only on the ratio 10 z 1 /z 2 . The R-matrix (2.68) corresponds to the case m = l = 1/2 and is equal to R (1/2,1/2) (z).
The fusion procedure consists in constructing the matrix R (m,l) (z), for general half-integers m and l, from R (1/2,1/2) (z), without having to deal with the universal R-matrix R which is a rather complicated object in practice.
For instance we would like to construct the representation π (1) z from the tensor product of the fundamental representations π (1/2) z 1 and π (1/2) z 2 which is defined with the help of the coproduct by (π The latter tensor product representation is irreducible for generic spectral parameters z 1 and z 2 but it becomes reducible when the ratio z 2 /z 1 = µ, with µ such that R (1/2,1/2) (µ) is a projector. The representation π (1) z is then obtained by projecting onto the invariant subspace, which is achieved by acting with the projector R (1/2,1/2) (µ). It yields for instance the formula Such kind of construction can be applied recursively to obtain R (m,l) (z), see for instance formulas (3.8) and (3.9) below. A similar fusion procedure also exists to deal with open boundary conditions and construct higher dimensional K-matrices [36,20].
We illustrate more precisely these procedures by applying it to the SSEP and ASEP.
The symmetric case. We are first interested in the SSEP case. To start with the fusion procedure for this model, we need to identify a particular point µ at which the R-matrix defined in (2.51) is a projector, i.e R(µ) 2 = R(µ). It is straightforward to check that µ = 1 is the only point that fulfills this condition. It is then easy to realize that R(1) is closely related to the rectangular matrices The rectangular matrices Q (l) and Q (r) will allow us to project on the invariant subspace of the tensor representation (going from a dimension 4 space, obtained as the tensor product of two fundamental representations of dimension 2, to a dimension 3 irreducible subspace). More precisely we can fuse the second space of the R-matrix as follows The resulting matrix R i,<jk> (z) acts on C 2 ⊗C 3 and satisfies some specific Yang-Baxter equation, see [20]. We are now left with the fusion of the first space The matrix R(z) acts on C 3 ⊗C 3 and will correspond to a process with at most 2 particles allowed on each site (i.e it corresponds to s = 2). It has the following explicit expression, given in the basis |0 ⊗ |0 , |0 ⊗ |1 , |0 ⊗ |2 , |1 ⊗ |0 , |1 ⊗ |1 , |1 ⊗ |2 , |2 ⊗ |0 , |2 ⊗ |1 , |2 ⊗ |2 (ordered this way) Using the fact that the matrix R(z) satisfies the Yang-Baxter equation (2.11) and also the properties (3.7), it is possible to show that R(z) satisfies the Yang-Baxter equation 11 It is also straightforward to show, using either the definition (3.9) or the explicit expression (3.10), that the matrix R(z) satisfies the Markovian, regularity and unitarity properties. A similar fusion procedure can also be applied to the K-matrices [20], defining and The matrices K(z) and K(z) both act on C 3 and take the following explicit expressions, written in the basis |0 , |1 and |2 It can also be checked directly using the explicit expression (3.10). and Note that in order to make these matrices easier to read we have replaced the diagonal entries (which do not enjoy a factorized form) by * . We also have used [·] to denote usual brackets. The value of the diagonal entries can be easily deduced from the rest of the matrix using the Markovian property: the sum of the entries on each column is equal to one. The regularity and unitarity properties are also fulfilled. These matrices satisfy the reflection equations (it can be checked from the explicit expressions (3.14) and (3.15)) and We are now equipped to describe the dynamics of the model, i.e the transition probabilities between the different configurations. We recall that the local stochastic dynamics is encoded in the following matrices As explained in subsection 2.2 the stochastic dynamics on the whole lattice is then defined by the Markov matrix ; 11 −→ 20 κ (2κ + 1)(κ + 1) ; 11 −→ 11 2κ 2 + κ + 1 (2κ + 1)(κ + 1) ; (3.20) Note that we can observe a left-right symmetry in the transition probabilities. It is explained by the fact that U = PU P, where P is the permutation operator in C 3 ⊗C 3 . We can also observe that the transition probabilities are not simply proportional to the number of particles lying on the departure site but depend in a non-trivial way on the composition of both departure and arrival sites. It thus describes non trivial interactions between the particles. On the left boundary the dynamics encoded by the matrix B reads 0 −→ 1 8aκ (2κ − 1)(c − a) + 2 (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; 0 −→ 2 8a 2 κ(2κ − 1) (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; 1 −→ 0 4cκ (2κ − 1)(c − a) + 2 (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; 1 −→ 2 4aκ (2κ − 1)(a − c) + 2 (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; 2 −→ 0 8c 2 κ(2κ − 1) (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; 2 −→ 1 8cκ (2κ − 1)(a − c) + 2 (2κ − 1)(a + c) + 2 (2κ + 1)(a + c) + 2 ; (3.21) It describes a non-trivial interaction with a particle reservoir. On the right boundary the dynamics encoded by the matrix B reads As far as we know, this model is the first example of integrable generalized exclusion process with open boundary conditions. It emphasizes the width of applications of the generic procedure presented in subsection 2.2 to construct discrete time models.
The asymmetric case. We are now interested in the ASEP case. To apply again the fusion procedure on this model, we need to identify the particular point µ at which the R-matrix defined in (2.68) is a projector. It is possible to establish that µ = t 2 is the only point that fulfills this condition. It is then easy to realize that R(t 2 ) is closely related to the rectangular matrices through similar relations as in the symmetric case The rectangular matrices Q (l) and Q (r) will again allow us to project on the invariant subspace of the tensor representation. Similarly to the symmetric case, we can fuse the second space of the R-matrix as follows The matrix R i,<jk> (z) acts on C 2 ⊗ C 3 and satisfies also a Yang-Baxter equation. We are then left with the fusion of the first space The matrix R(z) acts on C 3 ⊗ C 3 and will correspond again to a process with at most 2 particles allowed on each site (i.e it corresponds to s = 2). It has the following explicit expression, given in the basis |0 ⊗ |0 , |0 ⊗ |1 , |0 ⊗ |2 , |1 ⊗ |0 , |1 ⊗ |1 , |1 ⊗ |2 , |2 ⊗ |0 , |2 ⊗ |1 , |2 ⊗ |2 (ordered this way) The dynamics of the model can be considerably simplified by considering the particular value κ = t. The transitions involving the simultaneous jumps of two particles have in this case a vanishing probability.

Fused matrix ansatz.
We introduce in this subsection a new method to construct explicitly the stationary state of the (symmetric and asymmetric) generalized exclusion processes in a matrix product form, using the matrices involved in the solution of the strict exclusion case (when there is at most one particle per site, i.e for s = 1). We call this method, which is heavily inspired by the fusion procedure, "fused" matrix ansatz. We apply it to the generalized exclusion processes (where at most two particles are allowed on the same site, i.e where s = 2) introduced in subsection 3.1. We stress nevertheless that the fused matrix ansatz technique appears as quite general and may be adapted to solve also higher dimensional models (i.e for s > 2).
The symmetric case. Following the approach presented in subsection 2.3, the goal is to construct a vector A(z) with algebraic entries such that the ZF and GZ relations are fulfilled. The stationary state is then easily constructed from this vector. For such a purpose we take advantage of the particular construction of the fused matrices R, K and K and of the fact that we already know a vector that fulfill the ZF and GZ equations associated to the matrices R, K and K for the SSEP (i.e for s = 1). The idea is to mimic the fusion procedure by defining where A(z) is the vector, defined in (2.57), used to construct the stationary state of the SSEP (with at most one particle per site, i.e for s = 1). Writing it out explicitly it has the following expression From expression (3.40), a direct computation (using (3.7) and (2.40)) shows that Using (2.41) it is also possible to show that Following the procedure explained in subsection 2.3, we know that the stationary state of the symmetric generalized exclusion process is given by We will see below how we can use this algebraic structure to compute physical quantities.

Computation of physical quantities.
We stress that, despite the apparent complexity of the matrix product solution of the generalized exclusion processes which we introduced, it is possible to compute interesting physical quantities.
The symmetric case. Once again the first quantity we would like to evaluate is the normalization Z L , which is defined as Z L = W |C(κ)C(−κ) · · · C(κ)|V , with C(z) = σ|A(z). But looking at the expression (3.40) and at the definition of Q (l) , it is straightforward to realize that From that fact and the result (2.62) we deduce directly the following expression of the normalization which is identical to the normalization of the strict exclusion process (i.e with s = 1) with a system size doubled. We can also compute the mean particle density at site i (i.e the mean number of particles at site i). For such a purpose the first step is to evaluate, using the explicit expression Taking into account the averaging over the two steps of the Floquet dynamics (similarly to what we did in the s = 1 case), it follows that We can also compute the mean particle current. Analyzing carefully all the possible local transitions, we reach, after a long computation, the expression (3.57) The asymmetric case. We mention briefly some computations that can be done in the asymmetric case. Similarly to the s = 1 case, we will not enter the full details of the computations to avoid dealing with the explicit representation of the matrix ansatz algebra (which would require a separate study). We rather focus on showing that most of the computations are in fact very close to the s = 1 case. The first quantity we would like to compute is the normalization Z L , which is defined as Z L = W |C(κ)C(1/κ) · · · C(κ)|V , with C(z) = σ|A(z). But looking at the expression (3.45) and at the definition of Q (l) , we see that We can also provide a partial computation for the mean particle current. Analyzing again all the possible local transitions, we end up after a long computation to J = 1 κ − κ W | κ t + t κ + e + d L−1 2 κ + 1 κ + t + 1 t (e + d) κt + 1 κt + e + d L−1 |V W | κ t + t κ + e + d L κt + 1 κt + e + d L |V .
(3.63) Once again, using the explicit representation of the matrix ansatz algebra, we should be able to compute the asymptotic behavior of the current and draw the exact phase diagram of the model. The matrices e and d are indeed the same as the ones appearing in the solution of the continuous time ASEP and the asymptotic behavior of (3.63) should be derived through the study of contour integral expression similarly to the continuous time ASEP (see [7] for a review). This will be the subject of a future work.

Conclusion.
There remains a lot of questions to investigate, both from the physical and the mathematical point of view.
The most direct work is to analyse the physical behavior of the generalized exclusion processes which we introduced. We only gave a very brief flavor of the computation of physical quantities that could be performed using the matrix product solutions. As already mentioned, we need to perform the computation of the mean particle density and current in the asymmetric case, using the explicit representation of the matrices e and d. The asymptotic behavior of the current would yield the phase diagram of the model. In particular it would be interesting to study more precisely what are the precise physical effect of i) the discrete time update in comparison to the continuous time update ii) the small relaxation of the exclusion principle (i.e having s = 2) in comparison to the strict exclusion principle (s = 1).
We are also interested in computing more involved physical quantities such as the large deviation functional of the density profile in the stationary state. The starting point would be to write an additivity formula for the stationary weights using the matrix product structure, similarly to what was done in the strict exclusion (s = 1) case (see [14] for a review). The large deviation functional of the density profile has only been computed for a few different models including the SSEP [16], the ASEP [17], and more recently for a multi-species SSEP [47]. The computation of this functional in the generalized exclusion case could give more insight on the general structure of such functionals.
It would also be interesting to give a hydrodynamic description of the symmetric generalized exclusion process (introduced in subsection 3.1) in the Macroscopic Fluctuation Theory framework [6]. For such a purpose a challenging problem would be to compute the diffusion constant and the conductivity associated to the model. Some promising results have been already obtained in similar but different models [2,3,4]. The exact expression of the large deviation functional of the density profile could also be of great help to guess or check the expression of the transport coefficients (using the predictions of the Macroscopic Theory for the large deviation of the density profile). All these aforementioned problems will be the object of future works.
On a longer term perspective, we should investigate the construction of higher dimensional generalized exclusion processes, where more than two particles are allowed on the same site (i.e where s > 2). The fusion procedure can be again applied and interesting results have already been obtained concerning the exact expression of the fused R-matrix at a particular point [34]. The difficult part of the work would be to handle the fusion of the K-matrices (i.e deriving closed formulas for the entries). It should be possible to express the stationary state of these models using a fused matrix ansatz.
From a more mathematical perspective it would be interesting to explore the connexion of these new models with the theory of multivariate orthogonal polynomials. Some recent progress has been made in this direction in the strict exclusion case [9]. Another challenging task would be to investigate duality relations for these generalized exclusion processes. The pioneering work [44] led to many promising developments [27,21,8,31,10] which provide us with new tools for deriving duality relations.