Integrability of two-species partially asymmetric exclusion processes

We work towards the classification of all one-dimensional exclusion processes with two species of particles that can be solved by a nested coordinate Bethe Ansatz. Using the Yang-Baxter equations, we obtain conditions on the model parameters that ensure that the underlying system is integrable. Three classes of integrable models are thus found. Of these, two classes are well known in literature, but the third has not been studied until recently, and never in the context of the Bethe ansatz. The Bethe equations are derived for the latter model as well as for the associated dynamics encoding the large deviation of the currents.


Introduction
In statistical physics far from equilibrium, the one-dimensional asymmetric simple exclusion process (ASEP) is a paradigmatic model, playing a similar role to the Ising model for equilibrium systems [1,2,3]. Its simplicity means that it can be used as a generic model in many fields, including biophysics [4,5,6], traffic modelling [7] and microfluidics [8]. At the same time, a multitude of exact results have been derived, which is rare for interacting N-body problems. The two main approaches for obtaining exact results have been the matrix product ansatz [9] and the Bethe ansatz (BA) [10,11,12]. For general reviews of these methods in the ASEP context, see [2] for the matrix product approach (which will not be used in the present work) and [13] for the BA.
The BA was first developed to diagonalize one-dimensional quantum spin chain Hamiltonians [14], but has since been extended to many other fields, such as 2D vertex models [15] and ASEPs. It was first shown that it can be applied to the exclusion process on a ring [10,11,12]. Since then, it has also been used on systems with various boundary conditions [16,17,18,19], partially asymmetric (PASEP) models [20,21,22,23,24,25] and multi-species cases [26,21,22,27]. Using a deformation of the Markov operator, it was first shown in [28] that the BA can be used very effectively to derive higher order particle displacement statistics, such as the diffusion coefficient. For some particularly simple cases, it has even been possible to calculate the full statistics to all orders [28,25].
The question of which systems can be solved by the BA, i.e. are integrable, has a long history in the quantum and classical equilibrium contexts [29]. An answer to this problem is through the Yang-Baxter equations (YBE), which can be interpreted as consistency conditions for the BA; indeed, by checking whether the YBE hold on a general class of models, it is possible to select values of the parameters for which exact solutions exist. Though this approach was used for quantum [30] and vertex model [31], it has been applied less systematically to study classical interacting particles systems, see for example [32,27,33,34].
In this paper, we conduct such a systematic study for two-species exclusion processes on a ring. By postulating a very general coordinate BA and deriving the YBE, we are able to assess the integrability of a wide range of models, which suggests a classification of integrable two-species systems.
We find that there are three families of integrable models: two of them are wellknown in the literature, but there exists another, special, model that has only come to attention recently [35] and has not yet been studied from the viewpoint of integrability. There, its steady state was solved using the matrix product approach, and the phase diagram, currents and density profiles were calculated. This model has the peculiar feature that it is only integrable when there is only one particle of one of the species, which can be seen clearly in the YBE. We conjecture that our study entirely exhausts all integrable two-species exclusion processes.
The remainder of the paper is structured as follows. In section 2, we show how the time evolution problem of a two-species PASEP can be solved by a general form of the coordinate BA. In addition, we show that a very similar ansatz can be used to solve the conditioned (or deformed) time evolution problem, which allows direct calculation of the large deviations of the current. In section 3, we derive the YBE and present the three classes of solutions. Two of these are well-known in literature [36,37,28,21,27], but the third has only come to attention recently [35] and has not yet been studied in the context of the BA. In section 4, we focus on this new model and derive the Bethe equations for the deformed time evolution problem. In Appendix A, we explain the procedure that we used to find the solutions of the YBE. Finally, in Appendix B, we fill in some algebraic steps in the calculation of section 4.

