Symmetry, Integrability and Geometry: Methods and Applications Dynamical Critical Exponent for Two-Species Totally Asymmetric Diffusion on a Ring ⋆

We present a study of the two species totally asymmetric diffusion model using the Bethe ansatz. The Hamiltonian has $U_q(SU(3))$ symmetry. We derive the nested Bethe ansatz equations and obtain the dynamical critical exponent from the finite-size scaling properties of the eigenvalue with the smallest real part. The dynamical critical exponent is 3/2 which is the exponent corresponding to KPZ growth in the single species asymmetric diffusion model.


Introduction
The single species asymmetric diffusion process has attracted and continues to attract a lot of interest. It is a driven diffusive system describing a stochastic movement of hard-core particles in one dimension, where movements to the left and right take place with different probabilities. There are many interesting applications to shock formation [1,2,3], traffic flow [4], biopolymers [5,6,7] and other driven diffusive systems 1 . Since the underlying quantum spin chain is integrable, the powerful method of Bethe ansatz is applicable and was widely used in the past to obtain analytical results (see e.g. [8,9,10,11,12]). The two-species asymmetric diffusion problem however has been less investigated. While the stationary state has been studied in some detail [13,14,15,16,17], much less is known about the full dynamics of the system. Some results for the dynamical phase diagram were obtained from the matrix product ansatz [18] (however these are restricted to certain regions of the full parameter space) and from numerical studies of the model [19].
Recently, the multi-species asymmetric diffusion model with open boundaries has sparked interest [20,21,22]. This case can also be tackled using the matrix product ansatz.
In [23] the dynamics of a particular choice of the model is studied where two species hop in opposite directions on a ring, the diffusion constant being different from the passing rate of particles of different species [24]. This model corresponds directly to two coupled Burgers equations or to a Burgers equation coupled to a diffusion equation. The dynamical scaling properties are calculated from the time-evolution of the two-point correlation functions, based on a numerical study. The dynamical critical exponent is consistent with a subtle double Kardar-Parisi-Zhang type scaling (in the sense of factorization).
An analytical study of the dynamical scaling of the multi-species asymmetric exclusion process was described by [25] who found a dynamical critical exponent of 3 2 from a Bethe ansatz calculation if the diffusion was asymmetric and a critical exponent of 2 for symmetric diffusion.
In the present paper, we investigate the dynamical critical exponent for the totally asymmetric diffusion process with equal rates which is a related but dynamically different model that will be precisely defined in Section 2. We will derive the Bethe ansatz equations starting from the Rmatrix of the totally asymmetric diffusion model, calculate the energy gap and from a finite-size scaling analysis extract the critical exponent.

