Generalized hydrodynamics in complete box-ball system for $U_q(\widehat{sl}_n)$

We introduce the complete box-ball system (cBBS), which is an integrable cellular automaton on 1D lattice associated with the quantum group $U_q(\widehat{sl}_n)$. Compared with the conventional $(n-1)$-color BBS, it enjoys a remarkable simplification that scattering of solitons is totally diagonal. We also introduce a randomized cBBS and study its non-equilibrium behavior by thermodynamic Bethe ansatz and generalized hydrodynamics. Excellent agreement is demonstrated between theoretical predictions and numerical simulation on the density plateaux generated from domain wall initial conditions including their diffusive broadening.

The box-ball system (BBS) [1] and its generalizations are paradigm examples of integrable cellular automata on 1D lattice connected to Yang-Baxter solvable vertex models [2], Bethe ansatz [3], crystal base theory of quantum groups at q = 0 [4], ultradiscretization and tropical geometry etc. See for example the review [5] and the references therein. In this paper we introduce a new version of BBS associated with U q ( sl n ), which we call the complete box-ball system (cBBS). It possesses a number of distinguished features compared with the conventional BBS with (n−1) kinds of balls. We also introduce the randomized cBBS and study its non-equilibrium behavior by thermodynamic Bethe ansatz [6,7] and generalized hydrodynamics [8][9][10]. The results generalize the earlier work [11] for n = 2, where the cBBS itself reduces to the original one [1].
Let B k,l be the set of semistandard tableaux on the k × l rectangular Young diagram over the alphabet {1, 2, . . . , n} [12,13]. It is a labelling set of the basis of the irreducible sl n module corresponding to the mentioned Young diagram. By the definition only 1 ≤ k < n is relevant and the simplest B 1,1 is identified with the set {1, 2, . . . , n}. A typical or conventionally most studied BBS with (n − 1) kinds of balls [14] is a dynamical system on · · · ⊗ B 1,1 ⊗ B 1,1 ⊗ B 1,1 ⊗ · · · , (1.1) where ⊗ can just be regarded as the product of sets in this paper. One interprets it as an array of boxes which is empty for 1 ∈ B 1,1 or contains a color a ball for a = 2, 3, . . . , n ∈ B 1,1 .
All the distant boxes are assumed to be 1 = empty. By using the crystal theory of U q ( sl n ), one can formulate an integrable dynamics on such states yielding solitons [15,16]. 1 Here is an example with n = 4, where the t is the discrete time variable (the symbol ⊗ is omitted for simplicity): The cBBS we propose and study in this paper is obtained, among other things, by replacing (1.1) with Apparently it looks more involved than (1.1) since now each site has a larger internal structure B than B 1,1 for n ≥ 3. However, it turns out to enjoy much simpler and elegant features than the conventional BBS as we explain below. First, solitons in cBBS can be labeled just with color a ∈ {1, 2, . . . , n − 1} and amplitude i ∈ Z ≥1 as S j × S (a) i whose nontrivial effect is integrated into a phase shift on their asymptotic trajectories. This is a drastic simplification compared with the conventional BBS. For example in the notation explained in Sec. 2, scattering of S       The initial solitons S l (r ∈ {1, 2, . . . , n − 1}, l ∈ Z ≥1 ) naturally without introducing an artificial "barrier" which was necessary for the conventional BBS to prevent balls from escaping from the system for r ≥ 2. See a remark after [17, eq. (3)]. This is assured by the stability (2.23) which is another benefit of the choice of B in (1.2).
Third, the complete set of conserved quantities of cBBS are given as an (n − 1)-tuple of Young diagrams µ (1) , . . . , µ (n−1) , and they all admit most transparent meaning; µ (a) is nothing but the list of amplitude i of the color a solitons S (a) i . Such a direct interpretation of the conserved Young diagrams in terms of solitons was possible only for the first one µ (1) in the conventional BBS. Our cBBS puts all of µ (1) , . . . , µ (n−1) on an equal footing achieving the democracy of the conserved Young diagrams.
Given all these fascinating features which have escaped notice so far, we regard cBBS as the most natural as well as decent generalization of the original BBS [1] along U q ( sl n ) differing from the conventional ones for n ≥ 3. The nomenclature "complete" BBS is meant to indicate that the complete list B 1,1 , B 2,1 , . . . , B n−1,1 of (the labeling set of the basis of) the fundamental representations of sl n have been gathered into B (1.2) at each site. Put in the other way, the somewhat unsatisfactory aspects of the conventional higher rank BBS mentioned in the above may be attributed to the "incomplete" choice of B. We expect similar stories also in integrable cellular automata associated with the other quantum affine algebras.
After clarifying the nature of cBBS in Sec. 2, the rest of the paper is devoted to the study of a randomized version of cBBS by thermodynamic Bethe ansatz (TBA) and generalized hydrodynamics (GHD) extending the previous work on n = 2 case [11]. We focus on the i.i.d. (independent and identically distributed) randomness including (n−1) temperature generalized Gibbs ensemble (3.15). It corresponds to assigning (relative) fugacities z 1 , z 2 , . . . , z n to the letters 1, 2, . . . , n in the semistandard tableaux. We formulate the TBA and GHD equations on the so-called Y-function, the string/hole densities and the effective speed of solitons under any time evolution. They are fully solved in terms of the Schur functions with fugacity entries z 1 , z 2 , . . . , z n . These results provide a quantitative description of the cBBS soliton gas in a homogeneous system.
One of the central ideas in GHD is that the TBA Y-variable plays the role of normal mode in the Euler-scale hydrodynamics of an "integrable fluid". We apply it to the nonequilibrium dynamics of the cBBS soliton gas started from domain wall initial conditions. This is a typical setting in Riemann problem called partitioning protocol. See [10,18,19] and the references therein. As with the n = 2 case [11], density profile of solitons exhibits a rich plateaux structure depending on the type of solitons, time evolutions and the fugacity controlling the inhomogeneity of the system. We derive their position and height including the diffusive broadening by synthesizing all the ingredients in the preceding sections. They are shown to agree with extensive numerical simulations.
The outline of the paper is as follows. In Sec. 2 we introduce cBBS and explain its basic properties such as the commuting family of time evolutions, complete set of conserved quantities, solitons and their scattering rule, and the inverse scattering formalism. A key role is played by the Bethe ansatz structure which is realized as soliton/string correspondence. In Sec. 3 we proceed to the randomized version of cBBS. Explicit solutions are presented for the TBA equation (3.19) in (3.24), and string/hole densities in (3.52)-(3.53). They are expressed by the special combination (3.36) of the Schur function (3.6). These results can also be viewed as the solution to the variant of the limit shape problem of rigged configurations [17] adapted to cBBS. In Sec. 4 we apply GHD to the randomized cBBS. We present the speed equation of solitons for any time evolution (4.1) and its explicit solution in (4.7). The former coincides with [23, eq.(11.7)] for n = 2 and T (1) ∞ dynamics. The formulation of the GHD equations in matrix forms in Sec. 4.2 is useful to recognize the characteristic Bethe ansatz structure. In Sec. 5 we study the density plateaux generated from domain wall initial conditions by GHD. Each plateau is formed by particular subsets of solitons from two sides of the domain wall. These subsets, which we call soliton contents, show some intriguing patterns. Sec. 6 is devoted to summary and conclusions.
Appendix A recalls the algorithm for obtaining the image of the combinatorial R in the most general case. Appendix B includes the explicit piecewise linear formulas for the combinatorial R for n = 2 and 3. Appendix C is a brief exposition of the KSS bijection in the situation used in this paper.