Problem statement
We consider an L site ring with particles of two-species hopping stochastically in continuous time. Let there be M 1 particles of species 1 and M 2 particles of species 2, with M 1 + M 2 = M. For the most general two-species partially asymmetric problem with overtaking, the following processes take place with the given rates: where 0, 1, 2 represent empty sites, and particles of species 1 and 2, respectively.
Fixing the reference frame, we denote by Q j and x j respectively the species and position of particle j, counting from the left, where Q j ∈ {1, 2} and x j ∈ Z L . As a shorthand, we denote the configuration of the system as The time evolution of the probability of a configuration C at time t is given by the master equation where M is the continuous-time Markov operator containing the transition rates from configuration C ′ to C. We wish to solve the eigenvalue-eigenvector problem for the Markov operator. It can be written in the form where λ is an eigenvalue and ψ is the corresponding eigenvector. This equation holds in configurations in which no two particles are on adjacent sites. It can be assumed to hold in all configurations, provided an additional set of equations, called collision equations, also hold. Considering the j-th and (j + 1)-th particles, with x j+1 = x j + 1, these are given by the following, where Q j,j+1 = 1, 2, where b 22 = b 11 = 0. Finally, the periodic boundary conditions are imposed via the requirement Then equations (2), (3), (4) define the eigenvalue-eigenvector problem to be solved.

Bethe ansatz solution
We use the following modified nested coordinate BA, where σ indicates the permutations of the integers 1, . . . , M; the A's and z's are undetermined complex amplitudes and Bethe roots respectively; and the r's are species dependent pre-factors, to be determined. This is the most general known coordinate ansatz and is a straightforward nested generalization of the ansatz introduced in [26]. There, it was used for a totally asymmetric problem with one defect particle. The ansatz (5) is found to solve (2) provided the following condition holds: In this case, we can set for some constant α. Then the ansatz (5) works, with We assume α = 0, as otherwise the ansatz vanishes. Notably, this excludes the model studied by Arndt et al [38,39], which can be defined with our notation as p 1 = 0, q 2 = 0. However, as this is believed to be non-integrable [40], this is not of major concern. Putting (5) with (7) into (2) gives the following expression for the eigenvalue in terms of the Bethe roots, We remark that condition (6) comes from the form of the ansatz (5). However, a similar condition was derived from the the algebraic Bethe ansatz for the 15-vertex model [31]. These two cases are linked, as with certain choices of parameters, the 15vertex model can be mapped to a two-species PASEP with overtaking. This suggests that condition (6) is generic in the context of exactly solvable models and not a special feature of the ansatz (5).
In the subsequent calculations, it is convenient to use the following dimensionless parameters Putting the ansatz (5) into the collision equations (3) gives where and the notation (i ↔ j) indicates a term similar to the explicitly written one but with the indices i, j swapped.
The system of equations (12a)-(12b) can be written compactly in vector form where This can be formally rewritten as where division is to be understood as multiplication by the matrix inverse from the left. Performing this calculation explicitly, we get relations of the form Of these, only the following are non-zero where

Bethe ansatz for deformed Markov operator
Large deviations in interacting particle systems can be derived by conditioning the dynamics on a given, atypical value, of the observable under consideration. The resulting evolution is governed by a deformation of the original Markov operator. This method goes back to Donsker and Varadhan (see [41,42] and references therein). In the context of ASEP, this idea has been has been used systematically in conjunction with the BA, leading to exact results for the large deviation function of the current [28,26,24,43,44]. We now outline this calculation for the two-species case, merely for completeness, since much of the detail is similar to the undeformed case, analyzed in the previous section.
We define observables Y 1 t , Y 2 t , Y 12 t , which count the number of processes 10 → 01 , 20 → 02 , 12 → 21 minus the inverse processes respectively, up to time t. Thus, for example, is the net total displacement of particles of species 1 to the right. We can split the Markov operator M into where M 0 contains the diagonal entries and M i,+(−) contains the rates of the processes Its generating function is defined as Note that by marginalizing over configurations, we obtain the moment generating function of A key property is that this satisfies a large deviation principle as t → ∞, with some rate function µ(γ 1 , γ 2 , γ 12 ). From (1), we get the evolution equation where M γ is the deformed Markov operator The eigenvalue-eigenvector problem for M γ then becomes As with the undeformed case, this equation is strictly valid only for configurations with no two particles being on adjacent sites, which is why γ 12 does not appear here. The problem has to be supplemented with collision equations, which are given below.
A particularly important eigenvalue, is the one with the largest real part, which can be identified with the rate function µ(γ 1 , γ 2 , γ 12 ) in (25). From the Perron-Frobenius theorem it follows that this eigenvalue is unique. In the limit γ 1 → 0, γ 2 → 0, γ 12 → 0, it converges to the unique 0 eigenvalue of the (undeformed) Markov operator M, whose eigenvector is the stationary measure.