Discussion of the literature about the two-species asymmetric diffusion model
We will briefly survey the literature on the two-species asymmetric diffusion model. Although there seem to be several approaches building on previous work [25,26,27,28] by trying to extract an appropriate limit to totally asymmetric diffusion, the details of each such strategy are not obvious, however, and one runs into serious difficulties and subtleties as we will explain. We found that instead of trying to tweak the calculations in an ad-hoc manner, it is clearer and simpler to start from scratch and derive the Bethe ansatz equations directly for the totally asymmetric case. With hindsight, our calculation provides a firm ground for future comparison which might also facilitate to formulate the correct limiting procedure. The Bethe ansatz equations for the two-species asymmetric diffusion model were put forth in an article by Dahmen [26]. They depend on a parameter γ that is related to the asymmetric hopping rates by exp(γ) = q = Γ L Γ R , where Γ L and Γ R describe the hopping rate to the left and right, respectively. It is not at all clear how to take the limit of one of these hopping rates going to zero directly in the Bethe ansatz equations. Since in the limit q → 0 all the ratios tend to one, and even the usual rescaling of the roots by a factor of q does not work -the resulting equations either diverge or become trivial.
In the paper by Popkov at al. [27] Bethe ansatz equations for three-species models with different hopping rates for all six possible particles exchanges are derived from the matrix product ansatz. One of the cases they discuss is a particular asymmetric exclusion process where the hopping rates are coupled as g B0 = g A0 + g AB where g B0 , g A0 and g AB correspond to the processes A + 0 → 0 + A, B + 0 → 0 + B and A + B → B + A, respectively (and the three other rates are set to zero). In this paper we will choose all three non-zero rates to be the same, so the model considered here is different.
Although the way Cantini [28] derives the Bethe ansatz equations for the case of asymmetric diffusion seems to be the most amenable for taking the limit to our model, however, he has introduced phase factors e ν 10 , e ν 20 , e ν 12 that remain in the Bethe ansatz equations, and it is not clear which value they should take for the case of totally asymmetric diffusion.
In the article by Arita et al. [25] the Bethe ansatz equations for the multi-species asymmetric diffusion model were derived. However, the Bethe ansatz equations for the case of totally asymmetric diffusion are impossible to obtain from this calculation in a straightforward manner. Although it is possible to take the limit where one of the hopping rates is set equal to zero at the starting point, the R-matrix, in the subsequent calculations the authors often multiply by or divide by the parameter that would become zero in the totally asymmetric limit. So taking the limit in the final equation does not lead to reasonable equations, and would have to do a step by step analysis of the derivation and make the necessary alterations at each step to avoid division by or multiplication by zero. At this point, we found it simpler and less ambiguous to just derive the equations for our case directly. 2 Two-species totally asymmetric dif fusion

Asymmetric diffusion
The two-species asymmetric diffusion model consists of two species of particles, A and B, diffusing asymmetrically in one dimension. The following processes take place:

Totally asymmetric diffusion
In the totally asymmetric case, Γ L = 0, so particles A do not diffuse to the left since the interchange B + A → A + B is blocked. This leads to a different dynamics for the model that will be analyzed in this article. Assuming periodic boundary conditions, i.e. on the ring Z/LZ (where L is the number of discrete sites of that ring), the model is shown qualitatively in Fig. 1. Diffusion to the right corresponds to clockwise motion around the circle, diffusion to the left to counterclockwise motion. Black dots represent particles A, grey dots particles B and open dots vacancies. The arrows indicate processes that are still allowed of Γ L = 0, the blocked arrow indicates the process B + A → A + B that is forbidden in the totally asymmetric case.
Particles A do not see any difference between vacant sites and particles B, whereas particles B trying to diffuse to the right are blocked by particles A. Therefore, we call particles A "firstclass", and particles B "second-class" particles. There is an interesting connection between second-class particles and the study of current fluctuations [29].

Master equation and Hamiltonian
The dynamics of the asymmetric diffusion model is described by a master equation for the probability distribution p t (η) of the lattice configuration η(t) at time t. If we denote a configuration by η and the jump rates Γ R and Γ L by c(j, j + 1, η), when interchanging the particles and/or vacancies of the configuration η on sites j and j + 1, the master equation reads: where η jj+1 denotes the configuration obtained from η by interchanging the particles and/or vacancies at sites j and j + 1. Following [30,31], attaching a vector space C 3 at each discrete point j and using a vector (1, 0, 0) T for a particle of type A, a vector (0, 1, 0) T for a particle of type B and a vector (0, 0, 1) T for a vacancy, the operator of the master equation can be written as a quantum spin chain which is given by the following expression if one assumes periodic boundary conditions In the following, the standard vector space R 3 will just be called V . In this expression, the matrices E αβ j are 3 × 3 matrices with only one non-zero entry: (E αβ j ) γδ = δ αγ δ βδ and, as usual, The parameters D and q are real and depend on the diffusion rates: q = Γ R Γ L and D = √ Γ R Γ L . This Hamiltonian is integrable and the eigenvalues can be found by applying the Bethe ansatz. Since the Hamiltonian is non-Hermitian in general, we will encounter complex eigenvalues.

