Density and current profiles in $U_q(A^{(1)}_2)$ zero range process

The stochastic $R$ matrix for $U_q(A^{(1)}_n)$ introduced recently gives rise to an integrable zero range process of $n$ classes of particles in one dimension. For $n=2$ we investigate how finitely many first class particles fixed as defects influence the grand canonical ensemble of the second class particles. By using the matrix product stationary probabilities involving infinite products of $q$-bosons, exact formulas are derived for the local density and current of the second class particles in the large volume limit.


Introduction
Zero range processes (ZRPs) [19] are stochastic particle systems on lattice modeling various flows in granules, queuing networks, traffic and so forth. Their characteristic feature is that particles are allowed to share a site and hop over the lattice with the rate that only depends on the occupancy and the list of leaving particles at the departure site 1 . To describe their hydrodynamic limit and the rich behavior like condensation has been an important issue in non-equilibrium statistical mechanics. See for example [5,8,9] and references therein.
In the recent work [11], a stochastic R matrix for the quantum affine algebra U q (A (1) n ) was constructed. It gives rise to discrete and continuous time Markov processes associated with a commuting family of Markov transfer matrices. They are formulated as stochastic dynamics of n classes of particles in one dimension with zero range type interaction. Many integrable Markov processes studied earlier, e.g. [2,3,12,16,18,20,21] can be identified with their special cases as summarized in [10,Fig.1,2]. In this paper we will be concerned with the version of the model introduced in [11,Sec.3.3,3.4], which will be called the U q (A (1) n ) ZRP. When n = 1 it reduces to the zero range chipping model introduced in [16].
Stationary states of the U q (A (1) n ) ZRP were obtained in [13,14]. Let σ i = (σ i,1 , . . . , σ i,n ) ∈ Z n ≥0 be a local state, which means that there are σ i,k particles of class k at the lattice site i. For the length L periodic chain, the probability of finding the system in a given configuration (σ 1 , . . . , σ L ) ∈ (Z n ≥0 ) L is expressed, up to normalization, by the matrix product formula P(σ 1 , . . . , σ L ) = Tr(X σ1 · · · X σL ), (1.1) where X σi is an operator acting on the tensor product F ⊗ n 2 (n−1) of the q-boson Fock space F = ⊕ m≥0 C(q)|m . For n = 1, the X σi is just a scalar meaning that the stationary measure is factorized. However it is an exceptional feature limited to n = 1. For instance when n = 2, the operator X α1,α2 for the local state (α 1 , α 2 ) ∈ Z 2 ≥0 reads X α1,α2 = (µ; q) α1+α2 (q; q) α1 (q; q) α2 where µ is another model parameter and the symbol (z; q) m is the q-shifted factorial defined in the end of this section. The operators b, c and k are the q-boson creation, annihilation and the number operators acting on F as b|m = |m + 1 , c|m = (1 − q m )|m − 1 , k|m = q m |m .
For n general the operator X α1,...,αn for the local state (α 1 , . . . , α n ) ∈ Z n ≥0 possesses a nested structure with respect to the rank n [14]. Thus the U q (A (1) n ) ZRPs form the first systematic examples of multispecies (or multi-class) ZRPs whose stationary measure on the ring is not factorized. The relevant matrix product operators are also quite distinct from those in the exclusion type processes (cf. [4,6,15]) in that they involve quantum dilogarithm type infinite products of q-bosons, offering a challenge to extract physics of the model.
With this background in mind we present in this paper a modest analysis of stationary properties of the U q (A (1) 2 ) ZRP based on the matrix product formula (1.1)-(1.2). We introduce finitely many first class particles as defects and investigate their influence on the second class particles whose density is kept finite in the infinite volume limit. To motivate this setting, although somewhat technically, note that one must pick an equal number of b's and c's to get a nonzero contribution to the trace (1.1). Therefore the sum of the expansion index j in (1.2) coming from X σ1 , . . . , X σL in (1.1) should coincide with the total number of the first class particles. Gathering all such contributions in the limit L → ∞ is a feasible task at least if the first class particles are kept finite.
The second class particles will be treated in the grand canonical ensemble, which means that X α1,α2 is effectively replaced by the generating series with respect to α 2 in the fugacity y: See (4.1) and (3.3). We stay in the regime 0 < q, µ < 1 where there is no symptom of condensation. See the remarks around (5.9) concerning this point. The equivalence with the canonical ensemble treatment will be argued in Section 6.3. The basic quantity is the probability P (r, m) that exactly m second class particles are found at the site r under the condition that the sites 1, 2, . . . , s of the periodic lattice Z L of size L contain d 1 , . . . , d s first class particles. To avoid the ambiguity we assume d 1 , d s ≥ 1 but the choice d i = 0 is still allowed for 0 < i < s. So they form a cluster of defects of size s in general. Up to normalization the conditional probability P (r, m) is given by replacing the r th operator from the left in Tr A d1 · · · A ds A L−s 0 by y m X dr,m if 1 ≤ r ≤ s and by y m X 0,m elsewhere 2 . Our task is to evaluate it in the infinite volume limit L → ∞ with s and d 1 , . . . , d s kept fixed. The limit separates the periodic lattice Z L into the three distinct regions I, II and III, which are the inside (1 ≤ r ≤ s), the right (r > s), and the left (r ≤ 0) of the defect cluster, respectively. The Z L periodicity of the lattice implies that the probability P (r, m) for the right region r > s and the left region r ≤ 0 should match when r → ∞. This is indeed the case as seen from (6.15) and (6.19).
Once the conditional probability is determined, the local density and current of the second class particles are derived at any site r. They are physical quantities seen from the defects. The necessary calculations are elementary. The final results are summarized in Theorem 4, 7 and 9 for the regions I, II and III, respectively. They are expressed in terms of the q-digamma function and its derivative together with the functions G m,l (d 1 , . . . , d s ) (7.10) which incorporate the effect of defects. Curiously the latters are related to the monodromy matrices of the U q (A (1) 1 ) ZRP containing the fugacity y as a spectral parameter. In Proposition 10 we will also show that the defects decrease the second class particles in the entire system exactly by their number d 1 + · · · + d s compared from the defect-free situation.
We present the profiles of the local density and currents in a number of figures for various values of q, µ, ρ and the defect pattern (d 1 , . . . , d s ), where ρ denotes the average density of the second class particles. The density profiles exhibit a peak and a valley at the left and the right boundaries of the defect cluster, respectively. In the current profiles the peak is not observed but other behavior is more or less similar to the density. The detail is dependent on the pattern (d 1 , . . . , d s ) and it is not easy to provide an intuitive explanation in general. However our result captures the mode η j (5.3) controlling the correlation length, which shows that the influence of the defects reaches longer distance when the 2 Precise treatment involves a regularization as mentioned after (4.6).
average density ρ is lower. Moreover we provide especially simple profiles in the limits ρ → 0 and ρ → ∞ in the case of homogeneous defects in Figure 4, which makes the general case easy to infer.
We remark that analyses of density and current based on the grand canonical ensemble similar to the present paper have been done for a class of two species exclusion type processes. See for example [4,17] and references therein.
The layout of the paper is as follows. We recall the U q (A (1) n ) ZRP [11] in Section 2 and the matrix product formula for the stationary probabilities for n = 2 [13] in Section 3. The local density and current of the second class particles in the grand canonical ensemble are formulated in Section 4. We first deal with the defect-free single species case in Section 5 as a preparatory warm-up partly reproducing known facts in earlier works, e.g. [2,5,16]. Our main results in the presence of defects are given in Section 6 with a number of figures showing the density and current profiles. Their derivation are detailed in Section 7. Section 8 contains a brief summary and discussion. Technical lemmas are collected in Appendix A.
Henceforth we call the T (λ|µ 1 , . . . , µ L ) Markov transfer matrix assuming 0 < µ ǫ i < λ ǫ < 1, q ǫ < 1 always. The choice of ǫ = ±1 specifies one of the two physical regimes of the system. The evolution equation (2.11) describes a stochastic dynamics of n classes of particles hopping to the right periodically via an extra lane (horizontal arrows in (2.8)) which particles get on or get off when they leave or arrive at a site. The rate of such local processes is specified by (2.2), (2.3) and (2.4). For n = 1 and the homogeneous choice µ 1 = · · · = µ L , it reduces to the model introduced in [16].
The Markov processes (2.12) is naturally interpreted as the stochastic dynamics of n classes of particles on the ring of length L. The base vector |α 1 ⊗ · · · ⊗ |α L with α i = (α i,1 . . . , α i,n ) ∈ Z n ≥0 represents a state in which there are α i,a class a particles at the i th site. There is no constraint on the number of particles that occupy a site. The matrices H + and H − describe their stochastic hopping to the right and the left nearest neighbor sites, respectively. The transition rate can be read off the first terms on the RHS of (2.14) and (2.15), where the array γ = (γ 1 , . . . , γ n ) specifies the numbers of particles that are jumping out. The superposition H(a, b, ǫ, q, µ) corresponds to a mixture of such right and left moving dynamics. The rate is determined from the original occupancy (α for h + and β for h − ) and the list of leaving particles (γ for h ± ) at the departure site and it is independent of the status of the destination site. Thus it defines a ZRP of n classes of particles in a slightly generalized sense in that the rate is allowed to depend on γ. Here is a snapshot of the system for the n = 2 case.

