# Generalized hydrodynamics in complete box-ball system
for $U_q(\hat{s}l_n)$

Atsuo Kuniba, Grégoire Misguich, Vincent Pasquier

► To cite this version:

Atsuo Kuniba, Grégoire Misguich, Vincent Pasquier. Generalized hydrodynamics in complete box-ball system for $U_q(\hat{s}l_n)$. SciPost Phys., 2021, 10 (4), pp.095. 10.21468/SciPostPhys.10.4.095. hal-03047646

HAL Id: hal-03047646

https://hal.science/hal-03047646

Submitted on 18 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Generalized hydrodynamics in complete
box-ball system for $U_q(\hat{sl}_n)$

Atsuo Kuniba¹*, Grégoire Misguich²† and Vincent Pasquier²‡

1 Institute of Physics, University of Tokyo, Komaba, Tokyo 153-8902, Japan

2 Institut de Physique Théorique, Université Paris Saclay, CEA,
CNRS UMR 3681, 91191 Gif-sur-Yvette, France

★ atsuo.s.kuniba@gmail.com, † gregoire.misguich@ipht.fr,

‡ vincent.pasquier@ipht.fr

Abstract

We introduce the complete box-ball system (cBBS), which is an integrable cellular automaton on 1D lattice associated with the quantum group $U_q(\hat{sl}_n)$. Compared with the conventional ($n-1$)-color BBS, it enjoys a remarkable simplification that scattering of solitons is totally diagonal. We also submit the cBBS to randomized initial conditions 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.

<figure><img src="image_1.png" /></figure>

Copyright A. Kuniba et al.

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 17-11-2020

Accepted 12-04-2021

Published 30-04-2021

doi:10.21468/SciPostPhys.10.4.095

<figure><img src="image_2.png" /></figure>

Check for
updates

Contents

1 **Introduction** 2

2 **Complete box-ball system** 5

2.1 Preliminaries 5

2.2 Combinatorial $R$ 5

2.3 Vacuum state and stability 8

2.4 Time evolutions and conserved quantities 9

2.5 Solitons and their scattering 11

2.6 Inverse scattering method 15

3 **Randomized cBBS** 16

3.1 Generalized Gibbs ensemble 16

3.2 I.I.D. randomness 17

3.3 Thermodynamic Bethe ansatz 19

4 **Generalized hydrodynamics** 23

4.1 Effective speed in homogenous system 23

4.2 Matrix form 26
Generalized hydrodynamics in complete
box-ball system for $U_q(\hat{sl}_n)$

Atsuo Kuniba¹*, Grégoire Misguich²† and Vincent Pasquier²‡

1 Institute of Physics, University of Tokyo, Komaba, Tokyo 153-8902, Japan

2 Institut de Physique Théorique, Université Paris Saclay, CEA,
CNRS UMR 3681, 91191 Gif-sur-Yvette, France

★ atsuo.s.kuniba@gmail.com, † gregoire.misguich@ipht.fr,

‡ vincent.pasquier@ipht.fr

Abstract

We introduce the complete box-ball system (cBBS), which is an integrable cellular automaton on 1D lattice associated with the quantum group $U_q(\hat{sl}_n)$. Compared with the conventional ($n-1$)-color BBS, it enjoys a remarkable simplification that scattering of solitons is totally diagonal. We also submit the cBBS to randomized initial conditions 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.

<figure><img src="image_2.png" /></figure>

Copyright A. Kuniba et al.

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 17-11-2020

Accepted 12-04-2021

Published 30-04-2021

doi:10.21468/SciPostPhys.10.4.095

<figure><img src="image_3.png" /></figure>

Check for
updates

Contents

1 **Introduction** 2

2 **Complete box-ball system** 5

2.1 Preliminaries 5

2.2 Combinatorial $R$ 5

2.3 Vacuum state and stability 8

2.4 Time evolutions and conserved quantities 9

2.5 Solitons and their scattering 11

2.6 Inverse scattering method 15

3 **Randomized cBBS** 16

3.1 Generalized Gibbs ensemble 16

3.2 I.I.D. randomness 17

3.3 Thermodynamic Bethe ansatz 19

4 **Generalized hydrodynamics** 23

4.1 Effective speed in homogenous system 23

4.2 Matrix form 26
<table><tbody><tr><td>4.3</td><td>Stationary currents</td><td>27</td></tr><tr><td>4.4</td><td>Dynamics in an inhomogeneous system</td><td>27</td></tr><tr><td>5</td><td>Domain wall initial condition</td><td>28</td></tr><tr><td>5.1</td><td>Ballistic picture</td><td>28</td></tr><tr><td>5.2</td><td>Diffusive correction</td><td>31</td></tr><tr><td>5.3</td><td>Simulations</td><td>32</td></tr><tr><td>6</td><td>Summary and conclusions</td><td>40</td></tr><tr><td>A</td><td>Algorithm for the combinatorial R</td><td>42</td></tr><tr><td>B</td><td>Combinatorial R and energy H for n = 2, 3</td><td>43</td></tr><tr><td>C</td><td>KSS bijection</td><td>44</td></tr><tr><td>C.1</td><td>Paths</td><td>44</td></tr><tr><td>C.2</td><td>Rigged configurations</td><td>45</td></tr><tr><td>C.3</td><td>KSS bijection</td><td>47</td></tr><tr><td>References</td><td></td><td>50</td></tr></tbody></table>