Dynamical critical exponent
In non-equilibrium dynamics, the dynamical critical exponent describes a relation between the relaxation time towards equilibrium τ (or temporal correlation length) of a system and the spatial correlation length ξ, namely that τ ≃ ξ z with the dynamical critical exponent z. It can be shown that for one-dimensional quantum spin chains, τ ≃ L z . Since the relaxation time is dominated by the eigenvalue of the Hamiltonian with the smallest real part (the energy of the ground state being equal to zero), we can obtain the exponent z from a finite size analysis of the Hamiltonian of the system as .
For the single species asymmetric diffusion model, z = 3 2 was obtained in [11] from a Bethe ansatz calculation for the totally asymmetric case which therefore belongs to the KPZ [32] universality class, whereas z = 2 for the partially asymmetric case which describes the Edwards-Wilkinson universality class [33]. We will determine the exponent z from a careful study of (2). We will find the lowest lying eigenvalue of the totally asymmetric diffusion model by means of the Bethe ansatz.

Nested Bethe ansatz
We start with the totally asymmetric diffusion model, setting the parameters Γ L = 0 and Γ R = 1. This does not lead to any singularities in the Hamiltonian since its expression contains the products Dq = Γ R = 1 and Dq −1 = Γ L = 0. The new Hamiltonian reads This Hamiltonian is integrable, and we will use the algebraic Bethe ansatz to find its spectrum. Although the Bethe ansatz is well known for the Hamiltonian (1), the totally asymmetric case (given by the Hamiltonian (3)) cannot be obtained as a special case since setting Γ L = 0 causes the Bethe ansatz equations to become singular. Therefore, we will derive the Bethe ansatz equations for the totally asymmetric case following [34,35]. We use the following R-matrix elements where the indices denote the following tensor product in End(V ⊗ V ): or in the language of the associated vertex model this corresponds to the initial state denoted by the two indices mi scattering into the final state denoted by the indices kl.
This R-matrix satisfies the factorization equation We now define the 3 L × 3 L matrices T (θ) ab as The monodromy matrix satisfies the fundamental relation The transfer matrix is obtained as the trace of the monodromy matrix. It is a matrix acting on V ⊗L . In the above notation, it can be written as We recover the Hamiltonian (3) as the logarithmic derivative of the transfer matrix, if the spectral parameter θ is set to zero: Therefore, diagonalizing H will be the same problem as diagonalizing τ .