Stationary states
3.1. Definition and example. By definition a stationary state of the discrete time U q (A (1) n ) ZRP (2.11) is a vector |P ∈ W ⊗L such that The stationary state is unique within each sector m up to normalization. Apart from m, it depends on q and the inhomogeneity parameters µ 1 , . . . , µ L but not on λ thanks to the commutativity (2.10). Sectors m = (m 1 , . . . , m n ) such that ∀m a ≥ 1 are called basic. Non-basic sectors are equivalent to a basic sector of some n ′ < n models with a suitable relabeling of the classes. Henceforth we concentrate on the basic sectors. The coefficient appearing in the expansion |P = is the stationary probability if it is normalized as (σ1,...,σL)∈S(m) P(σ 1 , . . . , σ L ) = 1. In this paper we will abuse the terminology also for the unnormalized states and probabilities. The stationary states for the continuous Markov process (2.12) are those |P that satisfy H|P = 0. They are obtained from the discrete time ones just by the specialization to the homogeneous case µ 1 = · · · = µ L = µ. This is because (2.12) is an infinitesimal version of the commuting time evolutions (2.11) by the construction. In particular the stationary states are independent of a and b in (2.12). Example 1. Set (n, L) = (2, 3) and consider the homogeneous case µ 1 = µ 2 = µ 3 = −κ. The stationary state in the sector m = (1, 2) is given by where |∅, 2, 12 for example is the multiset representation of the base vector |(0, 0) ⊗ |(0, 1) ⊗ |(1, 1) in the multiplicity representation. The terms "cyclic" are those obtained by the shift |σ 1 , Similarly the stationary state in the sector m = (2, 1) is given by It has been conjectured [14,Ex.4] that for any sector m there is a normalization such that P(

3.2.
Matrix product formula. In [14] a matrix product formula for the stationary probability was obtained for general n and inhomogeneity µ 1 , . . . , µ L . In the rest of the paper we shall exclusively deal with the n = 2 case [13] with the homogeneous choice µ 1 = · · · = µ L = µ in the regime 0 < q, µ < 1. In the continuous time setting, it corresponds to ǫ = +1 in (2.13)-(2.15).
To describe the result we need the q-boson algebra B generated by b, c, k obeying the relations Let F = m≥0 C(q)|m be the Fock space and F * = m≥0 C(q) m| be its dual on which the q-boson operators b, c, k act as where |−1 = −1| = 0. They satisfy the defining relations (3.3). The bilinear pairing of F * and F is specified as m|m ′ = δ m,m ′ (q) m . Then m|(Q|m ′ ) = ( m|Q)|m ′ is valid and the trace is given by Tr Q = m≥0 m|Q|m (q)m . As a vector space, the q-boson algebra B has the direct sum decomposition The trace Tr Q is finite and nonzero only if We introduce the operator depending on µ, q which acts on (a completion of) F * and F : We have written X (α1,α2) as X α1,α2 for simplicity. The ratio of the infinite product of operators here is to be understood via the series expansion 2 ) ZRP (2.12) or the discrete time U q (A 2 ) ZRP (2.11) with homogeneous parameters µ 1 = · · · = µ L = µ is expressed in the matrix product form This result was obtained in [13, eq.(42)]. The formula (3.6) depends on α 1 and α 2 quite differently, indicating distinct features between the two classes of particles. Interestingly this is a reflection of a tiny asymmetry of the hopping rate (2.16) and (2.17) under the interchange (α 1 , β 1 , γ 1 ) ↔ (α 2 , β 2 , γ 2 ).
Example 3. Consider the first two terms in (3.2).
Setting µ = −κ and substituting 4 , we reproduce the first two terms in (3.2) after removing the common factor (q) 1