# 1 Introduction

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(\hat{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 consider a protocol where the cBBS is initially prepared in some random initial conditions corresponding to two different ball densities in the left and in the right halves. The associated non-equilibrium behavior is then studied by thermodynamic Bethe ansatz [6, 7] and generalized hydrodynamics [8–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 \times l$ rectangular Young diagram over the alphabet $\{1, 2, \dots, 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 \le k < n$ is relevant and the simplest $B^{1,1}$ is identified with the set $\{1, 2, \dots, n\}$. A typical or conventionally most studied BBS with $(n-1)$ kinds of balls [14] is a dynamical system on

$$ \dots \otimes B^{1,1} \otimes B^{1,1} \otimes B^{1,1} \otimes \dots, \qquad (1.1) $$

where $\otimes$ 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 \in B^{1,1}$ or contains a color a ball for $a = 2, 3, \dots, n \in B^{1,1}$. All the distant boxes are assumed to be 1 = empty. By using the crystal theory of $U_q(\hat{sl}_n)$, one can formulate an integrable dynamics on such states yielding solitons [15, 16].<fnref n="1" /> Here is

<footnote n="1">¹In this paper we will not use the crystal theory of $U_q(\hat{sl}_n)$ extensively since it is not the main theme, and only</footnote>
<table><tbody><tr><td>4.3</td><td>Stationary currents</td><td>27</td></tr><tr><td>4.4</td><td>Dynamics in an inhomogeneous system</td><td>27</td></tr><tr><td>5</td><td>Domain wall initial condition</td><td>28</td></tr><tr><td>5.1</td><td>Ballistic picture</td><td>28</td></tr><tr><td>5.2</td><td>Diffusive correction</td><td>31</td></tr><tr><td>5.3</td><td>Simulations</td><td>32</td></tr><tr><td>6</td><td>Summary and conclusions</td><td>40</td></tr><tr><td>A</td><td>Algorithm for the combinatorial R</td><td>42</td></tr><tr><td>B</td><td>Combinatorial R and energy H for n = 2, 3</td><td>43</td></tr><tr><td>C</td><td>KSS bijection</td><td>44</td></tr><tr><td>C.1</td><td>Paths</td><td>44</td></tr><tr><td>C.2</td><td>Rigged configurations</td><td>45</td></tr><tr><td>C.3</td><td>KSS bijection</td><td>47</td></tr><tr><td>References</td><td></td><td>50</td></tr></tbody></table>

# 1 Introduction

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(\hat{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 consider a protocol where the cBBS is initially prepared in some random initial conditions corresponding to two different ball densities in the left and in the right halves. The associated non-equilibrium behavior is then studied by thermodynamic Bethe ansatz [6, 7] and generalized hydrodynamics [8–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 \times l$ rectangular Young diagram over the alphabet $\{1, 2, \dots, 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 \le k < n$ is relevant and the simplest $B^{1,1}$ is identified with the set $\{1, 2, \dots, n\}$. A typical or conventionally most studied BBS with $(n-1)$ kinds of balls [14] is a dynamical system on

$$ \dots \otimes B^{1,1} \otimes B^{1,1} \otimes B^{1,1} \otimes \dots, \qquad (1.1) $$

where $\otimes$ 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 \in B^{1,1}$ or contains a color $a$ ball for $a = 2, 3, \dots, n \in B^{1,1}$. All the distant boxes are assumed to be 1 = empty. By using the crystal theory of $U_q(\hat{sl}_n)$, one can formulate an integrable dynamics on such states yielding solitons [15, 16].<fnref n="1" /> Here is

<footnote n="1">¹In this paper we will not use the crystal theory of $U_q(\hat{sl}_n)$ extensively since it is not the main theme, and only</footnote>
Third, the complete set of conserved quantities of cBBS are given as an $(n-1)$-tuple of Young diagrams $\mu^{(1)}, \dots, \mu^{(n-1)}$, and they all admit most transparent meaning; $\mu^{(a)}$ is nothing but the list of amplitude $i$ of the color $a$ solitons $S_i^{(a)}$. Such a direct interpretation of the conserved Young diagrams in terms of solitons was possible only for the first one $\mu^{(1)}$ in the conventional BBS. Our cBBS puts all of $\mu^{(1)}, \dots, \mu^{(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(\mathfrak{sl}_n)$ differing from the conventional ones for $n \ge 3$. The nomenclature “complete” BBS is meant to indicate that the complete list $B^{1,1}, B^{2,1}, \dots, B^{n-1,1}$ of (the labeling set of the basis of) the fundamental representations of $\mathfrak{sl}_n$ have been gathered into $\mathcal{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 $\mathcal{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]. Randomized means here that some measures over random initial conditions are considered, but we stress that for a given initial microscopic configuration the dynamics remains completely deterministic. 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, \dots, 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, \dots, 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 non-equilibrium 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.

Let us comment on some other extensions of the basic setting (1.1) than (1.2) in the literature. The best known example is $\cdots \otimes B^{1,l_i-1} \otimes B^{1,l_i} \otimes B^{1,l_{i+1}} \otimes \cdots$, which is the BBS with boxes having (possibly inhomogeneous) general capacities [15]. The case $\cdots \otimes B^{k,1} \otimes B^{k,1} \otimes B^{k,1} \otimes \cdots$ has also been investigated in [20]. The closest model to our cBBS is $\cdots \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes \cdots$, which yields, with an elaborate decoration at the boundary, the BBS with reflecting end [21]. Although the time evolution is designed differently there, the state space of the bulk part coincides with that of cBBS for the smallest choice $n=3$. In all these examples except the lowest rank situation, scatterings of solitons are non-diagonal, which testifies the novelty of cBBS. As for the recent progress on probabilistic and statistical aspects of the randomized BBS (non-complete BBS except $n=2$), see also [17, 22–27].

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.
Third, the complete set of conserved quantities of cBBS are given as an $(n-1)$-tuple of Young diagrams $\mu^{(1)}, \dots, \mu^{(n-1)}$, and they all admit most transparent meaning; $\mu^{(a)}$ is nothing but the list of amplitude $i$ of the color $a$ solitons $S_i^{(a)}$. Such a direct interpretation of the conserved Young diagrams in terms of solitons was possible only for the first one $\mu^{(1)}$ in the conventional BBS. Our cBBS puts all of $\mu^{(1)}, \dots, \mu^{(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(\mathfrak{sl}_n)$ differing from the conventional ones for $n \ge 3$. The nomenclature “complete” BBS is meant to indicate that the complete list $B^{1,1}, B^{2,1}, \dots, B^{n-1,1}$ of (the labeling set of the basis of) the fundamental representations of $\mathfrak{sl}_n$ have been gathered into $\mathcal{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 $\mathcal{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]. Randomized means here that some measures over random initial conditions are considered, but we stress that for a given initial microscopic configuration the dynamics remains completely deterministic. 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, \dots, 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, \dots, 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 non-equilibrium 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.

Let us comment on some other extensions of the basic setting (1.1) than (1.2) in the literature. The best known example is $\cdots \otimes B^{1,l_i-1} \otimes B^{1,l_i} \otimes B^{1,l_{i+1}} \otimes \cdots$, which is the BBS with boxes having (possibly inhomogeneous) general capacities [15]. The case $\cdots \otimes B^{k,1} \otimes B^{k,1} \otimes B^{k,1} \otimes \cdots$ has also been investigated in [20]. The closest model to our cBBS is $\cdots \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes (B^{1,1} \otimes B^{n-1,1}) \otimes \cdots$, which yields, with an elaborate decoration at the boundary, the BBS with reflecting end [21]. Although the time evolution is designed differently there, the state space of the bulk part coincides with that of cBBS for the smallest choice $n=3$. In all these examples except the lowest rank situation, scatterings of solitons are non-diagonal, which testifies the novelty of cBBS. As for the recent progress on probabilistic and statistical aspects of the randomized BBS (non-complete BBS except $n=2$), see also [17, 22–27].

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_\infty^{(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.

## 2 Complete box-ball system

Throughout the paper we use the notation

$$ [a, b] = \{a, a + 1, \dots, b\} \quad (a \le b \in \mathbb{Z}), \qquad \mathcal{I} = \{(a, i) \mid a \in [1, n-1], i \in \mathbb{Z}_{\ge 1}\}. \quad (2.1) $$

### 2.1 Preliminaries

Consider the classical simple Lie algebra $sl_n$ ($n \ge 2$). We denote its Cartan matrix by $(C_{ab})_{a,b=1}^{n-1}$, where

$$ C_{ab} = 2\delta_{a,b} - \delta_{|a-b|,1}. \qquad (2.2) $$

Let $\varpi_1, \dots, \varpi_n$ be the fundamental weights and $\alpha_1, \dots, \alpha_n$ be the simple roots. They are related by $\alpha_a = \sum_{b=1}^n C_{ab}\varpi_b$. Let $\widehat{sl}_n$ be the non-twisted affinization of $sl_n$ [28] and $U_q = U_q(\widehat{sl}_n)$ be the quantum affine algebra (without derivation operator) [29, 30]. There is a family of irreducible finite-dimensional representations $\{W_l^{(k)} | (k, l) \in [1, n-1] \times \mathbb{Z}_{\ge 0}\}$ of $U_q$ called Kirillov-Reshetikhin (KR) module<fnref n="2" /> named after the related work on the Yangian [31]. As a representation of $U_q(sl_n)$, $W_l^{(k)}$ is isomorphic to the type 1 irreducible highest weight module $V(l\varpi_k)$ with highest weight $l\varpi_k$. $W_l^{(k)}$ is known to have a crystal base $B^{k,l}$ [4, 32]. Roughly speaking, it is a set of basis vectors of the $W_l^{(k)}$ 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.

### 2.2 Combinatorial $R$

For $(k, l) \in \mathcal{I}$, let $B^{k,l}$ be the set of semistandard tableaux of $k \times l$ rectangular shape over the alphabet $\{1, 2, \dots, 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} < \dots < t_{kj}$ and $t_{i1} \le t_{i2} \le \dots \le t_{il}$ hold for any $i \in [1, k]$ and $j \in [1, l]$. The array row $(b) = t_k \dots t_2 t_1$ with $t_i = t_{i1} t_{i2} \dots t_{il}$ is

<footnote n="2">²The actual KR modules carry a spectral parameter. In this paper it is irrelevant and hence suppressed.</footnote>
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_{\infty}^{(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.

## 2 Complete box-ball system

Throughout the paper we use the notation

$$ [a, b] = \{a, a + 1, \dots, b\} \quad (a \le b \in \mathbb{Z}), \qquad \mathcal{I} = \{(a, i) \mid a \in [1, n-1], i \in \mathbb{Z}_{\ge 1}\}. \quad (2.1) $$

### 2.1 Preliminaries

Consider the classical simple Lie algebra $sl_n$ ($n \ge 2$). We denote its Cartan matrix by $(C_{ab})_{a,b=1}^{n-1}$, where

$$ C_{ab} = 2\delta_{a,b} - \delta_{|a-b|,1}. \qquad (2.2) $$

Let $\varpi_1, \dots, \varpi_n$ be the fundamental weights and $\alpha_1, \dots, \alpha_n$ be the simple roots. They are related by $\alpha_a = \sum_{b=1}^n C_{ab}\varpi_b$. Let $\hat{sl}_n$ be the non-twisted affinization of $sl_n$ [28] and $U_q = U_q(\hat{sl}_n)$ be the quantum affine algebra (without derivation operator) [29, 30]. There is a family of irreducible finite-dimensional representations $\{W_l^{(k)} | (k, l) \in [1, n-1] \times \mathbb{Z}_{\ge 0}\}$ of $U_q$ called Kirillov-Reshetikhin (KR) module<fnref n="2" /> named after the related work on the Yangian [31]. As a representation of $U_q(sl_n)$, $W_l^{(k)}$ is isomorphic to the type 1 irreducible highest weight module $V(l\varpi_k)$ with highest weight $l\varpi_k$. $W_l^{(k)}$ is known to have a crystal base $B^{k,l}$ [4, 32]. Roughly speaking, it is a set of basis vectors of the $W_l^{(k)}$ 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.

### 2.2 Combinatorial $R$

For $(k, l) \in \mathcal{I}$, let $B^{k,l}$ be the set of semistandard tableaux of $k \times l$ rectangular shape over the alphabet $\{1, 2, \dots, 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} < \dots < t_{kj}$ and $t_{i1} \le t_{i2} \le \dots \le t_{il}$ hold for any $i \in [1, k]$ and $j \in [1, l]$. The array row($b$) = $t_k \dots t_2 t_1$ with $t_i = t_{i1} t_{i2} \dots t_{il}$ is

<footnote n="2">²The actual KR modules carry a spectral parameter. In this paper it is irrelevant and hence suppressed.</footnote>
called the row word of b. For instance, we have

$$b = \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \end{bmatrix} \in B^{2,3}, \quad \text{row}(b) = 245123, \quad c = \begin{bmatrix} 1 & 2 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \in B^{3,2}, \quad \text{row}(c) = 452312. \quad (2.3)$$

For tableaux S and T, their product is defined in the following two ways which are known to be equivalent:

$$S \cdot T = (\cdots ((S \leftarrow u_1) \leftarrow u_2) \leftarrow \cdots) \leftarrow u_l \quad (\text{row}(T) = u_1 u_2 \cdots u_l) \quad (2.4)$$

$$= v_1 \rightarrow (v_2 \rightarrow (\cdots (v_m \rightarrow T) \cdots)) \quad (\text{row}(S) = v_1 v_2 \cdots v_m). \quad (2.5)$$

Here $\leftarrow$ denotes the row insertion and $\rightarrow$ does the column insertion [12].

### Row insertion $S \leftarrow 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 \rightarrow 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 \le 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 \cdot c = \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 3 \\ 3 & 4 \\ 4 & 5 \end{bmatrix}, \quad c \cdot b = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 \\ 3 & 5 \\ 4 \end{bmatrix}. \quad (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

$$R: B^{k,l} \otimes B^{k',l'} \rightarrow B^{k',l'} \otimes B^{k,l}; \quad b \otimes c \mapsto \tilde{c} \otimes \tilde{b}, \quad (2.7)$$

whose image $\tilde{c} \otimes \tilde{b} = R(b \otimes c)$ is characterized by the condition $c \cdot b = \tilde{b} \cdot \tilde{c}$. Here and in what follows, $\otimes$ is used to mean just the product of sets. It should not be confused with the product $\cdot$ of tableaux. Note that the dependence on $k, k', l, l'$ has been suppressed in $R$. For example, from

$$\begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix} \cdot \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 \\ 3 & 5 \\ 4 \end{bmatrix}, \quad (2.8)$$
called the row word of b. For instance, we have

$$b = \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \end{bmatrix} \in B^{2,3}, \quad \text{row}(b) = 245123, \quad c = \begin{bmatrix} 1 & 2 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \in B^{3,2}, \quad \text{row}(c) = 452312. \quad (2.3)$$

For tableaux S and T, their product is defined in the following two ways which are known to be equivalent:

$$S \cdot T = (\cdots ((S \leftarrow u_1) \leftarrow u_2) \leftarrow \cdots) \leftarrow u_l \quad (\text{row}(T) = u_1 u_2 \cdots u_l) \quad (2.4)$$

$$= v_1 \rightarrow (v_2 \rightarrow (\cdots (v_m \rightarrow T) \cdots)) \quad (\text{row}(S) = v_1 v_2 \cdots v_m). \quad (2.5)$$

Here $\leftarrow$ denotes the row insertion and $\rightarrow$ does the column insertion [12].

### Row insertion $S \leftarrow 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 \rightarrow 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 \le 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 \cdot c = \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 3 \\ 3 & 4 \\ 4 & 5 \end{bmatrix}, \quad c \cdot b = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 \\ 3 & 5 \\ 4 \end{bmatrix}. \quad (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

$$R: B^{k,l} \otimes B^{k',l'} \rightarrow B^{k',l'} \otimes B^{k,l}; \quad b \otimes c \mapsto \tilde{c} \otimes \tilde{b}, \quad (2.7)$$

whose image $\tilde{c} \otimes \tilde{b} = R(b \otimes c)$ is characterized by the condition $c \cdot b = \tilde{b} \cdot \tilde{c}$. Here and in what follows, $\otimes$ is used to mean just the product of sets. It should not be confused with the product $\cdot$ of tableaux. Note that the dependence on $k, k', l, l'$ has been suppressed in $R$. For example, from

$$\begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix} \cdot \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 \\ 3 & 5 \\ 4 \end{bmatrix}, \quad (2.8)$$
which agrees with the right tableau in (2.6), we find

$$R: \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \end{bmatrix} \otimes \begin{bmatrix} 1 & 2 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \to \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix} \otimes \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix}. \qquad (2.9)$$

This is a particular case of $R: B^{2,3} \otimes B^{3,2} \to B^{3,2} \otimes B^{2,3}$. The algorithm to find $\tilde{b}$ and $\tilde{c}$ satisfying the condition $\tilde{b} \cdot \tilde{c} = c \cdot b$ is described in Appendix A following [35, p55].

The local energy $H$ is the function defined by

$$H: B^{k,l} \otimes B^{k',l'} \to \mathbb{Z}_{\ge 0}; \quad b \otimes c \mapsto H(b \otimes c), \qquad (2.10)$$

$$H(b \otimes c) = \text{number of boxes strictly below the } \max(k, k')\text{-th row of the tableau } c \cdot b. \qquad (2.11)$$

In our working example (2.3), $\max(k, k') = \max(2, 3) = 3$, hence from (2.6) we have $H(b \otimes c) = 1$.

Let $u_{k,l} \in 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

$$u_{1,1} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \quad u_{2,1} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \quad u_{3,1} = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}, \quad u_{1,2} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad u_{2,2} = \begin{bmatrix} 1 \\ 1 \\ 2 \end{bmatrix}, \quad u_{3,2} = \begin{bmatrix} 1 \\ 1 \\ 2 \\ 3 \end{bmatrix}. \qquad (2.12)$$

It is easy to see

$$R(u_{k,l} \otimes u_{k',l'}) = u_{k',l'} \otimes u_{k,l}, \qquad (2.13)$$

$$H(b \otimes u_{k',l'}) = 0 \quad (\forall b \in B^{k,l}) \qquad (2.14)$$

for any $k, k', l, l'$.

If $R$ in (2.7) and $H$ in (2.10) are denoted by $R_{B^{k,l} \otimes B^{k',l'}}$ and $H_{B^{k,l} \otimes B^{k',l'}}$ respectively, they satisfy

$$R_{B^{k,l} \otimes B^{k',l'}} = \mathrm{id}_{B^{k,l} \otimes B^{k',l}}, \qquad (2.15)$$

$$R_{B^{k,l} \otimes B^{k',l'}} R_{B^{k',l'} \otimes B^{k,l}} = \mathrm{id}_{B^{k',l'} \otimes B^{k,l}}, \qquad (2.16)$$

$$H_{B^{k',l'} \otimes B^{k,l}} \circ R_{B^{k,l} \otimes B^{k',l'}} = H_{B^{k,l} \otimes B^{k',l'}} \qquad (2.17)$$

by the definition. In spite of the simplification (2.15), the local energy $H_{B^{k,l} \otimes B^{k',l}}(b \otimes c)$ still depends on $b \otimes c$ nontrivially. The relation in (2.16) is called the inversion relation, which implies that $R(b \otimes c) = \tilde{c} \otimes \tilde{b}$ is equivalent to $R(\tilde{c} \otimes \tilde{b}) = b \otimes c$. In view of this we write these relations also as $b \otimes c \simeq \tilde{c} \otimes \tilde{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 $\mathrm{Aff}(B^{k,l}) = \{b[\alpha] \mid b \in B^{k,l}, \alpha \in \mathbb{Z}\}$ and extend $R$ to the map $\hat{R}$ defined by

$$\hat{R}: \mathrm{Aff}(B^{k,l}) \otimes \mathrm{Aff}(B^{k',l'}) \to \mathrm{Aff}(B^{k',l'}) \otimes \mathrm{Aff}(B^{k,l}) \\ b[\alpha] \otimes c[\beta] \quad \mapsto \quad \tilde{c}[\beta + H(b \otimes c)] \otimes \tilde{b}[\alpha - H(b \otimes c)], \qquad (2.18)$$

where $\tilde{c} \otimes \tilde{b} = R(b \otimes c)$ as in (2.7). The integer $\alpha$ attached to $b$ in $b[\alpha]$ is called a *mode*. It is a reminiscent of the spectral parameter in quantum $R$ matrices. Note that the shift of the mode $H(b \otimes c)$ in (2.18) is equally presented as $H(\tilde{c} \otimes \tilde{b})$ since they coincide owing to
which agrees with the right tableau in (2.6), we find

$$ R: \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \end{bmatrix} \otimes \begin{bmatrix} 1 & 2 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} \mapsto \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix} \otimes \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix}. \qquad (2.9) $$

This is a particular case of $R: B^{2,3} \otimes B^{3,2} \to B^{3,2} \otimes B^{2,3}$. The algorithm to find $\tilde{b}$ and $\tilde{c}$ satisfying the condition $\tilde{b} \cdot \tilde{c} = c \cdot b$ is described in Appendix A following [35, p55].

The local energy $H$ is the function defined by

$$ H: B^{k,l} \otimes B^{k',l'} \to \mathbb{Z}_{\ge 0}; \quad b \otimes c \mapsto H(b \otimes c), \qquad (2.10) $$

$$ H(b \otimes c) = \text{number of boxes strictly below the } \max(k, k')\text{-th row of the tableau } c \cdot b. \qquad (2.11) $$

In our working example (2.3), $\max(k, k') = \max(2,3) = 3$, hence from (2.6) we have $H(b \otimes c) = 1$.

Let $u_{k,l} \in 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

$$ u_{1,1} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad u_{2,1} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \quad u_{3,1} = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}, \quad u_{1,2} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad u_{2,2} = \begin{bmatrix} 1 \\ 1 \\ 2 \end{bmatrix}, \quad u_{3,2} = \begin{bmatrix} 1 \\ 1 \\ 2 \\ 3 \end{bmatrix}. \qquad (2.12) $$

It is easy to see

$$ R(u_{k,l} \otimes u_{k',l'}) = u_{k',l'} \otimes u_{k,l}, \qquad (2.13) $$

$$ H(b \otimes u_{k',l'}) = 0 \quad (\forall b \in B^{k,l}) \qquad (2.14) $$

for any $k, k', l, l'$.

If $R$ in (2.7) and $H$ in (2.10) are denoted by $R_{B^{k,l} \otimes B^{k',l'}}$ and $H_{B^{k,l} \otimes B^{k',l'}}$ respectively, they satisfy

$$ R_{B^{k,l} \otimes B^{k',l'}} = \mathrm{id}_{B^{k,l} \otimes B^{k',l}}, \qquad (2.15) $$

$$ R_{B^{k,l} \otimes B^{k',l'}} R_{B^{k',l'} \otimes B^{k,l}} = \mathrm{id}_{B^{k',l'} \otimes B^{k,l}}, \qquad (2.16) $$

$$ H_{B^{k',l'} \otimes B^{k,l}} \circ R_{B^{k,l} \otimes B^{k',l'}} = H_{B^{k,l} \otimes B^{k',l'}} \qquad (2.17) $$

by the definition. In spite of the simplification (2.15), the local energy $H_{B^{k,l} \otimes B^{k',l}}(b \otimes c)$ still depends on $b \otimes c$ nontrivially. The relation in (2.16) is called the inversion relation, which implies that $R(b \otimes c) = \tilde{c} \otimes \tilde{b}$ is equivalent to $R(\tilde{c} \otimes \tilde{b}) = b \otimes c$. In view of this we write these relations also as $b \otimes c \simeq \tilde{c} \otimes \tilde{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 $\mathrm{Aff}(B^{k,l}) = \{b[\alpha] \mid b \in B^{k,l}, \alpha \in \mathbb{Z}\}$ and extend $R$ to the map $\hat{R}$ defined by

$$ \begin{array}{l} \hat{R}: \mathrm{Aff}(B^{k,l}) \otimes \mathrm{Aff}(B^{k',l'}) \to \mathrm{Aff}(B^{k',l'}) \otimes \mathrm{Aff}(B^{k,l}) \\ b[\alpha] \otimes c[\beta] \quad \mapsto \quad \tilde{c}[\beta + H(b \otimes c)] \otimes \tilde{b}[\alpha - H(b \otimes c)], \end{array} \qquad (2.18) $$

where $\tilde{c} \otimes \tilde{b} = R(b \otimes c)$ as in (2.7). The integer $\alpha$ attached to $b$ in $b[\alpha]$ is called a *mode*. It is a reminiscent of the spectral parameter in quantum $R$ matrices. Note that the shift of the mode $H(b \otimes c)$ in (2.18) is equally presented as $H(\tilde{c} \otimes \tilde{b})$ since they coincide owing to
(2.17). This guarantees that $\hat{R}$ also satisfies the inversion relation. Thus we also write (2.18) as $b[\alpha] \otimes c[\beta] \simeq \tilde{c}[\beta + H(b \otimes c)] \otimes \tilde{b}[\alpha - H(b \otimes c)]$ putting the two sides on a more equal footing.

Now the Yang-Baxter equation is presented as

$$ (\hat{R} \otimes \text{id})(\text{id} \otimes \hat{R})(\hat{R} \otimes \text{id}) = (\text{id} \otimes \hat{R})(\hat{R} \otimes \text{id})(\text{id} \otimes \hat{R}), \quad (2.19) $$

which is an equality of the maps $\text{Aff}(B^{k,l}) \otimes \text{Aff}(B^{k',l'}) \otimes \text{Aff}(B^{k'',l''}) \to \text{Aff}(B^{k'',l''}) \otimes \text{Aff}(B^{k',l'}) \otimes \text{Aff}(B^{k,l})$ for arbitrary $(k, l), (k', l'), (k'', l'') \in \mathcal{I}$. 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.

<figure><img src="image_4.png" /></figure>

In general, we will write $b_1 \otimes \cdots \otimes b_L \simeq b'_1 \otimes \cdots \otimes b'_L$ if the elements $b_1 \otimes \cdots \otimes b_L$ and $b'_1 \otimes \cdots \otimes b'_L$ are transformed to each other by successively applying the combinatorial $R$ to the neighboring components as above.

The relation $b[\alpha] \otimes c[\beta] \simeq \tilde{c}[\beta + h] \otimes \tilde{b}[\alpha - h]$ (2.18) or its classical part (2.7) without the mode will be depicted as $(h = H(b \otimes c))$

$$ \begin{tikzcd} & c[\beta] \arrow[r] \arrow[d, "b[\alpha]"'] & \tilde{c}[\beta + h] \arrow[d, "\tilde{c}"] \\ b[\alpha] \arrow[r] & \tilde{b}[\alpha - h] & b \arrow[r] & \tilde{b} \end{tikzcd} \quad (2.20) $$

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.

## 2.3 Vacuum state and stability

Collecting all the $l = 1$ cases in (2.12) for a given $n$, we set

$$ \text{vac} = u_{1,1} \otimes u_{2,1} \otimes \cdots \otimes u_{n-1,1} \in B, \quad B = B^{1,1} \otimes B^{2,1} \otimes \cdots \otimes B^{n-1,1}. \quad (2.21) $$

For instance, vac looks as

$$ \boxed{1} (n=2), \quad \boxed{1} \otimes \boxed{1/2} (n=3), \quad \boxed{1} \otimes \boxed{1/2} \otimes \boxed{1/2} \otimes \boxed{3} (n=4). \quad (2.22) $$
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

$$B^{k,l} \otimes B^{\otimes L} \ni b \otimes \mathrm{vac}^{\otimes L} \simeq v_1 \otimes \cdots \otimes v_m \otimes \mathrm{vac}^{\otimes L-m} \otimes u_{k,l} \in B^{\otimes L} \otimes B^{k,l} \quad (\forall b \in B^{k,l}), \qquad (2.23)$$

when L gets sufficiently large for some m and $v_1, \dots, v_m \in 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 vacuum 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$:

<figure><img src="image_3.png" /></figure>

Here $b = \begin{bmatrix} 2 & 3 & 3 \end{bmatrix}$ and $L = 3$. The diagram is constructed by concatenating (2.20). Another example for $(k, l) = (2, 3)$ with $n = 3$ is

<figure><img src="image_4.png" /></figure>

In the (right) distance, all the combinatorial R's tend to the situation (2.13).

## 2.4 Time evolutions and conserved quantities

Now we define the complete BBS (cBBS). It is a dynamical system on $B^{\otimes L}$. An element $b_1 \otimes \cdots \otimes b_L \in B^{\otimes L}$ is called a state where each $b_i \in B$ is a local state at site i. Thus there are $\#B = \binom{n}{1}\binom{n}{2}\cdots\binom{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) \in \mathcal{I}$ in (2.1), we define the time evolution $T_l^{(r)}: B^{\otimes L} \to B^{\otimes L}$ as follows:

$$T_l^{(r)}(b_1 \otimes \cdots \otimes b_L) = b'_1 \otimes \cdots \otimes b'_L \quad (b_i \in B), \qquad (2.24)$$

$$u_{r,l} \otimes b_1 \otimes \cdots \otimes b_L \simeq b'_1 \otimes \cdots \otimes b'_L \otimes u_{r,l}. \qquad (2.25)$$

The relation (2.25) is obtained by successively applying the combinatorial R sending $B^{r,l}$ through $B^{\otimes L}$ to the right. Like the previous examples, the procedure can be shown graphically as

$$u_{r,l} \begin{array}{c} \begin{array}{ccccccc} b_1 & b_2 & & & & & b_L \\ \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow \\ b'_1 & b'_2 & & & & & b'_L \end{array} \end{array} \qquad (2.26)$$
Here the boundary condition implies $b_j = b_{j+1} = \cdots = b_L = \text{vac for } L \gg L-j \gg 1$. Then, thanks to the stability (2.23), we always end up with $u_{r,l}$ in the right. In short, $T_l^{(r)}$ 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_l^{(r)})^{-1} : b'_1 \otimes \cdots \otimes b'_L \mapsto b_1 \otimes \cdots \otimes b_L$ by going from SE to NW under the boundary condition $b_1 = b_2 = \cdots = b_m = \text{vac for } 1 \ll m \ll L$.

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_l^{(r)}$ 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} \otimes B^{2,1}$ can be displayed as

$$ \text{vac} = \begin{array}{c} \boxed{1} \otimes \boxed{\frac{1}{2}} = \boxed{\phantom{00}} \\ \boxed{2} \otimes \boxed{\frac{1}{2}} = \boxed{\phantom{00}} \\ \boxed{3} \otimes \boxed{\frac{1}{2}} = \boxed{\phantom{00}} \end{array}, \quad \begin{array}{c} \boxed{1} \otimes \boxed{\frac{1}{3}} = \boxed{\phantom{00}} \\ \boxed{2} \otimes \boxed{\frac{1}{3}} = \boxed{\phantom{00}} \\ \boxed{3} \otimes \boxed{\frac{1}{3}} = \boxed{\phantom{00}} \end{array}, \quad \begin{array}{c} \boxed{1} \otimes \boxed{\frac{2}{3}} = \boxed{\phantom{00}} \\ \boxed{2} \otimes \boxed{\frac{2}{3}} = \boxed{\phantom{00}} \\ \boxed{3} \otimes \boxed{\frac{2}{3}} = \boxed{\phantom{00}} \end{array}, \qquad (2.27) $$

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, \dots, 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_l^{(r)}$ is associated with the energy $E_l^{(r)} : B^{\otimes L} \to \mathbb{Z}_{\ge 0}$. Its quickest definition is to declare that the mode of the carrier in (2.26) is shifted as follows:

$$ u_{r,l}[\alpha] \rightarrow \begin{array}{c} b_1 \quad b_2 \\ \downarrow \quad \downarrow \\ b'_1 \quad b'_2 \end{array} \rightarrow \begin{array}{c} b_L \\ \downarrow \\ b'_L \end{array} \rightarrow u_{r,l}[\alpha - E_l^{(r)}(b_1 \otimes \cdots \otimes b_L)] \qquad (2.28) $$

To be more detailed, suppose the local states have the form $b_i = c_{(n-1)(i-1)+1} \otimes \cdots \otimes c_{(n-1)i} \in B = B^{1,1} \otimes \cdots \otimes B^{n-1,1}$ so that $b_1 \otimes \cdots \otimes b_L = c_1 \otimes \cdots \otimes c_{(n-1)L}$. The carriers $u_1, \dots, u_{(n-1)L} \in B^{r,l}$ in the intermediate steps are determined from the initial condition $u_0 = u_{r,l}$ and the recursion

$$ u_{j-1} \otimes c_j \simeq c'_j \otimes u_j \qquad (1 \le j \le (n-1)L), \qquad (2.29) $$

for some $c'_j$. From the defining behavior of the mode in (2.18), the energy $E_l^{(r)}$ is expressed as the sum of local energies attached to (2.29) as

$$ E_l^{(r)}(b_1 \otimes \cdots \otimes b_L) = \sum_{j=1}^{(n-1)L} H(u_{j-1} \otimes c_j) \in \mathbb{Z}_{\ge 0}. \qquad (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

$$ \text{commutativity : } T_l^{(k)} T_{l'}^{(k')} = T_{l'}^{(k')} T_l^{(k)}, \qquad (2.31) $$

$$ \text{conservation : } E_l^{(k)} T_{l'}^{(k')} = E_l^{(k)} \qquad (2.32) $$
on any state and for arbitrary $(k, l), (k', l') \in \mathcal{I}$. This is a simple consequence of the Yang-Baxter equation. In fact, write $\mathfrak{s} = b_1 \otimes \cdots \otimes b_L$ for short and consider

$$ 
\begin{aligned}
u_{k,l}[\alpha] \otimes u_{k',l'}[\beta] \otimes \mathfrak{s} & \approx u_{k,l}[\alpha] \otimes T_{l'}^{(k')}(\mathfrak{s}) \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(\mathfrak{s})] \\
& \approx T_{l}^{(k)} T_{l'}^{(k')}(\mathfrak{s}) \otimes u_{k,l}[\alpha - E_{l}^{(k)}(T_{l'}^{(k')}(\mathfrak{s}))] \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(\mathfrak{s})] \quad (2.33)
\end{aligned}
$$

by using (2.28) successively. On the other hand, one may first apply $u_{k,l}[\alpha] \otimes u_{k',l'}[\beta] \approx u_{k',l'}[\beta] \otimes u_{k,l}[\alpha]$ which is a consequence of (2.13) and (2.14). It tells that (2.33) is also equal to

$$ 
\begin{aligned}
T_{l'}^{(k')} T_{l}^{(k)}(\mathfrak{s}) \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(T_{l}^{(k)}(\mathfrak{s}))] \otimes u_{k,l}[\alpha - E_{l}^{(k)}(\mathfrak{s})] \\
\approx T_{l'}^{(k')} T_{l}^{(k)}(\mathfrak{s}) \otimes u_{k,l}[\alpha - E_{l}^{(k)}(\mathfrak{s})] \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(T_{l}^{(k)}(\mathfrak{s}))].
\end{aligned}
\quad (2.34)
$$

Comparing (2.33) and (2.34), one obtains (2.31) and (2.32).

**Remark 2.1.** The composition of the combinatorial $R$'s achieving $B \otimes B \mapsto B \otimes B$ is the identity. Therefore if vac is put in place of $u_{r,l}$ at the left in the diagram (2.26), we have $b'_1 = \text{vac}, b'_2 = b_1$,

$b'_3 = b_2, \dots$. It follows that the composition $T_1^{(1)} T_1^{(2)} \cdots T_1^{(n-1)}$ is the translation to the right by one lattice unit.

**Remark 2.2.** Both the time evolution $T_l^{(k)}$ and the associated energy $E_l^{(k)}$ have the well-defined limit as $l \to \infty$. A simple explanation of this fact is provided by the inverse scattering scheme (2.48), which translates $T_l^{(k)}$ as in (2.49) and $E_l^{(k)}$ as in (2.50) and (C.10). In particular, if $\gamma_k$ denotes the maximal amplitude of color $k$ solitons that will be defined in the next subsection, one has $T_l^{(k)} = T_{\gamma_k}^{(k)}$ and $E_l^{(k)} = E_{\gamma_k}^{(k)}$ for all $l \ge \gamma_k$.

## 2.5 Solitons and their scattering

This subsection is the place where 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$^{\otimes 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_l^{(1)}, \dots, E_l^{(n-1)}$ ($l \in \mathbb{Z}_{\ge 1}$). In what follows we present a complete list of solitons together with their scattering rule.

First we introduce the elementary excitations $s_1, s_2, \dots, s_{n-1} \in B$ by

$$ n=2: \quad s_1 = [2], \qquad (2.35) $$

$$ n=3: \quad s_1 = [2] \otimes \frac{1}{2}, \qquad s_2 = [1] \otimes \frac{1}{3}, \qquad (2.36) $$

$$ n=4: \quad s_1 = [2] \otimes \frac{1}{2} \otimes \frac{1}{3}, \qquad s_2 = [1] \otimes \frac{1}{3} \otimes \frac{1}{2}, \qquad s_3 = [1] \otimes \frac{1}{2} \otimes \frac{1}{4}. \qquad (2.37) $$

They are different from vac $\in B$ in (2.22) by only one entry at the bottom of one tableau.
on any state and for arbitrary $(k, l), (k', l') \in \mathcal{I}$. This is a simple consequence of the Yang-Baxter equation. In fact, write $\mathfrak{s} = b_1 \otimes \cdots \otimes b_L$ for short and consider

$$ \begin{aligned} u_{k,l}[\alpha] \otimes u_{k',l'}[\beta] \otimes \mathfrak{s} &\approx u_{k,l}[\alpha] \otimes T_{l'}^{(k')}(\mathfrak{s}) \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(\mathfrak{s})] \\ &\approx T_{l}^{(k)} T_{l'}^{(k')}(\mathfrak{s}) \otimes u_{k,l}[\alpha - E_{l}^{(k)}(T_{l'}^{(k')}(\mathfrak{s}))] \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(\mathfrak{s})] \end{aligned} \quad (2.33) $$

by using (2.28) successively. On the other hand, one may first apply $u_{k,l}[\alpha] \otimes u_{k',l'}[\beta] \approx u_{k',l'}[\beta] \otimes u_{k,l}[\alpha]$ which is a consequence of (2.13) and (2.14). It tells that (2.33) is also equal to

$$ \begin{aligned} T_{l'}^{(k')} T_{l}^{(k)}(\mathfrak{s}) \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(T_{l}^{(k)}(\mathfrak{s}))] \otimes u_{k,l}[\alpha - E_{l}^{(k)}(\mathfrak{s})] \\ \approx T_{l'}^{(k')} T_{l}^{(k)}(\mathfrak{s}) \otimes u_{k,l}[\alpha - E_{l}^{(k)}(\mathfrak{s})] \otimes u_{k',l'}[\beta - E_{l'}^{(k')}(T_{l}^{(k)}(\mathfrak{s}))]. \end{aligned} \quad (2.34) $$

Comparing (2.33) and (2.34), one obtains (2.31) and (2.32).

**Remark 2.1.** The composition of the combinatorial $R$'s achieving $B \otimes B \mapsto B \otimes B$ is the identity. Therefore if vac is put in place of $u_{r,l}$ at the left in the diagram (2.26), we have $b'_1 = \text{vac}, b'_2 = b_1$,

$b'_3 = b_2, \dots$. It follows that the composition $T_1^{(1)} T_1^{(2)} \cdots T_1^{(n-1)}$ is the translation to the right by one lattice unit.

**Remark 2.2.** Both the time evolution $T_l^{(k)}$ and the associated energy $E_l^{(k)}$ have the well-defined limit as $l \to \infty$. A simple explanation of this fact is provided by the inverse scattering scheme (2.48), which translates $T_l^{(k)}$ as in (2.49) and $E_l^{(k)}$ as in (2.50) and (C.10). In particular, if $\gamma_k$ denotes the maximal amplitude of color $k$ solitons that will be defined in the next subsection, one has $T_l^{(k)} = T_{\gamma_k}^{(k)}$ and $E_l^{(k)} = E_{\gamma_k}^{(k)}$ for all $l \ge \gamma_k$.

## 2.5 Solitons and their scattering

This subsection is the place where 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$^{\otimes 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_l^{(1)}, \dots, E_l^{(n-1)}$ ($l \in \mathbb{Z}_{\ge 1}$). In what follows we present a complete list of solitons together with their scattering rule.

First we introduce the elementary excitations $s_1, s_2, \dots, s_{n-1} \in B$ by

$$ n=2: \quad s_1 = [2], \quad (2.35) $$

$$ n=3: \quad s_1 = [2] \otimes \frac{1}{2}, \quad s_2 = \frac{1}{3}, \quad (2.36) $$

$$ n=4: \quad s_1 = [2] \otimes \frac{1}{2} \otimes \frac{1}{3}, \quad s_2 = \frac{1}{3} \otimes \frac{1}{3} \otimes \frac{1}{2}, \quad s_3 = \frac{1}{2} \otimes \frac{1}{2} \otimes \frac{1}{4}. \quad (2.37) $$

They are different from vac $\in B$ in (2.22) by only one entry at the bottom of one tableau.
General definition is this:

$$s_a = u_{1,1} \otimes \cdots u_{a-1,1} \otimes \cdots \otimes u_{a+1,1} \otimes \cdots \otimes u_{n-1,1} \in B \quad (a \in [1, n-1]). \quad (2.38)$$

Using $s_a$ as the building block we next introduce

$$S_i^{(a)} = \overbrace{s_a \otimes \cdots \otimes s_a}^{i} \in B^{\otimes i} \quad ((a, i) \in \mathcal{I}) \quad (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.

(I) When isolated, a soliton $S_i^{(a)}$ proceeds to the right with the velocity $\delta_{ar} \min(i, l)$ under $T_l^{(r)}$.

(II) Solitons are stable under collisions. A collision of solitons of type (a, i) and (b, j) induces the displacement $\Delta$ in their asymptotic trajectories which is common in the magnitude and opposite in the direction.

<figure><img src="image_5.png" /></figure>

Here time grows from the top to the bottom. The quantity $\Delta$, we call it phase shift, is given by

$$\Delta = C_{ab} \min(i, j), \quad (2.40)$$

where $(C_{ab})_{1 \le a, b \le n-1}$ is the Cartan matrix (2.2). Note that $\Delta$ 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

$$\dots (\text{color 1 solitons}) \dots\dots (\text{color 2 solitons}) \dots\dots \dots (\text{color } n-1 \text{ solitons}) \dots\dots \quad (2.41)$$

where $\dots$ denotes vac's. Distance of the neighboring islands can be made as large as one wishes. Each (color a solitons) has the form

$$\dots (S_1^{(a)}, s) \dots\dots (S_2^{(a)}, s) \dots\dots \dots (S_{\gamma_a}^{(a)}, s) \dots\dots \quad (2.42)$$
General definition is this:

$$s_a = u_{1,1} \otimes \cdots u_{a-1,1} \otimes \cdots \otimes u_{a+1,1} \otimes \cdots \otimes u_{n-1,1} \in B \quad (a \in [1, n-1]). \quad (2.38)$$

Using $s_a$ as the building block we next introduce

$$S_i^{(a)} = \overbrace{s_a \otimes \cdots \otimes s_a}^{i} \in B^{\otimes i} \quad ((a, i) \in \mathcal{I}) \quad (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.

(I) When isolated, a soliton $S_i^{(a)}$ proceeds to the right with the velocity $\delta_{ar} \min(i, l)$ under $T_l^{(r)}$.

(II) Solitons are stable under collisions. A collision of solitons of type (a, i) and (b, j) induces the displacement $\Delta$ in their asymptotic trajectories which is common in the magnitude and opposite in the direction.

<figure><img src="image_5.png" /></figure>

Here time grows from the top to the bottom. The quantity $\Delta$, we call it phase shift, is given by

$$\Delta = C_{ab} \min(i, j), \quad (2.40)$$

where $(C_{ab})_{1 \le a, b \le n-1}$ is the Cartan matrix (2.2). Note that $\Delta$ 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

$$\dots (\text{color 1 solitons}) \dots\dots (\text{color 2 solitons}) \dots\dots \dots (\text{color } n-1 \text{ solitons}) \dots\dots \quad (2.41)$$

where $\dots$ denotes vac's. Distance of the neighboring islands can be made as large as one wishes. Each (color a solitons) has the form

$$\dots (S_1^{(a)}, s) \dots\dots (S_2^{(a)}, s) \dots\dots \dots (S_{\gamma_a}^{(a)}, s) \dots\dots \quad (2.42)$$
where $\gamma_a$ is the maximal length of the color a solitons as mentioned in Remark 2.2. Again the distance of $(S_k^{(a)})$ and $(S_{k+1}^{(a)})$ can be made arbitrarily large. Each $(S_i^{(a)})$ can be brought into the form

$$ \dots S_i^{(a)} \dots S_i^{(a)} \dots \dots S_i^{(a)} \dots \qquad (2.43) $$

In contrast to (2.41) and (2.42), the separation $\dots$ of the neighboring $S_i^{(a)}$'s remains constant under time evolutions, but it is at least $i$ everywhere, namely, $\text{vac}^{\otimes d}$ with $d \ge i$.

(IV) The energy $E_l^{(r)}$ of the asymptotic states (hence all those connected to it by time evolutions) described in (2.41) – (2.43) is given by

$$ E_l^{(r)} = \sum_{i=1}^{\gamma_r} \min(i, l) m_i^{(r)}, \qquad (2.44) $$

where $m_i^{(a)}$ is the number (multiplicity) of type $(a, i)$ solitons in (2.43).

By inverting (2.44) one can obtain the $m_i^{(r)}$, that is the soliton content of the state, from the conserved energies $E_l^{(r)}$. In addition, by using the local energies (2.30), one gets the local contribution to $m_i^{(r)}$, which is therefore a soliton density. This is how the soliton densities are computed in the simulations in Sec. 5.3.

**Remark 2.3.** Denote the Cartan matrix (2.2) by $C_n$ exhibiting the dependence on $n$ temporarily. In view of the inverse $(C_n^{-1})_{ab} = \min(a, b) - \frac{ab}{n}$, the phase shift $C_{ab} \min(i, j)$ in (2.40) is the matrix element of $C_n \otimes C_\infty^{-1}$. This originates in the $x = 0$ Fourier component of the TBA kernel $\hat{A}_{ab}^{ij}(x) \hat{M}_{ab}(x)$ for the level $\ell U_q(\bar{s}l_n)$ RSOS model in [37, eq.(14.14)] in the limit $\ell \to \infty$.

Reflecting the commutativity (2.31), the phase shift $\Delta$ is independent of the choice of the time evolution $T_l^{(r)}$ as long as the two solitons have different bare speeds in (I) under it so that they eventually collide.

Let us present a few examples of (II) concerning collisions. We take $n = 4$.

**Example 2.4.** Collision of $S_3^{(a)}$ and $S_1^{(a)}$ under the time evolution by successive applications of $T_l^{(a)}$ with any $a = 1, 2, 3$ and $l \ge 3$.

A blank, , signifies the vacuum vac. See (2.22). The larger (resp. smaller) soliton $S_3^{(a)}$ (resp. $S_1^{(a)}$) is pushed forward (resp. pulled backward) by 2 sites compared from its free motion. This agrees with the phase shift $C_{aa} \min(3, 1) = 2$ in (2.40). Note that we count B having the internal structure as one lattice unit. In this example, no local state other than vac and $s_a$ is generated by the collision. In general if solitons of only one color are present, cBBS is equivalent to the very original BBS [1] corresponding to $n = 2$.

**Example 2.5.** Successive collisions of $S_3^{(3)}$ with $S_1^{(1)}$ and $S_2^{(2)}$ under the time evolution $T_l^{(3)}$ with $l \ge 3$. We write $s_2$ in (2.37) for example as $\frac{1}{12}$ to save the space.
<figure><img src="image_6.png" /></figure>

The solitons $S_1^{(1)}$ and $S_2^{(2)}$ do not move by themselves as their bare speed is zero under $T_l^{(3)}$. See (I). One observes that the phase shift of $S_3^{(3)}$ vs $S_1^{(1)}$ collision is $C_{31} \min(3, 1) = 0$, whereas the one for $S_3^{(3)}$ vs $S_2^{(2)}$ collision is $C_{32} \min(3, 2) = -2$ in agreement with (2.40). The intermediate states contain the non-vacuum local states $\frac{1}{224}, \frac{1}{144}, \frac{1}{144}, \frac{1}{134}$ which are not included in the elementary excitations $s_1, s_2, s_3$ in (2.37).

**Example 2.6.** Let us demonstrate how a generic state is decomposed into solitons as indicated in (III). We pick the state on the top line and first apply $T_l^{(3)}$ with $l \ge 1$ repeatedly to get

<figure><img src="image_7.png" /></figure>

The two color 3 solitons $S_1^{(3)}$ are separated. Their distance is 1 which is indeed not less than their common length in agreement with the last claim in (III). The two solitons can be taken away to the right by further applying $T_l^{(3)}$ without changing their mutual distance. It allows us to focus on the remaining left part of the state. Applying $T_l^{(2)}$ with $l \ge 3$ to it successively, we find

<figure><img src="image_8.png" /></figure>

Thus the color 2 solitons $S_1^{(2)}$ and $S_3^{(2)}$ are separated to the right leaving the color 1 soliton $S_1^{(1)}$ in the left. These procedures achieve (2.41)–(2.43).

**Remark 2.7.** Instead of (2.41) one can achieve the asymptotic decomposition

$$ \dots (\text{color } a_1 \text{ solitons}) \dots\dots (\text{color } a_2 \text{ solitons}) \dots\dots (\text{color } a_{n-1} \text{ solitons}) \dots \tag{2.45} $$
for arbitrary permutation $a_1, a_2, \dots, a_{n-1}$ of $1, 2, \dots, n-1$. As indicated by Example 2.6, it is
done by applying $(T_{\infty}^{(a_1)})^{M_1} (T_{\infty}^{(a_2)})^{M_2} \dots (T_{\infty}^{(a_{n-1})})^{M_{n-1}}$ with $M_1, M_2, \dots, M_{n-1} \gg 1$.

**Remark 2.8.** The tableaux letter $a \in [1, n]$ carries the $sl_n$ weight $\varpi_1 - \alpha_1 - \cdots - \alpha_{a-1}$. In this counting, the elementary excitation $s_a$ in (2.38) has the extra weight $-\alpha_a$ compared from the background vac (2.21). Thus the soliton $S_i^{(a)}$ (2.39) carries the weight $-i\alpha_a$.

In terms of letters in the boxes of semistandard tableaux, the presence of a soliton $S_i^{(a)}$ changes the letter $a$ into $a+1$ for $i$ times compared with the vacuum background. On the other hand, vac in (2.21) contains a letter $a$ for $(n-a)$ times. Therefore the number $\lambda_a$ of the letter $a$ in the $L$ site system is

$$ \lambda_a = L(n-a) - \sum_{i=1}^{\infty} i m_i^{(a)} + \sum_{i=1}^{\infty} i m_i^{(a-1)} \quad (a \in [1, n]), \qquad (2.46) $$

where $m_i^{(a)}$ the number of solitons $S_i^{(a)}$ and $m_i^{(0)} = 0$ by convention.<fnref n="3" /> This formula is nothing but (C.14)<fnref n="4" /> with (C.10) under the special choice $N = (n-1)L$ and $k_i \equiv i \in \mathbb{Z}_{n-1}$ with $i \in [1, n-1]$. In terms of $\rho_i^{(a)}$ and $\varepsilon_i^{(a)}$ that will be introduced in (3.16) and (3.17), the density $\varrho_a := \lambda_a/L$ of the tableau letter $a$ is given by

$$ \varrho_a = n - a - \sum_{i=1}^{\infty} i \rho_i^{(a)} + \sum_{i=1}^{\infty} i \rho_i^{(a-1)} = n - a - \varepsilon_{\infty}^{(a)} + \varepsilon_{\infty}^{(a-1)} \quad (a \in [1, n]), \qquad (2.47) $$

where $\rho_i^{(0)} = \varepsilon_{\infty}^{(0)} = 0$.

## 2.6 Inverse scattering method

In Appendix C a brief exposition is given on the KSS bijection $\Phi$ between the set of highest paths $\mathcal{P}_+(B, \lambda)$ and rigged configurations $\text{RC}(B, \lambda)$. States of our cBBS satisfying the boundary condition are examples of the former with $B = B^{\otimes L}$. See the remark in the end of Appendix C. The array $\lambda = (\lambda_1, \dots, \lambda_n)$ specifies that the number of letter $a$ in a state is $\lambda_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_l^{(k)}$ of the BBS induces that of rigged configurations via the commutative diagram:<fnref n="5" />

$$ \begin{gathered} \begin{gathered} \mathcal{P}_+(B, \lambda) \xrightarrow{\Phi} \text{RC}(B, \lambda) \\ T_l^{(k)} \downarrow \qquad \downarrow T_l^{(k)}. \end{gathered} \qquad (2.48) \\ \begin{gathered} \mathcal{P}_+(B, \lambda) \xrightarrow{\Phi} \text{RC}(B, \lambda) \end{gathered} \end{gathered} $$

A remarkable feature of this scheme is that it achieves the linearization [38]:

$$ T_l^{(k)} : \{(a_i, j_i, r_i)\} \mapsto \{(a_i, j_i, r_i + \delta_{k,a_i} \min(l, j_i))\}, \qquad (2.49) $$

<footnote n="3">³The formula (2.46) is obvious for asymptotic states. General case follows from it and the fact that $m_i^{(a)}$'s are conserved quantities.</footnote>

<footnote n="4">⁴The formula (C.14) is given for rigged configurations, but it is known to agree with (C.1) under the KSS bijection.</footnote>

<footnote n="5">⁵We use the same notation $T_l^{(k)}$ to also represent the time evolution of the rigged configurations since they are to be identified via the bijection $\Phi$.</footnote>
for arbitrary permutation $a_1, a_2, \dots, a_{n-1}$ of $1, 2, \dots, n-1$. As indicated by Example 2.6, it is
done by applying $(T_{\infty}^{(a_1)})^{M_1} (T_{\infty}^{(a_2)})^{M_2} \dots (T_{\infty}^{(a_{n-1})})^{M_{n-1}}$ with $M_1, M_2, \dots, M_{n-1} \gg 1$.

**Remark 2.8.** The tableaux letter $a \in [1, n]$ carries the $sl_n$ weight $\varpi_1 - \alpha_1 - \cdots - \alpha_{a-1}$. In this counting, the elementary excitation $s_a$ in (2.38) has the extra weight $-\alpha_a$ compared from the background vac (2.21). Thus the soliton $S_i^{(a)}$ (2.39) carries the weight $-i\alpha_a$.

In terms of letters in the boxes of semistandard tableaux, the presence of a soliton $S_i^{(a)}$ changes the letter $a$ into $a+1$ for $i$ times compared with the vacuum background. On the other hand, vac in (2.21) contains a letter $a$ for $(n-a)$ times. Therefore the number $\lambda_a$ of the letter $a$ in the $L$ site system is

$$ \lambda_a = L(n-a) - \sum_{i=1}^{\infty} i m_i^{(a)} + \sum_{i=1}^{\infty} i m_i^{(a-1)} \quad (a \in [1, n]), \qquad (2.46) $$

where $m_i^{(a)}$ the number of solitons $S_i^{(a)}$ and $m_i^{(0)} = 0$ by convention.<fnref n="3" /> This formula is nothing but (C.14)<fnref n="4" /> with (C.10) under the special choice $N = (n-1)L$ and $k_i \equiv i \in \mathbb{Z}_{n-1}$ with $i \in [1, n-1]$. In terms of $\rho_i^{(a)}$ and $\varepsilon_i^{(a)}$ that will be introduced in (3.16) and (3.17), the density $\varrho_a := \lambda_a/L$ of the tableau letter $a$ is given by

$$ \varrho_a = n - a - \sum_{i=1}^{\infty} i \rho_i^{(a)} + \sum_{i=1}^{\infty} i \rho_i^{(a-1)} = n - a - \varepsilon_{\infty}^{(a)} + \varepsilon_{\infty}^{(a-1)} \quad (a \in [1, n]), \qquad (2.47) $$

where $\rho_i^{(0)} = \varepsilon_{\infty}^{(0)} = 0$.

## 2.6 Inverse scattering method

In Appendix C a brief exposition is given on the KSS bijection $\Phi$ between the set of highest paths $P_+(B, \lambda)$ and rigged configurations $RC(B, \lambda)$. States of our cBBS satisfying the boundary condition are examples of the former with $B = B^{\otimes L}$. See the remark in the end of Appendix C. The array $\lambda = (\lambda_1, \dots, \lambda_n)$ specifies that the number of letter $a$ in a state is $\lambda_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_l^{(k)}$ of the BBS induces that of rigged configurations via the commutative diagram:<fnref n="5" />

$$ \begin{gathered} \begin{gathered} \mathcal{P}_+(B, \lambda) \xrightarrow{\Phi} \mathrm{RC}(B, \lambda) \\ T_l^{(k)} \downarrow \qquad \qquad \qquad \qquad \downarrow T_l^{(k)}. \end{gathered} \qquad (2.48) \\ \begin{gathered} \mathcal{P}_+(B, \lambda) \xrightarrow{\Phi} \mathrm{RC}(B, \lambda) \end{gathered} \end{gathered} $$

A remarkable feature of this scheme is that it achieves the linearization [38]:

$$ T_l^{(k)} : \{(a_i, j_i, r_i)\} \mapsto \{(a_i, j_i, r_i + \delta_{k,a_i} \min(l, j_i))\}, \qquad (2.49) $$

<footnote n="3">³The formula (2.46) is obvious for asymptotic states. General case follows from it and the fact that $m_i^{(a)}$'s are conserved quantities.</footnote>

<footnote n="4">⁴The formula (C.14) is given for rigged configurations, but it is known to agree with (C.1) under the KSS bijection.</footnote>

<footnote n="5">⁵We use the same notation $T_l^{(k)}$ to also represent the time evolution of the rigged configurations since they are to be identified via the bijection $\Phi$.</footnote>
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 $\Phi$ and $\Phi^{-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_l^{(k)} = \Phi^{-1} \circ T_l^{(k)} \circ \Phi$.

Take a state $b_1 \otimes \cdots \otimes b_L$ of cBBS and let $m_j^{(a)}$ be the number of color $a$ length $j$ strings in the rigged configuration $\Phi(b_1 \otimes \cdots \otimes b_L)$. We set $\mathcal{E}_j^{(a)} = \mathcal{E}_j^{(a)}(b_1 \otimes \cdots \otimes b_L) = \sum_{k \ge 1} \min(j, k) m_k^{(a)}$ according to (C.10). Since $m_j^{(a)}$ is determined only by the action variables $\{(a_i, j_i)\}$, the quantity $\mathcal{E}_j^{(a)}$ is conserved under the time evolution. On the other hand recall the conserved energy of a state $E_l^{(k)} = E_l^{(k)}(b_1 \otimes \cdots \otimes b_L)$ introduced in (2.30). The two are known to coincide:

$$E_l^{(k)} = \mathcal{E}_l^{(k)}. \qquad (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_j^{(a)}$ in $E_l^{(k)}$ (2.44) was the number of color $a$ length $j$ solitons contained in a state, whereas the one in $\mathcal{E}_l^{(k)}$ is the number of color $a$ length $j$ strings. Therefore (2.50) implies the soliton/string correspondence [38]

$$\text{soliton} = \text{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:

$$\mathfrak{s} = \begin{array}{cccccccccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 123 & 123 & 123 & 134 & 134 & 123 & 123 & 123 & 123 \end{array} \in \mathcal{P}_+(B^{\otimes 8}, (23, 13, 10, 2)), \qquad (2.51)$$

where a few vacuums are supplied in the front and in the tail and displayed explicitly. Such operators do not influence the relevant rigged configurations essentially. See the explanation in the end of Appendix C. From the demonstration in the previous subsection, we know that $\mathfrak{s}$ “consists of” solitons $S_1^{(3)}$, $S_1^{(3)}$ of color 3, solitons $S_3^{(2)}$, $S_1^{(2)}$ of color 2 and a soliton $S_1^{(1)}$ of color 1 including multiplicity. This indeed matches the image of the map $\Phi$:

$$\Phi(\mathfrak{s}) = \begin{array}{c} \square 5 \\ \boxed{3} \end{array} \begin{array}{c} \boxed{0} \\ \boxed{3} \end{array}, \qquad (2.52)$$

where the precise meaning of the above graphical representation for the rigged configuration is detailed in Appendix C.2.

# 3 Randomized cBBS

## 3.1 Generalized Gibbs ensemble

Now we introduce the randomized version of cBBS which corresponds to the distribution of the initial states according to the generalized Gibbs ensemble (GGE). The GGE partition function for size $L$ system $B^{\otimes L}$ is the sum of the Boltzmann factor $\exp(-\sum_{r,l} \beta_l^{(r)} E_l^{(r)})$ over all the
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 $\Phi$ and $\Phi^{-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_l^{(k)} = \Phi^{-1} \circ T_l^{(k)} \circ \Phi$.

Take a state $b_1 \otimes \cdots \otimes b_L$ of cBBS and let $m_j^{(a)}$ be the number of color $a$ length $j$ strings in the rigged configuration $\Phi(b_1 \otimes \cdots \otimes b_L)$. We set $\mathcal{E}_j^{(a)} = \mathcal{E}_j^{(a)}(b_1 \otimes \cdots \otimes b_L) = \sum_{k \ge 1} \min(j, k) m_k^{(a)}$ according to (C.10). Since $m_j^{(a)}$ is determined only by the action variables $\{(a_i, j_i)\}$, the quantity $\mathcal{E}_j^{(a)}$ is conserved under the time evolution. On the other hand recall the conserved energy of a state $E_l^{(k)} = E_l^{(k)}(b_1 \otimes \cdots \otimes b_L)$ introduced in (2.30). The two are known to coincide:

$$E_l^{(k)} = \mathcal{E}_l^{(k)}. \qquad (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_j^{(a)}$ in $E_l^{(k)}$ (2.44) was the number of color $a$ length $j$ solitons contained in a state, whereas the one in $\mathcal{E}_l^{(k)}$ is the number of color $a$ length $j$ strings. Therefore (2.50) implies the soliton/string correspondence [38]

$$\text{soliton} = \text{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:

$$\mathfrak{s} = \begin{array}{cccccccccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 12 & \otimes & 12 & \otimes & 12 & \otimes & 13 & \otimes & 23 & \otimes & 12 & \otimes & 12 & \otimes & 12 \end{array} \in \mathcal{P}_+(B^{\otimes 8}, (23, 13, 10, 2)), \quad (2.51)$$

where a few vacuums $\mathfrak{s} = \begin{array}{c} 1 \\ 123 \end{array}$ are supplied in the front and in the tail and displayed explicitly. Such operators do not influence the relevant rigged configurations essentially. See the explanation in the end of Appendix C. From the demonstration in the previous subsection, we know that $\mathfrak{s}$ “consists of” solitons $S_1^{(3)}$, $S_1^{(3)}$ of color 3, solitons $S_3^{(2)}$, $S_1^{(2)}$ of color 2 and a soliton $S_1^{(1)}$ of color 1 including multiplicity. This indeed matches the image of the map $\Phi$:

$$\Phi(\mathfrak{s}) = \begin{array}{c} \square 5 \\ \boxed{3} \boxed{3} \boxed{0} \boxed{3} \boxed{3} \end{array}, \qquad (2.52)$$

where the precise meaning of the above graphical representation for the rigged configuration is detailed in Appendix C.2.

# 3 Randomized cBBS

## 3.1 Generalized Gibbs ensemble

Now we introduce the randomized version of cBBS which corresponds to the distribution of the initial states according to the generalized Gibbs ensemble (GGE). The GGE partition function for size $L$ system $B^{\otimes L}$ is the sum of the Boltzmann factor $\exp(-\sum_{r,l} \beta_l^{(r)} E_l^{(r)})$ over all the
states of the cBBS, where $E_l^{(k)}$ is the conserved energy. See (2.30), (2.44), (2.50) and (C.10) for the definition and the properties. Especially from (2.44), the system can really be regarded as a gas of interacting solitons with various colors and lengths. The parameter $\beta_l^{(r)}$ is the inverse temperature associated with $E_l^{(r)}$. We are interested in the asymptotic behavior as $L \to \infty$.

At this stage a nontrivial question arises as how to incorporate properly the boundary condition that all the distant local states are vac (2.20). 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

$$Z_L(\{\beta_l^{(r)}\}) = \sum_{\text{highest states } \in B^{\otimes L}} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} E_l^{(r)}). \quad (3.1)$$

By synthesizing various results it can be rewritten as follows:

$$Z_L(\{\beta_l^{(r)}\}) = \sum_{\lambda} \sum_{s \in \mathcal{P}_+(B^{\otimes L}, \lambda)} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} E_l^{(r)}(s)) \quad (\text{definition (C.5)}) \\ = \sum_{\lambda} \sum_{S \in \text{RC}(B^{\otimes L}, \lambda)} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \mathcal{E}_l^{(r)}(S))$$

(2.50) and 1:1 correspondence (C.17))

$$= \sum_{\lambda} \sum_{\{m_j^{(a)}\}}^{(\lambda)} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \sum_{j \ge 1} \min(j,l) m_j^{(r)}) \prod_{(a,j) \in \mathcal{I}} \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}} \quad (\text{from (C.10), (C.16)})$$

$$= \sum_{\{m_j^{(a)}\}} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \sum_{j \ge 1} \min(j,l) m_j^{(r)}) \prod_{(a,j) \in \mathcal{I}} \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}}. \quad (3.2)$$