Diagonalization of the transfer matrix
We start be deriving algebraic relations for the matrices A, B j , C j , j = 1, 2 and D ik , i, k = 1, 2 in equation (4).
Since we would like to diagonalize τ , we are looking for a reference state that is a simultaneous eigenstate of A and D ii and is annihilated by C i and D ij for i = j.
We choose to be our reference state. This is not necessarily the physical ground state. T [L] acting on this reference state is interpreted as follows: where a(θ) = exp(θ) and c(θ) = 2 sinh(θ). The operators B i (θ, α) will create new states when acting on the reference state. They will be called creation operators in the following.
In the next step we use the fundamental relation equation (5) to derive the bilinear algebra among the operators A(θ), B i (θ), C i (θ) and D ij (θ).
The relations we are going to use later on are: Here r ij pq are coefficients of a 4 × 4 matrix where only the following five entries are different from zero: r 11 11 = 1, r 12 12 = 1, r 21 12 = 2 sinh(θ)e −θ , r 21 21 = e −2θ , r 22 22 = 1 and the functions g(θ) and h(θ) are given by: .
Now we are ready to look for an eigenstate of the transfer matrix τ (θ). We will make the following ansatz for an eigenfunction of the transfer matrix having the creation operators B i (θ), i = 1, 2 act on the reference state Ω given by equation (7).
We know how A(θ) and D ii (θ) commute with the operators B i (λ j ) from the relations (9) and we know how A(θ) and D ii (θ) act on Ω from equation (8).
We will get two types of terms. In the first type of terms, we get the original combination of B σ(y 1 ) (λ 1 )B σ(y 2 ) (λ 2 ) · · · B σ(yp) (λ p )|Ω back. They originate from taking the first term in the commutation relation for B i (θ)B j (θ) in (9). Using standard terminology, we will call those terms wanted terms. In the second type of terms, one of the B σ(y i ) operators will depend on the parameter θ. These terms will be called unwanted terms and their sum has to be set to zero.
We will also use the following notation: We define a 2 p dimensional vector B containing all possible terms of the form B σ(y 1 ) (λ 1 )B σ(y 2 ) (λ 2 ) · · · B σ(yp) where, as above, σ denotes a permutation of the set of indices {y 1 , y 2 , . . . , y p } ∈ {1, 2}. We also define a 2 p dimensional vector X containing the corresponding x σ in the same order. This allows us to rewrite the sum as a scalar product of these two vectors: Then the action of the transfer matrix τ (θ) on Ψ(λ 1 , λ 2 , . . . , λ p ) leads to wanted terms that can be written as: The unwanted terms arise from taking the second term in the commutation relation for the B(θ) operators (9). E.g. one of the unwanted terms appears if the second term in the commutator is applied to A(θ)B(λ k ) and is of the form As observed by de Vega [36] the general unwanted term can easily written down if one uses the following fact about a cyclic permutation of p operators B(λ j ): aa (θ, {λ}) and where [t [2] ab (θ)] ij , a, b, i, j = 1, 2 are 2 × 2 matrices where only the following five entries are non-zero: [t [2] 11 (θ)] 11  Here r ib aj (θ) are the coefficients appearing in equation (9)  So the cyclic permutation B(λ i ) → B(λ i+1 ) followed by the multiplication of the matrix leaves the product B(λ 1 ) ⊗ B(λ 2 ) ⊗ · · · ⊗ B(λ p ) invariant. Therefore the general unwanted term in the result for A(θ)Ψ(λ 1 · · · λ p ) can be written as Similarly one can find the wanted and unwanted terms after acting with (D 11 (θ) + D 22 (θ)) on |Ψ . The wanted term is Using the same argument as before, the general unwanted term can be written as a sum of terms where B(θ) replaces B(λ k ): Altogether, the sum of wanted terms reads The sum of unwanted terms reads: Since we are looking for an eigenstate of the transfer matrix, the sum of all wanted terms has to be proportional to |Ψ and the sum of all unwanted terms has to be equal to zero. The first condition determines X to be an eigenvector of τ (p) (2) corresponding to the eigenvalue Λ (2) (θ, {λ}): The eigenvalue Λ (2) (θ, {λ}) is determined by requiring the sum of unwanted terms to become zero: We have reduced the original eigenvalue problem to equation (10) which is an eigenvalue equation for the transfer matrix τ (p) (2) of the six-vertex model. This is the crucial step in this paper. All we have to do now is repeat the same diagonalization procedure for equation (10). We repeat exactly the same steps as before, acting with the transfer matrix on a reference state with r creation matrices B, and deriving equations for the wanted and unwanted terms. This leads to another expression for the eigenvalue λ (2) that can be set equal to equation (11) and another consistency equation stemming from setting the unwanted terms equal to zero. In this way we obtain the following nested Bethe ansatz equations which have two types of unknowns, λ k , k = 1, . . . , p and Λ j , j = 1, . . . , r, The eigenvalues of H are obtained using the relation between the transfer matrix and the Hamiltonian given by equation (6):

Numerical solution of the Bethe ansatz equations
Since the Hamiltonian (3) is describing a reaction-diffusion system, the ground state is the steady state with energy zero and all excitation energies will have positive real parts. We fixed L, p and r and solved the coupled system of Bethe ansatz equations numerically. We also diagonalized the Hamiltonian (3) numerically for a small number of sites (typically up to L = 9). We compared the energies obtained from the Bethe ansatz with the ones obtained by numerical diagonalization of Hamiltonian H to single out the eigenvalue with the smallest non-zero real part which plays the role of the second lowest eigenvalue and determines the gap. This first excited state always lies in the sector with p = L 3 , r = 0. This state has an equal density of particles A, B and vacancies of 1 3 . We extrapolated the energy values of the first excited state for L → ∞.