Density and currents of second class particles:
Formulation of the problem Now we come to the main theme of the paper, the stationary quantities in a grand canonical ensemble with respect to the second class particles in the infinite volume limit. The first class particles are kept finite and regarded as defects. In this section we give a general formulation of the problem only deferring the results and their derivation to subsequent sections. We fix the parameters 0 < q, µ < 1 and consider the ǫ = +1 case of (2.13)-(2.15) for the continuous time ZRP (2.12).
Let us introduce the generating series of X m,n (3.6) with respect to the second class particles with fugacity y: The quantity g m introduced here will be used very frequently in the sequel. We consider the conditional probability in the stationary states supposing that there are d i first class particles at site i for i = 1, . . . , s, and no first class particle is present elsewhere. So they form a fixed cluster of defects whose spatial extension is s and the total number of defect particles is d 1 + · · · + d s . In terms of the site variable σ i = (σ i,1 , σ i,2 ), the condition is expressed as allowing σ i,2 still to fluctuate everywhere. We assume that d 1 ≥ 1 and d s ≥ 1 to fix the location of the defect cluster but allow the choice d i = 0 for 1 < i < s. The basic quantity is the conditional probability that a given site r contains exactly n second class particles. In the customary notation for the probability P (A|B) of the event A under the condition B, we set There are three distinguish regions in the infinite volume limit L → ∞ as where r ≤ 0 should be understood as the site r + L ∈ Z L in the limit L → ∞. They correspond to the inside, the right and the left side of the defect cluster, respectively. The quantity (4.2) will be denoted by P I (r, n), P II (r, n) and P III (r, n) accordingly.
From the matrix product formula (3.8) and the definition (4.1), we have where (· · · ) ′ is defined after (3.5). In view of the matrix product operators (3.6) and (4.1), the prime restricts the ensemble to those sectors containing at least one second class particle for which the traces become finite. The formulas (4.5)-(4.6) fix the relative weight of such sectors, thereby specify what is meant by the "grand canonical ensemble" with respect to the second class particles with fugacity y. Obviously for any r the normalization n≥0 P (r, n) = 1 should be fulfilled in each region. Once the probability P (r, n) is obtained, one can evaluate various physical quantities. In this paper we investigate the expectation number and currents of the second class particles at site r. The former is defined by According to the regions, it will be denoted by ρ I (r), ρ II (r) and ρ III (r). We call them density for simplicity. As for the current, there are two components J(r) + and J(r) − associated with the local Markov matrices h + (2.14) and h − (2.15) with ǫ = +1 respectively: (4.8) They sum up the contributions from the l hopping second class particles out of total n+θ(1 ≤ r ≤ s)d r occupants weighted by the probability P (r, n) and the hopping rate w ± in (2.16) and (2.17). From (2.12) and the picture before (2.16), the total current from the site r to r + 1 is obtained by the superposition Thus it suffices to investigate J(r) + and J(r) − separately. We will also write them as J I (r) ± , J II (r) ± and J III (r) ± according to the regions.