The collision equations read
where γ 21 = −γ 12 . This problem is solved by the ansatz Plugging the ansatz (30a)-(30b) into the collision equations (29) gives after some simplification an equation similar to (18), with the only difference being that now Then repeating the calculation of the form factors, we get that most of them are the same as in the undeformed case, with only the following modifications 3. Yang-Baxter equations

Derivation of YBE
The Yang-Baxter relations are obtained by analyzing collisions that involve 3 particles. In this case, two sequences of binary collisions can result in the same permutation: ijk → ikj → kij → kji and ijk → jik → jki → kji. In order for the nested coordinate BA to be consistent, the products of form factors obtained from these two sequences must be the same. This requirement yields the YBE (see [21] for a similar analysis). We shall follow this procedure both for the regular and deformed cases, with the only difference being the explicit forms of the form factors. Somewhat surprisingly, for the systems at hand, the results are insensitive to the deformation parameters. The consistency condition gives some relations between the amplitudes , which for the two-species case is in principle an equation for 8 × 8 matrices, called R-matrices. The YBE can be formulated as an algebraic relation for the Rmatrices. In the context of driven diffusive systems this was done for instance in [27]. However, due to particle conservation, the relations can be simplified by dividing them into disjoint blocks with a fixed number of particles of each species. It is not hard to .
Then the YBE are given by The YBE for the 122 block can be derived similarly.
We remark that all the entries of the matrices Ξ are rational functions of the Bethe roots. Thus, when we combine the two matrix products in (35), the result is a 3 × 3 matrix whose entries are rational functions of the Bethe roots. Their explicit forms are evidently very complicated but they can be analyzed with the help of a symbolic programming language, such as Mathematica.

Solutions of YBE
We wish to find conditions on the model parameters (α, x, β 12 , β 21 ) that ensure that (35) are satisfied identically, for all values of the Bethe roots.
Firstly, we note that the numerators of the entries are polynomials in z i , z j , z k , which must vanish identically. We may therefore examine individual coefficients of z i , z j , z k in the numerators and require them to vanish. This can be done in many ways, as there are many possible coefficients to look at. Fortunately, taking only a few is sufficient to obtain conditions that are strong enough to satisfy the full set of equations. One possible method is outlined in Appendix A.
We obtain three general solutions. For the sake of conciseness, we do not list onespecies models (which trivially form subsets of the models listed below and are known to be integrable), and models which are equivalent to the given ones up to permutation of species labels (including empty sites as a "zeroth species").
We remark that all the following models satisfy the relation which was not imposed a priori. This shows that, the condition (6) cannot be circumvented by a relabelling of the species, which further supports the claim that it might be a necessary condition.

Solution 1: TASEP case.
The first solution is defined by the conditions which corresponds to the process 10 αp 2 → 01 ; 20 This is a totally asymmetric process, with particles of species 1 hopping at rate αp 2 and overtaking particles of species 2 at rate b 12 . The matrix product solution for the steady state was first given in [36]. The case of a single particle of species 2 (referred to as a defect) has been much studied in the literature [37,26,27].