Here the vacancy $p_i^{(a)}$ defined in general by (C.10) takes the form

$$p_i^{(a)} = L - \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i,j) m_j^{(b)}, \quad (3.3)$$

for our cBBS reflecting the fact that the states belong to $B^{\otimes L}$ with $B$ given by (2.21). The constraint (C.8) is taken into account by setting $\binom{K}{J} = 0$ unless $J \in [0, K]$. In (3.2) we are summing over the number $m_j^{(a)}$ of solitons $S_j^{(a)}$. In this sense $Z_L(\{\beta_l^{(r)}\})$ is a GGE partition function for the grand canonical ensemble of solitons.

**Remark 3.1.** The appearence of the phase shift $\Delta = C_{ab} \min(i, j)$ (defined in (2.40)) in (3.3) 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).

## 3.2 I.I.D. randomness

From here we concentrate on the situation where local states $b \in B = B^{1,1} \otimes \dots \otimes B^{n-1,1}$ obey an i.i.d. randomness. More specifically we consider an i.i.d. measure $\mathbb{P}(b \in B)$ including
states of the cBBS, where $E_l^{(k)}$ is the conserved energy. See (2.30), (2.44), (2.50) and (C.10) for the definition and the properties. Especially from (2.44), the system can really be regarded as a gas of interacting solitons with various colors and lengths. The parameter $\beta_l^{(r)}$ is the inverse temperature associated with $E_l^{(r)}$. We are interested in the asymptotic behavior as $L \to \infty$.

At this stage a nontrivial question arises as how to incorporate properly the boundary condition that all the distant local states are vac (2.20). 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

$$Z_L(\{\beta_l^{(r)}\}) = \sum_{\text{highest states } \in B^{\otimes L}} \exp(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} E_l^{(r)}). \quad (3.1)$$

By synthesizing various results it can be rewritten as follows:

$$
\begin{align}
Z_L(\{\beta_l^{(r)}\}) &= \sum_{\lambda} \sum_{s \in \mathcal{P}_+(B^{\otimes L}, \lambda)} \exp\left(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} E_l^{(r)}(s)\right) \quad (\text{definition (C.5)}) \\
&= \sum_{\lambda} \sum_{S \in \text{RC}(B^{\otimes L}, \lambda)} \exp\left(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \mathcal{E}_l^{(r)}(S)\right) \tag{2.50 and 1:1 correspondence (C.17)} \\
&= \sum_{\lambda} \sum_{\{m_j^{(a)}\}}^{(\lambda)} \exp\left(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \sum_{j \ge 1} \min(j, l) m_j^{(r)}\right) \prod_{(a,j) \in \mathcal{I}} \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}} \tag{from (C.10), (C.16)} \\
&= \sum_{\{m_j^{(a)}\}} \exp\left(-\sum_{(r,l) \in \mathcal{I}} \beta_l^{(r)} \sum_{j \ge 1} \min(j, l) m_j^{(r)}\right) \prod_{(a,j) \in \mathcal{I}} \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}}. \tag{3.2}
\end{align}
$$

Here the vacancy $p_i^{(a)}$ defined in general by (C.10) takes the form

$$p_i^{(a)} = L - \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) m_j^{(b)}, \quad (3.3)$$

for our cBBS reflecting the fact that the states belong to $B^{\otimes L}$ with $B$ given by (2.21). The constraint (C.8) is taken into account by setting $\binom{K}{J} = 0$ unless $J \in [0, K]$. In (3.2) we are summing over the number $m_j^{(a)}$ of solitons $S_j^{(a)}$. In this sense $Z_L(\{\beta_l^{(r)}\})$ is a GGE partition function for the grand canonical ensemble of solitons.

**Remark 3.1.** The appearence of the phase shift $\Delta = C_{ab} \min(i, j)$ (defined in (2.40)) in (3.3) 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).

## 3.2 I.I.D. randomness

From here we concentrate on the situation where local states $b \in B = B^{1,1} \otimes \dots \otimes B^{n-1,1}$ obey an i.i.d. randomness. More specifically we consider an i.i.d. measure $\mathbb{P}(b \in B)$ including
the parameters $z_1 > z_2 > \cdots > z_n > 0$ as follows:

$$ \mathbb{P}(b = c_1 \otimes \cdots \otimes c_{n-1} \in B) = \mathbb{P}^{(1)}(c_1) \cdots \mathbb{P}^{(n-1)}(c_{n-1}), \quad (3.4) $$

$$ \mathbb{P}^{(k)}(c = (a_1, \dots, a_k) \in B^{k,1}) = \frac{z_{a_1} \cdots z_{a_k}}{\sum_{c \in B^{k,1}} z_{a_1} \cdots z_{a_k}}, \quad (3.5) $$

where $a_i$ denotes the entry of the tableau $c \in B^{k,1}$ in the $i$ th box from the top. The formula (3.5) implies that the probability of occurrence of $a \in [1, n]$ is proportional to $z_a$. In this sense, $z_1, \dots, z_n$ are tableau variables representing the fugacities of the letters $1, \dots, n$. Their ambiguity due to the invariance of (3.5) under the change $z_a \to 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_1^{(k)}$ in the notation (3.10).

For a partition $\nu = (\nu_1, \nu_2, \dots, \nu_n)$, let $s_\nu = s_\nu(z_1, \dots, z_n)$ denote the Schur polynomial<fnref n="6" /> [13]:

$$ s_{\nu}(z_1, \dots, z_n) = \frac{\det(z_j^{\nu_i + n - i})_{i,j=1}^n}{\det(z_j^{n-i})_{i,j=1}^n}. \quad (3.6) $$

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:

$$ z_a = w_{a-1}^{-1} w_a \quad (1 \le a \le n), \quad w_0 = w_n = 1, \quad z_1 z_2 \cdots z_n = 1, \quad (3.7) $$

$$ w_a = e^{\varpi_a} = z_1 z_2 \cdots z_a \quad (0 \le a \le n), \quad \varpi_0 = \varpi_n = 0, \quad (3.8) $$

