Randomized box-ball systems, limit shape of rigged configurations and Thermodynamic Bethe ansatz

We introduce a probability distribution on the set of states in a generalized box-ball system associated with Kirillov-Reshetikhin (KR) crystals of type $A^{(1)}_n$. Their conserved quantities induce $n$-tuple of random Young diagrams in the rigged configurations. We determine their limit shape as the system gets large by analyzing the Fermionic formula by thermodynamic Bethe ansatz. The result is expressed as a logarithmic derivative of a deformed character of the KR modules and agrees with the stationary local energy of the associated Markov process of carriers.


Background and main results
1.1. Box-ball systems. The box-ball system (BBS) [40] is an integrable cellular automaton in 1 + 1 dimension. By now it has been generalized widely, and numerous aspects have been explored connected to quantum groups, crystal base theory (theory of quantum groups at q = 0), solvable lattice models, Bethe ansatz, soliton equations, ultradiscretization, tropical geometry and so forth. See for example a review [14] and the references therein. Here is an example of time evolution T (1) ∞ in the 3-color BBS [38] in the notation specified later: A letter 1 denotes an empty box whereas a = 2, 3, 4 is the one filled with a ball with "color" a. Initially there are three solitons 2222, 332 and 43 with amplitude 4, 3, 2. They proceed to the right with the velocity equal to the amplitude and eventually undergo collisions with messy intermediate states. However after enough time steps they come back and separate in exactly the original amplitude 2, 3, 4 in the reverse order, not being smashed into pieces nor glued together 1 . This is a manifestation of the integrability, or put more practically, existence of conserved quantities, either explicit or hidden, governing the dynamics. The original time evolution t → t + 1 in the n-color BBS was defined by a ball moving algorithm [38] as K 2 • K 3 • · · · • K n+1 , where K a moves every ball with color a once starting from the leftmost one successively to its nearest right empty box. So the number of balls with each color is obviously preserved. With some effort it is also possible to show that the list of amplitude of solitons, if defined properly, is also conserved. In the above example it is a partition (4,3,2). However a quite nontrivial and essential question is; what is the complete set of conserved quantities for the general n-color BBS? 1 In n-color BBS in general, the balls are labeled with 2, 3, . . . , n + 1, and a consecutive array of balls n + 1 ≥ a 1 ≥ · · · ≥ am ≥ 2 separated sufficiently from other balls behaves as a soliton with amplitude = velocity = m. Various choices of a 1 ≥ · · · ≥ am yield internal degrees of freedom of solitons, like quarks in the hadrons uud (proton), udd (neutron), uds (Λ) etc. The example demonstrates essential features of the soliton scattering; interchange of internal degrees of freedom and phase shift of asymptotic trajectories. The final list of solitons is known to be independent of the order of collisions of the initial ones (the Yang-Baxter property).