Defect-free case
First we illustrate the analysis on the defect-free case s = 0 as a warm-up. It corresponds to the single species model, and some of the contents are well known by earlier works, e.g. [2,5,16].
In the absence of defects, the system acquires the Z L translational symmetry. As the result, the probabilities (4.5) and (4.6) are equal and independent of r. So we simply denote it by P (n). It is calculated as , where the prime is defined under (3.5) and we have set These quantities will be utilized frequently in the subsequent calculations. We will abbreviate η m (q k y) to η m if and only if k = 0. We suppose 0 < y < 1 whose consistency will be confirmed shortly. Then holds due to 0 < q, µ < 1. In Appendix A we show that the limit lim L→∞ and the infinite sum m≥0 in (5.1) may be interchanged, i.e., Since this is just m≥0 q mn δ m,0 = 1, we obtain the probability The correct normalization n≥0 P (n) = 1 has been achieved by virtue of (3.7). This formally agrees with the probability [2, eq.(53)] upon identification of the parameter α there with the fugacity y here. Next we relate the density ρ of the second class particles to the fugacity y. In the grand canonical ensemble under consideration, it is evaluated as As mentioned before (5.5), the "partition function" Tr(A L 0 ) ′ here may be replaced with Λ(y) L as L goes to infinity, leading to The f (ζ) is a version of the q-digamma function. It monotonously grows from 0 to ∞ as ζ changes from 0 to 1 behaving as f (ζ) ≃ ζ 1−q (ζ ց 0) and f (ζ) ≃ 1 1−ζ (ζ ր 1). Thus the fugacity and density are related asymptotically as The difference f (y) − f (µy) in (5.7) is similarly increasing monotonously from 0 to ∞ for y ∈ (0, 1). Thus ρ ∈ (0, ∞) is in one-to-one correspondence with y ∈ (0, 1), which confirms the consistency of the assumption made before (5.4). It implies that there is no phase transition typically recognized as condensation (cf. [5]) in ZRPs. To summarize, (5.7) determines the relation y = y(ρ) and ρ = ρ(y).
The density ρ(r) (4.7) is also independent of r. From (5.6) it is evaluated as in terms of the average density ρ in (5.7) confirming the consistency. Now we are ready to evaluate the currents of the second class particles (4.8). Since they are independent of the site r, we simply write it as For instance J − is computed from (2.17) and (5.6) as where the last equality is due to (3.7). A similar calculation for J + leads to current-density relation, a basic characteristic of the system, given as via y = y(ρ) (5.7). Here the function h(ζ) is the derivative of the q-digamma function f (ζ) in (5.8): (5.12) It obviously satisfies the relation The result (5.10) is restated as n≥l≥1 lw − ((0, l)|(0, n))P (n) = h(y), (5.14) which is µ-independent despite that w − ((0, l)|(0, n)) (2.17) and P (n) (5.6) depend on µ individually. An expression similar to the total current aJ + − bJ − (4.9) with J ± given by (5.11) has also been obtained in [2, Sec.4.1].
The function h(ζ) is monotonously increasing and tends to ∞ as ζ approaches 1 from below. Consequently the both currents in (5.11) grow monotonously with the density ρ via (5.7). Their leading asymptotic behavior is given by    We close the section with the description on the limiting cases q, µ → 0, 1. We shall only give the leading terms.