Complete box-ball system
Throughout the paper we use the notation

Preliminaries
Consider the classical simple Lie algebra sl n (n ≥ 2). We denote its Cartan matrix by (C ab ) n−1 a,b=1 , where Let 1 , . . . , n be the fundamental weights and α 1 , . . . , α n be the simple roots. They are related by α a = n b=1 C ab b . Let sl n be the non-twisted affinization of sl n [28] and U q = U q ( sl n ) be the quantum affine algebra (without derivation operator) [29,30]. There is a family of irreducible finite-dimensional representations {W Kirillov-Reshetikhin (KR) module 2 named after the related work on the Yangian [31]. As a representation of U q (sl n ), W (k) l is isomorphic to the type 1 irreducible highest weight module V (l k ) with highest weight l k . W (k) l is known to have a crystal base B k,l [4,32]. Roughly speaking, it is a set of basis vectors of the W (k) l at q = 0. Practically in this paper we only need its combinatorial definition as a set and the operation called combinatorial R. They are described in the next subsection.

Combinatorial R
For (k, l) ∈ I, let B k,l be the set of semistandard tableaux of k × l rectangular shape over the alphabet {1, 2, . . . , n}. Let b = (t ij ), where t ij is the entry at the i th row from the top and the j th column from the left. By the definition, t 1j < t 2j < · · · < t kj and t i1 ≤ t i2 ≤ · · · ≤ t il hold for any i ∈ [1, k] and j ∈ [1, l]. The array row For tableaux S and T , their product is defined in the following two ways which are known to be equivalent: Here ← denotes the row insertion and → does the column insertion [12].
Row insertion S ← x: 1. Start at the top row of S.
2. If x is larger than or equal to the rightmost number in the current row, add x to the right of the row, which is the end of the insertion.
3. Otherwise, replace the leftmost element y of the row such that x < y by x and go to step (1) starting at the next row, now with y to be inserted.
Column insertion x → S: 1. Start at the leftmost column of S.
2. If x is larger than the bottom number in the current column, add x to the bottom of the column, which is the end of the insertion.
3. Otherwise, replace the smallest element y of the column such that x ≤ y by x and go to step (1) starting at the next column, now with y to be inserted.
In the above example one has b · c = 1 1 2 2 5 2 2 3 3 4 4 5 , c · b = 1 1 2 2 3 2 2 4 5 3 5 4 . (2.6) Now we are ready to explain the combinatorial R and the local energy H, which are essential ingredients in our cBBS. These notions were introduced in the crystal base theory [33] as the proper analogue of the quantum R matrices at q = 0, motivated by the corner transfer matrix method [2]. The concrete description given below is due to [34]. The combinatorial R is the bijection Here and in what follows, ⊗ is used to mean just the product of sets. It should not be confused with the product · of tableaux. Note that the dependence on k, k , l, l has been suppressed in R. For example, from This is a particular case of R : 3 . The algorithm to findb andc satisfying the conditionb ·c = c · b is described in Appendix A following [35, p55]. The local energy H is the function defined by

10)
H(b ⊗ c) = number of boxes strictly below the max(k, k )-th row of the tableau c · b. (2.11) In our working example (2.3), max(k, k ) = max(2, 3) = 3, hence from (2.6) we have H(b⊗c) = 1. Let u k,l ∈ B k,l be the particular tableau whose entries in the i th row from the top are all i. It corresponds to the highest weight element in the representation of U q (sl n ). Some small examples are It is easy to see for any k, k , l, l . If R in (2.7) and H in (2.10) are denoted by R B k,l ⊗B k ,l and H B k,l ⊗B k ,l respectively, they satisfy by the definition. In spite of the simplification (2.15), the local energy H B k,l ⊗B k,l (b ⊗ c) still depends on b ⊗ c nontrivially . The relation in (2.16) is called the inversion relation, which implies that R(b ⊗ c) =c ⊗b is equivalent to R(c ⊗b) = b ⊗ c. In view of this we write these relations also as b ⊗ c c ⊗b.
The most important property of the combinatorial R is the Yang-Baxter equation [2]. It includes R as the classical part and H as the affine part. To unify them into a single equation, we introduce an infinite set Aff(B k,l ) = {b[α] | b ∈ B k,l , α ∈ Z} and extend R to the mapR defined byR It is a reminiscent of the spectral parameter in quantum R matrices. Note that the shift of the mode H(b ⊗ c) in (2.18) is equally presented as H(c ⊗b) since they coincide owing to (2.17). This guarantees thatR also satisfies the inversion relation. Thus we also write (2.18) putting the two sides on a more equal footing. Now the Yang-Baxter equation is presented as which is an equality of the maps Aff( Here is an example for (k, k , k , l, l , l ) = (1, 2, 3, 3, 2, 1), where the modes are indicated as the indices in the bottom right of the tableaux.
In general, we will write b 1 ⊗· · ·⊗b L b 1 ⊗· · ·⊗b L if the elements b 1 ⊗· · ·⊗b L and b 1 ⊗· · ·⊗b L are transformed to each other by successively applying the combinatorial R to the neighboring components as above.
As the examples shown so far indicate, one can always forget the modes without destroying the relations on the classical parts. For readers convenience, we include explicit formulas of R and H for n = 2, 3 cases in Appendix B.

Vacuum state and stability
Collecting all the l = 1 cases in (2.12) for a given n, we set For instance, vac looks as The special element vac is referred to as vacuum, which will be used in the background configuration in our complete BBS. The set B is the tensor product of (crystals of) the complete list of fundamental representations of U q (sl n ). Using the KSS bijection in Appendix C, one can show a stability when L gets sufficiently large for some m and v 1 , . . . , v m ∈ B. As we will discuss in detail in Sec. 2.4, in the equation above we may interpret b as the initial state of a 'carrier'. The system contains L boxes that are initially in the vaccum state. The equality describes the change of the state of the system after the carrier has passed across the system, from left to right. It indicates that for large enough L only a finite number (at most m) of boxes on the left side get modified, while the others remain in the vacuum state vac. The output state of the carrier is also the vacuum u k,l . Here is an example (k, l) = (1, 3) with n = 3: In the (right) distance, all the combinatorial R's tend to the situation (2.13).

Time evolutions and conserved quantities
Now we define the complete BBS (cBBS). It is a dynamical system on B ⊗L . An element b 1 ⊗ · · · ⊗ b L ∈ B ⊗L is called a state where each b i ∈ B is a local state at site i. Thus there are #B = n 1 n 2 · · · n n−1 local states which have an internal structure as in (2.21). We assume that L is sufficiently large and impose the boundary condition that the "distant" local states are all vac in (2.21). This is compatible with the dynamics that we are going to introduce below. For any (r, l) ∈ I in (2.1), we define the time evolution T (r) l : B ⊗L → B ⊗L as follows: The relation (2. Here the boundary condition implies b j = b j+1 = · · · = b L = vac for L L − j 1. Then, thanks to the stability (2.23), we always end up with u r,l in the right. In short, T (r) l is obtained by going from NW to SE in (2.26). Thanks to the inversion relation (2.16), one can reverse the procedure to define the inverse (T (r) The key to the above construction is B r,l attached to the horizontal arrow in the diagram (2.26). It induces the time evolution T (r) l via the interactions with the local states by the combinatorial R. This degrees of freedom is called carrier [36], and especially u r,l is referred to as the vacuum carrier.
A local state may be labeled with its deviation from vac. For instance when n = 3, the 9 local states in B = B 1,1 ⊗ B 2,1 can be displayed as where the entries common with vac are colored white. In general one may regard a local state as an arrangement of n − 1 kinds of balls in n − 1 columns each obeying the semistandard condition. Then the combinatorial R specifies the rule under which the balls are exchanged between the carrier and a box. The complete BBS is a nomenclature emphasizing that the complete list of the possible column length 1, 2, . . . , n − 1 have been built in each local state as in (2.21). See also a remark after (2.22). We will see that this will lead to a remarkable simplification of solitons and their scattering. The time evolution T (r) l is associated with the energy E (r) l : B ⊗L → Z ≥0 . Its quickest definition is to declare that the mode of the carrier in (2.26) is shifted as follows: To be more detailed, suppose the local states have the form b i = c (n−1)(i−1)+1 ⊗ · · · ⊗ c (n−1)i ∈ B = B 1,1 ⊗· · ·⊗B n−1,1 so that b 1 ⊗· · ·⊗b L = c 1 ⊗· · ·⊗c (n−1)L . The carriers u 1 , . . . , u (n−1)L ∈ B r,l in the intermediate steps are determined from the initial condition u 0 = u r,l and the recursion for some c j . From the defining behavior of the mode in (2.18), the energy E (r) l is expressed as the sum of local energies attached to (2.29) as (2.30) As j gets large, u j tends to u r,l by the stability (2.23) under the boundary condition. Then (2.14) assures that the sum (2.30) converges to a finite value which is independent of L if it gets large enough. The time evolution and the energy enjoy the properties on any state and for arbitrary (k, l), (k , l ) ∈ I. This is a simple consequence of the Yang-Baxter equation. In fact, write s = b 1 ⊗ · · · ⊗ b L for short and consider which is a consequence of (2.13) and (2.14). It tells that (2.33) is also equal to 1 · · · T (n−1) 1 is the translation to the right by one lattice unit. as in (2.50) and (C.10). In particular, if γ k denotes the maximal amplitude of color k solitons that will be defined in the next subsection, one has T γ k for all l ≥ γ k .

Solitons and their scattering
This subsection is the place where the things start to differ significantly from the conventional (non-complete) BBS. The claims can be proved by invoking the inverse scattering method explained in Sec. 2.6.
We observe the cBBS in terms of the deviation of the states from the background configuration vac ⊗L . The first question is to find the "collective modes" or "quasi particles" which are stable localized patterns having a constant speed under any time evolution when isolated from the other patterns different from vac. They deserve to be called solitons if the stability under multi-body collisions, a much more stringent postulate, is further obeyed. This turn out to be the case for the cBBS reflecting the existence of the n − 1 families of the conserved quantities E , (2.37) They are different from vac ∈ B in (2.22) by only one entry at the bottom of one tableau. General definition is this: Using s a as the building block we next introduce S (a) i = i s a ⊗ · · · ⊗ s a ∈ B ⊗i ((a, i) ∈ I) (2.39) and call it soliton of color a and length (or amplitude or size) i, or type (a, i) for short. They have the following properties.
Here time grows from the top to the bottom. The quantity ∆, we call it phase shift, is given by where (C ab ) 1≤a,b≤n−1 is the Cartan matrix (2.2). Note that ∆ can be either positive or negative or even zero.
(III) By applying time evolutions sufficiently many times, any state can be decomposed into isolated solitons. Such asymptotic states can be taken, for example, as   l . In addition, by using the local energies (2.30), one gets the local contribution to m (r) i , which is therefore a soliton density. This is how the soliton densities are computed in the simulations in Sec. 5.3. with l ≥ 3. We write s 2 in (2.37) for example as 1 12 133 to save the space.

Inverse scattering method
In Appendix C, a brief exposition is given on the KSS bijection Φ between the set of highest paths P + (B, λ) and rigged configurations RC(B, λ). States of our cBBS satisfying the boundary condition are examples of the former with B = B ⊗L . See the remark in the end of Appendix C. The array λ = (λ 1 , . . . , λ n ) specifies that the number of letter a in a state is λ a . On the other hand, a rigged configuration is a combinatorial object like (C.12) visualizing the collection of strings which are triplets (color, length, rigging) as in S 0 (C.11). The time evolution T (k) l of the BBS induces that of rigged configurations via the commutative diagram: A remarkable feature of this scheme is that it achieves the linearization [38]: where rigged configurations are represented as multisets as in (C.7). It implies that the partial data {(a i , j i )} forms the complete set of conserved quantities (action variables) and the riggings undergo a straight motion (angle variables). In short, rigged configurations are action-angle variables of our cBBS, and the maps Φ and Φ −1 are direct and inverse scattering transformations. The scheme (2.48) achieves the solution of the initial value problem of cBBS by the inverse scattering method as T 3 The formula (2.46) is obvious for asymptotic states. General case follows from it and the fact that m (a) i 's are conserved quantities. 4 The formula (C.14) is given for rigged configurations, but it is known to agree with (C.1) undedr the KSS bijection.
Take a state b 1 ⊗· · ·⊗b L of cBBS and let m (a) j be the number of color a length j strings in the The two are known to coincide: (2.50) Due to the invariance under the time evolution, the proof reduces to the doable case of asymptotic states where solitons are all far apart. See for example [39].
Recall that m is the number of color a length j strings. Therefore (2.50) implies the soliton/string correspondence [38] soliton = string.
It has played a key role in the Bethe ansatz study of the randomized BBS [11,17].
To illustrate, consider the state on the top line of Example 2.6 in the previous subsection: At this stage a nontrivial question arises as to how to incorporate the boundary condition that all the distant local states are vac (2.20) properly. Here we propose that it is done by restricting the state sum to the highest ones that are defined in Appendix C.1. Similar treatment has been done in [11,17,27], where more detailed justification was made especially in the first reference. Thus the GGE partition function of our concern is (3.1) By synthesizing various results it can be rewritten as follows: Here the vacancy p is not a coincidence. It is indeed a general property in Bethe ansatz integrable systems that such a shift appears in the relation between the hole and particle densities and is related to the TBA kernel (see also Remark 2.3).

I.I.D. randomness
From here we concentrate on the situation where local states b ∈ B = B 1,1 ⊗ · · · ⊗ B n−1,1 obey an i.i.d. randomness. More specifically we consider an i.i.d. measure P(b ∈ B) including the parameters z 1 > z 2 > · · · > z n > 0 as follows: where a i denotes the entry of the tableau c ∈ B k,1 in the i th box from the top. The formula (3.5) implies that the probability of occurrence of a ∈ [1, n] is proportional to z a . In this sense, z 1 , . . . , z n are tableau variables representing the fugacities of the letters 1, . . . , n. Their ambiguity due to the invariance of (3.5) under the change z a → uz a for a constant u will be fixed in the rightmost condition in (3.7) below. The denominator of (3.5) is the k th elementary symmetric function which is Q (k) 1 in the notation (3.10). For a partition ν = (ν 1 , ν 2 , . . . , ν n ), let s ν = s ν (z 1 , . . . , z n ) denote the Schur polynomial 5 [13]: Partitions are identified with Young diagrams as usual. Schur polynomials are irreducible characters of finite dimensional gl n modules which can be restricted to those for sl n modules.
Reflecting this fact we will use further sets of variables w a , x a related to z a as follows: The variables w a and x a are formal exponential of the fundamental weight a and the negated simple root −α a regarded as parameters. See Sec. 2.1. The last relation in (3.7) leads to s ν = sν whereν = (ν 1 − ν n , . . . , ν n−1 − ν n , 0). We use a special notation for rectangle ν's as It satisfies the Q-system [31,37] (Q (a) is the special (n − 1)-parameter reduction of the GGE in Sec. 3.1. Denote the weight of a state s ∈ B ⊗L in the sense of (C.1) by wt(s) = (λ 1 , . . . , λ n ). Then the formulas (3.4) and (3.5) mean that s is realized with the probability proportional to z λ 1 1 · · · z λn n . Since the KKS bijection is weight preserving (see (C.17)), the RHS of (C.1) and (C.14) may be identified. Thus the factor n a=1 z λa a is proportional to ∞ thanks to (2.50). Therefore the probability of the state s is proportional to Comparing this with (3.1) we see that the i.i.d. measure (3.4) corresponds to the special case of GGE in which the inverse temperatures β (a) l is zero except Thus the inverse temperatures (with index ∞) can be identified with the simple roots. We will nonetheless allow co-existence of the symbols β (a) ∞ and α a in what follows. The identification (3.13) will also be reconfirmed after (3.25). From (3.9), the regime of parameters we consider is also specified in the variables z a , α a , x a as In particular x a → 0 and x a → 1 correspond to the zero and infinite limit of the temperature ∞ , respectively.
Remark 3.2. The dilute limit where no soliton is allowed is given by z 1 z 2 · · · z n or equivalently x 1 , . . . , x n → 0. In fact, hold in this limit. One also has Q

Thermodynamic Bethe ansatz
In the large L limit, the dominant contribution in the sum (3.15) comes from those {m This is (−1/L) times logarithm of the summand in (3.15) to which Stirling's formula has been applied. The scaling (3.16) is consistent with the extensive property of the free energy, which has enabled us to remove the system size L as a common overall factor. We have introduced a temporary cut-off s for the length i of solitons. It will be taken to infinity properly later.
Accordingly the latter relation in (3.17) should be understood as ε The corresponding maximal value F of the free energy per site F [ρ] (3.18) is obtained by using (3.17), (3.19) and (3.20). The result reads The TBA equation (3.19) is equivalent to the constant Y-system with the boundary condition (y (a) The Y-system (3.22) follows from the Q-system (3.11) by the substitution (cf. [37, Prop. 14.1]) . (3.24) As for the boundary condition (3.23), the first relation is valid due to the convention Q (a) −1 = 0 mentioned after (3.11). On the other hand the second relation is translated into s ) = e a in the regime e α 1 , , . . . , e α n−1 > 1 under consideration. Thus the large s limit of (3.25) leads to e β (a) ∞ = n−1 b=1 (e b ) C ab = e αa by (3.9). This is consistent with (3.13). We also remark that lim s→∞ (Q in [17, eq.(52)].
Substituting the last expression in (3.24) into the free energy density (3.21), we get is the normalized character of the irreducible sl n module with highest weight i a .
From (3.17) one sees that ε (a) is the total density of color a solitons. Let where we have written β which agrees with the density of total number of "balls" 2. Similarly for n = 3, one has Q (1) (3.31) In this way one can relate the densities and the inverse temperatures.
To summarize so far, we have obtained the solution {y 1 · · · Q (n−1) 1 . (3.33) The equation (3.32) is just the first relation in (3.17), where the first term is a disguised form of σ (a) i . The last relation in (3.33) is the postulate from the equation of state (3.28). At this point we invoke [17,Th. 5.1]. In the setting of this paper, it states that for any (r, l) ∈ I, the unique solution h is given by where the outer sum in (3.36) runs over those Young diagrams ν labeling the irreducible gl n modules appearing in the decomposition of the tensor product of those corresponding to the rectangles (i a ) and (l r ). We note that the irreducible decomposition of two rectangles is multiplicity free.
Example 3.5. For small n, (3.46) reads as given either by ρ i ) due to (3.17) and (3.20). It is independent of T (r) l and the coincidence of the two expressions is assured by (3.32). For n = 2 and T i (r, l). (4.9) Example 4.1. Consider n = 2 case. The variables in (3.7)-(3.9) are related as for i ≤ l. In view of the symmetry ν (3.36)), general case is obtained by replacing i with j := min(i, l) and l with m := max(i, l) in this formula. In particular, holds for i ≥ 1. Now the effective speed (4.7) is evaluated as .

Matrix form
It is convenient to formulate the results in the previous subsection in a matrix form whose indices range over I in (2.1). We introduce the matrices where the component-wise (or Hadamard) product * is commutative and will be used also for other vectors. For example we have ρ = y * σ. 7 In view of (4.1), the vector κ(r, l) is the bare speed of solitons under the time evolution T        Similarly the current (4.33) is expressed as a quadratic form:

Dynamics in an inhomogeneous system
Let us study an inhomogeneous system, with the hypothesis that it can be locally described by using the above formalism, but with the densities that have acquired a dependence on space (x = j/L, j = lattice site number) and time (t = j/L, j= time step). The main assumption is that the soliton current J l . 9 In the matrix notation it reads ∂ t ρ + ∂ x (ρ * v(r, l)) = 0.  9 The associated time variable could also be denoted by t (r) l , which is however avoided for notational simplicity. This is a collection of separated equations for each (a, i) ∈ I meaning that y (a) i 's are the normal modes of the generalized hydrodynamics [10]. leads to the conservation of the hole current ∂ t σ + ∂ x (σ * v(r, l)) = 0.

Domain wall initial condition
Let us embark on the study of transient behaviors of the randomized cBBS started from the domain wall initial condition. This will give the possibility to check some of the previous results obtained by using TBA and GHD with direct numerical simulations. Recall that an i.i.d. distribution is specified by the fugacity (z 1 , . . . , z n ) as explained in (3.4)-(3.9). We prepare the system initially in the i.i.d. distribution with fugacity (z 1L , . . . , z nL ) in the left domain x < 0 and a different fugacity (z 1R , . . . , z nR ) in the right domain x > 0. Then at time t = 0, we begin running a dynamics T (r) l for some (r, l) ∈ I, and observe the long time non-equilibrium behavior in the mixture of the soliton gas. This kind of setting is called partitioning protocol, which has been studied in many recent literatures. See [10,18,19] and the references therein. We shall present our GHD analysis first in a ballistic picture (Sec. 5.1) and then its diffusive correction (Sec. 5.2) followed by numerical verifications by simulation (Sec. 5.3).

Ballistic picture
The first approximation we employ is the ballistic scaling which implies that the normal mode y   Obviously heights in such a plot depend on the local quantity to consider but the locations of the plateau edges are common for all of them.
In general the normal mode y iR which are determined from the fugacity (z 1L , . . . , z nL ) or (z 1R , . . . , z nR ) by (3.24) and (3.10). In the mixed inhomogeneous system under consideration, we have the boundary condition y iR . For the neighboring plateaux P and P as in (5.2), the relation (5.1) indicates that there is a subset J ⊂ I such that the following hold 11 : Here v The boundary condition means that P → ∅ as ζ → −∞ and P → I as ζ → ∞. Thus a plateau configuration can be viewed as specifying a monotonous growth pattern of the label P. One starts from P = ∅ in the far left and proceeds to the right inductively to eventually get P = I in the far right. In the inductive step, we assume that a plateau P and the associated data {v (a) i (r, l) P | (a, i) ∈ I} are at hand. Then its right edge ζ * as in (5.2) and the set J which specifies the next plateau P by (5.7) are determined as This inductive procedure to trace the plateau profile is convenient to arrange GHD based numerical calculations. It is worth mentioning that the two alternatives of the normal mode variable y In what follows we will signify the solitons in (5.10), (5.11) as S iL and refer to P also as soliton content. If the initial domain is empty 12 for x < 0 or x > 0, then (5.11) or (5.10) becomes void, respectively.
So far our description of the plateaux is quite general. Now we propose a simplifying conjecture based on numerical experiments.
Conjecture. For any fugacity (z 1L , . . . , z nL ) and (z 1R , . . . , z nR ) satisfying (3.14), the following properties are valid: 1. Actually realized soliton contents P ⊆ I are limited to the following form for some d:   iL . Consequently, the plateaux can be numbered as 0, 1, 2, . . . from the left to the right, where the leftmost 0 covers the entire ζ < 0 region at least.
Thanks to (5.12) a plateau can be labeled by a vector d ∈ (Z ≥0 ) n−1 , which corresponds to the soliton content  Example 5.1. Consider the cBBS with n = 2, which is nothing but the ordinary single color BBS studied in [11]. We have the time evolutions T (r=1) l with l = 1, 2, . . ., and the plateaux are associated with 0 = d 0 ≤ d 1 ≤ d 2 ≤ · · · specifying their solitons contents. A typical setting is to let the initial domain in x > 0 or x < 0 be empty. Examples of the ball density plot in such cases are available in [11, Figure 4]. In either case, there are exactly l +1 plateaux P(d 0 ), P(d 1 ), . . . , P(d l ) where d j is given by Analytic formulas have also been obtained for many quantities in [11,, which are compatible with the above conjecture.

Diffusive correction
As we shall demonstrate in the next subsection, the actual plateau edges exhibit broadening for finite t whose main cause is the diffusive effect [10,11,45,46]. To describe it quantitatively, consider the adjacent plateaux P and P as in (5.2). We assume that the soliton contents of P and P differ only by one kind of solitons S Then the same argument as [11] (see also [46]) can be applied to obtain the envelope of an averaged local quantity Q(x, t) including the diffusive correction. Here we just mention a few key formulas used in the derivation: in the vicinity of x/t = ζ. The symbols Q(P) and Q(P ) denote the (averaged) values of Q on the plateau P and P , respectively. The most essential quantity in (5.19) is the constant Σ P,P (> 0) which controls the width of the edge broadening. It is given by hence M dr are those on the plateau P. 13 The α in (5.16) should not be confused with the one in (3.28).

Simulations
In this section we present some numerical simulations of the domain wall problem in the n = 3 model, and we compare the results with the GHD predictions. The method we use to simulate the BBS dynamics is similar to that of Ref. [11]. We consider a large but finite system with L sites, each having a B 1,1 tableau and a B 2,1 semistandard tableau.
As already mentioned at the beginning of Sec. 5, at time t = 0 the left part is initialized in some random i.i.d. state characterized by fugacities z 1,L , z 2,L and z 3,L . The right part is initialized in some random i.i.d. state characterized by fugacities z 1,R , z 2,R and z 3,R . As for the buffer, it is initially in the vacuum state. The role of this buffer region is to ensure that no soliton emanating from the left or the right has had enough time to reach the last site (x = 2L/3 − 1) at the end of the simulation at t = t max (we typically take t max = 1000, but some simulations have bee carried out up to t = 3000)). It is also important to check that the fastest solitons starting from the leftmost sites of the system have not had enough time to reach the left/right boundary (x = 0) at the end of the simulation. These conditions ensure that the finite size does not influence the behavior of the system in the region x ∈ [0, L/3 − 1] at least up to time t max . To match these conditions we need to scale the system size L as t max times the largest effective velocity.
The time evolution is computed by successive applications of a transfer matrix T (a) l with some large l (unless specified otherwise we used l = 100), and a = 1 or a = 2. We also report some results obtained using the product T  k (x, t) of solitons of size k and color a (using the local versions of the energies defined in App. B) as well as the mean box occupancies ( i=1,2,3 ) are recorded. In practice these densities are estimated by some average over a large number N samples of initial conditions. We typically use N samples ≥ 300000. k -soliton densities (with sizes k = 1, 2, 3 and 'color' a = 1, 2) as a function of ζ = x/t for t = 600 and t = 1000 and for the T (1) l=100 dynamics. The red line indicates the GHD prediction for t = ∞. The initial state is defined by the following fugacities in the left half: z 1,L = 1 z 2,L = 0.7 and z 3,L = 0.3. The right half is initially in the vacuum state. Parameters of the simulations: system size L = 300000 and number of random initial conditions N samples = 400000.
The figures 1, 2, 3 and 4 display the six soliton densities ρ (a=1,2) k=1,2,3 (x, t) at time t = 600 and t = 1000 as a function of ζ = x/t. In all cases the curves associated with different times are almost on top of each other, which indicates the ballistic nature of the transport. The above figures correspond to different initial states: the right part is empty in Fig. 1 and 4, the left part is empty in Fig. 2, and both sides are in a nontrivial state in Fig. 3. In all cases we observe that the soliton densities are in good agreement with the GHD prediction (Sec. 5.1). It can also be observed that in the vicinity of each transition between two consecutive plateaux the numerical curves do not form steep steps and have a finite slope. In addition the results are not exactly on top of each other, and the curves appear to become steeper when increasing t. This is a manifestation of the diffusive broadening discussed in Sec. 5.2. We will come back to this point below when discussing Fig. 9.
The figure 1 offers an illustration of the evolution of the soliton content from one plateau to the next, as discussed in Sec. 5.1. In this particular example the right half is initially empty, and all the y (a) iR therefore vanish. So, if a soliton type (i, a) belongs to the set P of a given plateau, then this species is absent in this plateau, its density ρ (a) i vanishes. In Fig. 1 the first step is located at ζ(1) 0.23815, and beyond this value of ζ the solitons S To go further in the description of the next plateaux one needs to look at the densities of solitons with larger (> 3) amplitude, and at larger values of ζ (data not shown). Such an analysis shows that for this particular initial condition the soliton content follows a somewhat irregular pattern up to the 12 th plateau: d 6 = (2, 4), d 7 = (2, 5), d 8 = (2, 6), d 9 = (2, 7), d 10 = (2, 8), d 11 = (2, 9), d 12 = (3,9). And then the sequence becomes regular: d n≥12 = (3, n − 3). After this first infinite series of plateaux a second series begins. The plateaux of the second series are free of color-2 solitons and they are characterized by: d ∞+1 = (4, ∞), d ∞+2 = (5, ∞), · · · , d ∞+n = (n + 3, ∞). k -soliton densities (with sizes k = 1, 2, 3 and 'color' a = 1, 2) as a function of ζ = x/t for t = 600 and t = 1000 and for the T (1) l=100 dynamics. The red line indicates the GHD prediction for t = ∞. The initial state is defined by the following fugacities in the right half: z 1,R = 1 z 2,R = 0.7 and z 3,R = 0.3, and the left half is initially in the vacuum state. Parameters of the simulations: system size L = 300000 and number of random initial conditions N samples = 400000.
The evolution of the soliton content is qualitatively different in the example of Fig. 2, where the initial state is the vacuum in the left half. There, the leftmost plateau P 0 which starts at ζ = −∞ is empty, and in this region any test color-2 soliton would have a vanishing velocity for the T (r=1) dynamics considered in this example. We thus have a situation where several velocities are degenerate: ∀i v (2) i (r = 1, l) P 0 = 0. The edge of P 0 is thus at ζ = ζ(1) = 0 and all the color-2 soliton densities jump from 0 to some nonzero value at this point, as can be seen in the lower panels of Fig. 2. One can also notice some density spikes at ζ = 0, where the all S (2) k soliton densities jump from zero to some nonzero value. This phenomenon appears generically at the edge of a vacuum plateau and it also occurs in the n = 2 model [11].  In the figure 4 we go back to a situation where the initial state is the vacuum for the right part, but the dynamics is generated by the product T (1) l=100 T (2) l=100 . In the GHD framework this means that the bare velocity v l=100 dynamics. Note that by construction 1 + 2 + 3 = 3, since each site contains 3 boxes (one in the tableau B 1,1 and two in the tableau B 2,1 ). The vertical red dotted line indicate the locations ζ(k) of the steps between the different plateaux, as predicted by GHD. The orange, yellow and blue lines are 1 , 2 and 3 from GHD (at t = ∞). The initial state is defined by the following fugacities in the left half: z 1,L = 1 z 2,L = 0.7 and z 3,L = 0.3. The right half is initially in the vacuum state. A large number of very narrow plateaux are observed to accumulate in the vicinity of ζ c 1.88708. Note that additional plateaux are present for ζ > 4 (not shown) and the system is still not in the vacuum state at the right of the figure (contrary to Fig. 6). Parameters of the simulations: system size L = 300000 and number of random initial conditions N samples = 300000.
l=3 dynamics. We again have an infinite series of narrow plateaux which accumulate in the vicinity of ζ 0.81586. The last step is located at ζ = 3, beyond this point the system is in the vacuum state (i.e 1 = 2, 2 = 1 and 3 = 0). Notice the density spike just before ζ = 3. The soliton speeds and densities in the three four plateaux are displayed in Fig. 8.
The figure 5 represents the 'letter' densities of '1', '2' and '3', for the same parameters as Fig. 1. By specializing (2.47) to the case n = 3 we can express these letter densities using the soliton densities: Fig. 5 an infinite number of very narrow plateaux are observed to accumulate in the vicinity of ζ = ζ c 1.88708. When approaching this limiting value ζ c from below, color-2 solitons of increasing size gradually disappear. Beyond that limiting velocity only color-1 solitons are left. In terms of (5.12) the vectors describing the soliton content for ζ > ζ c take the form d = (k, ∞). A very similar infinite series of narrow plateaux are observed for the T (2) l=3 dynamic, as illustrated in Fig. 6. Fig. 10 illustrates this phenomenon for the same initial state but with a few other values of l. In other words, this phenomenon does not disappear when taking a small carrier capacity l. Such an infinite series is also present with the T (2) l=100 dynamics, as illustrated in Fig. 7. In that case the color-1 solitons are absent for ζ beyond the accumulation point (found at ζ c 0.80019).
It should finally be noted that such a phenomenon is absent in the simpler n = 2 model [11]. It is also absent if one considers the n = 3 complete BBS with the same initial state but with the T (1) l=100 T (2) l=100 dynamics (Fig. 4). i=1···10 in the four first plateaux of the domain wall problem where the dynamics is generated by T (r=1) l=3 , and the initial state is characterized by z 1,L = 1 z 2,L = 0.7 and z 3,L = 0.3 and is empty in the right (same as in Fig. 6). The upper panels illustrate the fact that the velocity is not necessarily a monotonous function of the soliton size. In this example the velocity of color-1 solitons is growing for i = 1, 2, 3 and then it decreases and approaches a finite limit when i → ∞. As for the color-2 solitons, their velocties turn out to be increasing functions of the size. In the lower panels (densities) on can check that i) the density of S  solitons in the vicinities of the three first steps (located at x/t = ζ(1), x/t = ζ(2) and x/t = ζ(3)) and at three different times. With the choice (x/t − ζ(k))t 1/2 for the horizontal scale the curves almost collapse onto each other. This illustrates the diffusive nature of the step broadening. The red curves have been obtained without any free parameter (no fit) and represent the GHD prediction. They take the form of a complementary error function (see Sec. 5.2 and (5.19)). Initial state fugacities: z 1,L = 1 z 2,L = 0.7 and z 3,L = 0.3, z 1,R = 1 z 2,R = 0.9 and z 3,R = 0.1. System size L = 1.210 6 . Number of random initial conditions N samples = 400000. This simulation represents (2 · L) · N samples · t max ∼ 0.3 · 10 16 applications of the combinatorial R.
Finally Fig. 9 shows that the widths of the steps between consecutive plateaux scale like t 1/2 , in agreement with the diffusive corrections calculated in Sec. 5.2. This figure represents some soliton densities in vicinity of the step k = 1, 2, 3, plotted as a function of (x/t−ζ(k))t 1/2 . The fact that the data measured at different times almost fall onto the same error function curve demonstrates the diffusive nature of the step broadening. Moreover, the data are in good agreement with the GHD prediction (5.19) calculated in Sec. 5.2, including the quantitative value of the parameter Σ P,P . It should in particular be noted that the red curves in Fig. 9 are the GHD predictions and contain no adjustable parameter. In the bottom right panel of Fig. 9 the numerical data appear to be slightly below the GHD curve. The reason for this mismatch is not clear to us but it should however be noticed that it tends to decrease when the time increases (compare t = 1000, 2000 and 3000). It is therefore likely due to some finite-time effect. l . The initial state is defined by the fugacities z 1,L = 1 z 2,L = 0.7, z 3,L = 0.3 and with the vacuum state in the right part. For l = 8 the results of simulations (at t = 1000 and N samples = 3.10 5 ) are also shown. The sum of the three densities is equal to 3 by construction (notice the different vertical scales in the different panels). Three regions can be distinguished: i) In the region ζ < ζ c (l) some soliton of color 1 and color 2 are present. An infinity of narrow plateaux develop when ζ → ζ c (l), where the color-2 solitons of increasing size gradually disappear. ii) For ζ c (l) < ζ < l only color-1 solitons are present. There is a finite number (of order l) of plateaux in this region. iii) For l < ζ the system is in the vacuum state. The critical values of ζ are: ζ c (2) 0.57484, ζ c (3) 0.81586, ζ c (4) 1.02025, ζ c (8) 1.55262, ζ c (100) 1.88708, and they can be identified in the right panel as the points where 3 vanishes.