$$ x_a = e^{-\alpha_a} = z_a^{-1} z_{a+1} = \prod_{b=1}^{n-1} w_b^{-C_{ab}} \quad (1 \le a \le n-1). \quad (3.9) $$

The variables $w_a$ and $x_a$ are formal exponential of the fundamental weight $\varpi_a$ and the negated simple root $-\alpha_a$ regarded as parameters. See Sec. 2.1. The last relation in (3.7) leads to $s_\nu = s_{\tilde{\nu}}$ where $\tilde{\nu} = (\nu_1 - \nu_n, \dots, \nu_{n-1} - \nu_n, 0)$. We use a special notation for rectangle $\nu$'s as

$$ Q_i^{(a)} = s_{(i^a)}(z_1, \dots, z_n). \quad (3.10) $$

It satisfies the Q-system [31,37]

$$ (Q_i^{(a)})^2 = Q_{i-1}^{(a)} Q_{i+1}^{(a)} + Q_i^{(a-1)} Q_i^{(a+1)}, \quad (3.11) $$

for $(a, i) \in \mathcal{I}$ with $Q_i^{(0)} = Q_i^{(n)} = 1$. To validate (3.11) also at $i=0$ with $Q_0^{(a)} = 1$, we employ the convention $Q_{-1}^{(a)} = 0$.

Let us show that the i.i.d. measure (3.4) is the special $(n-1)$-parameter reduction of the GGE in Sec. 3.1. Denote the weight of a state $s \in B^{\otimes L}$ in the sense of (C.1) by $\mathrm{wt}(s) = (\lambda_1, \dots, \lambda_n)$. Then the formulas (3.4) and (3.5) mean that $s$ is realized with the probability proportional to $z_1^{\lambda_1} \cdots z_n^{\lambda_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 $\prod_{a=1}^n z_a^{\lambda_a}$ is proportional to $\prod_{a=1}^n z_a^{-\mathcal{E}_\infty^{(a)} + \mathcal{E}_\infty^{(a-1)}}$, where we have dropped $\sum_{i=1}^N \theta(k_i \ge a)$ in (C.14) since it is the constant $(n-a)L$ for our states in $B^{\otimes L}$. (Note $\mathcal{E}_j^{(0)} = \mathcal{E}_j^{(n)} = 0$ as mentioned there.) We can further replace $\mathcal{E}_\infty^{(a)}$ by $E_\infty^{(a)}$ thanks to (2.50). Therefore the probability of the state $s$ is proportional to

$$ z_1^{-E_\infty^{(1)}} z_2^{-E_\infty^{(2)} + E_\infty^{(1)}} \cdots z_{n-1}^{-E_\infty^{(n-1)}} + E_\infty^{(n-2)} z_n^{E_\infty^{(n-1)}} = \prod_{a=1}^{n-1} \left( \frac{z_{a+1}}{z_a} \right)^{E_\infty^{(a)}} = \exp\left( - \sum_{a=1}^{n-1} \alpha_a E_\infty^{(a)} \right). \quad (3.12) $$

<footnote n="6">⁶The Schur polynomial $s_\nu$ should not be confused with $s_a$ in (2.38).</footnote>
the parameters $z_1 > z_2 > \dots > z_n > 0$ as follows:

$$ \mathbb{P}(b = c_1 \otimes \dots \otimes c_{n-1} \in B) = \mathbb{P}^{(1)}(c_1) \dots \mathbb{P}^{(n-1)}(c_{n-1}), \qquad (3.4) $$

$$ \mathbb{P}^{(k)}(c = (a_1, \dots, a_k) \in B^{k,1}) = \frac{z_{a_1} \cdots z_{a_k}}{\sum_{c \in B^{k,1}} z_{a_1} \cdots z_{a_k}}, \qquad (3.5) $$

where $a_i$ denotes the entry of the tableau $c \in B^{k,1}$ in the $i$ th box from the top. The formula (3.5) implies that the probability of occurrence of $a \in [1, n]$ is proportional to $z_a$. In this sense, $z_1, \dots, z_n$ are tableau variables representing the fugacities of the letters $1, \dots, n$. Their ambiguity due to the invariance of (3.5) under the change $z_a \to 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_1^{(k)}$ in the notation (3.10).

For a partition $\nu = (\nu_1, \nu_2, \dots, \nu_n)$, let $s_\nu = s_\nu(z_1, \dots, z_n)$ denote the Schur polynomial<fnref n="6" /> [13]:

$$ s_{\nu}(z_1, \dots, z_n) = \frac{\det(z_j^{\nu_i + n - i})_{i,j=1}^n}{\det(z_j^{n-i})_{i,j=1}^n}. \qquad (3.6) $$

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:

$$ z_a = w_{a-1}^{-1} w_a \qquad (1 \le a \le n), \qquad w_0 = w_n = 1, \qquad z_1 z_2 \cdots z_n = 1, \qquad (3.7) $$

$$ w_a = e^{\varpi_a} = z_1 z_2 \cdots z_a \qquad (0 \le a \le n), \qquad \varpi_0 = \varpi_n = 0, \qquad (3.8) $$

$$ x_a = e^{-\alpha_a} = z_a^{-1} z_{a+1} = \prod_{b=1}^{n-1} w_b^{-C_{ab}} \qquad (1 \le a \le n-1). \qquad (3.9) $$

The variables $w_a$ and $x_a$ are formal exponential of the fundamental weight $\varpi_a$ and the negated simple root $-\alpha_a$ regarded as parameters. See Sec. 2.1. The last relation in (3.7) leads to $s_\nu = s_{\tilde{\nu}}$ where $\tilde{\nu} = (\nu_1 - \nu_n, \dots, \nu_{n-1} - \nu_n, 0)$. We use a special notation for rectangle $\nu$'s as

$$ Q_i^{(a)} = s_{(i^a)}(z_1, \dots, z_n). \qquad (3.10) $$

It satisfies the Q-system [31,37]

$$ (Q_i^{(a)})^2 = Q_{i-1}^{(a)} Q_{i+1}^{(a)} + Q_i^{(a-1)} Q_i^{(a+1)}, \qquad (3.11) $$

for $(a, i) \in \mathcal{I}$ with $Q_i^{(0)} = Q_i^{(n)} = 1$. To validate (3.11) also at $i=0$ with $Q_0^{(a)} = 1$, we employ the convention $Q_{-1}^{(a)} = 0$.

Let us show that the i.i.d. measure (3.4) is the special $(n-1)$-parameter reduction of the GGE in Sec. 3.1. Denote the weight of a state $s \in B^{\otimes L}$ in the sense of (C.1) by $\mathrm{wt}(s) = (\lambda_1, \dots, \lambda_n)$. Then the formulas (3.4) and (3.5) mean that $s$ is realized with the probability proportional to $z_1^{\lambda_1} \cdots z_n^{\lambda_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 $\prod_{a=1}^n z_a^{\lambda_a}$ is proportional to $\prod_{a=1}^n z_a^{-\mathcal{E}_{\infty}^{(a)} + \mathcal{E}_{\infty}^{(a-1)}}$, where we have dropped $\sum_{i=1}^N \theta(k_i \ge a)$ in (C.14) since it is the constant $(n-a)L$ for our states in $B^{\otimes L}$. (Note $\mathcal{E}_j^{(0)} = \mathcal{E}_j^{(n)} = 0$ as mentioned there.) We can further replace $\mathcal{E}_{\infty}^{(a)}$ by $E_{\infty}^{(a)}$ thanks to (2.50). Therefore the probability of the state $s$ is proportional to

$$ z_1^{-E_{\infty}^{(1)}} z_2^{-E_{\infty}^{(2)} + E_{\infty}^{(1)}} \cdots z_{n-1}^{-E_{\infty}^{(n-1)} + E_{\infty}^{(n-2)}} z_n^{E_{\infty}^{(n-1)}} = \prod_{a=1}^{n-1} \left( \frac{z_{a+1}}{z_a} \right)^{E_{\infty}^{(a)}} = \exp\left( - \sum_{a=1}^{n-1} \alpha_a E_{\infty}^{(a)} \right). \qquad (3.12) $$

<footnote n="6">⁶The Schur polynomial $s_\nu$ should not be confused with $s_a$ in (2.38).</footnote>
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 $\beta_l^{(a)}$ is zero except

$$ \beta_{\infty}^{(a)} = \alpha_a \quad (a \in [1, n-1]). \qquad (3.13) $$

Thus the inverse temperatures (with index $\infty$) can be identified with the simple roots. We will nonetheless allow co-existence of the symbols $\beta_{\infty}^{(a)}$ and $\alpha_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, \alpha_a, x_a$ as

$$ z_1 > \cdots > z_n > 0 \quad (z_1 \cdots z_n = 1), \quad \alpha_1, \dots, \alpha_{n-1} > 0, \quad x_1, \dots, x_{n-1} < 1. \qquad (3.14) $$

In particular $x_a \to 0$ and $x_a \to 1$ correspond to the zero and infinite limit of the temperature $1/\beta_{\infty}^{(a)}$, respectively.

**Remark 3.2.** The dilute limit where no soliton is allowed is given by $z_1 \gg z_2 \gg \cdots \gg z_n$ or equivalently $x_1, \dots, x_n \to 0$. In fact, $\mathbb{P}^{(k)}(c) \to \delta_{c,u_{k,1}}$ in (3.5) hence $\mathbb{P}(b) \to \delta_{b,\text{vac}}$ in (3.4) hold in this limit. One also has $Q_i^{(a)}/z_a^i \to 1$. See (2.12) and (2.21) for the definition of $u_{k,1}$ and vac.

Now the partition function (3.2) is reduced to

$$ Z_L(\{\beta_{\infty}^{(r)}\}) = \sum_{\{m_j^{(a)}\}} \prod_{(a,j) \in \mathcal{I}} \exp(-\beta_{\infty}^{(a)} j m_j^{(a)}) \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}}. \qquad (3.15) $$

### 3.3 Thermodynamic Bethe ansatz

In the large $L$ limit, the dominant contribution in the sum (3.15) comes from those $\{m_j^{(a)}\}$ having the $L$-linear asymptotic behavior

$$ m_i^{(a)} \approx L \rho_i^{(a)}, \quad p_i^{(a)} \approx L \sigma_i^{(a)}, \quad \mathcal{E}_i^{(a)} \approx L \varepsilon_i^{(a)}, \qquad (3.16) $$

$$ \sigma_i^{(a)} = 1 - \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)}, \quad \varepsilon_i^{(a)} = \sum_{j \ge 1} \min(i, j) \rho_j^{(a)}, \qquad (3.17) $$

where the second line follows from the first by (3.3). The relation (3.17) is a spectral parameter free version of the Bethe equation in terms of the string density $\rho_i^{(a)}$ and the hole density $\sigma_i^{(a)}$. One should seek $\rho = (\rho_i^{(a)})$ that minimizes the free energy per site

$$ F[\rho] = \sum_{a=1}^{n-1} \beta_{\infty}^{(a)} \sum_{i=1}^{s} i \rho_i^{(a)} - \sum_{a=1}^{n-1} \sum_{i=1}^{s} \left( (\rho_i^{(a)} + \sigma_i^{(a)}) \log(\rho_i^{(a)} + \sigma_i^{(a)}) - \rho_i^{(a)} \log \rho_i^{(a)} - \sigma_i^{(a)} \log \sigma_i^{(a)} \right). \qquad (3.18) $$

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 $\varepsilon_i^{(a)} = \sum_{j=1}^{s} \min(i, j) \rho_j^{(a)}$.
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 $\beta_l^{(a)}$ is zero except

$$ \beta_{\infty}^{(a)} = \alpha_a \quad (a \in [1, n-1]). \qquad (3.13) $$

Thus the inverse temperatures (with index $\infty$) can be identified with the simple roots. We will nonetheless allow co-existence of the symbols $\beta_{\infty}^{(a)}$ and $\alpha_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, \alpha_a, x_a$ as

$$ z_1 > \cdots > z_n > 0 \quad (z_1 \cdots z_n = 1), \quad \alpha_1, \dots, \alpha_{n-1} > 0, \quad x_1, \dots, x_{n-1} < 1. \quad (3.14) $$

In particular $x_a \to 0$ and $x_a \to 1$ correspond to the zero and infinite limit of the temperature $1/\beta_{\infty}^{(a)}$, respectively.

**Remark 3.2.** The dilute limit where no soliton is allowed is given by $z_1 \gg z_2 \gg \cdots \gg z_n$ or equivalently $x_1, \dots, x_n \to 0$. In fact, $\mathbb{P}^{(k)}(c) \to \delta_{c,u_{k,1}}$ in (3.5) hence $\mathbb{P}(b) \to \delta_{b,\text{vac}}$ in (3.4) hold in this limit. One also has $Q_i^{(a)}/z_a^i \to 1$. See (2.12) and (2.21) for the definition of $u_{k,1}$ and vac.

Now the partition function (3.2) is reduced to

$$ Z_L(\{\beta_{\infty}^{(r)}\}) = \sum_{\{m_j^{(a)}\}} \prod_{(a,j) \in \mathcal{I}} \exp(-\beta_{\infty}^{(a)} j m_j^{(a)}) \binom{p_j^{(a)} + m_j^{(a)}}{m_j^{(a)}}. \qquad (3.15) $$

### 3.3 Thermodynamic Bethe ansatz

In the large $L$ limit, the dominant contribution in the sum (3.15) comes from those $\{m_j^{(a)}\}$ having the $L$-linear asymptotic behavior

$$ m_i^{(a)} \approx L \rho_i^{(a)}, \quad p_i^{(a)} \approx L \sigma_i^{(a)}, \quad \mathcal{E}_i^{(a)} \approx L \varepsilon_i^{(a)}, \qquad (3.16) $$

$$ \sigma_i^{(a)} = 1 - \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)}, \quad \varepsilon_i^{(a)} = \sum_{j \ge 1} \min(i, j) \rho_j^{(a)}, \qquad (3.17) $$

where the second line follows from the first by (3.3). The relation (3.17) is a spectral parameter free version of the Bethe equation in terms of the string density $\rho_i^{(a)}$ and the hole density $\sigma_i^{(a)}$. One should seek $\rho = (\rho_i^{(a)})$ that minimizes the free energy per site

$$ F[\rho] = \sum_{a=1}^{n-1} \beta_{\infty}^{(a)} \sum_{i=1}^{s} i \rho_i^{(a)} - \sum_{a=1}^{n-1} \sum_{i=1}^{s} \left( (\rho_i^{(a)} + \sigma_i^{(a)}) \log(\rho_i^{(a)} + \sigma_i^{(a)}) - \rho_i^{(a)} \log \rho_i^{(a)} - \sigma_i^{(a)} \log \sigma_i^{(a)} \right). \qquad (3.18) $$

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 $\varepsilon_i^{(a)} = \sum_{j=1}^{s} \min(i, j) \rho_j^{(a)}$.
From $\frac{\partial \sigma_j^{(b)}}{\partial \rho_i^{(a)}} = -C_{ab} \min(i, j)$, one finds that the condition $\frac{\partial F[\rho]}{\partial \rho_i^{(a)}} = 0$ is expressed as a TBA equation

$$ -i\beta_{\infty}^{(a)} + \log(1 + (y_i^{(a)})^{-1}) = \sum_{b=1}^{n-1} C_{ab} \sum_{j=1}^{s} \min(i, j) \log(1 + y_j^{(b)}), \quad (3.19) $$

in terms of the ratio<fnref n="7" />

$$ y_i^{(a)} = \frac{\rho_i^{(a)}}{\sigma_i^{(a)}}. \quad (3.20) $$

The corresponding maximal value $\mathcal{F}$ of the free energy per site $F[\rho]$ (3.18) is obtained by using (3.17), (3.19) and (3.20). The result reads

$$ \mathcal{F} = - \sum_{a=1}^{n-1} \sum_{i=1}^{s} \log(1 + y_i^{(a)}). \quad (3.21) $$

The TBA equation (3.19) is equivalent to the constant Y-system

$$ \frac{(1 + (y_i^{(a)})^{-1})^2}{(1 + (y_{i-1}^{(a)})^{-1})(1 + (y_{i+1}^{(a)})^{-1})} = \prod_{b=1}^{n-1} (1 + y_i^{(b)})^{C_{ab}} \quad (1 \le i \le s), \quad (3.22) $$

with the boundary condition

$$ (y_0^{(a)})^{-1} = 0, \quad 1 + (y_{s+1}^{(a)})^{-1} = e^{\beta_{\infty}^{(a)}} (1 + (y_s^{(a)})^{-1}). \quad (3.23) $$

The Y-system (3.22) follows from the Q-system (3.11) by the substitution (cf. [37, Prop. 14.1])

$$ y_i^{(a)} = \frac{Q_i^{(a-1)} Q_i^{(a+1)}}{Q_{i-1}^{(a)} Q_{i+1}^{(a)}}, \quad 1 + (y_i^{(a)})^{-1} = \prod_{b=1}^{n-1} (Q_i^{(b)})^{C_{ab}}, \quad 1 + y_i^{(a)} = \frac{(Q_i^{(a)})^2}{Q_{i-1}^{(a)} Q_{i+1}^{(a)}}. \quad (3.24) $$

As for the boundary condition (3.23), the first relation is valid due to the convention $Q_{-1}^{(a)} = 0$ mentioned after (3.11). On the other hand the second relation is translated into

$$ e^{\beta_{\infty}^{(a)}} = \prod_{b=1}^{n-1} \left( \frac{Q_{s+1}^{(b)}}{Q_s^{(b)}} \right)^{C_{ab}}. \quad (3.25) $$

At this point we let $s$ tend to infinity. The result [40, Th. 7.1 (C)] tells that $\lim_{s \to \infty} (Q_{s+1}^{(a)}/Q_s^{(a)}) = e^{\bar{\varpi}_a}$ in the regime $e^{\alpha_1}, \dots, e^{\alpha_{n-1}} > 1$ under consideration. Thus the large $s$ limit of (3.25) leads to $e^{\beta_{\infty}^{(a)}} = \prod_{b=1}^{n-1} (e^{\bar{\varpi}_b})^{C_{ab}} = e^{\alpha_a}$ by (3.9). This is consistent with (3.13). We also remark that $\lim_{s \to \infty} (Q_{s+1}^{(a)}/Q_s^{(a)}) = e^{\bar{\varpi}_a}$ and the last relation in (3.24) tells that $\lim_{i \to \infty} y_i^{(a)} = 0$, therefore (3.20) indicates $\lim_{i \to \infty} \rho_i^{(a)} = 0$ which is indeed necessary for the convergence of $\varepsilon_{\infty}^{(a)} = \sum_{i=1}^{\infty} i \rho_i^{(a)}$.

Substituting the last expression in (3.24) into the free energy density (3.21), we get

$$ \mathcal{F} = - \sum_{a=1}^{n-1} \log \left( \frac{Q_1^{(a)} Q_s^{(a)}}{Q_{s+1}^{(a)}} \right) \underset{s \to \infty}{\longrightarrow} - \sum_{a=1}^{n-1} \log \bar{Q}_1^{(a)}, \quad (3.26) $$

$$ \bar{Q}_i^{(a)} = e^{-i\bar{\varpi}_a} Q_i^{(a)} = 1 + \dots \in \mathbb{Z}_{\ge 0}[x_1, \dots, x_{n-1}]. \quad (3.27) $$

<footnote n="7">⁷This $y_i^{(a)}$ is the inverse of $Y_i^{(a)}$ in [17, eq.(52)].</footnote>
From $\frac{\partial \sigma_j^{(b)}}{\partial \rho_i^{(a)}} = -C_{ab} \min(i, j)$, one finds that the condition $\frac{\partial F[\rho]}{\partial \rho_i^{(a)}} = 0$ is expressed as a TBA equation

$$ -i\beta_{\infty}^{(a)} + \log(1 + (y_i^{(a)})^{-1}) = \sum_{b=1}^{n-1} C_{ab} \sum_{j=1}^{s} \min(i, j) \log(1 + y_j^{(b)}), \quad (3.19) $$

in terms of the ratio<fnref n="7" />

$$ y_i^{(a)} = \frac{\rho_i^{(a)}}{\sigma_i^{(a)}}. \quad (3.20) $$

The corresponding maximal value $\mathcal{F}$ of the free energy per site $F[\rho]$ (3.18) is obtained by using (3.17), (3.19) and (3.20). The result reads

$$ \mathcal{F} = - \sum_{a=1}^{n-1} \sum_{i=1}^{s} \log(1 + y_i^{(a)}). \quad (3.21) $$

The TBA equation (3.19) is equivalent to the constant Y-system

$$ \frac{(1 + (y_i^{(a)})^{-1})^2}{(1 + (y_{i-1}^{(a)})^{-1})(1 + (y_{i+1}^{(a)})^{-1})} = \prod_{b=1}^{n-1} (1 + y_i^{(b)})^{C_{ab}} \quad (1 \le i \le s), \quad (3.22) $$

with the boundary condition

$$ (y_0^{(a)})^{-1} = 0, \quad 1 + (y_{s+1}^{(a)})^{-1} = e^{\beta_{\infty}^{(a)}} (1 + (y_s^{(a)})^{-1}). \quad (3.23) $$

The Y-system (3.22) follows from the Q-system (3.11) by the substitution (cf. [37, Prop. 14.1])

$$ y_i^{(a)} = \frac{Q_i^{(a-1)} Q_i^{(a+1)}}{Q_{i-1}^{(a)} Q_{i+1}^{(a)}}, \quad 1 + (y_i^{(a)})^{-1} = \prod_{b=1}^{n-1} (Q_i^{(b)})^{C_{ab}}, \quad 1 + y_i^{(a)} = \frac{(Q_i^{(a)})^2}{Q_{i-1}^{(a)} Q_{i+1}^{(a)}}. \quad (3.24) $$

As for the boundary condition (3.23), the first relation is valid due to the convention $Q_{-1}^{(a)} = 0$ mentioned after (3.11). On the other hand the second relation is translated into

$$ e^{\beta_{\infty}^{(a)}} = \prod_{b=1}^{n-1} \left( \frac{Q_{s+1}^{(b)}}{Q_{s}^{(b)}} \right)^{C_{ab}}. \quad (3.25) $$

At this point we let $s$ tend to infinity. The result [40, Th. 7.1 (C)] tells that $\lim_{s \to \infty} (Q_{s+1}^{(a)}/Q_s^{(a)}) = e^{\bar{\varpi}_a}$ in the regime $e^{\alpha_1}, \dots, e^{\alpha_{n-1}} > 1$ under consideration. Thus the large $s$ limit of (3.25) leads to $e^{\beta_{\infty}^{(a)}} = \prod_{b=1}^{n-1} (e^{\bar{\varpi}_b})^{C_{ab}} = e^{\alpha_a}$ by (3.9). This is consistent with (3.13). We also remark that $\lim_{s \to \infty} (Q_{s+1}^{(a)}/Q_s^{(a)}) = e^{\bar{\varpi}_a}$ and the last relation in (3.24) tells that $\lim_{i \to \infty} y_i^{(a)} = 0$, therefore (3.20) indicates $\lim_{i \to \infty} \rho_i^{(a)} = 0$ which is indeed necessary for the convergence of $\varepsilon_{\infty}^{(a)} = \sum_{i=1}^{\infty} i \rho_i^{(a)}$.

Substituting the last expression in (3.24) into the free energy density (3.21), we get

$$ \mathcal{F} = - \sum_{a=1}^{n-1} \log \left( \frac{Q_1^{(a)} Q_s^{(a)}}{Q_{s+1}^{(a)}} \right) \underset{s \to \infty}{\longrightarrow} - \sum_{a=1}^{n-1} \log \bar{Q}_1^{(a)}, \quad (3.26) $$

$$ \bar{Q}_i^{(a)} = e^{-i\bar{\varpi}_a} Q_i^{(a)} = 1 + \dots \in \mathbb{Z}_{\ge 0}[x_1, \dots, x_{n-1}]. \quad (3.27) $$

<footnote n="7">⁷This $y_i^{(a)}$ is the inverse of $Y_i^{(a)}$ in [17, eq.(52)].</footnote>
The quantity $\overline{Q}_i^{(a)}$ is the normalized character of the irreducible sl$_n$ module with highest weight $i\overline{\omega}_a$.

From (3.17) one sees that $\epsilon_{\infty}^{(a)} = \sum_{i \ge 1} i \rho_i^{(a)}$ is the total density of color a solitons. Let $\epsilon_1, \dots, \epsilon_{n-1}$ be the expectation values of $\epsilon_{\infty}^{(1)}, \dots, \epsilon_{\infty}^{(n-1)}$. From (3.1), they are to be derived from the free energy density by the standard prescription $\epsilon_a = \frac{\partial F}{\partial \beta_{\infty}^{(a)}}$. Further substitution of (3.26) into this leads to the equation of state of our cBBS as follows:

$$ \epsilon_a = - \sum_{b=1}^{n-1} \frac{\partial \log \overline{Q}_1^{(b)}}{\partial \alpha_a} = x_a \frac{\partial}{\partial x_a} \log(\overline{Q}_1^{(1)} \cdots \overline{Q}_1^{(n-1)}) \quad (a \in [1, n-1]), \quad (3.28) $$

where we have written $\beta_{\infty}^{(a)}$ as $\alpha_a$ by (3.13) to express the RHS purely in terms of representation theoretical data. An analogous result on "non-complete BBS" was given in [17, eq. (66)].

**Example 3.3.** In the simplest case $n=2$, $z_2 = z_1^{-1}$, $x_1 = z_1^{-2}$. Thus we have

$$ Q_1^{(1)} = z_1 + z_2, \quad \overline{Q}_1^{(1)} = 1 + x_1, \quad \epsilon_1 = \frac{x_1}{1+x_1} = \frac{z_2}{z_1+z_2}, \quad (3.29) $$