Density and currents: Results in general case
Let us proceed to the general case in which d 1 , . . . , d s defect (first class) particles are present at the sites 1, . . . , s. We present the final results on the conditional probability P (r, n) (4.2), the local density ρ(r) (4.7) and the currents J(r) ± (4.8) in the regions I, II and III in Theorem 4, 7 and 9. Some of them look bit messy but the point is that they are always finite sums of appropriate building blocks, which allow accurate numerical evaluations. Another point is that they are expressed in the form that elucidates the difference from the defect-free case s = 0. We continue to stay in the range 0 < q, µ, y < 1 and use the following quantities that have already appeared in the previous section: ρ: average density of the second class particles in the entire system, y ∈ (0, 1): fugacity of the second class particles determined from ρ via (5.7), P (n): probability in the defect-free case (5.6), J ± : currents in the defect-free case (5.11), h(ζ): derivative of q-digamma function describing the currents (5.12), η j : quantity controlling the decay of correlations (5.3).
In addition to them the following functions, detailed in Section 7, will be the basic ingredients: φ(l|m): constituent of G m,l (d 1 , . . . , d s ) (7.4), G m,l (d 1 , . . . , d s ): building block incorporating the effect of defects (7.9).
The relation (5.7) between the fugacity y and the average density ρ of the second class particles was originally derived for the defect-free case. We will justify its use in the presence of defects in Proposition 10 and the comments following it.
6.1. Main results. First we consider the region inside the defect cluster.
Theorem 4. In the region I (1 ≤ r ≤ s) the following formulas are valid: 3)
The proof will be given in Section 7.2. By means of (3.7) one finds that the total probability is expressed as n≥0 P I (r, n) = m G 0,m (d 1 , . . . , d r ). This indeed gives 1 thanks to (7.14). In Section 7.5 we will show that (6.2) can also be expressed as where K(0) = 0 and the sum over m is taken in the same way as in Theorem 4. The difference structure K(r) − K(r − 1) in (6.5) matches the sum rule in Proposition 10.
Let us remark on the low and high density asymptotic behavior. They correspond to y ց 0 and y ր 1, respectively. See (5.9). From the properties of the unperturbed density and currents in (5.9), (5.15) and (5.16), one can derive the following behavior by utilizing (7.17) and (7.19). The result (6.7) with r = 1 agrees with the y → 0 case of Example 5 implied by (5.9).
Next we consider the right region II.
Theorem 7. In the region II (r > s) the following formulas are valid: 14) where the sum in (6.12) extends over i, j, m ∈ Z ≥0 such that 0 ≤ i ≤ j ≤ m and d s ≤ m ≤ d 1 +· · ·+d s . The sums in (6.13) and (6.14) are taken for 1 ≤ j ≤ m and d s ≤ m ≤ d 1 + · · · + d s .
Finally we consider the left region III.
Theorem 9. In the region III (r ≤ 0) the following formulas are valid: The proof will be given in Section 7.4. Let us introduce the quantity which represents the total excess of the number of the second class particles from the average value in the entire system. From (6.8) and (6.16) one can easily check lim ρ→∞ ∆ρ tot = −(d 1 + · · · + d s ). Our next result shows that this holds in general.
Proposition 10. The total excess of the second class particles is given by The proof will be given in Section 7.5. Since ∆ρ tot is finite, the parameter ρ indeed remains to be the average density of the second class particles in the infinite volume limit. It is curious that ∆ρ tot exactly cancels the total excess +(d 1 + · · · + d s ) coming from the first class particles. To find a simple explanation of this fact is an interesting open problem.
The results in Theorem 4 (resp. Theorem 9) are independent of d r+1 , . . . , d s (resp. d 1 , . . . , d s ). They imply that the defects only influence their right in the large volume limit despite that the component H − in the Markov matrix (2.12) represents the left moving particles. We do not have an intuitive explanation of this fact. One should however remember that the both H ± in (2.13) originate in the discrete time evolution T (λ|µ, . . . , µ) described by the right moving particles only as in (2.8). Moreover the derivative in (2.13) gives rise to the denominator T (λ|µ, . . . , µ) −1 which is 1 for H + whereas it is the left cyclic shift for H − due to S(µ, µ)(|α ⊗ |β ) = |β ⊗ |α . Thus we could have had the right moving dynamics only, if not so neatly described, just by taking the derivative rather than the logarithmic derivative at λ = µ. In other words, the curiosity is attributed to the commutativity [H + , H − ] = 0, a basic consequence of the integrability of the model, which implies that the left and the right moving dynamics should possess the same stationary states on the ring of finite size.