1.2.
Rigged configuration as action-angle variables. The answer is known to be an n-tuple of Young diagrams. It was derived from the solution of a more general problem of constructing the actionangle variables of the BBS [25,26]. By action variables we mean a set of conserved quantities and by the angle variables those linearizing the dynamics. The integrability of BBS allows us to transform the system bijectively into action-angle variables(!) For the BBS states in the above time evolution, they are combinatorial objects that look as follows: 4t 6 + 3t 11 + 2t For the n-color BBS in general, there are an n-tuple of Young diagrams µ 1 , . . . , µ n in which each row is assigned with an integer. The n-tuple of Young diagrams and the assigned integers are called configuration and rigging, respectively. Thus in short, the action-angle variables of the BBS are rigged configurations [26]. It is the configuration that is preserved and the riggings that change linearly in time. One indeed sees that the first Young diagram µ 1 = (4, 3,2) gives the list of amplitude of solitons which remains invariant under the time evolution. The other ones µ 2 , . . . , µ n are "higher" conserved quantities reflecting the internal degrees of freedom of solitons 2 . The n-color BBS is endowed with the higher time evolutions besides the simplest one K 2 • K 3 • · · · • K n+1 mentioned before. They are all commutative and change the riggings attached to µ 2 , . . . , µ n linearly.
Rigged configurations for type A n has been formulated most generally in [22] extending the invention [19,20] in 1980's. These works were devoted to a proof of the Fermionic formula for a generalized Kostka-Foulkes polynomials (cf. [30]) by establishing an elaborate bijection between rigged configurations and other combinatorial objects. Roughly speaking in the context of BBS, the bijection provides the direct and inverse scattering maps [25,26] {BBS states} ←→ {action-angle variables}, (1) which transform the nonlinear dynamics in BBS to a straight motion. The n-tuple of Young diagrams form a label of the iso-level sets of BBS. The Fermionic formula tells the multiplicity of a given iso-level set.
1.3. Randomized BBS and main result. Now let us embark on a randomized version of the story. We assume that some probability distribution on the set of BBS states has been introduced. Then it is natural to ask; (i) What is the probability measure on the n-tuple of Young diagrams induced by (1)? (ii) What is the limit shape of them when the system size L of the BBS tends to infinity?
In this paper we answer (i) and (ii) for the most general BBS associated with the Kirillov-Reshetikhin (KR) crystals [17] of the quantum affine algebra U q (A (1) n ). The randomness of the BBS states we will be concerned with is the product of the one at each site. The latter is the probability distribution on a single KR crystal just proportional to e wt . (See (16).) Under this simple choice, the answer to (i) is given by the Fermionic form itself multiplied with the Boltzmann factor accounting for the e wt contribution as a chemical potential term. (See (43).) The Fermionic form measure is quite distinct in nature from the well studied ones like the Plancherel measure for the symmetric group and/or its Poissonized versions. It fits an asymptotic analysis by the thermodynamic Bethe ansatz (TBA) [42]. The method employs the idea of the grand canonical ensemble and captures the equilibrium characteristics of the system by a variational principle. The equilibrium condition shows up as the so called TBA equation. It plays a central role together with the equation of state connecting the density of balls with fugacity. Our TBA analysis is essentially the spectral parameter free version of [24, sec.14]. In particular the Y-system and the Q-system (Q of the KR module come into the game naturally. It turns out that a proper scaling of the Young diagrams is to shrink them vertically by 1/L. This feature will be established in [23] by invoking the large deviation principle. The resulting limit shape is described by the logarithmic derivative of the deformed character of the KR modules as (2) See (72) and (24) for the definition of the deformed character Q The data (r, s) is specified according to the choice of the set of local states in the BBS 4 . The quantity (2) coincides with the stationary local energy of a carrier in the relevant KR crystal derived in (21). Independent variables in the deformed characters are linked to the prescribed fugacity of the BBS by the equation of state (63) or equivalently (66). This general and intrinsic answer to the above question (ii) is our main result in this paper. Further concrete formulas are available for the simplest n-color BBS [38] in terms of Schur functions in (93) and (94). It will be interesting to investigate the results in this paper further in the light of recent results on the BBS from probabilistic viewpoints [3,5,23,29].
1.4. Outline of the paper. In Section 2 we recall basic facts on generalized BBS necessary in this paper. In Section 3 we consider the BBS in a randomized setting. It amounts to introducing a Markov process of carriers associated to each time evolution T (a) i . We construct a stationary measure of the process quite generally by the character of the relevant KR modules (Proposition 3.2). It leads to the stationary local energy (21) or equivalently (22). In Section 4 we recall the Fermionic formula based on [10,22] as a preparation for subsequent sections. The deformed character in (37), (72) and its logarithmic derivative will be the building blocks in describing the limit shape. Section 5 is the main part of the paper. We identify the Fermionic form with the probability measure on the n-tuple of Young diagrams induced from the randomized BBS via its conserved quantities. By a TBA analysis, a difference equation characterizing the limit shape of the Young diagrams is derived. Our main result is Theorem 5.1, which identifies the solution to the difference equation with the stationary local energy obtained in Section 3. It reveals a new connection between TBA and crystal theory via the limit shape problem. In Section 6 we deal with the simplest example corresponding to the n-color BBS in [38]. The scaled column length of the Young diagrams are given explicitly in terms of the Schur function involving the ball densities. We check the result against the stationary local energy of a randomly generated BBS states numerically and confirm a good agreement. Section 7 contains a summary and discussion. We conjecturally describe the difference equation and its solution like Theorem 5.1 uniformly for the BBS associated with the simplylaced quantum affine algebras U q (ĝ) withĝ = A (1) n , D (1) n , E (1) 6,7,8 . We also suggest some future problems as concluding remarks.
2. Box-ball systems 2.1. KR crystals. Consider the classical simple Lie algebra of type A n . We denote its Cartan matrix by (C ab ) n a,b=1 , where C ab = C ba = 2δ ab − θ(a ∼ b) and a ∼ b means that the two nodes a and b are connected by a bond in the Dynkin diagram, i.e. |a − b| = 1. 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 . We use the set of positive roots ∆ + , the weight lattice P = n a=1 Z̟ a , the root lattice Q = n a=1 Zα a and their subsets P + = n a=1 Z ≥0 ̟ a , Q + = n a=1 Z ≥0 α a . Denote the irreducible module with highest weight λ ∈ P + by V (λ) and its character by ch V (λ). The latter belongs to n be the non-twisted affinization of A n [16] and U q = U q (A (1) n ) be the quantum affine algebra (without derivation operator) [4,15]. There is a family of irreducible finite-dimensional representations {W (r) s | (r, s) ∈ [1, n] × Z ≥0 } of U q called Kirillov-Reshetikhin (KR) module 5 named after the related work on the Yangian [21]. As a representation of A n , W (r) s is isomorphic to V (s̟ r ). W (r) s is known to have a crystal base B (r) s [18,17]. Roughly speaking, it is a set of basis vectors of a U q -module at q = 0. B (r) s is called a KR crystal. It is identified with the set of semistandard tableaux of rectangular shape (s r ) with letters from {1, 2, . . . , n + 1}. The highest weight element of B (r) s , which is the tableau whose lm . Representation theoretically, it is the character of a generalized Demazure module [31]. 4 The original n-color BBS [38] corresponds to the choice (r, s) = (1, 1). 5 The actual KR modules carries a spectral parameter. In this paper it is irrelevant and hence suppressed.
j-th row is all j, is denoted by u (r) s . For two crystals B 1 , B 2 their tensor product B 1 ⊗ B 2 is well defined, and as a set Before explaining necessary ingredients related to KR crystals, we review a notion of tableau product S · T for two tableaux S, T . Let row(T ) be a row word of a tableau T . It is obtained by reading letters from bottom to top, left to right in each row. Let row(T ) = u 1 u 2 · · · u l and we apply to S the row bumping algorithm [8] successively as The resulting tableau is nothing but S · T . Alternatively, let row(S) = v 1 v 2 · · · v m and apply to T the column bumping algorithm successively as The result also gives S · T .
We are ready to review the combinatorial R and the (local) energy H. Let b, c be elements of B (r) s determined by the following combinatorial rule [36]. The image of R is given in such a way that R(b ⊗ c) =c ⊗b is equivalent to c · b =b ·c. The fact that for c · b there is a unique such pair (b,c) is assured since the decomposition of the tensor product module V (i̟ a ) ⊗ V (s̟ r ) is multiplicity free. The value H(b ⊗ c)(= H(c ⊗b)) is defined to be the number of nodes strictly below the max(a, r)-th row of the tableau c · b. By definition, H is nonnegative and H(u

2.2.
Deterministic box-ball system. The original BBS was introduced in [40]. Since then it has been generalized from various viewpoints. One of such generalizations was done by using KR crystals as we describe below.
Take a sufficiently large integer L and consider B i . Graphically, it can be shown as below.
We call B    i . We call this extra tensor factor the barrier.
Next we recall the conserved quantity under the time evolution T (a) i introduced for a = 1 in [7] for the n-color BBS [38]. In order to make c in (3) to be u (a) i we attach the barrier if necessary and assume the number of the tensor factors in the quantum space to be L. Define c j (j = 1, 2, . . . , L) by The definition corresponds to setting the j-th vertex in the previous diagram as follows.
We introduce the (row transfer matrix) energy by i (b) following the same argument as [7] for a = a ′ = 1. Moreover, supplement of barriers does not change E (k) l for any k, l. We note that these features are valid even when the quantum space (B In general when the carrier is B (1) i and the quantum space is (B 1 ) ⊗L , the dynamics on the latter reproduces the ball-moving algorithm in the n-color BBS [38] as i → ∞.
The next example is the case when the carrier is again B  13 11 In general choosing the quantum space as (B (1) s ) ⊗L corresponds to the boxes with capacity s. The last example is the case when the carrier is B   This is the most general situation. Local states and carriers are no longer simple boxes but possess a structure of a shelf with a nontrivial constraint on the arrangement of balls from the semistandard condition of the tableaux.
Introduction of carriers [39] as a hidden dynamical variable of BBS was a corner stone in the development of the theory. It provided the apparently nonlocal ball moving algorithm with a local description encoded in a single vertex in the above diagrams. A further discovery that these vertices are nothing but the combinatorial R unveiled the nature of BBS as solvable vertex models [1] at q = 0, where time evolutions are naturally identified with commuting row transfer matrices [9,7]. As we will see in Section 3, carriers also play a fundamental role in the randomized BBS via their Markov processes.

Rigged configuration as action angle variables.
Here we review a combinatorial object called rigged configuration and see how it is used to linearize the BBS dynamics. Rigged configurations are defined based on data {(k j , l j )} 1≤j≤L such that (k j , l j ) ∈ [1, n] × Z ≥1 . Through the Kirillov-Schilling-Shimozono (KSS) bijection which we discuss later, it is related to the tensor product of KR crystals B lL . A rigged configuration consists of a configuration, an n-tuple of Young diagrams µ 1 , . . . , µ n , and riggings, sequence of nonnegative integers attached to each row of µ a for a ∈ [1, n]. Let m A configuration is required to satisfy p i . Among riggings of the rows of the same length in µ a , the order does not matter. So we label riggings in non increasing order when going downwards. From these definitions one can immediately write down the number of the rigged configurations with the prescribed configuration µ 1 , . . . , µ n as This is an ultimate generalization of the celebrated Bethe formula [2, eq.(45)] due to [19,20,22] for type A n . See [11, sec.1] for a historical account and [24, sec.13] for a concise review. We will come back to this Fermionic form as the main object of the TBA analysis in Section 5.
The KSS bijection [22] gives an algorithm to construct an element of the tensor product of KR crystals from a rigged configuration. The image of this bijection consists of special elements which we call highest states. Representation theoretically, they correspond to highest weight vectors of B. It is equivalent to saying that the tableau product b L · · · · · b 1 (b j ∈ B (kj ) lj ) is a tableau such that letters in the i-th row are all i.
The KSS bijection separates the BBS states into action and angle variables. It is known [26] that if b is a highest state, then the application of T (a) i causes, in the rigged configuration side, the increase of riggings by δ ac min(i, l) when they are attached to the length l row of µ c .
Identifying rigged configurations originating in the Bethe ansatz [2] with action-angle variables of BBS implies a correspondence between Bethe strings in the former and solitons in the latter. This is natural as we will also comment in the end of Section 5.1. As far as the action variables are concerned, this soliton/string correspondence [25,26] is quantified most generally as [35] Remember that the LHS is the row transfer matrix energy in (4), which was indeed known (for a = 1) to measure the amplitude of solitons [7] for the original n-color BBS [38]. The RHS is defined by (6) from the rigged configuration which is essentially an assembly of Bethe strings [2,19,20]. Thus the LHS and the RHS in (8) are referring to solitons and strings, respectively. Our main result Theorem 5.1 in this paper may be regarded as a generalization of (8) to a randomized situation.

Example 2.3. We give examples of the KSS bijection for A
3 . An element of (B is a highest state which corresponds to the rigged configuration below. The numbers left to the Young diagram are vacancies. p is a highest state of (B 2 ) ⊗10 which corresponds to the following one.
The condition on the sum is depicted as a vertex in (3) as In [23] it has been shown that this Markov process is irreducible and has the unique stationary measure for r = s = 1. We conjecture and assume the irreducibility for general r, s in what follows. Denote the resulting stationary measure bỹ π The combinatorial R is given by By the definition the transitions in the first, second and the third column happen with probabilities p 1 , p 2 , p 3 , respectively. Thus denotingπ (1) 2 (23) =π 23 etc for short, the stationary condition reads π 11 = p 1π11 + p 1π12 + p 1π13 , Under the assumption p 1 + p 2 + p 3 = 1, these equations admit a unique solution (π ij ) such that 1≤i≤j≤3π ij = 1. For instance let us parametrize p i as Thenπ ij is given bỹ is determined by the following stationary condition of the carrier process: where the latter equality follows from (9). The carrier and local states are taken from B For a partition λ = (λ 1 , λ 2 , . . . , λ n+1 ), let s λ (w 1 , . . . , w n+1 ) denote the associated Schur polynomial [30]: This is the well-known Weyl formula for the character chV (λ) under the identification of λ with n i=1 (λ i − λ i+1 )̟ i ∈ P + . We use a special notation when λ is a rectangle.
The following proposition gives an explicit expression for the stationary measureπ for any (a, i) ∈ [1, n] × Z ≥1 satisfies the stationary condition (12) and the normalization condition Proof. Consider the combinatorial R Take the exp wt of the both sides and sum over This yields (12) withπ s . The normalization condition is obvious from (14). Proposition 3.2 tells that as long as the randomness π (r) s of the local states are taken to be proportional to e wt as in (16), the stationary measureπ See (4) for the definition of E We call this stationary local energy for the carrier from B although it is suppressed in the notation. To write (21) more concretely, consider the irreducible decomposition of the tensor product V (i̟ a ) ⊗ V (s̟ r ). It is multiplicity free as noted in Section 2.1, and results in the identity of the Schur functions as Here P (a,r) i,s denotes the set of partitions (Young diagrams) labeling the irreducible components described by the Littlewood-Richardson rule. Concretely one has where ℓ(ν) denotes the length of the partition ν. From the description of the local energy H in Section 2.1, the result (21) is expressed as s . The local energy takes the value H(x ⊗ y) = k. Thus h .
We will use this formula with s = 1 in Section 6.
Here w is a parameter having nothing to do with w 1 , . . . , w n+1 in (15). The element b is the one occurring at the position r by sending b j ∈ B (kj ) lj to the left by successively applying the combinatorial R as In particular we set b (r) In contrast to the energy associated with the row transfer matrices (4), the quantity D in (25) corresponds to the energy of a corner transfer matrix 6 , which goes back to [1, chap.13]. In fact, using the Yang-Baxter equation for the combinatorial R it can be identified with the sum of the local energy associated to all the L(L − 1)/2 vertices in the following diagram (L = 3 example).
This quadrant structure is essentially a combinatorial counterpart of [1, Fig.13.1(b)]. By the definition we have due to (14). In this sense χ w (B) is a w-deformation of the character ch ⊗ L i=1 V (l i ̟ ki ) . See [31] for a representation theoretical study.

4.2.
Fermionic formula. Given a tensor product B = B (k1) l1 ⊗ · · · ⊗ B (kL) lL and λ ∈ P , we define the The quantity p By the definition the summand corresponding to m = (m The necessity for λ ∈ P + is seen by noting that (35) and (31)  i ) (a,i)∈[1,n]×Z ≥1 such that the above condition (i) and (ii) are satisfied for some λ obeying (36). Those m are called configurations. A configuration is equivalent to an n-tuple of Young diagrams via (44). They obey nontrivial constrains originating from the above (i) and (ii). To determine their asymptotic shape in the large L limit is a main theme of this paper. ⊗ · · · ⊗ B (kL) lL , the following equality is valid: where the sum extends over those λ satisfying (36).
The simplest case L = 1 of Theorem 4.2 gives (14). Namely one has which is actually independent of w 9 . Fermionic forms for general affine Lie algebra were introduced for non-twisted [10] and twisted [11] cases inspired by those for the Kostka-Foulkes polynomials [30] which correspond to A  ⊗ · · · ⊗ B (kL) lL be arbitrary. For any (a, i) ∈ [1, n] × Z ≥1 the following equality holds 10 : where φ = i + L j=1 δ a,kj min(i, l j ). Actually (39) was shown in [10] by substituting (37) to the three terms and using a decomposition of the Fermionic form. As a corollary of Proposition 4.3 and (27) with empty B we see that the classical character Q To validate this at i = 0 with Q in the remainder of this section and (73). The following result resembles the Wick theorem.
The equality is invalid without specialization to w = 1.
Proof. From (24) and (25) we have where the sum is taken over ⊗ · · · ⊗ B L in the notation of (26). By changing the summation variables, the summand corresponding to the pair i < j is expressed as  occurs with the probability proportional to e wt(u) = n a=1 z ua a for wt(u) = u 1 ̟ 1 + · · · + u n ̟ n ∈ P . We shall concentrate on the regime of the parameters z 1 , . . . , z n such that z a > 0 and n b=1 z C ab b > 1 for all a ∈ [1, n]. In view of n b=1 z C ab b = e n b=1 C ab ̟ b = e αa , it means α a > 0 for all the simple roots α a of A n . Thus the local states closer to the highest weight element u    s ) ⊗L whose probability distribution is proportional to e wt(b1⊗···⊗bL) . The randomized BBS in Section 3 corresponds to E L (r, s). 10 When i = 1, the factor B It slightly differs from E + L (r, s) in which the local states are not i.i.d. due to the nonlocal constraint of being highest. Both of them induce a probability distribution on the set of n-tuple of Young diagrams µ 1 , . . . , µ n by taking the conserved quantities. In the regime α 1 , . . . , α n > 0 under consideration, the highest condition on b 1 ⊗ · · · ⊗ b j ⊗ b j+1 ⊗ · · · ⊗ b L ∈ (B (r) s ) ⊗L becomes void almost surely for the right part b j+1 ⊗ · · · ⊗ b L in the limit L ≫ j → ∞. Since the large L asymptotics µ 1 , . . . , µ n does not depend on the left finite tail of b 1 ⊗ · · · ⊗ b L , we claim that those induced from E L (r, s) and E + L (r, s) coincide. (This "asymptotic equivalence" of E L (r, s) and E + L (r, s) is discussed in more detail for B For E + L (r, s), the conserved quantities µ 1 , . . . , µ n are the Young diagrams in the rigged configurations obtained by the KSS bijection. Therefore their (joint) probability distribution is explicitly given by In what follows we will identity the n-tuple of Young diagrams µ = (µ 1 , . . . , µ n ) with the data m = (m This will serve as the source of (L, r, s)-dependence of Prob(µ 1 , . . . , µ n ). The parameters β 1 , . . . , β n are chemical potentials or inverse temperatures in the context of the generalized Gibbs ensemble. As we will see in (57), they are actually the simple roots α 1 , . . . , α n . Therefore the factor e − n a=1 βa i≥1 im (a) i in (43) is just e −Ls̟r+λ due to (35). Besides the irrelevant constant e −Ls̟r , the factor e λ here indeed incorporates the relative probability e wt(b1⊗···⊗bL) adopted in E + L (r, s). Note that Prob(µ 1 , . . . , µ n ) = 0 unless m = (m (a) i ) (a,i)∈[1,n]×Z ≥1 is a configuration in the sense explained after (36). Finally Z L is given by (46).
Our aim is to determine the "equilibrium", i.e., most probable configuration under the probability distribution (43) when L tends to infinity. It will be done by the method of grand canonical ensemble with the partition function, namely the generating function of (43): The latter expression tells that Z L is a generating series of the branching coefficients [V (s̟ r ) ⊗L : V (λ)]. Numerous combinatorial objects labeling the irreducible components V (λ) and their counting formulas are known in combinatorial representation theory and algebraic combinatorics. In the original work by Bethe himself [2], a considerable effort was devoted to the completeness issue of his own string hypothesis. The succeeding development [20,21] assembled the Bethe strings and visualized them as rigged configurations. These works produced the Fermionic counting formula (46) for the representation theoretical quantity (47). A further insight, soliton/string correspondence (see Section 2.3) gained after entering this century, elucidated that the Bethe strings are nothing but the BBS solitons for which one can formulate an integrable dynamics based on KR crystals. It endowed the individual term in the sum (46) with a natural interpretation as the partition function of the BBS with a prescribed soliton content m [25,26]. In short the BBS provided the Fermionic formula with a refinement via a quasi-particle picture. Physically speaking the BBS solitons are bound states of magnons over a ferromagnetic ground state of an integrable A (1) n -symmetric spin chain deformed by U q=0 (A (1) n ). 5.2. TBA equation and Y-system. We are going to apply the idea of TBA [42] to the system governed by the grand canonical partition function (46). Similar problems have been studied in the context of ideal gas of Haldane exclusion statistics. See for example the original works [37,41,13] and a review from the viewpoint of a generalized Q-system [24, sec.13]. In fact our treatment here is a constant (spectral parameter free) version of the TBA analysis in [24, sec.14, 15]. In Theorem 5.1 it will be shown that the results coincide quite nontrivially with those obtained from the crystal theory consideration.
In the large L limit, the dominant contribution in (46) come from those m = (m (a) i ) exhibiting the L-linear asymptotic behavior where ρ . This fact will be justified by invoking the large deviation principle in [23]. From (45) and (32) the scaled variables are related as This is a constant version of the Bethe equation in terms of string density ρ i ) that minimizes the "free energy per site" This is (−1/L) times logarithm of the summand in (46) to which the Stirling formula has been applied. Note that (48) is consistent with the extensive property of the free energy, which enabled us to remove the system size L as a common overall factor. We have introduced a cut-off l for the index i, which will be sent to infinity later. Accordingly the latter relation in (49) should be understood as ε The TBA equation is equivalent to the Y-system The Y-system is known to follow from the Q-system (40) by the substitution (cf. [24,Prop. 14 where Q (a) i ∈ Z[z ±1 1 , . . . , z ±1 n ] is defined in (14). Now we take the boundary condition (54) into account. The left one Y The result [10, Th. 7.1 (C)] tells that lim l→∞ (Q > 1 under consideration. Thus the large l limit of (56) can be taken, giving In this way the chemical potentials β a are naturally identified with the simple roots α a . We shall keep using the both symbols although.
To summarize so far, we have determined the equilibrium configuration ρ eq of ρ = (ρ (a) i ) implicitly by (49), (52), (55) and (57) in terms of the chemical potentials β 1 , . . . , β n . The next task is to relate them to the canonically conjugate densities which are "physically more controllable". It amounts to formulating the equation of state. This we do in the next subsection.

Equation of state for randomized BBS.
From now on we will only treat the equilibrium values and frequently omit mentioning it. Let us calculate the equilibrium value of the free energy per site (50). First we use (52) to rewrite (50) as On the other hand taking the linear combination of the TBA equation as n Substituting this into the first term on the RHS of (58) and using σ (a) i from (49) we find where (55) is used and l ≫ s is assumed in the second equality. Now we resort to the general relation for 1 ≤ a ≤ n. In view of (57) it is convenient to take the linear combination of this as follows: Substituting (60) we arrive at the equation of state of the system: The LHS is an explicit rational function of z 1 , . . . , z n that can be calculated from (13) and (14). The variables z 1 , . . . , z n are simply related to the chemical potentials β 1 , . . . , β n or equivalently to the fugacity e −β1 , . . . , e −βn by (57) as Thus (63) relates the densities ν 1 , . . . , ν n with the fugacity e −β1 , . . . , e −βn , thereby enabling us to control either one by the other. Set where y 1 , . . . , y n are the fugacity mentioned just above since α a = β a according to (57). Then the equation of state (63) also admits a somewhat simpler presentation as To see this, note from y a = n b=1 z −C ab b and (65) that (63) is rewritten as In the language of the BBS, the relation ε where # a (b) for a ∈ [1, n + 1] denotes the number of the letter a in b ∈ B (r) s regarded as a semistandard tableau of shape (s r ). The empty space corresponds to the letter 1. Note that the weight wt(b) specifies an element b ∈ B The quantity in the parenthesis in the RHS is ρ In this way we have characterized the vertically 1/L-scaled equilibrium shape of the Young diagrams µ 1 , . . . , µ n under the prescribed densities (ν 1 , . . . , ν n ) = lim L→∞ L −1 (|µ 1 |, . . . , |µ n |) in terms of the variables ε According to (20) this is equal to the lim L→∞ i . On the other hand (for the highest state ensemble E + L (r, s) with L → ∞), the soliton/string correspondence (8) and the definition (48) indicate that the same quantity should also show up in the TBA analysis exactly as ε (a) i . Thus the asymptotic equivalence of the two ensembles indicates that they coincide. The next theorem, which is our main result in the paper, identifies them rigorously. Theorem 5.1 shows a nontrivial coincidence of the equilibrium configuration, i.e., the n-tuple of Young diagrams, determined by two different approaches: • stationary local energy for the Markov process of carriers in the randomized BBS, • difference equation arising from the TBA analysis of the Fermionic formula.
The result may be regarded as randomized version of the soliton/string correspondence (8). Being able to give an explicit formula for ε (a) i is a very rare event in the actual TBA analyses involving the spectral parameter.
For simplicity we temporarily write the w-deformed character (24) as In this notation, Lemma 4.4 reads as Then Theorem 5.1 is summarized in the following formula for the 1/L-scaled Young diagrams: From this and (70) the quantity η C ab H(x Here and in what follows x s , respectively. By removing the denominators, this is cast into In the derivation, we have used the Q-system (40) to cancel a factor Q (a) i in the first term of the RHS. In order to verify (77) we consider the two special cases of (39): where κ = i + δ a,r min(i, s) and the product over b means the one by * . Take the w-derivative of (78) at w = 1. By means of Lemma 4.4 or equivalently (73), it leads to The same calculation for (79) tells that the quantity in the big parenthesis of the last line of (80) is equal to i b∼a Q (b) i at w = 1. Therefore this term cancels the κ term on the second line of (80) partially. The resulting relation is nothing but (77).
Next we verify the boundary condition h Here the symmetry under the exchange of the indices is due to the invariance of weights and the local energy H by the combinatorial R.
Exchange the indices (a, i) ↔ (r, s) here and apply the symmetry h Proof. The Q-system (40) becomes (Q In the sequel we prove (85). From the w-derivative of (28) at w = 1, the LHS of (85) with fixed i(≥ s) is expressed as In the regime n b=1 z C ab b > 1 under consideration, the variables w a in (15) satisfy w 1 > w 2 > · · · > w n+1 . Then the large i limit of each summand in RHS of (86) is easily extracted from the determinantal formula in (13). It is decomposed into a product of Schur polynomials, which leads to LHS of (85) = z −s We have set w a = w a /z 1 = y 1 y 2 · · · y a−1 . See (15) and (65). As for RHS of (85), we invoke the formula for Q (1) s as the sum over semistandard tableaux on the Young diagram with length s single row shape. The entry b ∈ [1, n + 1] of the tableaux corresponds to w b . Therefore for any a ∈ [1, n] we have Since the first factor is free from y a whereas the latter contains it as the overall multiplier y k a , the derivative y a ∂Q (1) s ∂ya coincides with (87). This completes the proof of (85) hence that of Proposition 5.3. We have finished the proof of Theorem 5.1.

Example
In this section we focus on the simplest choice B  corresponding to the semistandard tableau containing a in the single box Young diagram. So p 1 is the density of empty sites and p a with a ∈ [2, n + 1] is the density of balls with color a. According to (16), one has π (1) 1 (a) = e ̟1−α1−···−αa−1 /Q (1) 1 . Therefore in the regime α 1 , . . . , α n > 0 under consideration, 1 > p 1 > p 2 > · · · > p n+1 > 0 holds. Of course p 1 + · · · + p n+1 = 1 should also be satisfied. For n = 2 this notation agrees with Example 3.1. According to (16) we set in terms of w j = z −1 j−1 z j given in (15). The denominator is Q 1 (14). Thus we find (cf. [23]) From (68), the ball densities p 1 , . . . , p n+1 are connected to the Young diagram densities ν 1 , . . . , ν n as The equation of state (63) reads where z 0 = z n+1 = 1 as in (15). One can easily check that (92) is satisfied by z a and ν a in (90) and (91) provided that p 1 + · · · + p n+1 = 1 is valid. This essentially achieves the step (i) in Section 5.4. For the remaining steps (ii) and (iii), we have already given the general solution in Theorem 5.1. In the present case the solution ε can be written down concretely by setting s = 1 in Example 3.3: For simplicity denote the Schur polynomial s λ (w 1 , . . . , w n+1 ) by s λ . Then the quantity (75) is given neatly as where we have used a bilinear identity among the Schur polynomials.
In the simplest case n = 1, the equation of state (92) becomes and s (i,i) (w 1 , w 2 ) = 1, the result (94) reduces to This agrees with a corresponding result in [23].
Applying this to (94) and using (90) we find that η (a) i tends to 0 as i → ∞ as up to exponentially small corrections. Thus the estimate (96) implies the logarithmic scaling in the leading order. For n = a = 1 and 1−p 1 = p 2 = p, (97) and (96) Therefore the estimate (96) gives This square root scaling behavior is a signal of criticality as observed in [29].
6.3. Numerical check. Here we deal with the n = 2 case, i.e., 2-color BBS. The relevant KR crystals are B Table 1. Notations and equation numbers except in the first column are those in [27].
We have generated a BBS state in {1, 2, 3} L with a prescribed ball densities 1 > p 1 > p 2 > p 3 > 0 and length L = 1000 by computer. Calculating the energy E    According to (96) we have truncated the scaled µ 1 and µ 2 at the width I 1 = 17 and I 2 = 16. The agreement of the numerical data from BBS and the TBA prediction is more or less satisfactory. 7. Discussion 7.1. Summary. We have elucidated a new interplay among the randomized BBS, Markov processes of carriers, KR modules/crystals, combinatorial R, local energy, deformed characters, Fermionic formulas, rigged configurations, Q and Y-systems, TBA equations and so forth. Our main result is Theorem 5.1 which identifies the stationary local energy of the KR crystal (71) as the explicit solution to the difference equation (69) originating from TBA. It determines the equilibrium shape of the Young diagrams µ 1 , . . . , µ n in the scaling limit as in (70), (74) and (75). These random Young diagrams arise as the conserved quantities (generalized soliton contents) of the randomized BBS and obey the probability distribution given by the Fermionic form (43). 7.2. Generalization to simply-laced case. Although the above results are concerned with the quantum affine algebra U q (ĝ) withĝ = A (1) n , all the essential ingredients are known or at least conjecturally/conceptually ready for general quantum affine algebras. In particular formulas for the simplylaced casesĝ = A  . Although it is no longer a character of an irreducible g module in general, it satisfies the Q-system (40). See for example [10] and [24, sec.13] and references therein.
The U q (ĝ) BBS is formulated in the same manner as Section 2. Take the set of local states to be B (r) s . We consider the randomized U q (ĝ) BBS where the local states and the stationary measure of the carrier for the time evolution T Concerning the deformed character (24), the corner transfer matrix energy D in (25) needs to be replaced by where b ♮ j ∈ B (kj ) lj is the unique element such that ϕ(b ♮ j ) = l j Λ 0 . See [32, sec.5.1] for a detailed account of this. The first sum on the RHS is referred to as the boundary energy. It is 0 for A (1) n but is nontrivial for the other types.
The Fermionic form M (B, λ, w) is defined by the same formulas as (30)- (35). Theorem 4.2 11 is valid for D (1) n [34,31] and conjecturally valid for E (1) 6,7,8 . Proposition 4.3 has been shown in [10]. Lemma 4.4 is influenced by the boundary energy and replaced by ∂ ∂w log χ w (B 1 ⊗ · · · ⊗ B L )| w=1 As for the TBA analysis, all the relations from (43) until (72) remain unchanged 12 . In particular, the property lim l→∞ (Q We conjecture that Theorem 5.1 is also valid for type D n and E 6,7,8 . In fact admitting Theorem 4.2, it can be shown that ε  13 . The proof uses 11 For general affine Lie algebra it is often called X = M conjecture [10,11] . 12 The unique exception is the last expression in (64) which is specific to type An. 13 This assertion is the analogue of Proposition 5.2, which was the "first half" of Theorem 5.1.
(102). In particular with the notation (72), its L = 2 case captures the stationary local energy h The point is that the effect of extra "boundary terms" containing b ♮ 1 ∈ B (a) is canceled by those in (102), leaving the difference equation unchanged from (69). As for the boundary condition for the difference equation, we conjecture that (81) or equivalently (85) holds universally for type D n and E 6,7,8 . It is an intriguing relation involving the local energy whose proof will shed new light into the KR crystals and the Q-system. 7.3. Further outlook. We expect the generalization to the non simply-laced cases and twisted affine Lie algebras is also feasible albeit with a slight technical complexity. Another obvious direction of a future research is periodic systems. The generalized BBS for A (1) n with the quantum space (B (1) 1 ) ⊗L has been studied under the periodic boundary condition [28]. It also has an n-tuple of Young diagrams as a label of iso-level sets for which a Fermionic formula [28, eq.(57), Th.3] for the multiplicity has been obtained under a technical assumption. It will be interesting to analyze it by TBA similarly to this paper. In the simplest case n = 1, the Fermionic formula has been fully justified and reduces to for system size L and M -ball sector with M < L 2 in the same notation as (31) 14 . So at least in this simplest situation, the scaled limit shape of the Young diagram remains the same as (95).
There are a number of further challenging problems to be investigated. We list a few of them as closing remarks.
(i) Study the limit shape problem when the BBS states are inhomogeneous as B (r1) s1 ⊗ · · · ⊗ B (rL) sL with a given statistical distribution of (r i , s i ). (ii) Can one architect a BBS like dynamical system whose Markov process of carriers has the stationary measure described by q-characters [6]? (iii) Can one extend the TBA analysis so as to include w-binomials in (43) with w = 1? What is the counterpart of the BBS corresponding to such a generalization? (iv) Our TBA analysis in this paper was spectral parameter free. See the remark after (49). Is there any Yang-Baxterization of Theorem 5.1?