Summary and conclusions
In this paper we have introduced a new family of box ball systems, dubbed complete BBS. The models are indexed by a parameter n ∈ Z >1 and are associated to the quantum group U q ( sl n ). The cBBS family generalizes the original single-color model [1], which is recovered for n = 2. Each 'site' of the model is made of the n − 1 column shape semistandard tableaux with length 1, 2, · · · , n − 1. The 'balls' can have n different 'colors' and occupy the boxes of the above tableaux, respecting the semistandard rules. These states are then equipped with some integrable dynamics, constructed from the so-called combinatorial R. As for the simplest BBS, the time evolutions can be viewed as generated by a 'carrier' which propagates through the system at each time step and which takes, moves and deposits 'balls' from one place to another. In the case of the time evolution T (k) l , the load of the carrier is itself a rectangular tableau of shape k × l with k ∈ {1, . . . , n − 1} and l ∈ Z >0 . The integrability manifests itself from the fact that all these time evolutions commute (whatever k and l) and from the existence of a family of conserved energies. 14 Another consequence of the integrability is the existence of solitons. Each soliton is characterized by a 'color' a and amplitude (or size) i. These solitons interact with each other in some nontrivial way during collisions. The local structure of each site of the cBBS may look complicated compared to the other versions of n-color BBS, but, as explained in Sec. 2, it offers a remarkable simplification: the soliton scattering is completely diagonal, namely each soliton regains its color and amplitude after collisions. The effect of the interactions is encoded in a simple phase shift. The end of Sec. 2 presents the inverse scattering scheme for this model, which is based on the KSS bijection, and which allows one to construct the action-angle variables.
The Sec. 3 discusses the randomized cBBS, where we have considered generalized Gibbs ensembles of cBBS configurations. We focused on a special family of GGE, called i.i.d., which are characterized by n fugacities (one for each color) and where all the tableaux are independently distributed. The properties of such i.i.d. GGE can be studied using TBA and we obtained the free energy, the mean value of the conserved energies and the associated soliton densities. In Sec. 4 we have presented the generalized hydrodynamics of the cBBS. The first ingredient is the effective speed of each soliton species in a homogenous state. This could be obtained thanks to the TBA approach. Next, the currents carried by the soliton of each color are assumed to be given by the product of their effective speed and their density. For a system which is inhomogeneous but where spatial gradients are sufficiently small so that it can locally be approximated by a GGE, the continuity equations for the soliton currents give the GHD equations. These equations take a particularly simple form when re-written in terms of the Y-variables of the TBA, which shows that these variables are the normal modes of the GHD. As a concrete application of GHD we studied in Sec. 5 the evolution of the n = 3 model starting from a domain wall initial state with two different i.i.d. states in the left and in the right halves. In very good quantitative agreement with the GHD prediction, the high-precision simulations showed that at sufficiently long times the state of the system becomes piecewise constant in the variable ζ = x/t. Inside each plateau, the soliton densities and the ball densities measured in the simulations match perfectly the GHD expectations. Compared to the simpler n = 2 model [11], the additional color degree of freedom leads to a complex evolution of the soliton content from one plateau to the next, and opens up the possibility to have an infinity of narrow plateaux which accumulate to some critical point ζ c . Finally, we have described the diffusion effect that is responsible for the broadening of the steps between consecutive plateaux. Here again we observed a good quantitative agreement between the simulations and the theoretical results. To our knowledge, this paper achieves the first systematic and successful application of GHD to an integrable system associated with a higher rank quantum algebra.
We wish to conclude this article by mentioning a few possible future directions. The equation (4.1) for the effective velocities is quite fundamental. It would be very interesting to have a proof of it, in the line of what was done in [11,Sec.4.2] for n = 2. See also [25] for n = 2. It would also be appealing to invent and explore other protocols to set the system out of equilibrium. This might be done, for example, with a system that is open at its boundaries.
One could also think of more general initial states, homogeneous or inhomogeneous, defined for instance by some nontrivial weights on the configurations, and look at the evolution toward equilibrium. This could be, for instance, a domain wall between two general (non-i.i.d.) GGE. Another interesting question would be to explore the possibility of some anomalous transport in these models. Investigating numerically other values of n, possibly large ones, might also reveal some new phenomena. It also seems worth comparing the dynamics of the present cBBS with non-complete n-color BBS [14,17,20]. Finally, a possible direction would be to include some integrability-breaking perturbations in the model. We may expect the disappearance of the density plateaux, and it should be possible to observe numerically and to analyze some possible crossover from the integrable regime to a chaotic one where solitons are no longer stable but only short-lived. Such a setup could give rise to some prethermalization [47], possibly described by some GGE.
A Algorithm for the combinatorial R We use the notation in Sec. 2.2. Given b ⊗ c ∈ B k,l ⊗ B k ,l , we present an algorithm for finding the imagec ⊗b = R(b ⊗ c) ∈ B k ,l ⊗ B k,l of the combinatorial R following [35, p55]. As noted after (2.7), the task is to construct the solution to the equationb ·c = c · b which is unique when the tableaux are rectangular.
A skew tableau θ is a part of a semistandard Young tableau having a skew Young diagram shape. We call θ an m-vertical strip if it contains at most one box on each row and the total number of the boxes is m.
Let P = c · b be the product tableau [12]. Its shape Young diagram includes a k × l rectangle. Denote by P the NW part of P corresponding to the k × l rectangle, and by P/P the skew tableau. The number of boxes of these tableaux (denoted by | · |) are given by |P | = kl + k l , |P | = kl, |P/P | = k l . (A.1) Step 1. Label each box of P/P with {1, 2, . . . , k l }. Let θ 1 be the rightmost vertical k -strip in P/P as upper as possible. Remove θ 1 from P/P and define the vertical k -strip θ 2 in a similar manner. This can be continued until P/P is decomposed into the disjoint union of the vertical k -strips θ 1 , θ 2 , . . . , θ l . Now the label is obtained by assigning the boxes of θ i with k (i − 1) + 1, k (i − 1) + 2, . . . , k (i − 1) + k from the bottom to the top.
Step 2. Reverse row bumping from P in the order of the label.
Find the semistandard Young tableaux P 1 , P 2 , . . . , P k l and w 1 , w 2 , . . . , w k l ∈ [1, n] such that P = (P 1 ← w 1 ); reverse bumping of the letter in the box 1, ; reverse bumping of the letter in the box 2, · · · P k l −1 = (P k l ← w k l ); reverse bumping of the letter in the box k l .