Case of homogeneous defect.
To grasp the results in Theorem 4 and 7 qualitatively, it is helpful to consider the homogeneous case d 1 = · · · = d s = d and the limits ρ → 0 and ρ → ∞. The behavior for general ρ can then be inferred as something in between. In this setting we still have the integer parameters d, s ∈ Z ≥1 and the continuous parameters 0 < q, µ < 1. Set They satisfy D > d and ν > 1 obviously. The results on the local density (6.7), (6.8), (6.16) and Schematic plots of them look as Figure 4 below. Figure 4. Density profiles under the homogeneous defect d 1 = · · · = d s = d in the two limits ρ → 0 and ρ → ∞ according to (6.21).
Actually in the limit ρ → 0, either ρ(s)/ρ ≷ 1 can happen according to νq (s−1)d ≷ 1. In the other limit ρ → ∞, the density completely resumes the average value for r ≥ s + 2. For general 0 < ρ < ∞, the local density ρ(r) gets larger than ρ at the left boundary r = 1 of the defect cluster forming a peak. It then decreases until the site r = s + 1 forming a valley at the left boundary of the region II. Then in the region II, it recovers toward the average value ρ exponentially as mentioned after (6.15). Such behavior and sensitivity to the defects are sharper in higher density case.
As for the local currents we are to compare J(r) + and J(r + 1) − in view of (4.9). The results given in (6.9) -(6.11) and (6.17) - (6.19) are summarized as Thus we see that J + (r) and J − (r+1) have similar asymptotics aside from the overall normalization for large ρ. As mentioned under (5.16) about J ± , the former is convergent whereas the latter is divergent as ρ tends to infinity. They do not exhibit a peak at r = 1 like the density ρ(r), but other behavior is more or less similar. In particular when ρ is sufficiently large, the both currents get smaller inside the defect cluster than the defect-free case J ± .
6.3. Comparison with numerical evaluation in canonical ensemble. One can evaluate the density numerically by the matrix product formula (3.8) in a fixed sector with system size L. It formally corresponds to the canonical ensemble average which will be denoted by ρ c L (r). In the actual calculation of ρ c L (r), we have truncated the operators (3.6) and (4.1) to the matrices acting on the finite dimensional subspace of F of the form d1+···+ds+t m=0 R|m and checked the convergence has been achieved sufficiently already for t = 1, 2, etc. Figure 5 compares ρ(r) and ρ c L (r) for 1 ≤ r ≤ L. As remarked after (6.15), the correlation length of the system is smaller for larger average density ρ. Thus the agreement of ρ(r) and ρ c L (r) is expected to be better for larger ρ, therefore is harder to observe in numerical calculations. Admittedly the agreement in Figure 5 is not quite excellent, but always exhibits the tendency to improve as L gets large. In general fluctuations in the density is of order L − 1 2 . So the grand canonical approach should coincide with the canonical one in the large volume limit for there is no symptom of phase transition in the range 0 < q, µ < 1. See the remarks following (5.9).

Density and current profiles.
Here we present the result of numerical evaluation of the formulas in Theorem 4, 7 and 9 in a number of figures. We first consider the profile of local density ρ(r) (6.2), (6.13) and (6.19) in the presence of various defects in Figure 6  Examples of the density profile for a finite range of ρ are given in Figure 7. To display the defect region I in the center, the lattice coordinate r has been shifted. A similar convention will be employed in the subsequent figures in this subsection except Figure 9. For a fixed average density ρ, dependence of ρ(r) on q and µ are shown in Figure 8. One sees that the inhomogeneous defects make the density profile in the region I irregular, but they still tend to produce a peak and a valley at their boundaries as observed in the homogeneous case in Figure 4.

Derivation of main results
In this section we derive the main results given in Theorem 4, 7, 9 and Proposition 10. 7.1. Preliminary. Our first step is to reduce the trace Tr(· · · ) ′ over the Fock space to the "vacuum expectation value" 0|(· · · )|0 in the infinite volume limit.
Proposition 11. The traces in the probabilities (4.4)-(4.6) are reduced to the following: This is a corollary of Lemma 18 in Appendix A. Let us prepare the functions that will serve as building blocks to express (7.1) and (7.2). For l, m ∈ Z ≥0 set where Φ q was defined in (2.3). Although the dependence on µ and y is suppressed in this notation, it deserves attention that the fugacity y plays the role of a spectral parameter here. From (2.6) we know l≥0 φ(l|m) = 1, (7.5) where the summand is nonzero only for 0 ≤ l ≤ m.

Lemma 12.
For any m ∈ Z ≥0 the following equality is valid:

Proof. Explicitly it reads as
where we have replaced j by m − j. We prove (7.6) by induction on m. The case m = 0 is obvious.
In what follows we assume (7.6) is valid when m is replaced by 0, 1, . . . , m − 1. Upon multiplication by (µy) m , the both sides of (7.6) become polynomials in µ of order m − 1. Thus it suffices to check the equality at m points, say µ = q −r with r = 0, 1, . . . , m − 1. So we set µ = q −r in (7.6) and further replace y by q r y to simplify the formula slightly. The result reads where A = m−1 i=r 1 1−yq i and the upper bound of the j sum has been reduced from m − 1 to r owing to the factor (q −r ) j . The coefficient of A is 1 due to (7.5) with (µ, y) replaced by (q −r , q r y). Let us subtract A from (7.7). The upper bound of the k sum in the RHS becomes r − 1. The second sum in the LHS restricts the upper bound of j to r − 1. Further applying in the process, we find that the result is equivalent to This coincides with (7.6) with (m, µ, y) replaced by (r, q −m , q m y). Since r ≤ m − 1, its validity is assured by the induction hypothesis.
We note that the limit y → 1 in Lemma 12 leads to the identity The following function plays the basic role in our working.