which agrees with the density of total number of "balls" 2. Similarly for $n=3$, one has

$$ Q_1^{(1)} = z_1 + z_2 + z_3, \quad \overline{Q}_1^{(1)} = 1 + x_1 + x_1 x_2, \quad \epsilon_1 = \frac{x_1(1+x_2)}{1+x_1+x_1 x_2} + \frac{x_1 x_2}{1+x_2+x_1 x_2}, \quad (3.30) $$

$$ Q_1^{(2)} = z_1 z_2 + z_1 z_3 + z_2 z_3, \quad \overline{Q}_1^{(2)} = 1 + x_2 + x_1 x_2, \quad \epsilon_2 = \frac{x_2(1+x_1)}{1+x_2+x_1 x_2} + \frac{x_1 x_2}{1+x_1+x_1 x_2}. \quad (3.31) $$

In this way one can relate the densities and the inverse temperatures.

To summarize so far, we have obtained the solution $\{y_i^{(a)}\}$ to the TBA equation (3.19)$_{|s=\infty}$ in (3.24), (3.10) and (3.6). They depend on $n-1$ independent parameters which may be taken either as the inverse temperatures $\beta_{\infty}^{(a)} = \alpha_a (1 \le a < n)$ or the fugacities $z_a (1 \le a \le n)$ as in (3.7)–(3.9). They can further be related to the prescribed expectation values $\epsilon_a$ of $\epsilon_{\infty}^{(a)} = \sum_{i \ge 1} i \rho_i^{(a)} (1 \le a < n)$ by the equation of state (3.28).

The remaining task is to express the string and hole densities $\rho_i^{(a)}$ and $\sigma_i^{(a)}$ in terms of these variables. From (3.20) and (3.17), it suffices to determine $\epsilon_i^{(a)}$. It is characterized by the second order difference equation and the boundary condition

$$ (y_i^{(a)})^{-1}(-\epsilon_{i-1}^{(a)} + 2\epsilon_i^{(a)} - \epsilon_{i+1}^{(a)}) + \sum_{b=1}^{n-1} C_{ab} \epsilon_i^{(b)} = 1, \quad (3.32) $$

$$ \epsilon_0^{(a)} = 0, \quad \epsilon_{\infty}^{(a)} = x_a \frac{\partial}{\partial x_a} \log(\overline{Q}_1^{(1)} \cdots \overline{Q}_1^{(n-1)}). \quad (3.33) $$

The equation (3.32) is just the first relation in (3.17), where the first term is a disguised form of $\sigma_i^{(a)}$. 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) \in \mathcal{I}$, the unique solution $h_i^{(a)} = h_i^{(a)}(r, l)$ to

$$ (y_i^{(a)})^{-1}(-h_{i-1}^{(a)} + 2h_i^{(a)} - h_{i+1}^{(a)}) + \sum_{b=1}^{n-1} C_{ab} h_i^{(b)} = \delta_{a,r} \min(i, l), \quad (3.34) $$

$$ h_0^{(a)} = 0, \quad h_{\infty}^{(a)} = x_a \frac{\partial}{\partial x_a} \log \overline{Q}_l^{(r)} \quad (3.35) $$
The quantity $\overline{Q}_i^{(a)}$ is the normalized character of the irreducible sl$_n$ module with highest weight $i\overline{\omega}_a$.

From (3.17) one sees that $\epsilon_{\infty}^{(a)} = \sum_{i \ge 1} i \rho_i^{(a)}$ is the total density of color a solitons. Let $\epsilon_1, \dots, \epsilon_{n-1}$ be the expectation values of $\epsilon_{\infty}^{(1)}, \dots, \epsilon_{\infty}^{(n-1)}$. From (3.1), they are to be derived from the free energy density by the standard prescription $\epsilon_a = \frac{\partial F}{\partial \beta_{\infty}^{(a)}}$. Further substitution of (3.26) into this leads to the equation of state of our cBBS as follows:

$$ \epsilon_a = - \sum_{b=1}^{n-1} \frac{\partial \log \overline{Q}_1^{(b)}}{\partial \alpha_a} = x_a \frac{\partial}{\partial x_a} \log(\overline{Q}_1^{(1)} \cdots \overline{Q}_1^{(n-1)}) \quad (a \in [1, n-1]), \quad (3.28) $$

where we have written $\beta_{\infty}^{(a)}$ as $\alpha_a$ by (3.13) to express the RHS purely in terms of representation theoretical data. An analogous result on "non-complete BBS" was given in [17, eq. (66)].

**Example 3.3.** In the simplest case $n=2$, $z_2 = z_1^{-1}$, $x_1 = z_1^{-2}$. Thus we have

$$ Q_1^{(1)} = z_1 + z_2, \quad \overline{Q}_1^{(1)} = 1 + x_1, \quad \epsilon_1 = \frac{x_1}{1+x_1} = \frac{z_2}{z_1+z_2}, \quad (3.29) $$

which agrees with the density of total number of "balls" 2. Similarly for $n=3$, one has

$$ Q_1^{(1)} = z_1 + z_2 + z_3, \quad \overline{Q}_1^{(1)} = 1 + x_1 + x_1 x_2, \quad \epsilon_1 = \frac{x_1(1+x_2)}{1+x_1+x_1 x_2} + \frac{x_1 x_2}{1+x_2+x_1 x_2}, \quad (3.30) $$

$$ Q_1^{(2)} = z_1 z_2 + z_1 z_3 + z_2 z_3, \quad \overline{Q}_1^{(2)} = 1 + x_2 + x_1 x_2, \quad \epsilon_2 = \frac{x_2(1+x_1)}{1+x_2+x_1 x_2} + \frac{x_1 x_2}{1+x_1+x_1 x_2}. \quad (3.31) $$

In this way one can relate the densities and the inverse temperatures.

To summarize so far, we have obtained the solution $\{y_i^{(a)}\}$ to the TBA equation (3.19)$_{|s=\infty}$ in (3.24), (3.10) and (3.6). They depend on $n-1$ independent parameters which may be taken either as the inverse temperatures $\beta_{\infty}^{(a)} = \alpha_a (1 \le a < n)$ or the fugacities $z_a (1 \le a \le n)$ as in (3.7)-(3.9). They can further be related to the prescribed expectation values $\epsilon_a$ of $\epsilon_{\infty}^{(a)} = \sum_{i \ge 1} i \rho_i^{(a)} (1 \le a < n)$ by the equation of state (3.28).

The remaining task is to express the string and hole densities $\rho_i^{(a)}$ and $\sigma_i^{(a)}$ in terms of these variables. From (3.20) and (3.17), it suffices to determine $\epsilon_i^{(a)}$. It is characterized by the second order difference equation and the boundary condition

$$ (y_i^{(a)})^{-1}(-\epsilon_{i-1}^{(a)} + 2\epsilon_i^{(a)} - \epsilon_{i+1}^{(a)}) + \sum_{b=1}^{n-1} C_{ab} \epsilon_i^{(b)} = 1, \quad (3.32) $$

$$ \epsilon_0^{(a)} = 0, \quad \epsilon_{\infty}^{(a)} = x_a \frac{\partial}{\partial x_a} \log(\overline{Q}_1^{(1)} \cdots \overline{Q}_1^{(n-1)}). \quad (3.33) $$

The equation (3.32) is just the first relation in (3.17), where the first term is a disguised form of $\sigma_i^{(a)}$. 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) \in \mathcal{I}$, the unique solution $h_i^{(a)} = h_i^{(a)}(r, l)$ to

$$ (y_i^{(a)})^{-1}(-h_{i-1}^{(a)} + 2h_i^{(a)} - h_{i+1}^{(a)}) + \sum_{b=1}^{n-1} C_{ab} h_i^{(b)} = \delta_{a,r} \min(i, l), \quad (3.34) $$

$$ h_0^{(a)} = 0, \quad h_{\infty}^{(a)} = x_a \frac{\partial}{\partial x_a} \log \overline{Q}_l^{(r)} \quad (3.35) $$
is given by

$$h_i^{(a)}(r, l) = \frac{\sum_{\nu=(\nu_1, \dots, \nu_n)} (\sum_{j=\max(a,r)+1}^{n} \nu_j) s_\nu}{s_{(i^a)} s_{(l^r)}} (= h_l^{(r)}(a, i)), \quad (3.36)$$

where the outer sum in (3.36) runs over those Young diagrams $\nu$ 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.4.** For generic $n$, the Littlewood-Richardson rule [13] leads to

$$h_i^{(1)}(1, l) = \frac{1}{s_{(i)} s_{(l)}} \sum_{j=1}^{\min(i,l)} j s_{(i+l-j,j)}, \quad (3.37)$$

$$h_i^{(2)}(1, l) = \frac{1}{s_{(i,i)} s_{(l)}} \sum_{j=1}^{\min(i,l)} j s_{(i+l-j,i,j)}, \quad (3.38)$$

$$h_i^{(2)}(2, l) = \frac{1}{s_{(i,i)} s_{(l,l)}} \sum_{j,j'} (2l - 2j - j') s_{(i+j+j',i+j,l-j,l-j-j')}, \quad (3.39)$$

$$h_i^{(a)}(r, 1) = \frac{1}{s_{(i^a)} s_{(1^r)}} \sum_{j=1}^{\min(r,n-a)} j s_{((i+1)^{r-j}, i^{a-r+j}, 1^j)}, \quad (1 \le r \le a), \quad (3.40)$$

$$= \frac{1}{s_{(i^a)} s_{(1^r)}} \sum_{j=1}^{\min(a,n-r)} j s_{(i+1)^{a-j}, ij, 1^{r-a+j}}, \quad (a \le r \le n). \quad (3.41)$$

The sum (3.39) is taken over $j, j' \ge 0$ such that $j + j' \le l$ and $\max(l - i, 0) \le j \le l$. Given $n$, the sums (3.37)-(3.39) should be truncated so that the length of the partitions appearing in the indices of the Schur functions does not exceed $n$. For instance when $n = 3$, (3.39) is reduced to $(s_{(i,i)}s_{(l,l)})^{-1} \sum_{0 \le j < l} (l-j)s_{(i+l,i+j,l-j)}$. For $n = 4$ (3.40) and (3.41) read

$$h_i^{(1)}(1, 1) = \frac{s_{(i,1)}}{s_{(i)} s_{(1)}}, \quad h_i^{(1)}(2, 1) = \frac{s_{(i,1,1)}}{s_{(i)} s_{(1,1)}}, \quad h_i^{(1)}(3, 1) = \frac{s_{(i,1,1,1)}}{s_{(i)} s_{(1,1,1)}}, \quad (3.42)$$

$$h_i^{(2)}(1, 1) = \frac{s_{(i,i,1)}}{s_{(i,i)} s_{(1)}}, \quad h_i^{(2)}(2, 1) = \frac{s_{(i+1,i,1)} + 2s_{(i,i,1,1)}}{s_{(i,i)} s_{(1,1)}}, \quad h_i^{(2)}(3, 1) = \frac{s_{(i+1,i,1,1)}}{s_{(i,i)} s_{(1,1,1)}}, \quad (3.43)$$

$$h_i^{(3)}(1, 1) = \frac{s_{(i,i,i,1)}}{s_{(i,i,i)} s_{(1)}}, \quad h_i^{(3)}(2, 1) = \frac{s_{(i+1,i,i,1)}}{s_{(i,i,i)} s_{(1,1)}}, \quad h_i^{(3)}(3, 1) = \frac{s_{(i+1,i+1,i,1)}}{s_{(i,i,i)} s_{(1,1,1)}}, \quad (3.44)$$

Comparing the two systems (3.32), (3.33) and (3.34), (3.35), we find that the solution to the former is constructed as the superposition of the latter as

$$\varepsilon_i^{(a)} = \sum_{r=1}^{n-1} h_i^{(a)}(r, 1). \quad (3.45)$$

From (3.40)-(3.41) we obtain the solution to (3.32), (3.33) as

$$\varepsilon_i^{(a)} = \frac{1}{s_{(i^a)}} \sum_{r=1}^{n-1} \frac{1}{s_{(1^r)}} \sum_{k=1}^{\min(a,r,n-a,n-r)} k s_{((i+1)^{\min(a,r)-k}, i^{a+k-\min(a,r)}, 1^{r+k-\min(a,r)})}, \quad (3.46)$$

where $z_1 \cdots z_n = 1$ (3.7) is imposed on the Schur funciton (3.6).
is given by

$$h_i^{(a)}(r, l) = \frac{\sum_{\nu=(\nu_1, \dots, \nu_n)} (\sum_{j=\max(a,r)+1}^{n} \nu_j) s_\nu}{s_{(i^a)} s_{(l^r)}} (= h_l^{(r)}(a, i)), \quad (3.36)$$

where the outer sum in (3.36) runs over those Young diagrams $\nu$ 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.4.** For generic $n$, the Littlewood-Richardson rule [13] leads to

$$h_i^{(1)}(1, l) = \frac{1}{s_{(i)} s_{(l)}} \sum_{j=1}^{\min(i,l)} j s_{(i+l-j,j)}, \quad (3.37)$$

$$h_i^{(2)}(1, l) = \frac{1}{s_{(i,i)} s_{(l)}} \sum_{j=1}^{\min(i,l)} j s_{(i+l-j,i,j)}, \quad (3.38)$$

$$h_i^{(2)}(2, l) = \frac{1}{s_{(i,i)} s_{(l,l)}} \sum_{j,j'} (2l - 2j - j') s_{(i+j+j',i+j,l-j,l-j-j')}, \quad (3.39)$$

$$h_i^{(a)}(r, 1) = \frac{1}{s_{(i^a)} s_{(1^r)}} \sum_{j=1}^{\min(r,n-a)} j s_{((i+1)^{r-j}, i^{a-r+j}, 1^j)}, \quad (1 \le r \le a), \quad (3.40)$$

$$= \frac{1}{s_{(i^a)} s_{(1^r)}} \sum_{j=1}^{\min(a,n-r)} j s_{(i+1)^{a-j}, ij, 1^{r-a+j}}, \quad (a \le r \le n). \quad (3.41)$$

The sum (3.39) is taken over $j, j' \ge 0$ such that $j + j' \le l$ and $\max(l - i, 0) \le j \le l$. Given $n$, the sums (3.37)-(3.39) should be truncated so that the length of the partitions appearing in the indices of the Schur functions does not exceed $n$. For instance when $n = 3$, (3.39) is reduced to $(s_{(i,i)} s_{(l,l)})^{-1} \sum_{0 \le j < l} (l - j) s_{(i+l,i+j,l-j)}$. For $n = 4$ (3.40) and (3.41) read

$$h_i^{(1)}(1, 1) = \frac{s_{(i,1)}}{s_{(i)} s_{(1)}}, \quad h_i^{(1)}(2, 1) = \frac{s_{(i,1,1)}}{s_{(i)} s_{(1,1)}}, \quad h_i^{(1)}(3, 1) = \frac{s_{(i,1,1,1)}}{s_{(i)} s_{(1,1,1)}}, \quad (3.42)$$

$$h_i^{(2)}(1, 1) = \frac{s_{(i,i,1)}}{s_{(i,i)} s_{(1)}}, \quad h_i^{(2)}(2, 1) = \frac{s_{(i+1,i,1)} + 2s_{(i,i,1,1)}}{s_{(i,i)} s_{(1,1)}}, \quad h_i^{(2)}(3, 1) = \frac{s_{(i+1,i,1,1)}}{s_{(i,i)} s_{(1,1,1)}}, \quad (3.43)$$

$$h_i^{(3)}(1, 1) = \frac{s_{(i,i,i,1)}}{s_{(i,i,i)} s_{(1)}}, \quad h_i^{(3)}(2, 1) = \frac{s_{(i+1,i,i,1)}}{s_{(i,i,i)} s_{(1,1)}}, \quad h_i^{(3)}(3, 1) = \frac{s_{(i+1,i+1,i,1)}}{s_{(i,i,i)} s_{(1,1,1)}}, \quad (3.44)$$

Comparing the two systems (3.32), (3.33) and (3.34), (3.35), we find that the solution to the former is constructed as the superposition of the latter as

$$\varepsilon_i^{(a)} = \sum_{r=1}^{n-1} h_i^{(a)}(r, 1). \quad (3.45)$$

From (3.40)-(3.41) we obtain the solution to (3.32), (3.33) as

$$\varepsilon_i^{(a)} = \frac{1}{s_{(i^a)}} \sum_{r=1}^{n-1} \frac{1}{s_{(1^r)}} \sum_{k=1}^{\min(a,r,n-a,n-r)} k s_{((i+1)^{\min(a,r)-k}, i^{a+k-\min(a,r)}, 1^{r+k-\min(a,r)})}, \quad (3.46)$$

where $z_1 \cdots z_n = 1$ (3.7) is imposed on the Schur funciton (3.6).
**Example 3.5.** For small $n$, (3.46) reads as

$$n=2; \quad \varepsilon_i^{(1)} = \frac{s_{(i,1)}}{s_{(i)}s_{(1)}} = \frac{x_1(1-x_1^i)}{(1+x_1)(1-x_1^{i+1})}, \qquad (3.47)$$

$$n=3; \quad \varepsilon_i^{(1)} = \frac{1}{s_{(i)}} \left( \frac{s_{(i,1)}}{s_{(1)}} + \frac{s_{(i,1,1)}}{s_{(1,1)}} \right), \quad \varepsilon_i^{(2)} = \frac{1}{s_{(i,i)}} \left( \frac{s_{(i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,1)}}{s_{(1,1)}} \right), \quad (3.48)$$

$$n=4; \quad \varepsilon_i^{(1)} = \frac{1}{s_{(i)}} \left( \frac{s_{(i,1)}}{s_{(1)}} + \frac{s_{(i,1,1)}}{s_{(1,1)}} + \frac{s_{(i,1,1,1)}}{s_{(1,1,1)}} \right), \qquad (3.49)$$

$$\varepsilon_i^{(2)} = \frac{1}{s_{(i,i)}} \left( \frac{s_{(i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,1)}}{s_{(1,1)}} + \frac{2s_{(i,i,1,1)}}{s_{(1,1)}} + \frac{s_{(i+1,i,1,1)}}{s_{(1,1,1)}} \right), \qquad (3.50)$$

$$\varepsilon_i^{(3)} = \frac{1}{s_{(i,i,i)}} \left( \frac{s_{(i,i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,i,1)}}{s_{(1,1)}} + \frac{s_{(i+1,i+1,i,1)}}{s_{(1,1,1)}} \right). \qquad (3.51)$$

The result (3.47) reproduces [11, eq. (3.24)] with $a = z = x_1$.

Finally the densities $\rho_i^{(a)}$ and $\sigma_i^{(a)}$ are obtained from (3.17), (3.20) and (3.45) as

$$\rho_i^{(a)} = \sum_{r=1}^{n-1} (-h_{i-1}^{(a)}(r,1) + 2h_i^{(a)}(r,1) - h_{i+1}^{(a)}(r,1)) = (y_i^{(a)})^{-1} \sigma_i^{(a)}, \qquad (3.52)$$

$$\sigma_i^{(a)} = 1 - \sum_{b,r=1}^{n-1} C_{ab} h_i^{(b)}(r,1). \qquad (3.53)$$

# 4 Generalized hydrodynamics

## 4.1 Effective speed in homogenous system

According to Sec. 2.5, the configurations $S_i^{(a)}$ defined by (2.39) and (2.38) with $(a, i) \in \mathcal{I}$ provide the complete list of solitons. Under the time evolution $T_l^{(r)}$, a soliton $S_i^{(a)}$ possesses the bare speed $\delta_{ar} \min(i, l)$ and acquires the phase shift (2.40) in a collision with $S_j^{(b)}$. Let $v_i^{(a)} = v_i^{(a)}(r, l)$ denote the effective speed of $S_i^{(a)}$ under $T_l^{(r)}$. Then the speed equation in the sense of [41–43] takes the form

$$v_i^{(a)} = \delta_{ar} \min(i, l) + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) (v_i^{(a)} - v_j^{(b)}) \rho_j^{(b)}. \qquad (4.1)$$

Here $\rho_i^{(a)}$ is density of $S_i^{(a)}$ given either by $\rho_i^{(a)} = -\varepsilon_{i-1}^{(a)} + 2\varepsilon_i^{(a)} - \varepsilon_{i+1}^{(a)}$ or $\rho_i^{(a)} = y_i^{(a)}(1 - \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)})$ due to (3.17) and (3.20). It is independent of $T_l^{(r)}$ and the coincidence of the two expressions is assured by (3.32). For $n=2$ and $T_\infty^{(1)}$ dynamics, the equation (4.1) reduces to [23, eq.(11.7)].

From (3.17), the first term in the sum in (4.1) is equal to $v_i^{(a)} \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)} = v_i^{(a)}(1-\sigma_i^{(a)})$. Substitution of this into the speed equation (4.1) causes a cancellation of $v_i^{(a)}$, after which the result is expressed neatly in terms of the combination $\nu_i^{(a)}$ as follows:

$$\nu_i^{(a)} + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) y_j^{(b)} \nu_j^{(b)} = \delta_{ar} \min(i, l), \qquad (4.2)$$

$$\nu_i^{(a)} = \nu_i^{(a)}(r, l) = \sigma_i^{(a)} v_i^{(a)}(r, l). \qquad (4.3)$$
**Example 3.5.** For small $n$, (3.46) reads as

$$n = 2; \quad \varepsilon_i^{(1)} = \frac{s_{(i,1)}}{s_{(i)}s_{(1)}} = \frac{x_1(1 - x_1^i)}{(1 + x_1)(1 - x_1^{i+1})}, \qquad (3.47)$$

$$n = 3; \quad \varepsilon_i^{(1)} = \frac{1}{s_{(i)}} \left( \frac{s_{(i,1)}}{s_{(1)}} + \frac{s_{(i,1,1)}}{s_{(1,1)}} \right), \quad \varepsilon_i^{(2)} = \frac{1}{s_{(i,i)}} \left( \frac{s_{(i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,1)}}{s_{(1,1)}} \right), \qquad (3.48)$$

$$n = 4; \quad \varepsilon_i^{(1)} = \frac{1}{s_{(i)}} \left( \frac{s_{(i,1)}}{s_{(1)}} + \frac{s_{(i,1,1)}}{s_{(1,1)}} + \frac{s_{(i,1,1,1)}}{s_{(1,1,1)}} \right), \qquad (3.49)$$

$$\varepsilon_i^{(2)} = \frac{1}{s_{(i,i)}} \left( \frac{s_{(i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,1)}}{s_{(1,1)}} + \frac{2s_{(i,i,1,1)}}{s_{(1,1)}} + \frac{s_{(i+1,i,1,1)}}{s_{(1,1,1)}} \right), \qquad (3.50)$$

$$\varepsilon_i^{(3)} = \frac{1}{s_{(i,i,i)}} \left( \frac{s_{(i,i,i,1)}}{s_{(1)}} + \frac{s_{(i+1,i,i,1)}}{s_{(1,1)}} + \frac{s_{(i+1,i+1,i,1)}}{s_{(1,1,1)}} \right). \qquad (3.51)$$

The result (3.47) reproduces [11, eq. (3.24)] with $a = z = x_1$.

Finally the densities $\rho_i^{(a)}$ and $\sigma_i^{(a)}$ are obtained from (3.17), (3.20) and (3.45) as

$$\rho_i^{(a)} = \sum_{r=1}^{n-1} (-h_{i-1}^{(a)}(r, 1) + 2h_i^{(a)}(r, 1) - h_{i+1}^{(a)}(r, 1)) = (y_i^{(a)})^{-1} \sigma_i^{(a)}, \qquad (3.52)$$

$$\sigma_i^{(a)} = 1 - \sum_{b,r=1}^{n-1} C_{ab} h_i^{(b)}(r, 1). \qquad (3.53)$$

# 4 Generalized hydrodynamics

## 4.1 Effective speed in homogenous system

According to Sec. 2.5, the configurations $S_i^{(a)}$ defined by (2.39) and (2.38) with $(a, i) \in \mathcal{I}$ provide the complete list of solitons. Under the time evolution $T_l^{(r)}$, a soliton $S_i^{(a)}$ possesses the bare speed $\delta_{ar} \min(i, l)$ and acquires the phase shift (2.40) in a collision with $S_j^{(b)}$. Let $v_i^{(a)} = v_i^{(a)}(r, l)$ denote the effective speed of $S_i^{(a)}$ under $T_l^{(r)}$. Then the speed equation in the sense of [41–43] takes the form

$$v_i^{(a)} = \delta_{ar} \min(i, l) + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) (v_i^{(a)} - v_j^{(b)}) \rho_j^{(b)}. \qquad (4.1)$$

Here $\rho_i^{(a)}$ is density of $S_i^{(a)}$ given either by $\rho_i^{(a)} = -\varepsilon_{i-1}^{(a)} + 2\varepsilon_i^{(a)} - \varepsilon_{i+1}^{(a)}$ or $\rho_i^{(a)} = y_i^{(a)}(1 - \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)})$ due to (3.17) and (3.20). It is independent of $T_l^{(r)}$ and the coincidence of the two expressions is assured by (3.32). For $n = 2$ and $T_\infty^{(1)}$ dynamics, the equation (4.1) reduces to [23, eq.(11.7)].

From (3.17), the first term in the sum in (4.1) is equal to $v_i^{(a)} \sum_{b=1}^{n-1} C_{ab} \varepsilon_i^{(b)} = v_i^{(a)}(1 - \sigma_i^{(a)})$. Substitution of this into the speed equation (4.1) causes a cancellation of $v_i^{(a)}$, after which the result is expressed neatly in terms of the combination $\nu_i^{(a)}$ as follows:

$$\nu_i^{(a)} + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) y_j^{(b)} \nu_j^{(b)} = \delta_{ar} \min(i, l), \qquad (4.2)$$

$$\nu_i^{(a)} = \nu_i^{(a)}(r, l) = \sigma_i^{(a)} v_i^{(a)}(r, l). \qquad (4.3)$$
The equation (4.2) coincides with (3.34) under the identification $h_i^{(a)} = \sum_{j \ge 1} \min(i, j) y_j^{(a)} \nu_j^{(a)}$. Therefore it has a solution

$$ \nu_i^{(a)}(r, l) = \delta_{ar} \min(i, l) - \sum_{b=1}^{n-1} C_{ab} h_i^{(b)}(r, l) \quad (4.4) $$

$$ = (y_i^{(a)})^{-1}(-h_{i-1}^{(a)}(r, l) + 2h_i^{(a)}(r, l) - h_{i+1}^{(a)}(r, l)). \quad (4.5) $$

From the first formula (4.4) with $l = 1$ and (3.53) we have

$$ \sigma_i^{(a)} = \sum_{k=1}^{n-1} \nu_i^{(a)}(k, 1). \quad (4.6) $$

Combining (4.3) and (4.6) we obtain the effective speed

$$ v_i^{(a)}(r, l) = \frac{\nu_i^{(a)}(r, l)}{\sum_{k=1}^{n-1} \nu_i^{(a)}(k, 1)}, \quad (4.7) $$

in terms of $\nu_i^{(a)}(r, l)$ in (4.4)–(4.5) with $h_i^{(a)}(r, l)$ given by (3.36). The result leads to the relation

$$ \sum_{r=1}^{n-1} v_i^{(a)}(r, 1) = 1, \quad (4.8) $$

which is consistent with Remark 2.1.

Among the family of time evolutions $T_l^{(r)}$, typical ones are $T_\infty^{(1)}, \dots, T_\infty^{(n-1)}$. The corresponding limit of the numerator $\lim_{l \to \infty} \nu_i^{(a)}(r, l)$ in (4.7) can be evaluated compactly as

$$ \nu_i^{(a)}(r, \infty) = \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}h_i^{(b)}(r, \infty) = \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}h_\infty^{(r)}(b, i) \\ \stackrel{(3.35)}{=} \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}x_r \frac{\partial}{\partial x_r} \log(e^{-i\varpi_b} Q_i^{(b)}) = x_r \frac{\partial}{\partial x_r} \log \left( \frac{Q_i^{(a-1)} Q_i^{(a+1)}}{(Q_i^{(a)})^2} \right). \quad (4.9) $$

**Example 4.1.** Consider $n = 2$ case. The variables in (3.7)–(3.9) are related as $z_1 = z_2^{-1} = w_1 = x_1^{-1/2}$. The Schur function is given by $s_{(\nu_1, \nu_2)}(z_1, z_2) = (z_1^{d+1} - z_1^{-d-1})/(z_1 - z_1^{-1})$ with $d = \nu_1 - \nu_2$. From (4.4) and (3.37) we have

$$ \nu_i^{(1)}(1, l) = i - \frac{2}{s_{(i)}s_{(l)}} \sum_{k=1}^{i} k s_{(i+l-k,k)} = i \frac{(1+x_1^{1+i})(1+x_1^{1+l})}{(1-x_1^{1+i})(1-x_1^{1+l})} - \frac{2x_1(1-x_1^i)(1+x_1^{1+l})}{(1-x_1)(1-x_1^{1+i})(1-x_1^{1+l})} \quad (4.10) $$

for $i \le l$. In view of the symmetry $\nu_i^{(1)}(1, l) = \nu_l^{(1)}(1, i)$ (see (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,

$$ \nu_i^{(1)}(1, 1) = \frac{(1 - x_1)(1 + x_1^{1+i})}{(1 + x_1)(1 - x_1^{1+i})} \quad (4.11) $$
The equation (4.2) coincides with (3.34) under the identification $h_i^{(a)} = \sum_{j \ge 1} \min(i, j) y_j^{(a)} \nu_j^{(a)}$. Therefore it has a solution

$$ \nu_i^{(a)}(r, l) = \delta_{ar} \min(i, l) - \sum_{b=1}^{n-1} C_{ab} h_i^{(b)}(r, l) \quad (4.4) $$

$$ = (y_i^{(a)})^{-1}(-h_{i-1}^{(a)}(r, l) + 2h_i^{(a)}(r, l) - h_{i+1}^{(a)}(r, l)). \quad (4.5) $$

From the first formula (4.4) with $l = 1$ and (3.53) we have

$$ \sigma_i^{(a)} = \sum_{k=1}^{n-1} \nu_i^{(a)}(k, 1). \quad (4.6) $$

Combining (4.3) and (4.6) we obtain the effective speed

$$ v_i^{(a)}(r, l) = \frac{\nu_i^{(a)}(r, l)}{\sum_{k=1}^{n-1} \nu_i^{(a)}(k, 1)}, \quad (4.7) $$

in terms of $\nu_i^{(a)}(r, l)$ in (4.4)–(4.5) with $h_i^{(a)}(r, l)$ given by (3.36). The result leads to the relation

$$ \sum_{r=1}^{n-1} v_i^{(a)}(r, 1) = 1, \quad (4.8) $$

which is consistent with Remark 2.1.

Among the family of time evolutions $T_l^{(r)}$, typical ones are $T_{\infty}^{(1)}, \dots, T_{\infty}^{(n-1)}$. The corresponding limit of the numerator $\lim_{l \to \infty} \nu_i^{(a)}(r, l)$ in (4.7) can be evaluated compactly as

$$ \begin{aligned} \nu_i^{(a)}(r, \infty) &= \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}h_i^{(b)}(r, \infty) = \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}h_{\infty}^{(r)}(b, i) \\ &\stackrel{(3.35)}{=} \delta_{ar}i - \sum_{b=1}^{n-1} C_{ab}x_r \frac{\partial}{\partial x_r} \log(e^{-i\varpi_b} Q_i^{(b)}) = x_r \frac{\partial}{\partial x_r} \log \left( \frac{Q_i^{(a-1)} Q_i^{(a+1)}}{(Q_i^{(a)})^2} \right). \end{aligned} \quad (4.9) $$

**Example 4.1.** Consider $n = 2$ case. The variables in (3.7)–(3.9) are related as $z_1 = z_2^{-1} = w_1 = x_1^{-1/2}$. The Schur function is given by $s_{(\nu_1, \nu_2)}(z_1, z_2) = (z_1^{d+1} - z_1^{-d-1})/(z_1 - z_1^{-1})$ with $d = \nu_1 - \nu_2$. From (4.4) and (3.37) we have

$$ \nu_i^{(1)}(1, l) = i - \frac{2}{s_{(i)}s_{(l)}} \sum_{k=1}^{i} k s_{(i+l-k,k)} = i \frac{(1+x_1^{1+i})(1+x_1^{1+l})}{(1-x_1^{1+i})(1-x_1^{1+l})} - \frac{2x_1(1-x_1^i)(1+x_1^{1+l})}{(1-x_1)(1-x_1^{1+i})(1-x_1^{1+l})} \quad (4.10) $$