Solution 2: PASEP with first-and second-class particles.
The second solution is defined by the conditions which corresponds to the process This is a two-species partially asymmetric process, with both species hopping at the same rates and the species 1 particles overtaking species 2 particles with their usual hopping rates. Since species 1 does not distinguish between species 2 and a vacancy, species 1 can be seen as having priority in the dynamics and is therefore referred to in the literature as first-class and species 2 as second-class. This case has been studied in the literature as well, with the matrix product solution for the steady state given in [36] and BA equations derived in [21,22]. We remark that although the Bethe equations have been derived, the statistics of the current fluctuations have not yet been calculated.

Solution 3: PASEP with a single first-class defect.
The third solution is defined by the conditions This is a partially asymmetric process, with a single first-class defect particle with hopping rates αp 2 , q 2 /α that overtakes the second-class particles with its usual hopping rates. This model has not been analyzed in the literature until recently, when its steady state was solved using a matrix product approach [35]. We remark that the constraint on particle number M 1 = 1 is unusual in the context of YBE. It emerges here as the given conditions on the rates are sufficient to satisfy (35) for the 122 block but not the 112 block. This means that the ansatz works for this general case only if there is exactly one particle of species 1 in the system. To satisfy (35) for the 112 block as well, it is necessary to impose either α = 1 or α = x or x = 0. In all these cases, this solution reduces down to solution 1 or 2.

Derivation of Bethe equations for Solution 3
As solution 3 has not appeared much in the literature, it presents a case of particular interest. As mentioned previously, it was studied using a matrix product approach in [35]. The phase diagram was shown to consist of localized phases, in which the defect particle has no macroscopic effect on the system, and a shock phase, in which the defect creates two extended regions with different bulk densities. The density profiles and currents were also calculated. The BA approach presented here complements these results, as it should allow one to calculate the fluctuations of the current as well.
To that end, we now derive the Bethe equations for the deformed eigenvalue problem from section 2.3. As we have only one particle of species 1, it is not necessary to use the full nested ansatz. Instead, we can fix the particle indexing so that the species 1 particle is always the leftmost one. We can always do this using the boundary condition (4), which gives the following relation for amplitudes Now that the order of species is always fixed to be 12 . . . 2, we can drop the upper indices without ambiguity. Then we get for 22, 12 and 21 collisions respectively, We will denote the form factor in (41a) as S 22 22 (z i , z j ) for conciseness. Applying (41a) M 2 − 1 times to (41b) and (41c) gives where In order for the eigenvector not to vanish, the system of equations (42a)-(42b) for A ij... , A ji... must be degenerate. The equation resulting from this requirement can be simplified, although this procedure is algebraically somewhat cumbersome and not very illuminating. We give the details in Appendix B for completeness. In this calculation, it becomes convenient to consider the following transformation of the Bethe roots [20,24], Ultimately, one obtains a set of equations, called the Bethe equations, for the transformed Bethe roots y i , where B is a constant. The constant B can be fixed by multiplying the Bethe equations for all i, where we have implemented the periodic boundary condition (4) as It is interesting to compare (45) to the corresponding equations for the well-known cases of a one-species PASEP [24] and the TASEP with a defect [26,27]. The one-species PASEP can also be simplified using a change of variable similar to (44). The resulting equation involves a product of fractions of linear differences of the transformed Bethe roots y i , like (45).
Meanwhile, the TASEP with a defect also features an additional constant, like B, that has to be fixed as part of solving the Bethe equations. As will be shown in a forthcoming publication, this leads to one of the Bethe roots converging to a different point in the complex plane to the other roots in the γ 1 → 0, γ 2 → 0, γ 12 → 0 limit.
We now briefly comment on how one would proceed to solve the Bethe equations. (45) is a system of polynomial equations of degree L + M + 1 for each Bethe root y i , in which the coefficients in turn depend on all the roots. One would first need to solve (45) consistently with (46), then choose M roots to satisfy (47) for some L-th root of unity. This would then give an eigenvalue through (9). In practice this is very challenging both analytically and numerically. However, special functional methods have been developed that allow to calculate perturbative expansions of the largest eigenvalue in terms of γ 1 , γ 2 , γ 12 [28,26,24,25]. In a forthcoming publication, we will show how these methods can be applied to the problem at hand.