Proof. Substitute the expansion
where the sum is over l 1 , . . . , l s ∈ Z ≥0 and the constraint is imposed to pick the non-vanishing bracket.
the claimed formula follows.
Remark 15. By Lemma 18, the result (7.20) with m = 0 is equivalent to The factor g d1 · · · g ds here is proportional to the stationary probability of the configuration (d 1 , . . . , d s ) of the single species model [16] of the first class particles on the length s ring. In this sense the effect of the second class particles in the grand canonical picture is "renormalized" into the extra factor y −1 for the first class particles' fugacity, reducing the U q (A 2 ) ZRP effectively to the U q (A 1 ) ZRP. We expect a similar recursive feature in the U q (A (1) n ) ZRP with general n [11], which may be viewed as a reminiscent of the nested Bethe ansatz. The result (7.22) is valid including d i = 0. It implies that separated clusters of the first class particles do not attract nor repel depending on their distance under the grand canonical background of the second class particles. They feel each other only when they merge together at the same site or split from a common site.
The following result reduces the multiple sum in G m,l (d 1 , . . . , d s ) (7.9) into a single sum when d 1 = · · · = d s = 0 and elucidates the s-dependence explicitly. It will be utilized in the region II to extract the large distance behavior of the density and currents away from the defects.
Lemma 16. For 0 ≤ i ≤ m, the following formula is valid: Proof. Introduce the function where the sum is over l 1 , . . . , l r ∈ Z ≥0 under the specified condition. It is related to the LHS by where F m−i,r (q i y) is given by (7.23) with η j = η j (y) (5.3) replaced by η j (q i y). The relation (7.24) is easily checked by means of (7.16) and η i+j = η i+j (y) = η i η j (q i y). Replacing y, m and the summation variable j with q −i y, m + i, j + i and using (y) m−l = (−q m y) −l q l(l+1)/2 (y) m (q 1−m y −1 ) l (7.25) with y = q, the formula to show becomes We prove this by induction on r. When r = 1, the RHS is expressed in terms of the q-hypergeometric q −n , a c ; q, q = a n (c/a) n (c) n , (7.27) we obtain (µ)m (q)m , which agrees with (7.23) with r = 1. We now prove the r + 1 case assuming the r case holds. From (7.23) we have Using the induction hypothesis and (7.25), exchanging the orders of the summations with respect to l and j, it can be written as Using (7.27) and (q 1−m y −1 ) m = (−y) −m q −m(m−1)/2 (y) m , we arrive at the expression (7.26) with r replaced by r + 1, which finishes the proof.
Lemma 16 and (7.24) with i = 0 yield the formula We remark that by the definition (7.9) and (7.4) the LHS in Lemma 16 is regular at q = 0 and manifestly positive in the range 0 < q, µ, y < 1 under consideration, whereas such features are highly nontrivial in the RHS. In particular, individual terms therein can become O(q − 1 2 m 2 ) which needs care in numerical evaluation for small q.