Results
We found the solution for p = L 3 , r = 0 corresponding to the first excited state of Hamiltonian. To make sure that this really is the lowest lying excited state, we compared to the results of a numerical diagonalization of the Hamiltonian for finite lattice lengths. This state has an equal density of particles A, B and vacancies of 1 3 . The roots all lie in the complex plane. A typical pattern of roots for the first excited state of the Hamiltonian is shown for L = 36 in Fig. 2. The solid line is the unit circle that is given as a guide to the eye.
The data was extrapolated using the Bulirsch-Stoer algorithm [37]. Since we expect a scaling of the form (remember the ground state energy is always zero) The result of the extrapolation with BST-algorithm is −z = −1.50000009 with error 0.00000323. This clearly shows that the exponent z is 3 2 . It is interesting to note that this state corresponds to the Bethe ansatz equations with only one type of roots. The vanishing of the second type of roots does not correspond to the simple inclusion of the one-species model into the two-species model given by making the density of the second type of particles zero but rather for the state we consider corresponds to equal densities of these in-equivalent particles. The fact that the densities are coupled may suggest that there is some underlying quasi-particle formalism which may give a theoretical way to explain the connection to the single species exclusion model and the occurrence of the exponent 3 2 well known from the totally asymmetric single species exclusion model.

Analytical solution of the Bethe ansatz equations
First we would like to rewrite the Bethe ansatz equations equations (12) and (13) in integral form. To that end we use the following two changes of variables: e λ k = z k , z 2 k = Z k , k = 1, . . . , p and e Λ j = y j , y 2 j = Y j , j = 1, . . . , r. Equation (12) becomes: , k = 1, . . . , p, The new equation for the energies reads: Our numerical work described above suggests that the first excited state lies in the sector with p = L 3 , r = 0. After applying the above transformations, the transformed roots Z k , k = 1, . . . , L 3 will lie on a curve that is shown for L = 36 in Fig. 3. In order to analyze the Bethe ansatz equations in the limit of large lattice lengths L, it is convenient to introduce a function g(z) = ln z z − 1 and a function K(z l , z) = ln z z l .
For both definitions the branch cut of the ln-function is taken along the negative real axis. The Bethe ansatz equations for p = L 3 , r = 0 can now be written as with a so-called counting function iY L (z) = g(z) + 1 L L/3 l=1 K(z l , z).
Each excited state of the Hamiltonian in the sector with p = L 3 , r = 0 will correspond to one set of integers {I j | j = 1, . . . , L 3 }. The situation is very similar to the one described in [38] for the one-species asymmetric exclusion model. Therefore we can adopt the technique to transform the discrete Bethe ansatz equations (14) to integral equations by using an identity that follows from the residue theorem: In this equation, C is the contour enclosing all the roots Z j . We can view this contour C as the union of two contours C 1 and C 2 as shown in Fig. 4.
Rewriting this integral by separating the contributions coming from two contours C 1 and C 2 , this becomes: The formula for the energy reads: where the function ǫ(z) is given by It should be possible to analyze these equations along the lines of de Gier and Essler [38] by an expansion in inverse powers of L; this will be left for a future publication.

Conclusion and outlook
We have shown that the dynamical critical exponent of the totally asymmetric exclusion model is 3 2 which in the single-species asymmetric diffusion model is the exponent for the KPZ universality class. This is a new result that cannot be deducted from the recent paper by Arita et al. [25].
The Bethe ansatz equations and their solutions are qualitatively very different from the case of a single species asymmetric exclusion model. It would be very interesting to derive this result analytically once the patterns of the solutions of the Bethe ansatz equations are fully understood. Building on the results reported here, the next step would be to solve the integral equations representing the Bethe ansatz equations in the continuum along the lines of de Gier and Essler [38]. Since the integration is performed over the curve formed by the roots in the thermodynamic limit, a qualitative understanding of this curve is a prerequisite for this calculation. We hope to report on this soon.
Another new direction will be to generalize the Bethe ansatz to the case of open boundaries.