Conclusion
We have presented a nested coordinate BA that solves the two-species PASEP problem with overtaking, in both the deformed and simple (undeformed) cases, provided the condition (6) is satisfied.
The YBE were derived, which allowed us to assess the BA integrability of all models in this class. Although the derivation presented here is not entirely general, as a particular form of the nested coordinate Bethe wave function (5) was used, it is conjectured that the list of integrable models presented is exhaustive.
It was shown that the solutions fall into three cases, namely a TASEP with two species that move at different rates, with one species overtaking the other with an arbitrary rate; a PASEP with two species moving at the same rate and one species treating the other as holes; and a PASEP with a single defect particle that moves with a different rate and does not distinguish between particles and holes. These results are summarized in table 1. Two of these cases are well-known in literature. Namely, the TASEP with defect particles, and the second-class PASEP.

Model
Parameter restrictions Free parameters TASEP with defect species Table 1. All Bethe ansatz integrable sub-classes of the two-species partially asymmetric exclusion process with overtaking on a ring. The rate p 2 is not included as a free parameter as it can be seen as just setting the timescale. The particle numbers M 1 , M 2 are included as free parameters to highlight that M 1 is fixed in the third case.
Solution 3 is a PASEP with a single first-class defect particle and overtaking. This model has only been investigated recently, when its steady state was solved using the matrix product approach [35]. It is peculiar among the solutions listed here in that it is only integrable for a single defect particle. This follows from the fact that the YBE are satisfied for three-body collisions in which the particle species are 122 and 222 but not 112.
It is interesting to note that that among the models examined here, all the models that are integrable in the normal (undeformed) case remain integrable after a deformation of the Markov operator. This is a priori not obvious and it would be of interest to investigate whether integrability is preserved under deformation in other contexts as well.
Finally, we derived the Bethe equations for the deformed time evolution problem of solution 3 of the YBE. This should allow one to calculate the long-time current fluctuations in this model. This would complement the result for the mean current, which has already been calculated using a matrix product ansatz [35]. The Bethe equations have similar features to the well-known one-species PASEP and two-species TASEP, which makes it amenable to the techniques developed for those cases.
Another possible direction for future work would be to consider some extensions of the traditional YBE, which have been proposed in the literature. These include the dynamical YBE [45] and the braided YBE [46,47]. It would be interesting to investigate whether those generalizations can provide further insight in the context of asymmetric exclusion processes. Setting α = 1 is enough to satisfy all YBE but it gives a essentially a one-species PASEP, which is not of interest here, as we are looking for two-species solutions.
Choosing x = 0 gives a special case of solution 1 with β 12 = 0.
Appendix B. Algebraic steps in derivation of Bethe equations for solution 3 Requiring the system of equations (42a)-(42b) to be degenerate, we get D ij − xz j + xz i S 22 22 (z i , z j )E j (D ij − z j )S 22 22 (z j , z i )E i + z i − (i ↔ j) = 0 . (B.1) Expanding and simplifying this, we get Now plugging in the explicit forms of S 22 22 , this becomes −(D ij − (1 + x)z i )E i + (D ij − (1 + x)z j )E j + (z i − z j )(1 + xE i E j ) = 0 . (B.3) With the change of variable (44), we get, after some simplification, where now The terms with index i and j can now be separated to give the following equation, As this procedure can be carried out for any two indices i, j, this means that this expression must equal some constant B for all i, j. Then equating the left-hand side with B and rearranging to isolate E i , we finally obtain (45).