7.2.
Proof of Theorem 4. Theorem 4 is concerned with the region I (4.3). Let us derive (6.1) from (7.1). Substituting (3.6) into it and applying (7.20) to the denominator we get |i i| (q)i at the two places so as to use By Lemma 13 and (7.20) the result reads where φ(l − m + d r |l) is defined in (7.4). Thus (6.1) follows by taking the l sum using (7.15) with (s, t) = (r, r − 1) and (7.11). Next we show (6.2) for the local density. After substituting (6.1) into (4.7), the sum is taken by n≥0 ny n q n(m−dr) (q dr µ) n (q) n where ρ is the average density (5.7). The resulting expression for ρ I (r) agrees with (6.2) except that ρ comes with the coefficient m G 0,m (d 1 , . . . , d r ). But this equals 1 by (7.14).
Substitution of (2.16), (2.17) and (6.1) leads to the calculation similar to (5.10). Its essential part is where h(y) is the derivative of the q-digamma function defined in (5.12). Now we have Expand h(q m µy) and h(q m−dr y) into the sum of h(y) and the remainder terms by (5.13). From (5.11) the contribution from the h(y) is expressed as m G 0,m (d 1 , . . . , d r )J ± = J ± by (7.14). The remainder terms give the RHS of (6.3) and (6.4).
where P (n) is the probability in the defect-free case (5.6) and the sum extends over 0 ≤ i ≤ d 1 +· · ·+d s . On the other hand the decomposition (7.15) and Lemma 16 lead to = (−1) j q mj−j(j−1)/2 (q −m )j (q)j into this, we get (6.12). The normalization n≥0 P II (r, n) = 1 is attributed to that of P I (r, n). In fact it can also be directly verified from (7.30) by noting n≥0 P (n)q ni = η i and (7.14).
To derive the density (6.13), we again set (d 1 , . . . , d r ) = (d 1 , . . . , d s , r−s 0, . . . , 0) in the region I result (6.2) and use (7.31) to get where the sum extends over i, j, m ∈ Z ≥0 such that 0 ≤ i ≤ j ≤ m and d s ≤ m ≤ d 1 + · · · + d s . We may restrict j to 1 ≤ j ≤ m. Then the sums over i and k are taken by applying the identity where the q-binomial theorem and (3.7) are used. After some simplification the result becomes (6.13).
To derive the current (6.14), we set (d 1 , . . . , d r ) = (d 1 , . . . , d s , r−s 0, . . . , 0) in the region I result (6.3), (6.4) and use (7.31) to get where the sum i,j,m is taken in the same way as (7.32). This time we apply the identity which follows from (7.33) by differentiation. The result yields (6.14). 7.4. Proof of Theorem 9. Theorem 9 is concerned with the region III (4.3). From the definition of X 0,n in (3.6) and 0| (µb)∞ (b)∞ = 0|k n = 0| by (3.4), the conditional probability (7.3) in the region III becomes = y n g n Λ(y) −1 due to (7.20). This is independent of r and coincides with the probability P (n) (5.6) for the defect-free case. Consequently the local density and current also reduce to the average density ρ and J ± in (5.11).
Thus Proposition 10 follows from K(1) = 0, K(r) =K(r + 1) (1 ≤ r < s), The first equality is obvious from (7.11) and the last one has already been derived in (7.34). So we are left to show K(r) =K(r + 1) only. By using (7.15) with (t, s) = (r, r + 1) and (7.11) we find which tells thatK(r + 1) is actually independent of d r+1 as with K(r). Now the proof is finished by Lemma 12.

Discussion
In this paper we have investigated the stationary properties of the U q (A 2 ) ZRP [11] on a ring whose stationary probability is not factorized but expressed in a nontrivial matrix product form. The second class particles are treated in the grand canonical ensemble subject to the influence of the d 1 , . . . , d s first class particles fixed as defects at the sites 1, . . . , s. The local density and current of the second class particles in the infinite volume limit are obtained in Theorem 4, 7 and 9. The density profile is shown to exhibit a peak and a valley at the boundaries of the defect cluster. The currents are locally suppressed by the defects. The detail is dependent on the inhomogeneity of the data (d 1 , . . . , d s ). Outside the defect cluster, its influence reaches longer distance when the average density of the second class particles is lower. It was shown that the second class particles are decreased by the number of the defect particles in Proposition 10.
The effect of the defects is expressed through the function G l,m (d 1 , . . . , d s ), which is essentially a column monodromy matrix of the U q (A (1) 1 ) ZRP containing the fugacity y as the spectral parameter. Compare the diagrams (7.10) and (2.8). We regard it as another reminiscent of the nested Bethe ansatz structure in addition to Remark 15. It is natural to generalize the analysis in this paper to the grand canonical ensemble for the both classes of particles simultaneously. Denoting the fugacity of the first class particles by x, one is led to the analogue of the grand canonical partition function in two variables as where we have used the cyclicity of the trace and the commutativity [Γ(x, y), Γ(µx, µy)] = 0 which follows from (3.6) by applying (µ) m+n = (µ) m (q m µ) n = (µ) n (q n µ) m . See also [14,Rem.3]. We leave the large L asymptotics of Z L (x, y) and related issues for a future study.
Appendix A. Proof of (5.5) and Proposition 11 We invoke the interchangeability of limits and infinite sums in elementary calculus adapted in the following form.
From (A.4), there exists l 0 ≥ 0 such that Γ l > 0 holds for all l ≥ l 0 3 . Then the above estimation leads to l≥0 V L,l = C + l≥l0 Γ l v L,l ≤ C + Lη ∞ l1+···+lL=λ g l1 · · · g lL j≥1 f j g j y j , where C = 0≤l<l0 V L,l is a finite constant and we have set f j = l≥l0 Γ l q jl . This series f j is convergent because of lim l→∞ Γ l = 1 as noted in (a). The series j≥1 f j g j y j in y is also convergent due to 0 < y < 1 and lim j→∞ fj gj fj+1gj+1 = q −l0 lim j→∞ 1−q j+1 1−q j µ = q −l0 > 1. Thus l≥0 V L,l is convergent for any fixed L. On the other hand (b) 1 tells that l≥0 V L,l is monotonously decreasing with respect to L for sufficiently large L. This proves (b) 2 .