for $i \le l$. In view of the symmetry $\nu_i^{(1)}(1, l) = \nu_l^{(1)}(1, i)$ (see (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,

$$ \nu_i^{(1)}(1, 1) = \frac{(1 - x_1)(1 + x_1^{1+i})}{(1 + x_1)(1 - x_1^{1+i})} \quad (4.11) $$
holds for $i \ge 1$. Now the effective speed (4.7) is evaluated as

$$v_i^{(1)}(1, l) = \frac{\nu_i^{(1)}(1, l)}{\nu_i^{(1)}(1, 1)} = \frac{1 + x_1^{1+m}}{1 - x_1^{1+m}} \left( j \frac{1 + x_1}{1 - x_1} - \frac{2x_1(1 + x_1)(1 - x^j)}{(1 - x_1)^2(1 + x_1^{1+j})} \right). \quad (4.12)$$

This reproduces [11, eq.(D.2)] with $a = z = x_1$. Let us also check (4.9) which says

$$\nu_i^{(1)}(1, \infty) = -2x_1 \frac{\partial}{\partial x_1} \log \frac{x_1^{\frac{i+1}{2}} - x_1^{-\frac{i+1}{2}}}{x_1^{\frac{1}{2}} - x_1^{-\frac{1}{2}}} = i \frac{1 + x_1^{i+1}}{1 - x_1^{i+1}} - \frac{2x_1(1 - x_1^i)}{(1 - x_1)(1 - x_1^{i+1})}. \quad (4.13)$$

In the regime (3.14) under consideration, this indeed coincides $\lim_{l \to \infty} \nu_i^{(1)}(1, l)$ derived from (4.10).

**Example 4.2.** Consider $n = 3$ case. Explicit formulas for the effective speed $v_i^{(a)}(r, l)$ are messy for generic $z_1, z_2, z_3$ and $l$. Here we treat $v_i^{(a)}(r, \infty)$ for the one parameter specialization $(z_1, z_2, z_3) = (q^{-1}, 1, q)$ with $0 < q < 1$, which satisfies (3.14). This is the so called *principal specialization* reducing the Schur function to the $q$-dimension as $s_{(\lambda_1, \lambda_2, \lambda_3)} = \frac{q^{\lambda_3 - \lambda_1}(1 - q^{\lambda_1 - \lambda_2 + 1})(1 - q^{\lambda_1 - \lambda_3 + 2})(1 - q^{\lambda_2 - \lambda_3 + 1})}{(1 - q)^2(1 - q^2)}$. To calculate the denominator of (4.7) we substitute

$$\begin{aligned} h_i^{(1)}(1, 1) &= \frac{s_{(i,1)}}{s_{(i)}s_{(1)}}, \quad h_i^{(2)}(1, 1) = \frac{s_{(i-1, i-1)}}{s_{(i,i)}s_{(1)}}, \quad h_i^{(1)}(2, 1) = \frac{s_{(i-1)}}{s_{(1,1)}s_{(i)}}, \\ h_i^{(2)}(2, 1) &= \frac{s_{(i,i-1)}}{s_{(i,i)}s_{(1,1)}} \end{aligned} \quad (4.14)$$

into (4.4), i.e., $\nu_i^{(1)}(1, 1) = 1 - 2h_i^{(1)}(1, 1) + h_i^{(2)}(1, 1)$, $\nu_i^{(2)}(1, 1) = h_i^{(1)}(1, 1) - 2h_i^{(1)}(1, 1)$, $\nu_i^{(1)}(2, 1) = -2h_i^{(1)}(2, 1) + h_i^{(2)}(2, 1)$, $\nu_i^{(2)}(2, 1) = 1 + h_i^{(1)}(2, 1) - 2h_i^{(1)}(2, 1)$ to find

$$\nu_i^{(1)}(1, 1) = \nu_i^{(2)}(2, 1) = \frac{(1 - q)^2(1 + q^{i+1})(1 - q^{i+3})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}, \quad (4.15)$$

$$\nu_i^{(1)}(2, 1) = \nu_i^{(2)}(1, 1) = \frac{q(1 - q)^2(1 - q^i)(1 + q^{i+2})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.16)$$

Thus the denominator of (4.7) is factorized as

$$\nu_i^{(1)}(1, 1) + \nu_i^{(1)}(2, 1) = \nu_i^{(2)}(1, 1) + \nu_i^{(2)}(2, 1) = \frac{(1 - q)(1 - q^2)(1 - q^{2i+3})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.17)$$

Similar calculations of (4.4)–(4.5) lead to

$$\nu_i^{(1)}(1, \infty) = \nu_i^{(2)}(2, \infty) = i \frac{1 + q^{i+1} + q^{i+2}}{(1 - q^{i+1})(1 - q^{i+2})} - \frac{q(1 - q^i)(2 + 3q + q^{i+3})}{(1 - q^2)(1 - q^{i+1})(1 - q^{i+2})}, \quad (4.18)$$

$$\nu_i^{(1)}(2, \infty) = \nu_i^{(2)}(1, \infty) = -i \frac{q^{i+1}(1 + q + q^{i+2})}{(1 - q^{i+1})(1 - q^{i+2})} + \frac{q(1 - q^i)(1 + 3q^{i+2} + 2q^{i+3})}{(1 - q^2)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.19)$$

These are related by $\nu_i^{(1)}(2, \infty) = -\nu_i^{(1)}(1, \infty) \Big|_{q \to q^{-1}}$. Finally the effective speed is given by

$$v_i^{(1)}(1, \infty) = v_i^{(2)}(2, \infty) = i \frac{(1 - q^3)(1 + q^{i+1} + q^{i+2})}{(1 - q)(1 - q^2)(1 - q^{2i+3})} - \frac{q(1 - q^3)(1 - q^i)(2 + 3q + q^{i+3})}{(1 - q)(1 - q^2)^2(1 - q^{2i+3})}, \quad (4.20)$$

$$\begin{aligned} v_i^{(1)}(2, \infty) = v_i^{(2)}(1, \infty) = & -i \frac{q^{i+1}(1 - q^3)(1 + q + q^{i+2})}{(1 - q)(1 - q^2)(1 - q^{2i+3})} \\ & + \frac{q(1 - q^3)(1 - q^i)(1 + 3q^{i+2} + 2q^{i+3})}{(1 - q)(1 - q^2)^2(1 - q^{2i+3})}. \end{aligned} \quad (4.21)$$
holds for $i \ge 1$. Now the effective speed (4.7) is evaluated as

$$v_i^{(1)}(1, l) = \frac{\nu_i^{(1)}(1, l)}{\nu_i^{(1)}(1, 1)} = \frac{1 + x_1^{1+m}}{1 - x_1^{1+m}} \left( j \frac{1 + x_1}{1 - x_1} - \frac{2x_1(1 + x_1)(1 - x^j)}{(1 - x_1)^2(1 + x_1^{1+j})} \right). \quad (4.12)$$

This reproduces [11, eq.(D.2)] with $a = z = x_1$. Let us also check (4.9) which says

$$\nu_i^{(1)}(1, \infty) = -2x_1 \frac{\partial}{\partial x_1} \log \frac{x_1^{\frac{i+1}{2}} - x_1^{-\frac{i+1}{2}}}{x_1^{\frac{1}{2}} - x_1^{-\frac{1}{2}}} = i \frac{1 + x_1^{i+1}}{1 - x_1^{i+1}} - \frac{2x_1(1 - x_1^i)}{(1 - x_1)(1 - x_1^{i+1})}. \quad (4.13)$$

In the regime (3.14) under consideration, this indeed coincides $\lim_{l \to \infty} \nu_i^{(1)}(1, l)$ derived from (4.10).

**Example 4.2.** Consider $n = 3$ case. Explicit formulas for the effective speed $v_i^{(a)}(r, l)$ are messy for generic $z_1, z_2, z_3$ and $l$. Here we treat $v_i^{(a)}(r, \infty)$ for the one parameter specialization $(z_1, z_2, z_3) = (q^{-1}, 1, q)$ with $0 < q < 1$, which satisfies (3.14). This is the so called *principal specialization* reducing the Schur function to the $q$-dimension as $s_{(\lambda_1, \lambda_2, \lambda_3)} = \frac{q^{\lambda_3 - \lambda_1}(1 - q^{\lambda_1 - \lambda_2 + 1})(1 - q^{\lambda_1 - \lambda_3 + 2})(1 - q^{\lambda_2 - \lambda_3 + 1})}{(1 - q)^2(1 - q^2)}$. To calculate the denominator of (4.7) we substitute

$$\begin{aligned} h_i^{(1)}(1, 1) &= \frac{s_{(i,1)}}{s_{(i)}s_{(1)}}, \quad h_i^{(2)}(1, 1) = \frac{s_{(i-1, i-1)}}{s_{(i,i)}s_{(1)}}, \quad h_i^{(1)}(2, 1) = \frac{s_{(i-1)}}{s_{(1,1)}s_{(i)}}, \\ h_i^{(2)}(2, 1) &= \frac{s_{(i,i-1)}}{s_{(i,i)}s_{(1,1)}} \end{aligned} \quad (4.14)$$

into (4.4), i.e., $\nu_i^{(1)}(1, 1) = 1 - 2h_i^{(1)}(1, 1) + h_i^{(2)}(1, 1)$, $\nu_i^{(2)}(1, 1) = h_i^{(1)}(1, 1) - 2h_i^{(1)}(1, 1)$, $\nu_i^{(1)}(2, 1) = -2h_i^{(1)}(2, 1) + h_i^{(2)}(2, 1)$, $\nu_i^{(2)}(2, 1) = 1 + h_i^{(1)}(2, 1) - 2h_i^{(1)}(2, 1)$ to find

$$\nu_i^{(1)}(1, 1) = \nu_i^{(2)}(2, 1) = \frac{(1 - q)^2(1 + q^{i+1})(1 - q^{i+3})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}, \quad (4.15)$$

$$\nu_i^{(1)}(2, 1) = \nu_i^{(2)}(1, 1) = \frac{q(1 - q)^2(1 - q^i)(1 + q^{i+2})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.16)$$

Thus the denominator of (4.7) is factorized as

$$\nu_i^{(1)}(1, 1) + \nu_i^{(1)}(2, 1) = \nu_i^{(2)}(1, 1) + \nu_i^{(2)}(2, 1) = \frac{(1 - q)(1 - q^2)(1 - q^{2i+3})}{(1 - q^3)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.17)$$

Similar calculations of (4.4)–(4.5) lead to

$$\nu_i^{(1)}(1, \infty) = \nu_i^{(2)}(2, \infty) = i \frac{1 + q^{i+1} + q^{i+2}}{(1 - q^{i+1})(1 - q^{i+2})} - \frac{q(1 - q^i)(2 + 3q + q^{i+3})}{(1 - q^2)(1 - q^{i+1})(1 - q^{i+2})}, \quad (4.18)$$

$$\nu_i^{(1)}(2, \infty) = \nu_i^{(2)}(1, \infty) = -i \frac{q^{i+1}(1 + q + q^{i+2})}{(1 - q^{i+1})(1 - q^{i+2})} + \frac{q(1 - q^i)(1 + 3q^{i+2} + 2q^{i+3})}{(1 - q^2)(1 - q^{i+1})(1 - q^{i+2})}. \quad (4.19)$$

These are related by $\nu_i^{(1)}(2, \infty) = -\nu_i^{(1)}(1, \infty) \rvert_{q \to q^{-1}}$. Finally the effective speed is given by

$$v_i^{(1)}(1, \infty) = v_i^{(2)}(2, \infty) = i \frac{(1 - q^3)(1 + q^{i+1} + q^{i+2})}{(1 - q)(1 - q^2)(1 - q^{2i+3})} - \frac{q(1 - q^3)(1 - q^i)(2 + 3q + q^{i+3})}{(1 - q)(1 - q^2)^2(1 - q^{2i+3})}, \quad (4.20)$$

$$\begin{aligned} v_i^{(1)}(2, \infty) = v_i^{(2)}(1, \infty) = -i \frac{q^{i+1}(1 - q^3)(1 + q + q^{i+2})}{(1 - q)(1 - q^2)(1 - q^{2i+3})} &+ \frac{q(1 - q^3)(1 - q^i)(1 + 3q^{i+2} + 2q^{i+3})}{(1 - q)(1 - q^2)^2(1 - q^{2i+3})}. \end{aligned} \quad (4.21)$$
Again they are related by $v_i^{(1)}(2, \infty) = -v_i^{(1)}(1, \infty) \big|_{q \to q^{-1}}$. The speeds (4.20) and (4.21) are positive for $0 < q < 1$.

## 4.2 Matrix form

It is convenient to formulate the results in the previous subsection in a matrix form whose indices range over $\mathcal{I}$ in (2.1). We introduce the matrices

$$I = (\delta_{ab}\delta_{ij})_{(a,i),(b,j)\in\mathcal{I}}, \qquad M = (C_{ab}\min(i,j))_{(a,i),(b,j)\in\mathcal{I}}, \qquad (4.22)$$

$$\mathbf{y} = (\delta_{ab}\delta_{ij}y_i^{(a)})_{(a,i),(b,j)\in\mathcal{I}}, \qquad \mathbf{v} = \mathbf{v}(r,l) = (\delta_{ab}\delta_{ij}v_i^{(a)}(r,l))_{(a,i),(b,j)\in\mathcal{I}} \qquad (4.23)$$

and the vectors

$$\mathbb{I} = (1)_{(a,i)\in\mathcal{I}}, \qquad \kappa = \kappa(r,l) = (\delta_{ar}\min(i,l))_{(a,i)\in\mathcal{I}}, \qquad (4.24)$$

$$\rho = (\rho_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad \sigma = (\sigma_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad y = (y_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad (4.25)$$

$$v = v(r, l) = (v_i^{(a)}(r, l))_{(a,i)\in\mathcal{I}}, \qquad \nu = \nu(r, l) = \sigma * v := (\sigma_i^{(a)} v_i^{(a)}(r, l))_{(a,i)\in\mathcal{I}}, \qquad (4.26)$$

where the component-wise (or Hadamard) product $*$ is commutative and will be used also for other vectors. For example we have $\rho = y * \sigma$.<fnref n="8" /> In view of (4.1), the vector $\kappa(r, l)$ is the bare speed of solitons under the time evolution $T_l^{(r)}$. Note that $I$ and $\mathbb{I}$ are different.

Now the equations (3.17) and (4.2) regarded as the ones for $\sigma$ and $\nu$ are presented as follows:

$$\sigma + M\rho = (I + My)\sigma = \mathbb{I}, \qquad (4.27)$$

$$\sigma * v(r, l) + M(\rho * v(r, l)) = (I + My)\nu(r, l) = \kappa(r, l). \qquad (4.28)$$

In general we define the “dressing” $\eta^{\text{dr}}$ of a vector $\eta$ by

$$\eta^{\text{dr}} = \eta - My\eta^{\text{dr}}, \qquad (4.29)$$

namely, $\eta^{\text{dr}} = (I + My)^{-1}\eta$. Then (4.27) and (4.28) are rephrased as

$$\sigma = \mathbb{I}^{\text{dr}}, \qquad \nu(r, l) = \mathbb{I}^{\text{dr}} * v(r, l) = \kappa(r, l)^{\text{dr}}. \qquad (4.30)$$

In particular the effective speed is given by

$$v(r, l) = \frac{\kappa(r, l)^{\text{dr}}}{\mathbb{I}^{\text{dr}}} = \frac{(I + My)^{-1}\kappa(r, l)}{(I + My)^{-1}\mathbb{I}}, \qquad (4.31)$$

where the division is component-wise. From $\rho = y * \sigma$ we also have

$$\rho = (\mathbf{y}^{-1} + M)^{-1}\mathbb{I}, \qquad \rho * v(r, l) = (\mathbf{y}^{-1} + M)^{-1}\kappa(r, l). \qquad (4.32)$$

**Remark 4.3.** Regard the effective speed (4.31) as a function of $\{y_i^{(a)} | (a, i) \in \mathcal{I}\}$. Then from Cramer's formula and the fact that $\mathbf{y}$ (4.23) is diagonal, it follows that $v_s^{(p)}(r, l)$ does not depend on the variable $y_s^{(p)}$ for any $(p, s) \in \mathcal{I}$. This property will be used to ensure (5.5).

**Remark 4.4.** In the dilute limit explained in Remark 3.2, one has $y_i^{(a)} = 0$ and $\rho_i^{(a)} = 0$ from (3.24) and (3.20) for all $(a, i) \in \mathcal{I}$. Therefore the dressing (4.29) trivializes into $\eta^{\text{dr}} = \eta$. This brings the effective speed back to the bare speed, i.e., $v(r, l) = \kappa(r, l)$ in (4.31).

<footnote n="8">⁸The $y$ in (4.25) is a vector which should be distinguished from the diagonal matrix $\mathbf{y}$ in (4.23).</footnote>
Again they are related by $v_i^{(1)}(2, \infty) = -v_i^{(1)}(1, \infty) \big|_{q \to q^{-1}}$. The speeds (4.20) and (4.21) are positive for $0 < q < 1$.

## 4.2 Matrix form

It is convenient to formulate the results in the previous subsection in a matrix form whose indices range over $\mathcal{I}$ in (2.1). We introduce the matrices

$$I = (\delta_{ab}\delta_{ij})_{(a,i),(b,j)\in\mathcal{I}}, \qquad M = (C_{ab}\min(i,j))_{(a,i),(b,j)\in\mathcal{I}}, \qquad (4.22)$$

$$\mathbf{y} = (\delta_{ab}\delta_{ij}y_i^{(a)})_{(a,i),(b,j)\in\mathcal{I}}, \qquad \mathbf{v} = \mathbf{v}(r,l) = (\delta_{ab}\delta_{ij}v_i^{(a)}(r,l))_{(a,i),(b,j)\in\mathcal{I}}, \qquad (4.23)$$

and the vectors

$$\mathbb{I} = (1)_{(a,i)\in\mathcal{I}}, \qquad \kappa = \kappa(r,l) = (\delta_{ar}\min(i,l))_{(a,i)\in\mathcal{I}}, \qquad (4.24)$$

$$\rho = (\rho_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad \sigma = (\sigma_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad y = (y_i^{(a)})_{(a,i)\in\mathcal{I}}, \qquad (4.25)$$

$$v = v(r, l) = (v_i^{(a)}(r, l))_{(a,i)\in\mathcal{I}}, \qquad \nu = \nu(r, l) = \sigma * v := (\sigma_i^{(a)} v_i^{(a)}(r, l))_{(a,i)\in\mathcal{I}}, \qquad (4.26)$$

where the component-wise (or Hadamard) product $*$ is commutative and will be used also for other vectors. For example we have $\rho = y * \sigma$.<fnref n="8" /> In view of (4.1), the vector $\kappa(r, l)$ is the bare speed of solitons under the time evolution $T_l^{(r)}$. Note that $I$ and $\mathbb{I}$ are different.

Now the equations (3.17) and (4.2) regarded as the ones for $\sigma$ and $\nu$ are presented as follows:

$$\sigma + M\rho = (I + My)\sigma = \mathbb{I}, \qquad (4.27)$$

$$\sigma * v(r, l) + M(\rho * v(r, l)) = (I + My)\nu(r, l) = \kappa(r, l). \qquad (4.28)$$

In general we define the “dressing” $\eta^{\text{dr}}$ of a vector $\eta$ by

$$\eta^{\text{dr}} = \eta - My\eta^{\text{dr}}, \qquad (4.29)$$

namely, $\eta^{\text{dr}} = (I + My)^{-1}\eta$. Then (4.27) and (4.28) are rephrased as

$$\sigma = \mathbb{I}^{\text{dr}}, \qquad \nu(r, l) = \mathbb{I}^{\text{dr}} * v(r, l) = \kappa(r, l)^{\text{dr}}. \qquad (4.30)$$

In particular the effective speed is given by

$$v(r, l) = \frac{\kappa(r, l)^{\text{dr}}}{\mathbb{I}^{\text{dr}}} = \frac{(I + My)^{-1}\kappa(r, l)}{(I + My)^{-1}\mathbb{I}}, \qquad (4.31)$$

where the division is component-wise. From $\rho = y * \sigma$ we also have

$$\rho = (\mathbf{y}^{-1} + M)^{-1}\mathbb{I}, \qquad \rho * v(r, l) = (\mathbf{y}^{-1} + M)^{-1}\kappa(r, l). \qquad (4.32)$$

**Remark 4.3.** Regard the effective speed (4.31) as a function of $\{y_i^{(a)} | (a, i) \in \mathcal{I}\}$. Then from Cramer's formula and the fact that $\mathbf{y}$ (4.23) is diagonal, it follows that $v_s^{(p)}(r, l)$ does not depend on the variable $y_s^{(p)}$ for any $(p, s) \in \mathcal{I}$. This property will be used to ensure (5.5).

**Remark 4.4.** In the dilute limit explained in Remark 3.2, one has $y_i^{(a)} = 0$ and $\rho_i^{(a)} = 0$ from (3.24) and (3.20) for all $(a, i) \in \mathcal{I}$. Therefore the dressing (4.29) trivializes into $\eta^{\text{dr}} = \eta$. This brings the effective speed back to the bare speed, i.e., $v(r, l) = \kappa(r, l)$ in (4.31).

<footnote n="8">⁸The $y$ in (4.25) is a vector which should be distinguished from the diagonal matrix $\mathbf{y}$ in (4.23).</footnote>
## 4.3 Stationary currents

Let us define the (density of) color a stationary current under the time evolution $T_l^{(r)}$ as

$$J^{(a)}(r, l) = \sum_{i \ge 1} i J_i^{(a)}(r, l), \quad J_i^{(a)}(r, l) = \rho_i^{(a)} v_i^{(a)}(r, l). \qquad (4.33)$$

In view of Remark 2.8, this is a component of the weight current $-\sum_{a=1}^{n-1} J^{(a)}(r, l)\alpha_a$ transported by solitons. It is expressed as a quadratic form of the bare velocities as

$$J^{(a)}(r, l) = \kappa(a, \infty)^{\mathrm{T}} \rho * v(r, l) = \kappa(a, \infty)^{\mathrm{T}} (\mathbf{y}^{-1} + M)^{-1} \kappa(r, l) \qquad (4.34)$$

due to (4.32), where $\eta^{\mathrm{T}}$ signifies a row vector obtained by the transpose of a column vector $\eta$.

As an additional remark, one can also express $J^{(a)}(r, l)$ in a form that generalizes [11, eq.(4.4)] naturally. Let $m_i^{(a)}$ be the number of soliton $S_i^{(a)}$ (2.39) in the system. We assume $m_i^{(a)} = 0$ for $i$ sufficiently large. Using (3.3) one can rewrite the speed equation (4.2) as

$$p_i^{(a)} v_i^{(a)} + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) m_j^{(b)} v_j^{(b)} = L \delta_{ar} \min(i, l). \qquad (4.35)$$

Let $\hat{\mathcal{I}} = \{(a, i, \alpha) \mid (a, i) \in \mathcal{I}, \alpha \in [1, m_i^{(a)}]\}$ be an extension of the index set $\mathcal{I}$ in (2.1). The extra component $\alpha$ may be viewed as labeling individual solitons $S_i^{(a)}$ of color a and length $i$. The cardinality of $\hat{\mathcal{I}}$ is $\sum_{(a,i) \in \mathcal{I}} m_i^{(a)}$ which is finite by the assumption. With a vector $\eta = (\eta_i^{(a)})_{(a,i) \in \mathcal{I}}$ we associate a similar extension $\hat{\eta} = (\hat{\eta}_i^{(a)})_{(a,i,\alpha) \in \hat{\mathcal{I}}}$. It repeats the same element $\eta_i^{(a)}$ for $m_i^{(a)}$ times which is the size of the $(a, i)$ block. Introduce the matrix<fnref n="9" />

$$B = (p_i^{(a)} \delta_{ab} \delta_{ij} \delta_{\alpha\beta} + C_{ab} \min(i, j))_{(a,i,\alpha),(b,j,\beta) \in \hat{\mathcal{I}}}. \qquad (4.36)$$

These definitions enable us to present the speed equation (4.35) as

$$B \hat{v}(r, l) = L \hat{\kappa}(r, l). \qquad (4.37)$$

Similarly the current (4.33) is expressed as a quadratic form:

$$J^{(a)}(r, l) = \frac{1}{L} \sum_{i \ge 1} im_i^{(a)} v_i^{(a)}(r, l) = \frac{1}{L} \hat{\kappa}(a, \infty)^{\mathrm{T}} \hat{v}(r, l) = \hat{\kappa}(a, \infty)^{\mathrm{T}} B^{-1} \hat{\kappa}(r, l). \qquad (4.38)$$

This formula is a natural generalization of the $n=2$ case in [11, eq.(4.4)], where a further connection to the period matrix of the tropical Riemann theta function was explored. At present, time averaged current of the $n$-color cBBS on a circle is yet to be obtained. However we expect the result should coincide with (4.38) supported by the $n=2$ case established in [11].

## 4.4 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

<footnote n="9">⁹B appeared in [44, eq.(4.6)] for a non-complete multi-color BBS with a periodic boundary condition with the vacancy $p_i^{(a)}$ different from (3.3) reflecting the “non-completeness”. Of course this B should not be confused with B in (2.21).</footnote>
## 4.3 Stationary currents

Let us define the (density of) color a stationary current under the time evolution $T_l^{(r)}$ as

$$J^{(a)}(r, l) = \sum_{i \ge 1} i J_i^{(a)}(r, l), \qquad J_i^{(a)}(r, l) = \rho_i^{(a)} v_i^{(a)}(r, l). \qquad (4.33)$$

In view of Remark 2.8, this is a component of the weight current $-\sum_{a=1}^{n-1} J^{(a)}(r, l)\alpha_a$ transported by solitons. It is expressed as a quadratic form of the bare velocities as

$$J^{(a)}(r, l) = \kappa(a, \infty)^{\mathrm{T}} \rho * v(r, l) = \kappa(a, \infty)^{\mathrm{T}} (\mathbf{y}^{-1} + M)^{-1} \kappa(r, l) \qquad (4.34)$$

due to (4.32), where $\eta^{\mathrm{T}}$ signifies a row vector obtained by the transpose of a column vector $\eta$.

As an additional remark, one can also express $J^{(a)}(r, l)$ in a form that generalizes [11, eq.(4.4)] naturally. Let $m_i^{(a)}$ be the number of soliton $S_i^{(a)}$ (2.39) in the system. We assume $m_i^{(a)} = 0$ for $i$ sufficiently large. Using (3.3) one can rewrite the speed equation (4.2) as

$$p_i^{(a)} v_i^{(a)} + \sum_{(b,j) \in \mathcal{I}} C_{ab} \min(i, j) m_j^{(b)} v_j^{(b)} = L \delta_{ar} \min(i, l). \qquad (4.35)$$

Let $\hat{\mathcal{I}} = \{(a, i, \alpha) \mid (a, i) \in \mathcal{I}, \alpha \in [1, m_i^{(a)}]\}$ be an extension of the index set $\mathcal{I}$ in (2.1). The extra component $\alpha$ may be viewed as labeling individual solitons $S_i^{(a)}$ of color a and length $i$. The cardinality of $\hat{\mathcal{I}}$ is $\sum_{(a,i) \in \mathcal{I}} m_i^{(a)}$ which is finite by the assumption. With a vector $\eta = (\eta_i^{(a)})_{(a,i) \in \mathcal{I}}$ we associate a similar extension $\hat{\eta} = (\hat{\eta}_i^{(a)})_{(a,i,\alpha) \in \hat{\mathcal{I}}}$. It repeats the same element $\eta_i^{(a)}$ for $m_i^{(a)}$ times which is the size of the $(a, i)$ block. Introduce the matrix<fnref n="9" />

$$B = (p_i^{(a)} \delta_{ab} \delta_{ij} \delta_{\alpha\beta} + C_{ab} \min(i, j))_{(a,i,\alpha),(b,j,\beta) \in \hat{\mathcal{I}}} \qquad (4.36)$$

These definitions enable us to present the speed equation (4.35) as

$$B \hat{v}(r, l) = L \hat{\kappa}(r, l). \qquad (4.37)$$

Similarly the current (4.33) is expressed as a quadratic form:

$$J^{(a)}(r, l) = \frac{1}{L} \sum_{i \ge 1} im_i^{(a)} v_i^{(a)}(r, l) = \frac{1}{L} \hat{\kappa}(a, \infty)^{\mathrm{T}} \hat{v}(r, l) = \hat{\kappa}(a, \infty)^{\mathrm{T}} B^{-1} \hat{\kappa}(r, l). \qquad (4.38)$$

This formula is a natural generalization of the $n=2$ case in [11, eq.(4.4)], where a further connection to the period matrix of the tropical Riemann theta function was explored. At present, time averaged current of the $n$-color cBBS on a circle is yet to be obtained. However we expect the result should coincide with (4.38) supported by the $n=2$ case established in [11].

## 4.4 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

<footnote n="9">⁹B appeared in [44, eq.(4.6)] for a non-complete multi-color BBS with a periodic boundary condition with the vacancy $p_i^{(a)}$ different from (3.3) reflecting the “non-completeness”. Of course this B should not be confused with B in (2.21).</footnote>
that the soliton current $J_i^{(a)}(r, l)$ in (4.33) associated with $S_i^{(a)}$ obeys the continuity equation $\partial_t \rho_i^{(a)} + \partial_x J_i^{(a)}(r, l) = 0$ for any time evolution $T_l^{(r)}$.<fnref n="10" /> In the matrix notation it reads

$$ \partial_t \rho + \partial_x (\rho * v(r, l)) = 0. \qquad (4.39) $$

Substituting (4.32) and using $\partial_\alpha (\mathbf{y}^{-1} + M)^{-1} = \beta \partial_\alpha \mathbf{y}(\mathbf{y}^{-1} + M)^{-1}$ with $\beta = (\mathbf{y}^{-1} + M)^{-1} \mathbf{y}^{-2}$ for $\alpha = t, x$, we get $\partial_t \mathbf{y}(\mathbf{y}^{-1} + M)^{-1} \mathbf{I} + \partial_x \mathbf{y}(\mathbf{y}^{-1} + M)^{-1} \kappa(r, l) = 0$. By utilizing (4.32) again and noting that $\mathbf{y}$ and $\mathbf{v}$ in (4.23) are diagonal, the result becomes

$$ \partial_t y + \mathbf{v}(r, l) \partial_x y = 0. \qquad (4.40) $$

This is a collection of separated equations for each $(a, i) \in \mathcal{I}$ meaning that $y_i^{(a)}$'s are the normal modes of the generalized hydrodynamics [10].

**Remark 4.5.** Replacing $\rho = (\rho_i^{(a)})_{(a,i)\in\mathcal{I}}$ in (4.39) by $(f_i^{(a)}(y_i^{(a)})\rho_i^{(a)})_{(a,i)\in\mathcal{I}}$ with any function $f_i^{(a)}$ leads to the same equation (4.40). In particular, taking $f_i^{(a)}(y_i^{(a)}) = y_i^{(a)}$ leads to the conservation of the hole current $\partial_t \sigma + \partial_x (\sigma * v(r, l)) = 0$.

# 5 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, \dots, z_n)$ as explained in (3.4)–(3.9). We prepare the system initially in the i.i.d. distribution with fugacity $(z_{1L}, \dots, z_{nL})$ in the left domain $x < 0$ and a different fugacity $(z_{1R}, \dots, z_{nR})$ in the right domain $x > 0$. Then at time $t = 0$, we begin running a dynamics $T_l^{(r)}$ for some $(r, l) \in \mathcal{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 the recent literature. 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 using simulations (Sec. 5.3).

## 5.1 Ballistic picture

The first approximation we employ is the ballistic scaling which implies that the normal mode $y_i^{(a)}(x, t)$ (see Sec. 4.4) depends on $x, t$ only through the ray variable $\zeta = x/t$ as $y_i^{(a)}(\zeta)$. Then (4.40) is reduced to

$$ (\zeta - v_i^{(a)}(r, l)) \partial_\zeta y_i^{(a)}(\zeta) = 0 \quad ((a, i) \in \mathcal{I}). \qquad (5.1) $$

Since the set of velocities $v_i^{(a)}(r, l)$ is discrete in our setting, it follows that $y_i^{(a)}(\zeta)$, hence the local state of the system itself, is piecewise constant and uniform when observed via the variable $\zeta$.¹¹¹² Thus the plot of any local quantity, typically density of solitons $S_i^{(a)}$, against $\zeta$

<footnote n="10">¹⁰The associated time variable could also be denoted by $t_i^{(r)}$, which is however avoided for notational simplicity.</footnote>

<footnote n="11">¹¹One can show that $y_i^{(a)}(\zeta)$ can change discontinuously across $\zeta = v_i^{(a)}(r, l) \pm 0$ without violating the current conservation. The proof is formally identical with [11, Sec. 5.3].</footnote>

<footnote n="12">¹²A short comment is here in order about the word state. So far we have been using state to represent a single microscopic ball configuration. However, in the present GHD context, the (local) state of the system must be understood in a thermodynamic or probabilistic sense, as a measure over microscopic configurations.</footnote>
that the soliton current $J_i^{(a)}(r, l)$ in (4.33) associated with $S_i^{(a)}$ obeys the continuity equation $\partial_t \rho_i^{(a)} + \partial_x J_i^{(a)}(r, l) = 0$ for any time evolution $T_l^{(r)}$.<fnref n="10" /> In the matrix notation it reads

$$ \partial_t \rho + \partial_x (\rho * v(r, l)) = 0. \qquad (4.39) $$

Substituting (4.32) and using $\partial_\alpha (\mathbf{y}^{-1} + M)^{-1} = \beta \partial_\alpha \mathbf{y}(\mathbf{y}^{-1} + M)^{-1}$ with $\beta = (\mathbf{y}^{-1} + M)^{-1} \mathbf{y}^{-2}$ for $\alpha = t, x$, we get $\partial_t \mathbf{y}(\mathbf{y}^{-1} + M)^{-1} \mathbf{I} + \partial_x \mathbf{y}(\mathbf{y}^{-1} + M)^{-1} \kappa(r, l) = 0$. By utilizing (4.32) again and noting that $\mathbf{y}$ and $\mathbf{v}$ in (4.23) are diagonal, the result becomes

$$ \partial_t y + \mathbf{v}(r, l) \partial_x y = 0. \qquad (4.40) $$

This is a collection of separated equations for each $(a, i) \in \mathcal{I}$ meaning that $y_i^{(a)}$'s are the normal modes of the generalized hydrodynamics [10].

**Remark 4.5.** Replacing $\rho = (\rho_i^{(a)})_{(a,i)\in\mathcal{I}}$ in (4.39) by $(f_i^{(a)}(y_i^{(a)})\rho_i^{(a)})_{(a,i)\in\mathcal{I}}$ with any function $f_i^{(a)}$ leads to the same equation (4.40). In particular, taking $f_i^{(a)}(y_i^{(a)}) = y_i^{(a)}$ leads to the conservation of the hole current $\partial_t \sigma + \partial_x (\sigma * v(r, l)) = 0$.

# 5 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, \dots, z_n)$ as explained in (3.4)–(3.9). We prepare the system initially in the i.i.d. distribution with fugacity $(z_{1L}, \dots, z_{nL})$ in the left domain $x < 0$ and a different fugacity $(z_{1R}, \dots, z_{nR})$ in the right domain $x > 0$. Then at time $t = 0$, we begin running a dynamics $T_l^{(r)}$ for some $(r, l) \in \mathcal{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 the recent literature. 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 using simulations (Sec. 5.3).

## 5.1 Ballistic picture

The first approximation we employ is the ballistic scaling which implies that the normal mode $y_i^{(a)}(x, t)$ (see Sec. 4.4) depends on $x, t$ only through the ray variable $\zeta = x/t$ as $y_i^{(a)}(\zeta)$. Then (4.40) is reduced to

$$ (\zeta - v_i^{(a)}(r, l)) \partial_\zeta y_i^{(a)}(\zeta) = 0 \quad ((a, i) \in \mathcal{I}). \qquad (5.1) $$

Since the set of velocities $v_i^{(a)}(r, l)$ is discrete in our setting, it follows that $y_i^{(a)}(\zeta)$, hence the local state of the system itself, is piecewise constant and uniform when observed via the variable $\zeta$.¹¹¹² Thus the plot of any local quantity, typically density of solitons $S_i^{(a)}$, against $\zeta$

<footnote n="10">¹⁰The associated time variable could also be denoted by $t_i^{(r)}$, which is however avoided for notational simplicity.</footnote>

<footnote n="11">¹¹One can show that $y_i^{(a)}(\zeta)$ can change discontinuously across $\zeta = v_i^{(a)}(r, l) \pm 0$ without violating the current conservation. The proof is formally identical with [11, Sec. 5.3].</footnote>

<footnote n="12">¹²A short comment is here in order about the word state. So far we have been using state to represent a single microscopic ball configuration. However, in the present GHD context, the (local) state of the system must be understood in a thermodynamic or probabilistic sense, as a measure over microscopic configurations.</footnote>
collapses onto a single (and “static”) curve as $t \to \infty$, which is a collection of plateaux having sharp edges as follows:

<figure><img src="image_9.png" /></figure>

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_i^{(a)}(\zeta)$ can take two values $y_{iL}^{(a)}$ or $y_{iR}^{(a)}$ which are determined from the fugacity ($z_{1L}, \dots, z_{nL}$) or ($z_{1R}, \dots, z_{nR}$) by (3.24) and (3.10). In the mixed inhomogeneous system under consideration, we have the boundary condition $y_i^{(a)}(-\infty) = y_{iL}^{(a)}$ and $y_i^{(a)}(\infty) = y_{iR}^{(a)}$. For the neighboring plateaux $\mathcal{P}$ and $\mathcal{P}'$ as in (5.2), the relation (5.1) indicates that there is a subset $\mathcal{J} \subset \mathcal{I}$ such that the following hold<fnref n="13" />:

$$y_i^{(a)}(\zeta^* - 0) = y_i^{(a)}(\zeta^* + 0) \quad (= y_{iL}^{(a)} \text{ or } y_{iR}^{(a)}) \quad \text{for } (a, i) \in \bar{\mathcal{J}}, \qquad (5.3)$$

$$y_i^{(a)}(\zeta) = \begin{cases} y_{iL}^{(a)} & \zeta < \zeta^* \\ y_{iR}^{(a)} & \zeta > \zeta^* \end{cases} \quad \text{for } (a, i) \in \mathcal{J}, \qquad (5.4)$$

$$v_i^{(a)}(r, l)_{\mathcal{P}} = v_i^{(a)}(r, l)_{\mathcal{P}'} \quad \text{for } (a, i) \in \mathcal{J}. \qquad (5.5)$$

Here $v_i^{(a)}(r, l)_{\mathcal{P}}$, $v_i^{(a)}(r, l)_{\mathcal{P}'}$, denote the effective speed of $S_i^{(a)}$ in $\mathcal{P}$, $\mathcal{P}'$ which can be computed from the respective normal mode $y$ by (4.31). The equality (5.5) can be derived from (5.3) by using Remark 4.3. Put plainly, $\mathcal{J}$ is the set such that (5.1) holds as $0 \times \infty = 0$ if $(a, i) \in \mathcal{J}$ and as (finite nonzero) $\times 0 = 0$ if $(a, i) \in \bar{\mathcal{J}}$. From the consistency with this picture, we expect that the value (5.5) coincides with $\zeta^*$ for all $(a, i) \in \mathcal{J}$.

As we will see, $\mathcal{J}$ is not necessarily a single element set. In fact the point $\zeta^* = 0$ in the figure 2 corresponds to the infinite set $\mathcal{J} = \{(2, i) \mid i \in \mathbb{Z}_{\ge 1}\}$.

Let us label a plateau $\mathcal{P}$ with a subset $\mathcal{P} \subseteq \mathcal{I}$ (denoted by the same symbol) such that

$$y_i^{(a)}(\zeta) = \begin{cases} y_{iL}^{(a)} & (a, i) \in \bar{\mathcal{P}}, \\ y_{iR}^{(a)} & (a, i) \in \mathcal{P}. \end{cases} \qquad (5.6)$$

Then the transition rule (5.3)–(5.4) from $\mathcal{P}$ to $\mathcal{P}'$ across the edge $\zeta = \zeta^*$ is stated as

$$\mathcal{P} \sqcup \mathcal{J} = \mathcal{P}'. \qquad (5.7)$$

The boundary condition means that $\mathcal{P} \to \emptyset$ as $\zeta \to -\infty$ and $\mathcal{P} \to \mathcal{I}$ as $\zeta \to \infty$. Thus a plateau configuration can be viewed as specifying a monotonic growth pattern of the label $\mathcal{P}$. One starts from $\mathcal{P} = \emptyset$ in the far left and proceeds to the right inductively to eventually get $\mathcal{P} = \mathcal{I}$ in the far right. In the inductive step, we assume that a plateau $\mathcal{P}$ and the associated data $\{v_i^{(a)}(r, l)_{\mathcal{P}} \mid (a, i) \in \mathcal{I}\}$ are at hand. Then its right edge $\zeta^*$ as in (5.2) and the set $\mathcal{J}$ which specifies the next plateau $\mathcal{P}'$ by (5.7) are determined as

$$\zeta^* = \min_{(a,i) \in \bar{\mathcal{P}}} \{v_i^{(a)}(r,l)_{\mathcal{P}}\}, \qquad (5.8)$$

$$\mathcal{J} = \{(a, i) \in \bar{\mathcal{P}} \mid v_i^{(a)}(r, l)_{\mathcal{P}} = \zeta^*\}. \qquad (5.9)$$

<footnote n="13">¹³For any subset $\mathcal{K} \subseteq \mathcal{I}$, we let $\bar{\mathcal{K}}$ denote the complement $\mathcal{I} \setminus \mathcal{K}$.</footnote>
collapses onto a single (and “static”) curve as $t \to \infty$, which is a collection of plateaux having sharp edges as follows:

<figure><img src="image_6.png" /></figure>

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_i^{(a)}(\zeta)$ can take two values $y_{iL}^{(a)}$ or $y_{iR}^{(a)}$ which are determined from the fugacity ($z_{1L}, \dots, z_{nL}$) or ($z_{1R}, \dots, z_{nR}$) by (3.24) and (3.10). In the mixed inhomogeneous system under consideration, we have the boundary condition $y_i^{(a)}(-\infty) = y_{iL}^{(a)}$ and $y_i^{(a)}(\infty) = y_{iR}^{(a)}$. For the neighboring plateaux $\mathcal{P}$ and $\mathcal{P}'$ as in (5.2), the relation (5.1) indicates that there is a subset $\mathcal{J} \subset \mathcal{I}$ such that the following hold<fnref n="13" />:

$$y_i^{(a)}(\zeta^* - 0) = y_i^{(a)}(\zeta^* + 0) \quad (= y_{iL}^{(a)} \text{ or } y_{iR}^{(a)}) \quad \text{for } (a, i) \in \bar{\mathcal{J}}, \qquad (5.3)$$

$$y_i^{(a)}(\zeta) = \begin{cases} y_{iL}^{(a)} & \zeta < \zeta^* \\ y_{iR}^{(a)} & \zeta > \zeta^* \end{cases} \quad \text{for } (a, i) \in \mathcal{J}, \qquad (5.4)$$

$$v_i^{(a)}(r, l)_{\mathcal{P}} = v_i^{(a)}(r, l)_{\mathcal{P}'} \quad \text{for } (a, i) \in \mathcal{J}. \qquad (5.5)$$

Here $v_i^{(a)}(r, l)_{\mathcal{P}}$, $v_i^{(a)}(r, l)_{\mathcal{P}'}$, denote the effective speed of $S_i^{(a)}$ in $\mathcal{P}$, $\mathcal{P}'$ which can be computed from the respective normal mode $y$ by (4.31). The equality (5.5) can be derived from (5.3) by using Remark 4.3. Put plainly, $\mathcal{J}$ is the set such that (5.1) holds as $0 \times \infty = 0$ if $(a, i) \in \mathcal{J}$ and as (finite nonzero) $\times 0 = 0$ if $(a, i) \in \bar{\mathcal{J}}$. From the consistency with this picture, we expect that the value (5.5) coincides with $\zeta^*$ for all $(a, i) \in \mathcal{J}$.

As we will see, $\mathcal{J}$ is not necessarily a single element set. In fact the point $\zeta^* = 0$ in the figure 2 corresponds to the infinite set $\mathcal{J} = \{(2, i) \mid i \in \mathbb{Z}_{\ge 1}\}$.

Let us label a plateau $\mathcal{P}$ with a subset $\mathcal{P} \subseteq \mathcal{I}$ (denoted by the same symbol) such that

$$y_i^{(a)}(\zeta) = \begin{cases} y_{iL}^{(a)} & (a, i) \in \bar{\mathcal{P}}, \\ y_{iR}^{(a)} & (a, i) \in \mathcal{P}. \end{cases} \qquad (5.6)$$

Then the transition rule (5.3)–(5.4) from $\mathcal{P}$ to $\mathcal{P}'$ across the edge $\zeta = \zeta^*$ is stated as

$$\mathcal{P} \sqcup \mathcal{J} = \mathcal{P}'. \qquad (5.7)$$

The boundary condition means that $\mathcal{P} \to \emptyset$ as $\zeta \to -\infty$ and $\mathcal{P} \to \mathcal{I}$ as $\zeta \to \infty$. Thus a plateau configuration can be viewed as specifying a monotonic growth pattern of the label $\mathcal{P}$. One starts from $\mathcal{P} = \emptyset$ in the far left and proceeds to the right inductively to eventually get $\mathcal{P} = \mathcal{I}$ in the far right. In the inductive step, we assume that a plateau $\mathcal{P}$ and the associated data $\{v_i^{(a)}(r, l)_{\mathcal{P}} \mid (a, i) \in \mathcal{I}\}$ are at hand. Then its right edge $\zeta^*$ as in (5.2) and the set $\mathcal{J}$ which specifies the next plateau $\mathcal{P}'$ by (5.7) are determined as

$$\zeta^* = \min_{(a,i) \in \bar{\mathcal{P}}} \{v_i^{(a)}(r, l)_{\mathcal{P}}\}, \qquad (5.8)$$

$$\mathcal{J} = \{(a, i) \in \bar{\mathcal{P}} \mid v_i^{(a)}(r, l)_{\mathcal{P}} = \zeta^*\}. \qquad (5.9)$$

<footnote n="13">¹³For any subset $\mathcal{K} \subseteq \mathcal{I}$, we let $\bar{\mathcal{K}}$ denote the complement $\mathcal{I} \setminus \mathcal{K}$.</footnote>
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_i^{(a)} = y_{iL}^{(a)}$ or $y_{iR}^{(a)}$ on a plateau can be interpreted as whether the solitons $S_i^{(a)}$ there originate in the domain $x < 0$ or $x > 0$. This is indicated in (3.20) which equates $y_i^{(a)}$ to the density $\rho_i^{(a)}$ of $S_i^{(a)}$ up to the factor $\sigma_i^{(a)}$ reflecting the interaction effect. Thus we may interpret physically a plateau label $\mathcal{P}$ as the soliton content meaning that it is the region of $\zeta$ where

solitons $\{S_i^{(a)} | (a, i) \in \mathcal{P}\}$ originating in the right domain have not gone away yet, (5.10)

solitons $\{S_i^{(a)} | (a, i) \in \bar{\mathcal{P}}\}$ originating in the left domain have already reached. (5.11)

In what follows we will signify the solitons in (5.10), (5.11) as $S_{iR}^{(a)}$, $S_{iL}^{(a)}$ and refer to $\mathcal{P}$ also as soliton content. If the initial domain is empty<fnref n="14" /> 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}, \dots, z_{nL}$) and ($z_{1R}, \dots, z_{nR}$) satisfying (3.14), the following properties are valid:

1. Actually realized soliton contents $\mathcal{P} \subseteq \mathcal{I}$ are limited to the following form for some $\mathbf{d}$:

$$ \mathcal{P}(\mathbf{d}) = \{(a, i) | 1 \le i \le d_a, a \in [1, n-1]\}, \quad \mathbf{d} = (d_1, \dots, d_{n-1}) \in (\mathbb{Z}_{\ge 0})^{n-1}. \tag{5.12} $$

2. The effective speeds are nonnegative for any time evolution and on any plateau, i.e.,

$$ v_i^{(a)}(r, l)_{\mathcal{P}(\mathbf{d})} \ge 0 \quad (\forall (a, i), (r, l) \in \mathcal{I}). \tag{5.13} $$

Nonnegativity of the velocities is not obvious since the solitons $S_i^{(a)}$ with $a \ne r$ have the vanishing bare speed under $T_l^{(r)}$ (see Sec. 2.5 (I)), and they may receive a “recoil effect” from the march of color $r$ solitons. For instance there is a minus sign in (4.21) although these velocities are indeed positive for $0 < q < 1$. Due to (5.13) all the solitons are right movers, therefore the domain $\zeta < 0$ is not influenced by the domain $\zeta > 0$. This feature allows us to replace one of the boundary condition $y_i^{(a)}(-\infty) = y_{iL}^{(a)}$ by $y_i^{(a)}(-0) = y_{iL}^{(a)}$. Consequently, the plateaux can be numbered as 0, 1, 2, ... from the left to the right, where the leftmost 0 covers the entire $\zeta < 0$ region at least.

Thanks to (5.12) a plateau can be labeled by a vector $\mathbf{d} \in (\mathbb{Z}_{\ge 0})^{n-1}$, which corresponds to the soliton content

$$ \begin{gathered} S_{1R}^{(1)}, \dots, S_{d_1 R}^{(1)}, S_{d_1+1 L}^{(1)}, S_{d_1+2 L}^{(1)}, \dots, \\ S_{1R}^{(2)}, \dots, S_{d_2 R}^{(2)}, S_{d_2+1 L}^{(2)}, S_{d_2+2 L}^{(2)}, \dots, \\ \dots \\ S_{1R}^{(n-1)}, \dots, S_{d_{n-1} R}^{(n-1)}, S_{d_{n-1}+1 L}^{(n-1)}, S_{d_{n-1}+2 L}^{(n-1)}, \dots \end{gathered} \tag{5.14} $$

For each color $a$, it is like a Dirac sea of level $i = d_a$ without a hole. Obviously $\mathcal{P}(\mathbf{d}) \subseteq \mathcal{P}(\mathbf{d}')$ holds if and only if $\mathbf{d} \le \mathbf{d}'$, i.e. $\mathbf{d}' - \mathbf{d} \in (\mathbb{Z}_{\ge 0})^{n-1}$. Thus one can describe a plateau configuration just by an increasing sequence as $(0, \dots, 0) = \mathbf{d}_0 \le \mathbf{d}_1 \le \mathbf{d}_2 \le \dots \in (\mathbb{Z}_{\ge 0})^{n-1}$ which should tend to $(\infty, \dots, \infty)$.

<footnote n="14">¹⁴Only vac (2.21) without a soliton, which corresponds to the dilute limit as mentioned in Remark 3.2 and Remark 4.4.</footnote>
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_i^{(a)} = y_{iL}^{(a)}$ or $y_{iR}^{(a)}$ on a plateau can be interpreted as whether the solitons $S_i^{(a)}$ there originate in the domain $x < 0$ or $x > 0$. This is indicated in (3.20) which equates $y_i^{(a)}$ to the density $\rho_i^{(a)}$ of $S_i^{(a)}$ up to the factor $\sigma_i^{(a)}$ reflecting the interaction effect. Thus we may interpret physically a plateau label $\mathcal{P}$ as the soliton content meaning that it is the region of $\zeta$ where

solitons $\{S_i^{(a)} | (a, i) \in \mathcal{P}\}$ originating in the right domain have not gone away yet, (5.10)

solitons $\{S_i^{(a)} | (a, i) \in \bar{\mathcal{P}}\}$ originating in the left domain have already reached. (5.11)

In what follows we will signify the solitons in (5.10), (5.11) as $S_{iR}^{(a)}$, $S_{iL}^{(a)}$ and refer to $\mathcal{P}$ also as soliton content. If the initial domain is empty<fnref n="14" /> 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}, \dots, z_{nL}$) and ($z_{1R}, \dots, z_{nR}$) satisfying (3.14), the following properties are valid:

1. Actually realized soliton contents $\mathcal{P} \subseteq \mathcal{I}$ are limited to the following form for some $\mathbf{d}$:

$$ \mathcal{P}(\mathbf{d}) = \{(a, i) | 1 \le i \le d_a, a \in [1, n-1]\}, \quad \mathbf{d} = (d_1, \dots, d_{n-1}) \in (\mathbb{Z}_{\ge 0})^{n-1}. \tag{5.12} $$

2. The effective speeds are nonnegative for any time evolution and on any plateau, i.e.,

$$ v_i^{(a)}(r, l)_{\mathcal{P}(\mathbf{d})} \ge 0 \quad (\forall (a, i), (r, l) \in \mathcal{I}). \tag{5.13} $$

Nonnegativity of the velocities is not obvious since the solitons $S_i^{(a)}$ with $a \ne r$ have the vanishing bare speed under $T_l^{(r)}$ (see Sec. 2.5 (I)), and they may receive a “recoil effect” from the march of color $r$ solitons. For instance there is a minus sign in (4.21) although these velocities are indeed positive for $0 < q < 1$. Due to (5.13) all the solitons are right movers, therefore the domain $\zeta < 0$ is not influenced by the domain $\zeta > 0$. This feature allows us to replace one of the boundary condition $y_i^{(a)}(-\infty) = y_{iL}^{(a)}$ by $y_i^{(a)}(-0) = y_{iL}^{(a)}$. Consequently, the plateaux can be numbered as 0, 1, 2, ... from the left to the right, where the leftmost 0 covers the entire $\zeta < 0$ region at least.

Thanks to (5.12) a plateau can be labeled by a vector $\mathbf{d} \in (\mathbb{Z}_{\ge 0})^{n-1}$, which corresponds to the soliton content

$$ \begin{gathered} S_{1R}^{(1)}, \dots, S_{d_1 R}^{(1)}, S_{d_1+1 L}^{(1)}, S_{d_1+2 L}^{(1)}, \dots, \\ S_{1R}^{(2)}, \dots, S_{d_2 R}^{(2)}, S_{d_2+1 L}^{(2)}, S_{d_2+2 L}^{(2)}, \dots, \\ \dots \\ S_{1R}^{(n-1)}, \dots, S_{d_{n-1} R}^{(n-1)}, S_{d_{n-1}+1 L}^{(n-1)}, S_{d_{n-1}+2 L}^{(n-1)}, \dots \end{gathered} \tag{5.14} $$

For each color $a$, it is like a Dirac sea of level $i = d_a$ without a hole. Obviously $\mathcal{P}(\mathbf{d}) \subseteq \mathcal{P}(\mathbf{d}')$ holds if and only if $\mathbf{d} \le \mathbf{d}'$, i.e. $\mathbf{d}' - \mathbf{d} \in (\mathbb{Z}_{\ge 0})^{n-1}$. Thus one can describe a plateau configuration just by an increasing sequence as $(0, \dots, 0) = \mathbf{d}_0 \le \mathbf{d}_1 \le \mathbf{d}_2 \le \dots \in (\mathbb{Z}_{\ge 0})^{n-1}$ which should tend to $(\infty, \dots, \infty)$.

<footnote n="14">¹⁴Only vac (2.21) without a soliton, which corresponds to the dilute limit as mentioned in Remark 3.2 and Remark 4.4.</footnote>
## 5.3 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. The system has open boundary conditions and is initially divided into three regions: the ‘left’ correspond to $x \in [-L/3, -1]$, the ‘right’ corresponds to $x \in [0, L/3 - 1]$, and the ‘buffer’ to $x \in [L/3, 2L/3 - 1]$.

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 been 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 \in [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_l^{(a)}$ 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_l^{(1)}T_l^{(2)}$. The application of such a transfer matrix on a given configuration is done using the combinatorial $R$ defined in App. B. At a few specific times the local densities $\rho_k^{(a)}(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 ($\rho_{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} \ge 300000$.

The Figs. 1, 2, 3 and 4 display the six soliton densities $\rho_{k=1,2,3}^{(a=1,2)}(x, t)$ at time $t = 600$ and $t = 1000$ as a function of $\zeta = 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. 10.

Fig. 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_{iR}^{(a)}$ therefore vanish. So, if a soliton type $(i, a)$ belongs to the set $\mathcal{P}$ of a given plateau, then this species is absent in this plateau, its density $\rho_i^{(a)}$ vanishes. In Fig. 1 the first step is located at $\zeta(1) \approx 0.23815$, and beyond this value of $\zeta$ the solitons $S_1^{(2)}$ are absent. At this transition we have $\mathcal{J} = \{(a = 2, i = 1)\}$. The second transition ($\zeta(2) \approx 0.30285$) is the place where the $S_1^{(1)}$ solitons disappear, hence $\mathcal{J} = \{(a = 1, i = 1)\}$ for this step. The third step is located at $\zeta(3) \approx 0.53924$. Beyond this value of $\zeta$ the $S_2^{(2)}$ soliton are absent and this step is characterized by $\mathcal{J} = \{(a = 2, i = 2)\}$. At the fourth step, $\zeta(4) \approx 0.80940$, the $S_3^{(2)}$ soliton
## 5.3 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. The system has open boundary conditions and is initially divided into three regions: the ‘left’ correspond to $x \in [-L/3, -1]$, the ‘right’ corresponds to $x \in [0, L/3 - 1]$, and the ‘buffer’ to $x \in [L/3, 2L/3 - 1]$.

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 been 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 \in [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_l^{(a)}$ 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_l^{(1)}T_l^{(2)}$. The application of such a transfer matrix on a given configuration is done using the combinatorial $R$ defined in App. B. At a few specific times the local densities $\rho_k^{(a)}(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 ($\rho_{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} \ge 300000$.

The Figs. 1, 2, 3 and 4 display the six soliton densities $\rho_{k=1,2,3}^{(a=1,2)}(x, t)$ at time $t = 600$ and $t = 1000$ as a function of $\zeta = 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. 10.

Fig. 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_{iR}^{(a)}$ therefore vanish. So, if a soliton type $(i, a)$ belongs to the set $\mathcal{P}$ of a given plateau, then this species is absent in this plateau, its density $\rho_i^{(a)}$ vanishes. In Fig. 1 the first step is located at $\zeta(1) \approx 0.23815$, and beyond this value of $\zeta$ the solitons $S_1^{(2)}$ are absent. At this transition we have $\mathcal{J} = \{(a = 2, i = 1)\}$. The second transition ($\zeta(2) \approx 0.30285$) is the place where the $S_1^{(1)}$ solitons disappear, hence $\mathcal{J} = \{(a = 1, i = 1)\}$ for this step. The third step is located at $\zeta(3) \approx 0.53924$. Beyond this value of $\zeta$ the $S_2^{(2)}$ soliton are absent and this step is characterized by $\mathcal{J} = \{(a = 2, i = 2)\}$. At the fourth step, $\zeta(4) \approx 0.80940$, the $S_3^{(2)}$ soliton
<figure><img src="image_10.png" /><figcaption>Figure 1: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$.</figcaption></figure>

density drops to zero and $\mathcal{J} = \{(a = 2, i = 3)\}$. Next we have $\zeta(5) \simeq 0.84452$, where the $S_2^{(1)}$ soliton disappear ($\mathcal{J} = \{(a = 1, i = 2)\}$). So far all this information can be read out from Fig. 1, by looking at the values of $\zeta$ where one soliton density drops to zero. Using the notations of (5.12) this translates into: $\mathbf{d}_0 = (0,0)$, $\mathbf{d}_1 = (0,1)$, $\mathbf{d}_2 = (1,1)$, $\mathbf{d}_3 = (1,2)$, $\mathbf{d}_4 = (1,3)$, $\mathbf{d}_5 = (2,3)$.

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 $\zeta$ (data not shown). Such an analysis shows that for this particular initial condition the soliton content follows a somewhat irregular pattern up to the 12th plateau: $\mathbf{d}_6 = (2,4)$, $\mathbf{d}_7 = (2,5)$, $\mathbf{d}_8 = (2,6)$, $\mathbf{d}_9 = (2,7)$, $\mathbf{d}_{10} = (2,8)$, $\mathbf{d}_{11} = (2,9)$, $\mathbf{d}_{12} = (3,9)$. And then the sequence becomes regular: $\mathbf{d}_{n \ge 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: $\mathbf{d}_{\infty+1} = (4, \infty)$, $\mathbf{d}_{\infty+2} = (5, \infty)$, ..., $\mathbf{d}_{\infty+n} = (n+3, \infty)$.

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
<figure><img src="image_7.png" /><figcaption>Figure 1: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$.</figcaption></figure>

density drops to zero and $\mathcal{J} = \{(a = 2, i = 3)\}$. Next we have $\zeta(5) \simeq 0.84452$, where the $S_2^{(1)}$ soliton disappears ($\mathcal{J} = \{(a = 1, i = 2)\}$). So far all this information can be read out from Fig. 1, by looking at the values of $\zeta$ where one soliton density drops to zero. Using the notations of (5.12) this translates into: $\mathbf{d}_0 = (0,0)$, $\mathbf{d}_1 = (0,1)$, $\mathbf{d}_2 = (1,1)$, $\mathbf{d}_3 = (1,2)$, $\mathbf{d}_4 = (1,3)$, $\mathbf{d}_5 = (2,3)$.

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 $\zeta$ (data not shown). Such an analysis shows that for this particular initial condition the soliton content follows a somewhat irregular pattern up to the 12th plateau: $\mathbf{d}_6 = (2,4)$, $\mathbf{d}_7 = (2,5)$, $\mathbf{d}_8 = (2,6)$, $\mathbf{d}_9 = (2,7)$, $\mathbf{d}_{10} = (2,8)$, $\mathbf{d}_{11} = (2,9)$, $\mathbf{d}_{12} = (3,9)$. And then the sequence becomes regular: $\mathbf{d}_{n \ge 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: $\mathbf{d}_{\infty+1} = (4, \infty)$, $\mathbf{d}_{\infty+2} = (5, \infty)$, ..., $\mathbf{d}_{\infty+n} = (n+3, \infty)$.

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
<figure><img src="image_11.png" /><figcaption>Figure 2: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$.</figcaption></figure>

$\zeta = -\infty$ 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: $\forall i \ v_i^{(2)}(r = 1, l)_{P_0} = 0$. The edge of $P_0$ is thus at $\zeta = \zeta(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 $\zeta = 0$, where the all $S_k^{(2)}$ 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].

Fig. 3 shows the soliton densities in a situation where both the left and the right are in a nontrivial i.i.d. state at $t = 0$ (the left and right fugacities for this example are given in the figure caption). Here again the comparison between the simulations and the densities obtained from GHD are in very good agreement.

In Fig. 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_{l=100}^{(1)}T_{l=100}^{(2)}$. In the GHD framework this means that the bare velocity $v_{i,bare}^{(a)}$ of the $S_i^{(a)}$ solitons is the sum of its value $\delta_{1,a}\min(i, l)$ (see (4.1)) for
<figure><img src="image_8.png" /><figcaption>Figure 2: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$.</figcaption></figure>

$\zeta = -\infty$ 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: $\forall i \ v_i^{(2)}(r = 1, l)_{P_0} = 0$. The edge of $P_0$ is thus at $\zeta = \zeta(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 $\zeta = 0$, where the all $S_k^{(2)}$ 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].

Fig. 3 shows the soliton densities in a situation where both the left and the right are in a nontrivial i.i.d. state at $t = 0$ (the left and right fugacities for this example are given in the figure caption). Here again the comparison between the simulations and the densities obtained from GHD are in very good agreement.

In Fig. 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_{l=100}^{(1)}T_{l=100}^{(2)}$. In the GHD framework this means that the bare velocity $v_{i,bare}^{(a)}$ of the $S_i^{(a)}$ solitons is the sum of its value $\delta_{1,a}\min(i, l)$ (see (4.1)) for
<figure><img src="image_12.png" /><figcaption>Figure 3: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 500$ and $t = 2000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$, and the following fugacities in the right: $z_{1,R} = 1$, $z_{2,R} = 0.9$ and $z_{3,R} = 0.1$. Parameters of the simulations: system size $L = 640000$ and number of random initial conditions $N_{\text{samples}} = 200000$.</figcaption></figure>

$T_l^{(r=1)}$ and its value $\delta_{2,a}\min(i, l)$ for $T_l^{(r=2)}$. We thus have $v_{i,bare}^{(a)} = \min(i, l)$ for $T_{l=100}^{(1)}T_{l=100}^{(2)}$. Contrary to the case where the dynamics is generated by $T_l^{(1)}$ or $T_l^{(2)}$, there is no infinite series of narrow plateaux which accumulate in the vicinity of some finite value of $\zeta$.

Fig. 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: $\varrho_1 = 2-A$, $\varrho_2 = 1+A-B$ and $\varrho_3 = B$ with $A = \sum_{i=1}^{\infty} i\rho_i^{(a=1)}$ and $B = \sum_{i=1}^{\infty} i\rho_i^{(a=2)}$. In Fig. 5 an infinite number of very narrow plateaux are observed to accumulate in the vicinity of $\zeta = \zeta_c \simeq 1.88708$. When approaching this limiting value $\zeta_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 $\zeta > \zeta_c$ take the form $\mathbf{d} = (k, \infty)$. A very similar infinite series of narrow plateaux are observed for the $T_{l=3}^{(2)}$ dynamic, as illustrated in Fig. 6. Fig. 7 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_{l=100}^{(2)}$ dynamics, as
<figure><img src="image_9.png" /><figcaption>Figure 3: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 500$ and $t = 2000$ and for the $T_{l=100}^{(1)}$ dynamics. The red line indicates the GHD prediction for $t = \infty$. 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$, and the following fugacities in the right: $z_{1,R} = 1$, $z_{2,R} = 0.9$ and $z_{3,R} = 0.1$. Parameters of the simulations: system size $L = 640000$ and number of random initial conditions $N_{\text{samples}} = 200000$.</figcaption></figure>

$T_l^{(r=1)}$ and its value $\delta_{2,a}\min(i, l)$ for $T_l^{(r=2)}$. We thus have $v_{i,bare}^{(a)} = \min(i, l)$ for $T_{l=100}^{(1)}T_{l=100}^{(2)}$. Contrary to the case where the dynamics is generated by $T_l^{(1)}$ or $T_l^{(2)}$, there is no infinite series of narrow plateaux which accumulate in the vicinity of some finite value of $\zeta$.

Fig. 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: $\varrho_1 = 2-A$, $\varrho_2 = 1+A-B$ and $\varrho_3 = B$ with $A = \sum_{i=1}^{\infty} i\rho_i^{(a=1)}$ and $B = \sum_{i=1}^{\infty} i\rho_i^{(a=2)}$. In Fig. 5 an infinite number of very narrow plateaux are observed to accumulate in the vicinity of $\zeta = \zeta_c \simeq 1.88708$. When approaching this limiting value $\zeta_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 $\zeta > \zeta_c$ take the form $\mathbf{d} = (k, \infty)$. A very similar infinite series of narrow plateaux are observed for the $T_{l=3}^{(2)}$ dynamic, as illustrated in Fig. 6. Fig. 7 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_{l=100}^{(2)}$ dynamics, as
<figure><img src="image_13.png" /><figcaption>Figure 4: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and ‘color’ $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}T_{l=100}^{(2)}$ dynamics. The dotted lines indicate the GHD prediction for the step positions. 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_{\text{samples}} = 300000$.</figcaption></figure>

illustrated in Fig. 8. In that case the color-1 solitons are absent for $\zeta$ beyond the accumulation point (found at $\zeta_c \simeq 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_{l=100}^{(1)}T_{l=100}^{(2)}$ dynamics (Fig. 4).

Fig. 9 shows the velocities $v_{i=1,...,10}^{(a=1,2)}$ and soliton densities $\rho_{i=1,...,10}^{(a=1,2)}$ in the four first plateaux of the domain wall problem where the dynamics is generated by $T_{l=3}^{(r=1)}$, 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 monotonic 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 \to \infty$. As for the color-2 solitons, their velocities turn out to be increasing functions of the size. In the lower panels (densities) one can check that i) the density of $S_1^{(2)}$ solitons vanishes in the plateau $k \ge 1$ (since they are the slowest species in the plateau 0), ii) the density of $S_2^{(2)}$ solitons vanishes in the
<figure><img src="image_10.png" /><figcaption>Figure 4: Numerical results for several $S_k^{(a)}$-soliton densities (with sizes $k = 1, 2, 3$ and 'color' $a = 1, 2$) as a function of $\zeta = x/t$ for $t = 600$ and $t = 1000$ and for the $T_{l=100}^{(1)}T_{l=100}^{(2)}$ dynamics. The dotted lines indicate the GHD prediction for the step positions. 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_{\text{samples}} = 300000$.</figcaption></figure>

illustrated in Fig. 8. In that case the color-1 solitons are absent for $\zeta$ beyond the accumulation point (found at $\zeta_c \simeq 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_{l=100}^{(1)}T_{l=100}^{(2)}$ dynamics (Fig. 4).

Fig. 9 shows the velocities $v_{i=1...10}^{(a=1,2)}$ and soliton densities $\rho_{i=1...10}^{(a=1,2)}$ in the four first plateaux of the domain wall problem where the dynamics is generated by $T_{l=3}^{(r=1)}$, 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 monotonic 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 \to \infty$. As for the color-2 solitons, their velocities turn out to be increasing functions of the size. In the lower panels (densities) one can check that i) the density of $S_1^{(2)}$ solitons vanishes in the plateau $k \ge 1$ (since they are the slowest species in the plateau 0), ii) the density of $S_2^{(2)}$ solitons vanishes in the