(A.2)
Here ← stands of the row insertion as in Sec. 2.2. The splitting of P i−1 into (P i ← w i ) is done by the reverse row bumping. Starting from the box labeled i in Step 1 and the letter in it, the bumping goes upwards row by row as follows. Given a letter α in a row, find the box in the adjacent upper row filled with the largest β such that β < α and the rightmost one in case more than one β is contained. Then let α occupy the box by bumping out β. Repeating this procedure one eventually bumps out a letter w i from the top row of P i−1 thereby getting also P i as the remnant tableau. Note that |P k l | = kl.

C KSS bijection C.1 Paths
Consider the product set of the form B = B k 1 ,1 ⊗ B k 2 ,1 ⊗ · · · ⊗ B k N ,1 with k i ∈ [1, n − 1] containing column shape tableaux only. For the definition of B r,s , see Sec. 2.2. States of our cBBS correspond to the choice k i ≡ i mod n − 1 and N = (n − 1)L. In this appendix, we do not assume the boundary condition of the BBS, and call the elements of B paths. Our aim is to explain the special case of the KSS bijection [48] corresponding to the above B along the convention adapted to this paper.

C.2 Rigged configurations
We keep B = B k 1 ,1 ⊗ · · · ⊗ B k N ,1 as above. Consider a multiset (a set allowing multiplicity of elements) where M ≥ 0 is arbitrary. Each triplet s = (a, j, r) in S is called a string. It carries color, length and rigging denoted by cl(s) = a, lg(s) = j and rg(s) = r, respectively. S is a rigged configuration for B (or just a rigged configuration for short) if 0 ≤ r i ≤ p lg(s) ≥ 0 for all s ∈ S, which is already a nontrivial constraint on the multiset {(a i , j i ) | i = 1, . . . , M } depending on k 1 , . . . , k N . This is the reason why the rigged configurations are defined with respect to a given B = B k 1 ,1 ⊗B k 2 ,1 ⊗· · ·⊗B k N ,1 . The dependence on k 1 , . . . , k N comes from (C.8) via the first terms in the RHS of (C.9) or (C.10).
We regard a string s = (a, j, r) as a length j row with the rigging r attached to its right side. Collecting such color a strings yields a Young diagram whose rows are assigned with riggings. Further collecting them for the colors a = 1, 2, . . . , n − 1 leads to a combinatorial object, which was the original rigged configuration in the literatures [49,50]. For instance for n = 3, 1 = 2. In general, strings with the same color and length form a rectangular block. Objects obtained by permuting the riggings attached to them should not be distinguished since originally it is just the multiset (C.7). Practically, one needs to keep track of the vacancy (C.10) to check (C.8) and this requires the data B as well. So it is customary and convenient to also include these information in the graphical representation. For instance if B = B ⊗7 with B = B 1,1 ⊗ B 2,1 , (C.12) is more detailed as The number assigned in the left of each rectangular block (not a row) is the vacancy. In this example they are p 1 = 3, and the condition (C.8) is certainly satisfied.
In general we have the Young diagrams µ (1) , . . . , µ (n−1) encoding the color and length of strings, which form the partial data (µ (1) , . . . , µ (n−1) ) called configuration. It corresponds to the multiset {(a i , j i ) | i = 1, . . . , M } obtained from (C.7). The symbol r (a) in (C.13) stands for the rigging attached to µ (a) . To compute the vacancy p Note that the same symbol λ has been used to denote the weight of a path in (C.1) and that of a rigged configuration in (C.14). This we did intentionally since they are going to be identified under the KSS bijection (C.17) in the next subsection. The number of rigged configurations is given by the so called

C.3 KSS bijection
There is a one to one correspondence between the finite sets P + (B, λ) and RC(B, λ). It was shown in a more general case of B = B k 1 ,s 1 ⊗ · · · ⊗ B k N ,s N in [48] generalizing the pioneering works [49,50]. Their original motivation was to provide a bijective (combinatorial) proof of the so called Fermionic character formula for the finite dimensional representations of quantum affine algebras. It originates in the string hypothesis in the Bethe ansatz. See the introductions of [40,51] for the rich history going back to Bethe himself [3]. Here we describe the algorithm for the bijection Φ : P + (B, λ) lg(s) ≥ 0 for the existing strings s only. 16 By Φ we actually mean the totality of the maps (C.17) for all (B, λ).