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

Figure 5: Densities $\rho_{i=1,2,3}$ of boxes with the entry $i$, as a function of $\zeta = x/t$ for $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. Note that by construction $\rho_1 + \rho_2 + \rho_3 = 3$, since each site contains 3 boxes (one in the tableau $B^{1,1}$ and two in the tableau $B^{2,1}$). The vertical red dotted line indicate the locations $\zeta(k)$ of the steps between the different plateaux, as predicted by GHD. The orange, yellow and blue lines are $\rho_1$, $\rho_2$ and $\rho_3$ from GHD (at $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. A large number of very narrow plateaux are observed to accumulate in the vicinity of $\zeta_c \simeq 1.88708$. Note that additional plateaux are present for $\zeta > 4$ (not shown) and the system is still not in the vacuum state at the right of the figure (contrary to Fig. 6). Parameters of the simulations: system size $L = 300000$ and number of random initial conditions $N_{\text{samples}} = 300000$.

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

Figure 6: Same as Fig. 5 but for the $T_{l=3}^{(1)}$ dynamics. We again have an infinite series of narrow plateaux which accumulate in the vicinity of $\zeta \simeq 0.81586$. The last step is located at $\zeta = 3$, beyond this point the system is in the vacuum state (i.e $\rho_1 = 2$, $\rho_2 = 1$ and $\rho_3 = 0$). Notice the density spike just before $\zeta = 3$. The soliton speeds and densities in the three four plateaux are displayed in Fig. 9.

plateau $k \ge 2$ (since they are the slowest species with non-zero density in the plateau 1), and
iii) the density of $S_1^{(1)}$ solitons vanishes in the plateau $k \ge 3$ (since they are the slowest species
<figure><img src="image_11.png" /></figure>

Figure 5: Densities $\rho_{i=1,2,3}$ of boxes with the entry $i$, as a function of $\zeta = x/t$ for $t = 1000$ and for the $T_{l=100}^{(1)}$ dynamics. Note that by construction $\rho_1 + \rho_2 + \rho_3 = 3$, since each site contains 3 boxes (one in the tableau $B^{1,1}$ and two in the tableau $B^{2,1}$). The vertical red dotted line indicate the locations $\zeta(k)$ of the steps between the different plateaux, as predicted by GHD. The orange, yellow and blue lines are $\rho_1$, $\rho_2$ and $\rho_3$ from GHD (at $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. A large number of very narrow plateaux are observed to accumulate in the vicinity of $\zeta_c \simeq 1.88708$. Note that additional plateaux are present for $\zeta > 4$ (not shown) and the system is still not in the vacuum state at the right of the figure (contrary to Fig. 6). Parameters of the simulations: system size $L = 300000$ and number of random initial conditions $N_{\text{samples}} = 300000$.

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

Figure 6: Same as Fig. 5 but for the $T_{l=3}^{(1)}$ dynamics. We again have an infinite series of narrow plateaux which accumulate in the vicinity of $\zeta \simeq 0.81586$. The last step is located at $\zeta = 3$, beyond this point the system is in the vacuum state (i.e $\rho_1 = 2$, $\rho_2 = 1$ and $\rho_3 = 0$). Notice the density spike just before $\zeta = 3$. The soliton speeds and densities in the three four plateaux are displayed in Fig. 9.

plateau $k \ge 2$ (since they are the slowest species with non-zero density in the plateau 1), and
iii) the density of $S_1^{(1)}$ solitons vanishes in the plateau $k \ge 3$ (since they are the slowest species
<figure><img src="image_16.png" /><figcaption>Figure 7: Ball densities $\rho_{i=1,2,3}$ obtained from GHD calculations and plotted as a function of $\zeta$. The different colors correspond to different values of the parameter $l$ in $T_l^{(1)}$. The initial state is defined by the fugacities $z_{1,L} = 1$, $z_{2,L} = 0.7$, $z_{3,L} = 0.3$ and with the vacuum state in the right part. For $l=8$ the results of simulations (at $t=1000$ and $N_{samples} = 3.10^5$) are also shown. The sum of the three densities is equal to 3 by construction (notice the different vertical scales in the different panels). Three regions can be distinguished: i) In the region $\zeta < \zeta_c(l)$ some soliton of color 1 and color 2 are present. An infinity of narrow plateaux develop when $\zeta \to \zeta_c(l)$, where the color-2 solitons of increasing size gradually disappear. ii) For $\zeta_c(l) < \zeta < l$ only color-1 solitons are present. There is a finite number (of order $l$) of plateaux in this region. iii) For $l < \zeta$ the system is in the vacuum state. The critical values of $\zeta$ are: $\zeta_c(2) \approx 0.57484$, $\zeta_c(3) \approx 0.81586$, $\zeta_c(4) \approx 1.02025$, $\zeta_c(8) \approx 1.55262$, $\zeta_c(100) \approx 1.88708$, and they can be identified in the right panel as the points where $\rho_3$ vanishes.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 8: Same as Fig. 5 but for the $T_{l=100}^{(2)}$ dynamics. We again have an infinite series of narrow plateaux, they accumulate here in the vicinity of $\zeta_c \approx 0.80019$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 7: Ball densities $\rho_{i=1,2,3}$ obtained from GHD calculations and plotted as a function of $\zeta$. The different colors correspond to different values of the parameter $l$ in $T_l^{(1)}$. The initial state is defined by the fugacities $z_{1,L} = 1$, $z_{2,L} = 0.7$, $z_{3,L} = 0.3$ and with the vacuum state in the right part. For $l=8$ the results of simulations (at $t=1000$ and $N_{samples} = 3.10^5$) are also shown. The sum of the three densities is equal to 3 by construction (notice the different vertical scales in the different panels). Three regions can be distinguished: i) In the region $\zeta < \zeta_c(l)$ some soliton of color 1 and color 2 are present. An infinity of narrow plateaux develop when $\zeta \to \zeta_c(l)$, where the color-2 solitons of increasing size gradually disappear. ii) For $\zeta_c(l) < \zeta < l$ only color-1 solitons are present. There is a finite number (of order $l$) of plateaux in this region. iii) For $l < \zeta$ the system is in the vacuum state. The critical values of $\zeta$ are: $\zeta_c(2) \approx 0.57484$, $\zeta_c(3) \approx 0.81586$, $\zeta_c(4) \approx 1.02025$, $\zeta_c(8) \approx 1.55262$, $\zeta_c(100) \approx 1.88708$, and they can be identified in the right panel as the points where $\rho_3$ vanishes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 8: Same as Fig. 5 but for the $T_{l=100}^{(2)}$ dynamics. We again have an infinite series of narrow plateaux, they accumulate here in the vicinity of $\zeta_c \approx 0.80019$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 9: Soliton velocities $v_{i=1,...,10}^{(a=1,2)}$ (top panels) and soliton densities $\rho_{i=1,...,10}^{(a=1,2)}$ (bottom panels) plotted as a function of soliton size $i$ in the four first plateaux. These curves are obtained from GHD calculations and for the same parameters as in Fig. 6, that is with the $T_{l=3}^{(r=1)}$ dynamics, fugacities $z_{1,L} = z_{2,L} = 0.7$, $z_{3,L} = 0.3$ and with the vacuum state in the right part at $t = 0$.</figcaption></figure>

with non-zero density in the plateau 2). For this particular initial state we have checked for the $T_l^{(r=1)}$ dynamics with $l = 2, 3, ..., 6$ that the fastest color-1 soliton has size $l$.

Finally Fig. 10 shows that the widths of the steps between consecutive plateaux scale like $t^{1/2}$, in agreement with the diffusive corrections calculated in Sec. 5.2. This figure represents some soliton densities in vicinity of the step $k = 1, 2, 3$, plotted as a function of $(x/t - \zeta(k))t^{1/2}$. The fact that the data measured at different times almost fall onto the same error function curve demonstrates the diffusive nature of the step broadening. Moreover, the data are in good agreement with the GHD prediction (5.19) calculated in Sec. 5.2, including the quantitative value of the parameter $\Sigma_{P,P'}$. It should in particular be noted that the red curves in Fig. 10 are the GHD predictions and contain no adjustable parameter. In the bottom right panel of Fig. 10 the numerical data appear to be slightly below the GHD curve. The reason for this mismatch is not clear to us but it should however be noticed that it tends to *decrease* when the time increases (compare $t = 1000, 2000$ and $3000$). It is therefore likely due to some finite-time effect.
<figure><img src="image_15.png" /><figcaption>Figure 9: Soliton velocities $v_{i=1,...,10}^{(a=1,2)}$ (top panels) and soliton densities $\rho_{i=1,...,10}^{(a=1,2)}$ (bottom panels) plotted as a function of soliton size $i$ in the four first plateaux. These curves are obtained from GHD calculations and for the same parameters as in Fig. 6, that is with the $T_{l=3}^{(r=1)}$ dynamics, fugacities $z_{1,L} = z_{2,L} = 0.7$, $z_{3,L} = 0.3$ and with the vacuum state in the right part at $t = 0$.</figcaption></figure>

with non-zero density in the plateau 2). For this particular initial state we have checked for the $T_l^{(r=1)}$ dynamics with $l = 2, 3, ..., 6$ that the fastest color-1 soliton has size $l$.

Finally Fig. 10 shows that the widths of the steps between consecutive plateaux scale like $t^{1/2}$, in agreement with the diffusive corrections calculated in Sec. 5.2. This figure represents some soliton densities in vicinity of the step $k = 1, 2, 3$, plotted as a function of $(x/t - \zeta(k))t^{1/2}$. The fact that the data measured at different times almost fall onto the same error function curve demonstrates the diffusive nature of the step broadening. Moreover, the data are in good agreement with the GHD prediction (5.19) calculated in Sec. 5.2, including the quantitative value of the parameter $\Sigma_{P,P'}$. It should in particular be noted that the red curves in Fig. 10 are the GHD predictions and contain no adjustable parameter. In the bottom right panel of Fig. 10 the numerical data appear to be slightly below the GHD curve. The reason for this mismatch is not clear to us but it should however be noticed that it tends to *decrease* when the time increases (compare $t = 1000, 2000$ and $3000$). It is therefore likely due to some finite-time effect.
<figure><img src="image_19.png" /><figcaption>Figure 10: Densities $\rho_{i=1}^{(a)}$ of $S_1^{(a=1,2)}$ solitons in the vicinities of the three first steps (located at $x/t = \zeta(1)$, $x/t = \zeta(2)$ and $x/t = \zeta(3)$) and at three different times. With the choice $(x/t - \zeta(k))t^{1/2}$ for the horizontal scale the curves almost collapse onto each other. This illustrates the diffusive nature of the step broadening. The red curves have been obtained without any free parameter (no fit) and represent the GHD prediction. They take the form of a complementary error function (see Sec. 5.2 and (5.19)). Initial state fugacities: $z_{1,L} = 1$, $z_{2,L} = 0.7$ and $z_{3,L} = 0.3$, $z_{1,R} = 1$, $z_{2,R} = 0.9$ and $z_{3,R} = 0.1$. System size $L = 1.210^6$. Number of random initial conditions $N_{\text{samples}} = 400000$. This simulation represents $(2 \cdot L) \cdot N_{\text{samples}} \cdot t_{\text{max}} \sim 0.3 \cdot 10^{16}$ applications of the combinatorial R.</figcaption></figure>

# 6 Summary and conclusions

In this paper we have introduced a new family of box ball systems, dubbed *complete BBS*. The models are indexed by a parameter $n \in \mathbb{Z}_{>1}$ and are associated to the quantum group $U_q(\sl{sl}_n)$. The cBBS family generalizes the original single-color model [1], which is recovered for $n=2$. Each ‘site’ of the model is made of the $n-1$ column shape semistandard tableaux with length 1, 2, ..., $n-1$. The ‘balls’ can have $n$ different ‘colors’ and occupy the boxes of the above tableaux, respecting the semistandard rules. These states are then equipped with some integrable dynamics, constructed from the so-called combinatorial $R$ (Sec. 2.2). As for the simplest BBS, the time evolutions can be viewed as generated by a ‘carrier’ which propagates through the system at each time step and which takes, moves and deposits ‘balls’ from one place to another. In the case of the time evolution $T_l^{(k)}$ (2.26), the load of the carrier is itself a rectangular tableau of shape $k \times l$ with $k \in \{1, \dots, n-1\}$ and $l \in \mathbb{Z}_{>0}$ (Sec. 2.4). The integrability manifests itself from the fact that all these time evolutions commute (whatever $k$
<figure><img src="image_16.png" /><figcaption>Figure 10: Densities $\rho_{i=1}^{(a)}$ of $S_1^{(a=1,2)}$ solitons in the vicinities of the three first steps (located at $x/t = \zeta(1)$, $x/t = \zeta(2)$ and $x/t = \zeta(3)$) and at three different times. With the choice $(x/t - \zeta(k))t^{1/2}$ for the horizontal scale the curves almost collapse onto each other. This illustrates the diffusive nature of the step broadening. The red curves have been obtained without any free parameter (no fit) and represent the GHD prediction. They take the form of a complementary error function (see Sec. 5.2 and (5.19)). Initial state fugacities: $z_{1,L} = 1$, $z_{2,L} = 0.7$ and $z_{3,L} = 0.3$, $z_{1,R} = 1$, $z_{2,R} = 0.9$ and $z_{3,R} = 0.1$. System size $L = 1.210^6$. Number of random initial conditions $N_{\text{samples}} = 400000$. This simulation represents $(2 \cdot L) \cdot N_{\text{samples}} \cdot t_{\text{max}} \sim 0.3 \cdot 10^{16}$ applications of the combinatorial R.</figcaption></figure>

# 6 Summary and conclusions

In this paper we have introduced a new family of box ball systems, dubbed *complete BBS*. The models are indexed by a parameter $n \in \mathbb{Z}_{>1}$ and are associated to the quantum group $U_q(\sl{n})$. The cBBS family generalizes the original single-color model [1], which is recovered for $n=2$. Each ‘site’ of the model is made of the $n-1$ column shape semistandard tableaux with length 1, 2, ..., $n-1$. The ‘balls’ can have $n$ different ‘colors’ and occupy the boxes of the above tableaux, respecting the semistandard rules. These states are then equipped with some integrable dynamics, constructed from the so-called combinatorial $R$ (Sec. 2.2). As for the simplest BBS, the time evolutions can be viewed as generated by a ‘carrier’ which propagates through the system at each time step and which takes, moves and deposits ‘balls’ from one place to another. In the case of the time evolution $T_l^{(k)}$ (2.26), the load of the carrier is itself a rectangular tableau of shape $k \times l$ with $k \in \{1, \dots, n-1\}$ and $l \in \mathbb{Z}_{>0}$ (Sec. 2.4). The integrability manifests itself from the fact that all these time evolutions commute (whatever $k$
and $l$) and from the existence of a family of conserved energies (2.30).<fnref n="16" /> Another consequence of the integrability is the existence of solitons (Sec. 2.5). Each soliton is characterized by a 'color' $a$ and amplitude (or size) $i$. These solitons interact with each other in some nontrivial way during collisions. The local structure of each site of the cBBS may look complicated compared to the other versions of $n$-color BBS, but, as explained in Sec. 2, it offers a remarkable simplification: the soliton scattering is completely *diagonal*, namely each soliton regains its color and amplitude after collisions. The effect of the interactions is encoded in a simple phase shift. The end of Sec. 2 presents the inverse scattering scheme for this model, which is based on the KSS bijection (details in Appendix C), and which allows one to construct the action-angle variables.

The Sec. 3 discusses the randomized cBBS, where we have considered generalized Gibbs ensembles of cBBS configurations. We focused on a special family of GGE, called i.i.d., which are characterized by $n$ fugacities (one for each color) and where all the tableaux are independently distributed (Sec. 3.2). The properties of such i.i.d. GGE can be studied using TBA (Sec. 3.3) and we obtained the free energy, the mean value of the conserved energies and the associated soliton densities. In Sec. 4 we have presented the generalized hydrodynamics of the cBBS. The first ingredient is the effective speed of each soliton species in a homogenous state (see (4.1) and its solution in matrix form (4.31)). This could be obtained thanks to the TBA approach. Next, the currents carried by the soliton of each color are assumed to be given by the product of their effective speed and their density (4.33). For a system which is inhomogeneous but where spatial gradients are sufficiently small so that it can locally be approximated by a GGE, the continuity equations for the soliton currents give the GHD equations (4.39). These equations take a particularly simple form when re-written in terms of the $Y$-variables of the TBA (4.40), which shows that these variables are the normal modes of the GHD. As a concrete application of GHD we studied in Sec. 5 the evolution of the $n=3$ model starting from a domain wall initial state with two different i.i.d. states in the left and in the right halves. In good quantitative agreement with the GHD prediction, the high-precision simulations showed that at sufficiently long times the state of the system becomes piecewise constant in the variable $\zeta = x/t$ (5.2). Inside each plateau, the soliton densities and the ball densities measured in the simulations are in good agreement with the GHD expectations. Compared to the simpler $n=2$ model [11], the additional color degree of freedom leads to a complex evolution of the soliton content from one plateau to the next, and opens up the possibility to have an infinity of narrow plateaux which accumulate to some critical point $\zeta_c$ (Figs. 5, 6 and 8). Finally, we have computed in Sec. 5.2 the diffusion effect that is responsible for the broadening of the steps between consecutive plateaux. Here again we observed a good quantitative agreement between the simulations (Fig. 10) and the theoretical results. To our knowledge, this paper achieves the first systematic and successful application of GHD to an integrable system associated with a higher rank quantum algebra.

We wish to conclude this article by mentioning a few possible future directions. The equation (4.1) for the effective velocities is quite fundamental. It would be very interesting to have a proof of it, in the line of what was done in [11, Sec. 4.2] for $n=2$. See also [25] for $n=2$. It would also be appealing to invent and explore other protocols to set the system out of equilibrium. This might be done, for example, with a system that is open at its boundaries. One could also think of more general initial states, homogeneous or inhomogeneous, defined for instance by some nontrivial weights on the configurations, and look at the evolution toward equilibrium. This could be, for instance, a domain wall between two general (non-i.i.d.) GGE. In

<footnote n="16">¹⁶As a comparison, the classical hard rod model [19] (for which the GHD has also been worked out) has an infinite number of conserved quantities (the numbers of rods of each size) but it lacks the commuting dynamics. A similar remark also holds, for instance, for the classical cellular automata studied in Refs. [47,48], where domain-wall problems similar to those studied here have been solved exactly but where commuting transfer matrices are not known.</footnote>
and *l*) and from the existence of a family of conserved energies (2.30).<fnref n="16" /> Another consequence of the integrability is the existence of solitons (Sec. 2.5). Each soliton is characterized by a 'color' *a* and amplitude (or size) *i*. These solitons interact with each other in some nontrivial way during collisions. The local structure of each site of the cBBS may look complicated compared to the other versions of *n*-color BBS, but, as explained in Sec. 2, it offers a remarkable simplification: the soliton scattering is completely *diagonal*, namely each soliton regains its color and amplitude after collisions. The effect of the interactions is encoded in a simple phase shift. The end of Sec. 2 presents the inverse scattering scheme for this model, which is based on the KSS bijection (details in Appendix C), and which allows one to construct the action-angle variables.

The Sec. 3 discusses the randomized cBBS, where we have considered generalized Gibbs ensembles of cBBS configurations. We focused on a special family of GGE, called i.i.d., which are characterized by *n* fugacities (one for each color) and where all the tableaux are independently distributed (Sec. 3.2). The properties of such i.i.d. GGE can be studied using TBA (Sec. 3.3) and we obtained the free energy, the mean value of the conserved energies and the associated soliton densities. In Sec. 4 we have presented the generalized hydrodynamics of the cBBS. The first ingredient is the effective speed of each soliton species in a homogenous state (see (4.1) and its solution in matrix form (4.31)). This could be obtained thanks to the TBA approach. Next, the currents carried by the soliton of each color are assumed to be given by the product of their effective speed and their density (4.33). For a system which is inhomogeneous but where spatial gradients are sufficiently small so that it can locally be approximated by a GGE, the continuity equations for the soliton currents give the GHD equations (4.39). These equations take a particularly simple form when re-written in terms of the *Y*-variables of the TBA (4.40), which shows that these variables are the normal modes of the GHD. As a concrete application of GHD we studied in Sec. 5 the evolution of the *n* = 3 model starting from a domain wall initial state with two different i.i.d. states in the left and in the right halves. In good quantitative agreement with the GHD prediction, the high-precision simulations showed that at sufficiently long times the state of the system becomes piecewise constant in the variable $\zeta = x/t$ (5.2). Inside each plateau, the soliton densities and the ball densities measured in the simulations are in good agreement with the GHD expectations. Compared to the simpler *n* = 2 model [11], the additional color degree of freedom leads to a complex evolution of the soliton content from one plateau to the next, and opens up the possibility to have an infinity of narrow plateaux which accumulate to some critical point $\zeta_c$ (Figs. 5, 6 and 8). Finally, we have computed in Sec. 5.2 the diffusion effect that is responsible for the broadening of the steps between consecutive plateaux. Here again we observed a good quantitative agreement between the simulations (Fig. 10) and the theoretical results. To our knowledge, this paper achieves the first systematic and successful application of GHD to an integrable system associated with a higher rank quantum algebra.

We wish to conclude this article by mentioning a few possible future directions. The equation (4.1) for the effective velocities is quite fundamental. It would be very interesting to have a proof of it, in the line of what was done in [11, Sec. 4.2] for *n* = 2. See also [25] for *n* = 2. It would also be appealing to invent and explore other protocols to set the system out of equilibrium. This might be done, for example, with a system that is open at its boundaries. One could also think of more general initial states, homogeneous or inhomogeneous, defined for instance by some nontrivial weights on the configurations, and look at the evolution toward equilibrium. This could be, for instance, a domain wall between two general (non-i.i.d.) GGE. In

<footnote n="16">¹⁶As a comparison, the classical hard rod model [19] (for which the GHD has also been worked out) has an infinite number of conserved quantities (the numbers of rods of each size) but it lacks the commuting dynamics. A similar remark also holds, for instance, for the classical cellular automata studied in Refs. [47, 48], where domain-wall problems similar to those studied here have been solved exactly but where commuting transfer matrices are not known.</footnote>
the spirit of what have been done for some simple one-dimensional cellular automata [47,48] it would be interesting to compute some two-time correlation functions and to relate these correlations to transport coefficients. Another interesting question would be to explore the possibility of some anomalous transport in these models. Investigating numerically other values of $n$, possibly large ones, might also reveal some new phenomena. It also seems worth comparing the dynamics of the present cBBS with non-complete $n$-color BBS [14, 17, 20]. Finally, a possible direction would be to include some integrability-breaking perturbations in the model. We may expect the disappearance of the density plateaux, and it should be possible to observe numerically and to analyze some possible crossover from the integrable regime to a chaotic one where solitons are no longer stable but only short-lived. Such a setup could give rise to some prethermalization [49], possibly described by some GGE.

## A Algorithm for the combinatorial $R$

We use the notation in Sec. 2.2. Given $b \otimes c \in B^{k,l} \otimes B^{k',l'}$, we present an algorithm for finding the image $\tilde{c} \otimes \tilde{b} = R(b \otimes c) \in B^{k',l'} \otimes B^{k,l}$ of the combinatorial $R$ following [35, p55]. As noted after (2.7), the task is to construct the solution to the equation $\tilde{b} \cdot \tilde{c} = c \cdot b$ which is unique when the tableaux are rectangular.

A skew tableau $\theta$ is part of a semistandard Young tableau having a skew Young diagram shape. We call $\theta$ an $m$-vertical strip if it contains at most one box on each row and the total number of the boxes is $m$.

Let $P = c \cdot b$ be the product tableau [12]. Its shape Young diagram includes a $k \times l$ rectangle. Denote by $P'$ the NW part of $P$ corresponding to the $k \times l$ rectangle, and by $P/P'$ the skew tableau. The number of boxes of these tableaux (denoted by $|\cdot|$) are given by

$$|P| = kl + k'l', \quad |P'| = kl, \quad |P/P'| = k'l'. \qquad (A.1)$$

**Step 1.** Label each box of $P/P'$ with $\{1, 2, \dots, k'l'\}$.

Let $\theta_1$ be the rightmost vertical $k'$-strip in $P/P'$ as upper as possible. Remove $\theta_1$ from $P/P'$ and define the vertical $k'$-strip $\theta_2$ in a similar manner. This can be continued until $P/P'$ is decomposed into the disjoint union of the vertical $k'$-strips $\theta_1, \theta_2, \dots, \theta_{l'}$. Now the label is obtained by assigning the boxes of $\theta_i$ with $k'(i-1)+1, k'(i-1)+2, \dots, k'(i-1)+k'$ from the bottom to the top.

**Step 2.** Reverse row bumping from $P$ in the order of the label.

Find the semistandard Young tableaux $P_1, P_2, \dots, P_{k'l'}$ and $w_1, w_2, \dots, w_{k'l'} \in [1, n]$ such that

$$ \begin{array}{ll} P = (P_1 \leftarrow w_1); & \text{reverse bumping of the letter in the box 1,} \\ P_1 = (P_2 \leftarrow w_2); & \text{reverse bumping of the letter in the box 2,} \\ \dots & \\ P_{k'l' - 1} = (P_{k'l'} \leftarrow w_{k'l'}); & \text{reverse bumping of the letter in the box } k'l'. \end{array} \qquad (A.2) $$

Here $\leftarrow$ stands of the row insertion as in Sec. 2.2. The splitting of $P_{i-1}$ into $(P_i \leftarrow w_i)$ is done by the reverse row bumping. Starting from the box labeled $i$ in Step 1 and the letter in it, the bumping goes upwards row by row as follows. Given a letter $\alpha$ in a row, find the box in the adjacent upper row filled with the largest $\beta$ such that $\beta < \alpha$ and the rightmost one in case more than one $\beta$ is contained. Then let $\alpha$ occupy the box by bumping out $\beta$. Repeating this procedure one eventually bumps out a letter $w_i$ from the top row of $P_{i-1}$ thereby getting also $P_i$ as the remnant tableau. Note that $|P_{k'l'}| = kl$.
the spirit of what have been done for some simple one-dimensional cellular automata [47,48] it would be interesting to compute some two-time correlation functions and to relate these correlations to transport coefficients. Another interesting question would be to explore the possibility of some anomalous transport in these models. Investigating numerically other values of $n$, possibly large ones, might also reveal some new phenomena. It also seems worth comparing the dynamics of the present cBBS with non-complete $n$-color BBS [14, 17, 20]. Finally, a possible direction would be to include some integrability-breaking perturbations in the model. We may expect the disappearance of the density plateaux, and it should be possible to observe numerically and to analyze some possible crossover from the integrable regime to a chaotic one where solitons are no longer stable but only short-lived. Such a setup could give rise to some prethermalization [49], possibly described by some GGE.

## A Algorithm for the combinatorial $R$

We use the notation in Sec. 2.2. Given $b \otimes c \in B^{k,l} \otimes B^{k',l'}$, we present an algorithm for finding the image $\tilde{c} \otimes \tilde{b} = R(b \otimes c) \in B^{k',l'} \otimes B^{k,l}$ of the combinatorial $R$ following [35, p55]. As noted after (2.7), the task is to construct the solution to the equation $\tilde{b} \cdot \tilde{c} = c \cdot b$ which is unique when the tableaux are rectangular.

A skew tableau $\theta$ is part of a semistandard Young tableau having a skew Young diagram shape. We call $\theta$ an $m$-vertical strip if it contains at most one box on each row and the total number of the boxes is $m$.

Let $P = c \cdot b$ be the product tableau [12]. Its shape Young diagram includes a $k \times l$ rectangle. Denote by $P'$ the NW part of $P$ corresponding to the $k \times l$ rectangle, and by $P/P'$ the skew tableau. The number of boxes of these tableaux (denoted by $|\cdot|$) are given by

$$|P| = kl + k'l', \quad |P'| = kl, \quad |P/P'| = k'l'. \qquad (A.1)$$

**Step 1.** Label each box of $P/P'$ with $\{1, 2, \dots, k'l'\}$.

Let $\theta_1$ be the rightmost vertical $k'$-strip in $P/P'$ as upper as possible. Remove $\theta_1$ from $P/P'$ and define the vertical $k'$-strip $\theta_2$ in a similar manner. This can be continued until $P/P'$ is decomposed into the disjoint union of the vertical $k'$-strips $\theta_1, \theta_2, \dots, \theta_{l'}$. Now the label is obtained by assigning the boxes of $\theta_i$ with $k'(i-1)+1, k'(i-1)+2, \dots, k'(i-1)+k'$ from the bottom to the top.

**Step 2.** Reverse row bumping from $P$ in the order of the label.

Find the semistandard Young tableaux $P_1, P_2, \dots, P_{k'l'}$ and $w_1, w_2, \dots, w_{k'l'} \in [1, n]$ such that

$$ \begin{array}{ll} P = (P_1 \leftarrow w_1); & \text{reverse bumping of the letter in the box 1,} \\ P_1 = (P_2 \leftarrow w_2); & \text{reverse bumping of the letter in the box 2,} \\ \dots & \\ P_{k'l' - 1} = (P_{k'l'} \leftarrow w_{k'l'}); & \text{reverse bumping of the letter in the box } k'l'. \end{array} \qquad (A.2) $$

Here $\leftarrow$ stands of the row insertion as in Sec. 2.2. The splitting of $P_{i-1}$ into $(P_i \leftarrow w_i)$ is done by the reverse row bumping. Starting from the box labeled $i$ in Step 1 and the letter in it, the bumping goes upwards row by row as follows. Given a letter $\alpha$ in a row, find the box in the adjacent upper row filled with the largest $\beta$ such that $\beta < \alpha$ and the rightmost one in case more than one $\beta$ is contained. Then let $\alpha$ occupy the box by bumping out $\beta$. Repeating this procedure one eventually bumps out a letter $w_i$ from the top row of $P_{i-1}$ thereby getting also $P_i$ as the remnant tableau. Note that $|P_{k'l'}| = kl$.
Step 3. The image $\tilde{c} \otimes \tilde{b} \in B^{k,l} \otimes B^{k',l'}$ is constructed as

$$ \tilde{b} = P_{k'l'}, \qquad (A.3) $$

$$ \tilde{c} = ((\cdots (\emptyset \leftarrow w_{k'l'} \leftarrow \cdots) \leftarrow w_2) \leftarrow w_1). \qquad (A.4) $$

The semistandard tableau (A.4) is known as the P-symbol of the word $w_{k'l'} \dots w_2 w_1$.

**Example A.1.** Consider $b \in B^{2,3}$ and $c \in B^{3,2}$ in (2.3). Their product $c \cdot b$ is the right tableau in (2.6). The labeling of the boxes of $P/P'$ with $\{1, 2, \dots, 6\}$ according to Step 1 is shown by the indices in

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 & 2 \\ 3 & 5 & 5 & 1 & 4 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix}. \qquad (A.5) $$

The reverse bumping in Step 2 is done along the order of these indices as follows:

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 & 3 \\ 3 & 5 & 4 & 4 & 4 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 4 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 3 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 4 \end{array} \right), \qquad (A.6) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 4 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 4 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \end{array} \right), \qquad (A.7) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \end{array} \right), \qquad (A.8) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 \\ 2 & 2 & 5 & 5 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 5 & 4 \end{bmatrix} \leftarrow 1 \\ \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 5 & 4 \end{bmatrix} \leftarrow 1 \end{array} \right), \qquad (A.9) $$

$$ \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 5 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 5 & 4 \end{bmatrix} \leftarrow 2 \\ \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 5 & 4 \end{bmatrix} \leftarrow 2 \end{array} \right), \qquad (A.10) $$

$$ \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 5 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 5 & 4 \end{bmatrix} \leftarrow 3 \\ \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 5 & 4 \end{bmatrix} \leftarrow 3 \end{array} \right). \qquad (A.11) $$

Therefore from Step 3 we obtain

$$ \tilde{b} = \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix}, \qquad (A.12) $$

$$ \tilde{c} = (((((\emptyset \leftarrow 3) \leftarrow 2) \leftarrow 1) \leftarrow 5) \leftarrow 4) \leftarrow 3 = \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix}. \qquad (A.13) $$

This yields (2.9).

## B Combinatorial $R$ and energy $H$ for $n = 2, 3$

We consider $n = 3$ case first. The sets $B^{1,l}, B^{2,l}$ of semistandard tableaux are parameterized as

$$ B^{1,l} = \{x = (x_1, \dots, x_n) \in (\mathbb{Z}_{\ge 0})^n \mid x_1 + \dots + x_n = l\}, \quad x_i = \# \text{ of } i \text{ in tableau}, \qquad (B.1) $$

$$ B^{2,l} = \{x = (x_1, \dots, x_n) \in (\mathbb{Z}_{\ge 0})^n \mid x_1 + \dots + x_n = l\}, \quad x_i = \# \text{ of columns without } i \text{ in tableau}. \qquad (B.2) $$
Step 3. The image $\tilde{c} \otimes \tilde{b} \in B^{k,l} \otimes B^{k',l'}$ is constructed as

$$ \tilde{b} = P_{k'l'}, \qquad (A.3) $$

$$ \tilde{c} = ((\cdots (\emptyset \leftarrow w_{k'l'} \leftarrow \cdots) \leftarrow w_2) \leftarrow w_1). \qquad (A.4) $$

The semistandard tableau (A.4) is known as the P-symbol of the word $w_{k'l'} \dots w_2 w_1$.

**Example A.1.** Consider $b \in B^{2,3}$ and $c \in B^{3,2}$ in (2.3). Their product $c \cdot b$ is the right tableau in (2.6). The labeling of the boxes of $P/P'$ with $\{1, 2, \dots, 6\}$ according to Step 1 is shown by the indices in

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 & 2 \\ 3 & 5 & 5 & 1 & 4 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix}. \qquad (A.5) $$

The reverse bumping in Step 2 is done along the order of these indices as follows:

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 3 \\ 2 & 2 & 4 & 5 & 3 \\ 3 & 5 & 4 & 4 & 4 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 4 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 3 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 4 \end{array} \right), \qquad (A.6) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 4 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 4 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \end{array} \right), \qquad (A.7) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 3 & 4 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \\ \begin{bmatrix} 1 & 1 & 2 & 2 & 5 \\ 2 & 2 & 5 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 5 \end{array} \right), \qquad (A.8) $$

$$ \begin{bmatrix} 1 & 1 & 2 & 2 \\ 2 & 2 & 5 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 4 & 3 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 1 \\ \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 4 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 1 \end{array} \right), \qquad (A.9) $$

$$ \begin{bmatrix} 1 & 2 & 2 & 2 \\ 2 & 3 & 5 & 4 & 3 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 4 & 3 \\ 4 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 2 \\ \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 4 & 3 \\ 3 & 4 & 4 & 4 & 4 \end{bmatrix} \leftarrow 2 \end{array} \right), \qquad (A.10) $$

$$ \begin{bmatrix} 1 & 2 & 2 & 3 \\ 2 & 4 & 5 & 4 & 3 \end{bmatrix} = \left( \begin{array}{c} \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 & 3 \\ 4 & 4 & 4 & 4 & 3 \end{bmatrix} \leftarrow 3 \\ \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 & 3 \\ 3 & 4 & 4 & 4 & 3 \end{bmatrix} \leftarrow 3 \end{array} \right). \qquad (A.11) $$

Therefore from Step 3 we obtain

$$ \tilde{b} = \begin{bmatrix} 1 & 2 & 2 \\ 2 & 4 & 5 \end{bmatrix}, \qquad (A.12) $$

$$ \tilde{c} = (((((\emptyset \leftarrow 3) \leftarrow 2) \leftarrow 1) \leftarrow 5) \leftarrow 4) \leftarrow 3 = \begin{bmatrix} 1 & 3 \\ 2 & 4 \\ 3 & 5 \end{bmatrix}. \qquad (A.13) $$

This yields (2.9).

## B Combinatorial R and energy H for $n = 2, 3$

We consider $n = 3$ case first. The sets $B^{1,l}, B^{2,l}$ of semistandard tableaux are parameterized as

$$ B^{1,l} = \{x = (x_1, \dots, x_n) \in (\mathbb{Z}_{\ge 0})^n \mid x_1 + \dots + x_n = l\}, \quad x_i = \# \text{ of } i \text{ in tableau}, \qquad (B.1) $$

$$ B^{2,l} = \{x = (x_1, \dots, x_n) \in (\mathbb{Z}_{\ge 0})^n \mid x_1 + \dots + x_n = l\}, \quad x_i = \# \text{ of columns without } i \text{ in tableau}. \qquad (B.2) $$
We extend the indices to $\mathbb{Z}$ by $x_i = x_{i+n}$ and similarly also for $Q_i(x, y)$ and $P_i(x, y)$ introduced below. Then the combinatorial $R$ and the energy $H$ are given by the piecewise linear formulas as follows ($n = 3$) [21, eqs. (2.1)-(2.4)]:

$$R_{B^{1,l} \otimes B^{1,l'}} : \quad B^{1,l} \otimes B^{1,l'} \to B^{1,l'} \otimes B^{1,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.3)$$

$$\tilde{x}_i = x_i + Q_i(x, y) - Q_{i-1}(x, y), \quad \tilde{y}_i = y_i + Q_{i-1}(x, y) - Q_i(x, y),$$

$$Q_i(x, y) = \min\{\sum_{j=1}^{k-1} x_{i+j} + \sum_{j=k+1}^{n} y_{i+j} \mid 1 \le k \le n\},$$

$$H_{B^{1,l} \otimes B^{1,l'}}(x \otimes y) = Q_0(x, y). \qquad (B.4)$$

$$R_{B^{1,l} \otimes B^{2,l'}} : \quad B^{1,l} \otimes B^{2,l'} \to B^{2,l'} \otimes B^{1,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.5)$$

$$\tilde{x}_i = x_i + P_{i+1}(x, y) - P_i(x, y), \quad \tilde{y}_i = y_i + P_{i+1}(x, y) - P_i(x, y),$$

$$P_i(x, y) = \min(x_i, y_i),$$

$$H_{B^{1,l} \otimes B^{2,l'}}(x \otimes y) = P_1(x, y). \qquad (B.6)$$

$$R_{B^{2,l} \otimes B^{1,l'}} : \quad B^{2,l} \otimes B^{1,l'} \to B^{1,l'} \otimes B^{2,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.7)$$

$$\tilde{x}_i = x_i + P_{i-1}(x, y) - P_i(x, y), \quad \tilde{y}_i = y_i + P_{i-1}(x, y) - P_i(x, y),$$

$$H_{B^{2,l} \otimes B^{1,l'}}(x \otimes y) = P_0(x, y). \qquad (B.8)$$

$$R_{B^{2,l} \otimes B^{2,l'}} : \quad B^{2,l} \otimes B^{2,l'} \to B^{2,l'} \otimes B^{2,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.9)$$

$$\tilde{x}_i = x_i + Q_{i-1}(y, x) - Q_i(y, x), \quad \tilde{y}_i = y_i - Q_{i-1}(y, x) + Q_i(y, x),$$

$$H_{B^{2,l} \otimes B^{2,l'}}(x \otimes y) = Q_0(y, x). \qquad (B.10)$$

These formulas do not depend on $l, l'$ explicitly. When $n = 2$, the relevant objects are (B.1), (B.3) and (B.4) only, and one can just set $n = 2$ there.

For instance $Q_0(x, y) = \min\{y_2+y_3, x_1+y_3, x_1+x_2\}$ for $n = 3$ and $Q_0(x, y) = \min\{y_2, x_1\}$ for $n = 2$.

# C KSS bijection

## C.1 Paths

Consider the product set of the form $\mathcal{B} = B^{k_1,1} \otimes B^{k_2,1} \otimes \dots \otimes B^{k_N,1}$ with $k_i \in [1, n-1]$ containing column shape tableaux only. For the definition of $B^{r,s}$, see Sec. 2.2. States of our cBBS correspond to the choice $k_i \equiv i \mod n-1$ and $N = (n-1)L$. In this appendix, we do not assume the boundary condition of the BBS, and call the elements of $\mathcal{B}$ paths. Our aim is to explain the special case of the KSS bijection [50] corresponding to the above $\mathcal{B}$ along the convention adapted to this paper.

For a path $\mathfrak{p} = c_1 \otimes \dots \otimes c_N \in \mathcal{B}$ with $c_j \in B^{k_j,1}$, the $n$-array $\lambda = (\lambda_1, \dots, \lambda_n)$ defined by

$$\mathrm{wt}(\mathfrak{p}) = (\lambda_1, \dots, \lambda_n) \in (\mathbb{Z}_{\ge 0})^n, \quad \lambda_a = \sum_{j=1}^{N} (\text{number of letter } a \text{ contained in } c_j) \qquad (C.1)$$

is called the weight of $\mathfrak{p}$. Let $c_{jl}$ be the $l$ th entry of the tableau $c_j \in B^{k_j,1}$ from the top. Consider the word

$$w_1 w_2 \dots w_{N'} := c_{11} c_{12} \dots c_{1k_1} c_{21} c_{22} \dots c_{2k_2} \dots c_{N1} c_{N2} \dots c_{Nk_N}, \qquad (C.2)$$
We extend the indices to $\mathbb{Z}$ by $x_i = x_{i+n}$ and similarly also for $Q_i(x, y)$ and $P_i(x, y)$ introduced below. Then the combinatorial $R$ and the energy $H$ are given by the piecewise linear formulas as follows ($n = 3$) [21, eqs. (2.1)-(2.4)]:

$$R_{B^{1,l} \otimes B^{1,l'}} : \quad B^{1,l} \otimes B^{1,l'} \to B^{1,l'} \otimes B^{1,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.3)$$
$$\tilde{x}_i = x_i + Q_i(x, y) - Q_{i-1}(x, y), \quad \tilde{y}_i = y_i + Q_{i-1}(x, y) - Q_i(x, y),$$

$$Q_i(x, y) = \min\{\sum_{j=1}^{k-1} x_{i+j} + \sum_{j=k+1}^{n} y_{i+j} \mid 1 \le k \le n\},$$

$$H_{B^{1,l} \otimes B^{1,l'}}(x \otimes y) = Q_0(x, y). \qquad (B.4)$$

$$R_{B^{1,l} \otimes B^{2,l'}} : \quad B^{1,l} \otimes B^{2,l'} \to B^{2,l'} \otimes B^{1,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.5)$$
$$\tilde{x}_i = x_i + P_{i+1}(x, y) - P_i(x, y), \quad \tilde{y}_i = y_i + P_{i+1}(x, y) - P_i(x, y),$$

$$P_i(x, y) = \min(x_i, y_i), \qquad (B.6)$$

$$R_{B^{2,l} \otimes B^{1,l'}} : \quad B^{2,l} \otimes B^{1,l'} \to B^{1,l'} \otimes B^{2,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.7)$$
$$\tilde{x}_i = x_i + P_{i-1}(x, y) - P_i(x, y), \quad \tilde{y}_i = y_i + P_{i-1}(x, y) - P_i(x, y),$$

$$H_{B^{2,l} \otimes B^{1,l'}}(x \otimes y) = P_0(x, y). \qquad (B.8)$$

$$R_{B^{2,l} \otimes B^{2,l'}} : \quad B^{2,l} \otimes B^{2,l'} \to B^{2,l'} \otimes B^{2,l}; \quad x \otimes y \mapsto \tilde{y} \otimes \tilde{x} \qquad (B.9)$$
$$\tilde{x}_i = x_i + Q_{i-1}(y, x) - Q_i(y, x), \quad \tilde{y}_i = y_i - Q_{i-1}(y, x) + Q_i(y, x),$$

$$H_{B^{2,l} \otimes B^{2,l'}}(x \otimes y) = Q_0(y, x). \qquad (B.10)$$

These formulas do not depend on $l, l'$ explicitly. When $n = 2$, the relevant objects are (B.1), (B.3) and (B.4) only, and one can just set $n = 2$ there.

For instance $Q_0(x, y) = \min\{y_2+y_3, x_1+y_3, x_1+x_2\}$ for $n = 3$ and $Q_0(x, y) = \min\{y_2, x_1\}$ for $n = 2$.

# C KSS bijection

## C.1 Paths

Consider the product set of the form $\mathcal{B} = B^{k_1,1} \otimes B^{k_2,1} \otimes \dots \otimes B^{k_N,1}$ with $k_i \in [1, n-1]$ containing column shape tableaux only. For the definition of $B^{r,s}$, see Sec. 2.2. States of our cBBS correspond to the choice $k_i \equiv i \mod n-1$ and $N = (n-1)L$. In this appendix, we do not assume the boundary condition of the BBS, and call the elements of $\mathcal{B}$ paths. Our aim is to explain the special case of the KSS bijection [50] corresponding to the above $\mathcal{B}$ along the convention adapted to this paper.

For a path $\mathfrak{p} = c_1 \otimes \dots \otimes c_N \in \mathcal{B}$ with $c_j \in B^{k_j,1}$, the $n$-array $\lambda = (\lambda_1, \dots, \lambda_n)$ defined by

$$\mathrm{wt}(\mathfrak{p}) = (\lambda_1, \dots, \lambda_n) \in (\mathbb{Z}_{\ge 0})^n, \quad \lambda_a = \sum_{j=1}^{N} \text{(number of letter } a \text{ contained in } c_j) \qquad (C.1)$$

is called the weight of $\mathfrak{p}$. Let $c_{jl}$ be the $l$ th entry of the tableau $c_j \in B^{k_j,1}$ from the top. Consider the word

$$w_1 w_2 \dots w_{N'} := c_{11} c_{12} \dots c_{1k_1} c_{21} c_{22} \dots c_{2k_2} \dots c_{N1} c_{N2} \dots c_{Nk_N}, \qquad (C.2)$$
where the superscript of the sum means the weight constraint (C.14) on {$m_j^{(a)}$} via (C.10). This follows simply from (C.8) since the number of choices of the riggings attached to the color a length j strings are equal to the number of integers $r_1, \dots, r_{m_j^{(a)}}$ satisfying $0 \le r_1 \le \dots \le r_{m_j^{(a)}} \le$

$p_j^{(a)}$, which is the binomial coefficient in (C.16).

### C.3 KSS bijection

There is a one to one correspondence between the finite sets $P_+(\mathcal{B}, \lambda)$ and $\text{RC}(\mathcal{B}, \lambda)$. It was shown in a more general case of $\mathcal{B} = B^{k_1, s_1} \otimes \dots \otimes B^{k_N, s_N}$ in [50] generalizing the pioneering works [51, 52]. Their original motivation was to provide a bijective (combinatorial) proof of the so called Fermionic character formula for the finite dimensional representations of quantum affine algebras. It originates in the string hypothesis in the Bethe ansatz. See the introductions of [40, 53] for the rich history going back to Bethe himself [3]. Here we describe the algorithm for the bijection

$$ \Phi: P_+(\mathcal{B}, \lambda) \xrightarrow{1:1} \text{RC}(\mathcal{B}, \lambda). \qquad (C.17) $$

It works recursively with respect to $\mathcal{B}$ and $\lambda$ along the commutative diagram<fnref n="18" />:

$$ \begin{tikzcd}[column sep=2.8em, row sep=2.8em]
    \mathcal{P}_+(\mathcal{B} \otimes B^{k,1}, \lambda) \arrow[r, "\Phi"] \arrow[d, "\text{rb}"] & \text{RC}(\mathcal{B} \otimes B^{k,1}, \lambda) \arrow[d, "\delta"] \\
    \bigcup_{\nu \in \lambda^{-}} \mathcal{P}_+(\mathcal{B} \otimes B^{k-1,1}, \nu) \arrow[r, "\Phi"] & \bigcup_{\nu \in \lambda^{-}} \text{RC}(\mathcal{B} \otimes B^{k-1,1}, \nu)
\end{tikzcd} \qquad (C.18) $$

Here $\lambda^{-} = \{\lambda - \mathbf{e}_1, \lambda - \mathbf{e}_2, \dots, \lambda - \mathbf{e}_n\} \cap (\mathbb{Z}_{\ge 0})^n$ with $\mathbf{e}_a$ being the elementary vector whose only nonvanishing component is 1 at the a th position. The map rb (right box removal) is defined by

$$ \text{rb}(c_1 \otimes \cdots \otimes c_N) = \begin{cases} c_1 \otimes \cdots \otimes c_{N-1} \otimes c'_N & (k > 1), \\ c_1 \otimes \cdots \otimes c_{N-1} & (k = 1), \end{cases} \qquad (C.19) $$

where $c'_N \in B^{k-1,1}$ is obtained from the tableau $c_N \in B^{k,1}$ just by removing its bottom box and letter. If the letter removed by rb in (C.18) is $x$, then $\nu = \lambda - \mathbf{e}_x$ by the definition (C.1). From (C.3) it is clear that if $p$ is highest, $\text{rb}(p)$ is also highest. In the example (C.6), $\text{rb}(p_0)$ differs from $p_0$ only by its rightmost component which becomes [1]. The main part of the algorithm lies in $\delta$ in (C.18), which we shall explain below in two ways which fit the calculation of $\Phi^{-1}$ and $\Phi$. A string $(a_i, j_i, r_i)$ in a rigged configuration is singular if the rigging attains the allowed maximum, i.e., $r_i = p_{j_i}^{(a_i)}$.

**Algorithm for $\Phi^{-1}$**. Given $k \in [1, n-1]$ and a rigged configuration $(\mu, r) = ((\mu^{(1)}, r^{(1)}), \dots, (\mu^{(n-1)}, r^{(n-1)})) \in \text{RC}(\mathcal{B} \otimes B^{k,1}, \lambda)$ on the top right corner of (C.18), we are going to define $\text{rk}(\mu, r) \in [1, n]$ called *rank* and a new rigged configuration $\delta(\mu, r) \in \text{RC}(\mathcal{B} \otimes B^{k-1,1}, \lambda - \mathbf{e}_{\text{rk}(\mu, r)})$ which is 'smaller' than $(\mu, r)$. The rank $\text{rk}(\mu, r)$ is the letter inscribed in the removed box in (C.19) in the corresponding path $\Phi^{-1}(\mu, r)$. Therefore, once we know $\text{rk}(\mu, r)$ and $\delta(\mu, r)$, the calculation of the image $\Phi^{-1}(\mu, r)$ is reduced to $\Phi^{-1}(\delta(\mu, r))$. This reduction can be iterated until $\emptyset = \Phi^{-1}(\emptyset)$ is reached, producing the image path $p \in \Phi^{-1}(\mu, r)$ letter by letter from the right in the notation of (C.2).

<footnote n="18">¹⁸By $\Phi$ we actually mean the totality of the maps (C.17) for all $(\mathcal{B}, \lambda)$.</footnote>
where the superscript of the sum means the weight constraint (C.14) on {$m_j^{(a)}$} via (C.10). This follows simply from (C.8) since the number of choices of the riggings attached to the color a length j strings are equal to the number of integers $r_1, \dots, r_{m_j^{(a)}}$ satisfying $0 \le r_1 \le \dots \le r_{m_j^{(a)}} \le$

$p_j^{(a)}$, which is the binomial coefficient in (C.16).

### C.3 KSS bijection

There is a one to one correspondence between the finite sets $P_+(\mathcal{B}, \lambda)$ and $\text{RC}(\mathcal{B}, \lambda)$. It was shown in a more general case of $\mathcal{B} = B^{k_1, s_1} \otimes \dots \otimes B^{k_N, s_N}$ in [50] generalizing the pioneering works [51, 52]. Their original motivation was to provide a bijective (combinatorial) proof of the so called Fermionic character formula for the finite dimensional representations of quantum affine algebras. It originates in the string hypothesis in the Bethe ansatz. See the introductions of [40, 53] for the rich history going back to Bethe himself [3]. Here we describe the algorithm for the bijection

$$ \Phi: P_+(\mathcal{B}, \lambda) \xrightarrow{1:1} \text{RC}(\mathcal{B}, \lambda). \qquad (C.17) $$

It works recursively with respect to $\mathcal{B}$ and $\lambda$ along the commutative diagram<fnref n="18" />:

$$ \begin{tikzcd}[column sep=2.8em, row sep=2.8em]
    \mathcal{P}_+(\mathcal{B} \otimes B^{k,1}, \lambda) \arrow[r, "\Phi"] \arrow[d, "\text{rb}"] & \text{RC}(\mathcal{B} \otimes B^{k,1}, \lambda) \arrow[d, "\delta"] \\
    \bigcup_{\nu \in \lambda^{-}} \mathcal{P}_+(\mathcal{B} \otimes B^{k-1,1}, \nu) \arrow[r, "\Phi"] & \bigcup_{\nu \in \lambda^{-}} \text{RC}(\mathcal{B} \otimes B^{k-1,1}, \nu)
\end{tikzcd} \qquad (C.18) $$

Here $\lambda^{-} = \{\lambda - \mathbf{e}_1, \lambda - \mathbf{e}_2, \dots, \lambda - \mathbf{e}_n\} \cap (\mathbb{Z}_{\ge 0})^n$ with $\mathbf{e}_a$ being the elementary vector whose only nonvanishing component is 1 at the a th position. The map rb (right box removal) is defined by

$$ \text{rb}(c_1 \otimes \cdots \otimes c_N) = \begin{cases} c_1 \otimes \cdots \otimes c_{N-1} \otimes c'_N & (k > 1), \\ c_1 \otimes \cdots \otimes c_{N-1} & (k = 1), \end{cases} \qquad (C.19) $$

where $c'_N \in B^{k-1,1}$ is obtained from the tableau $c_N \in B^{k,1}$ just by removing its bottom box and letter. If the letter removed by rb in (C.18) is $x$, then $\nu = \lambda - \mathbf{e}_x$ by the definition (C.1). From (C.3) it is clear that if $p$ is highest, $\text{rb}(p)$ is also highest. In the example (C.6), $\text{rb}(p_0)$ differs from $p_0$ only by its rightmost component which becomes [1]. The main part of the algorithm lies in $\delta$ in (C.18), which we shall explain below in two ways which fit the calculation of $\Phi^{-1}$ and $\Phi$. A string $(a_i, j_i, r_i)$ in a rigged configuration is singular if the rigging attains the allowed maximum, i.e., $r_i = p_{j_i}^{(a_i)}$.

**Algorithm for $\Phi^{-1}$**. Given $k \in [1, n-1]$ and a rigged configuration $(\mu, r) = ((\mu^{(1)}, r^{(1)}), \dots, (\mu^{(n-1)}, r^{(n-1)})) \in \text{RC}(\mathcal{B} \otimes B^{k,1}, \lambda)$ on the top right corner of (C.18), we are going to define $\text{rk}(\mu, r) \in [1, n]$ called *rank* and a new rigged configuration $\delta(\mu, r) \in \text{RC}(\mathcal{B} \otimes B^{k-1,1}, \lambda - \mathbf{e}_{\text{rk}(\mu, r)})$ which is 'smaller' than $(\mu, r)$. The rank $\text{rk}(\mu, r)$ is the letter inscribed in the removed box in (C.19) in the corresponding path $\Phi^{-1}(\mu, r)$. Therefore, once we know $\text{rk}(\mu, r)$ and $\delta(\mu, r)$, the calculation of the image $\Phi^{-1}(\mu, r)$ is reduced to $\Phi^{-1}(\delta(\mu, r))$. This reduction can be iterated until $\emptyset = \Phi^{-1}(\emptyset)$ is reached, producing the image path $p \in \Phi^{-1}(\mu, r)$ letter by letter from the right in the notation of (C.2).

<footnote n="18">¹⁸By $\Phi$ we actually mean the totality of the maps (C.17) for all $(\mathcal{B}, \lambda)$.</footnote>
(i) Set $\ell^{(k-1)} = 1$. Do the following procedure for $a = k, k + 1, \dots, n$ in this order until stopped. Find the shortest color $a$ singular string $(a, l, p_l^{(a)})$ in $(\mu^{(a)}, r^{(a)})$ such that $\ell^{(a-1)} \le l$.<fnref n="19" /> If no such $l$ exists, set $rk(\mu, r) = a$ and stop. Otherwise set $\ell^{(a)} = l$ and continue with $a+1$. This procedure certainly ends at most with $rk(\mu, r) = n$ since there is no $(\mu^{(n)}, r^{(n)})$.

(ii) Suppose $rk(\mu, r) = b$ and let $(k, \ell^{(k)}, *)$, $(k + 1, \ell^{(k+1)}, *)$, \dots, $(b - 1, \ell^{(b-1)}, *)$ be the so found singular strings, where $*$ denotes the rigging equal to the respective vacancy.

(iii) The new rigged configuration $\delta(\mu, r)$ is obtained by replacing them with $(k, \ell^{(k)} - 1, \sharp)$, $(k+1, \ell^{(k+1)} - 1, \sharp)$, \dots, $(b-1, \ell^{(b-1)} - 1, \sharp)$ keeping all the other strings unchanged.<fnref n="20" /> Here the riggings $\sharp$ are taken so that these strings again become singular, namely, they are chosen to be the respective vacancy in the new environment $B \otimes B^{k-1,1}$. When $rk(\mu, r) = k$ in particular, we have no string to replace, hence $\delta(\mu, r) = (\mu, r)$. It is known that $\delta(\mu, r) \in RC(B \otimes B^{k-1,1}, \lambda - e_{rk(\mu,r)})$.

**Algorithm for $\Phi$.** Let $\mathfrak{p} \in \mathcal{P}_+(B \otimes B^{k-1,1}, \lambda - e_d)$ be a highest path in the bottom left corner of (C.18) for some $d \in [1, n]$. Let further $\mathfrak{p}' \in \mathcal{P}_+(B \otimes B^{k,1}, \lambda)$ be the highest path such that $rb(\mathfrak{p}') = \mathfrak{p}$. Then (C.19) tells that $\mathfrak{p}'$ is obtained from $\mathfrak{p}$ by adding a box containing $d$ to its 'bottom right' component. (Note that $d \ge k$ holds due to the semistandardness of the rightmost component tableau of $\mathfrak{p}'$.) Given such $k \in [1, n-1]$, $d \in [1, n]$ and a rigged configuration $(\mu, r) = ((\mu^{(1)}, r^{(1)}), \dots, (\mu^{(n-1)}, r^{(n-1)})) = \Phi(\mathfrak{p})$, we are going define a new one $(\mu', r') = \Phi(\mathfrak{p}')$ on the top right corner of (C.18). This enables us to go backward vertically in (C.18) and to translate the growth of the paths into that of rigged configurations starting from $\Phi(\emptyset) = \emptyset$. Here the growth means $w_1, w_1w_2, w_1w_2w_3, \dots$ in the notation of (C.2).

(i)' Set $\ell^{(d)} = \infty$. Do the following procedure for $a = d-1, d-2, \dots, k$ in this order. Find the longest color $a$ singular string $(a, l, p_l^{(a)})$ in $(\mu^{(a)}, r^{(a)})$ such that $l \le \ell^{(a+1)}$ and continue with $a-1$.<fnref n="21" /> If there is no such string, set $\ell^{(a)} = 0$ and continue with $a-1$.<fnref n="22" />

(ii)' Let $(d-1, \ell^{(d-1)}, *), (d-2, \ell^{(d-2)}, *), \dots, (k, \ell^{(k)}, *)$ be the so found singular, or 'length 0' strings, where $*$ denotes the rigging equal to the respective vacancy.<fnref n="23" />

(iii)' The new rigged configuration $(\mu', r')$ such that $\delta(\mu', r') = (\mu, r)$ is obtained by replacing them with $(d-1, \ell^{(d-1)} + 1, \sharp), (d-2, \ell^{(d-2)} + 1, \sharp), \dots, (k, \ell^{(k)} + 1, \sharp)$ keeping all the other strings unchanged.<fnref n="24" /> Here riggings $\sharp$ are taken so that these strings become singular, namely, they are chosen to be respective vacancy in the new environment $B \otimes B^{k,1}$. When $k=d$ in particular, we do not have any string to replace, hence $(\mu', r') = (\mu, r)$.

**Example C.1.** Let us show $\Phi(\mathfrak{p}_0) = S_0$, where $\mathfrak{p}_0$ is given by (C.6) and $S_0$ is by (C.11) or graphically (C.12). We utilize the detailed representation (C.13) and display how $\mathfrak{p}_0 = \Phi^{-1}(S_0)$ is constructed by applying $\delta$ successively. In each line, $rk(\mu, r)$ is calculated and it is newly inscribed in the box of the path in the next line which is just the removed one in $B$ by rb. The leftmost data $B$ corresponds to the blank boxes of the path which are yet to be determined. The image of $\Phi^{-1}$ of the rigged configuration on each line is the part of $\mathfrak{p}_0$ (C.6) that corresponds

<footnote n="19">¹⁹If there are more than one such strings, pick any one of it.</footnote>

<footnote n="20">²⁰The selected strings with length $\ell^{(c)} = 1$ (if any) are to be just removed.</footnote>

<footnote n="21">²¹If there are more than one such strings, pick any one of it.</footnote>

<footnote n="22">²²Once this happens, the rest becomes $\ell^{(a)} = \ell^{(a-1)} = \dots = \ell^{(k)} = 0$ by the definition.</footnote>

<footnote n="23">²³In case $\ell^{(a)} = 0$, $*$ is undefined but it does not matter.</footnote>

<footnote n="24">²⁴When $\ell^{(a)} = 0$, replacing $(a, \ell^{(a)}, *)$ by $(a, \ell^{(a)} + 1, \sharp)$ is to create a color $a$ length 1 singular string.</footnote>
(i) Set $\ell^{(k-1)} = 1$. Do the following procedure for $a = k, k + 1, \dots, n$ in this order until stopped. Find the shortest color $a$ singular string $(a, l, p_l^{(a)})$ in $(\mu^{(a)}, r^{(a)})$ such that $\ell^{(a-1)} \le l$.<fnref n="19" /> If no such $l$ exists, set $rk(\mu, r) = a$ and stop. Otherwise set $\ell^{(a)} = l$ and continue with $a+1$. This procedure certainly ends at most with $rk(\mu, r) = n$ since there is no $(\mu^{(n)}, r^{(n)})$.

(ii) Suppose $rk(\mu, r) = b$ and let $(k, \ell^{(k)}, *), (k+1, \ell^{(k+1)}, *), \dots, (b-1, \ell^{(b-1)}, *)$ be the so found singular strings, where * denotes the rigging equal to the respective vacancy.

(iii) The new rigged configuration $\delta(\mu, r)$ is obtained by replacing them with $(k, \ell^{(k)}-1, \sharp), (k+1, \ell^{(k+1)}-1, \sharp), \dots, (b-1, \ell^{(b-1)}-1, \sharp)$ keeping all the other strings unchanged.<fnref n="20" /> Here the riggings $\sharp$ are taken so that these strings again become singular, namely, they are chosen to be the respective vacancy in the new environment $\mathcal{B} \otimes B^{k-1,1}$. When $rk(\mu, r) = k$ in particular, we have no string to replace, hence $\delta(\mu, r) = (\mu, r)$. It is known that $\delta(\mu, r) \in \text{RC}(\mathcal{B} \otimes B^{k-1,1}, \lambda - \mathbf{e}_{rk(\mu,r)})$.

**Algorithm for $\Phi$.** Let $\mathfrak{p} \in \mathcal{P}_+(\mathcal{B} \otimes B^{k-1,1}, \lambda - \mathbf{e}_d)$ be a highest path in the bottom left corner of (C.18) for some $d \in [1, n]$. Let further $\mathfrak{p}' \in \mathcal{P}_+(\mathcal{B} \otimes B^{k,1}, \lambda)$ be the highest path such that $rb(\mathfrak{p}') = \mathfrak{p}$. Then (C.19) tells that $\mathfrak{p}'$ is obtained from $\mathfrak{p}$ by adding a box containing $d$ to its 'bottom right' component. (Note that $d \ge k$ holds due to the semistandardness of the rightmost component tableau of $\mathfrak{p}'$.) Given such $k \in [1, n-1]$, $d \in [1, n]$ and a rigged configuration $(\mu, r) = ((\mu^{(1)}, r^{(1)}), \dots, (\mu^{(n-1)}, r^{(n-1)})) = \Phi(\mathfrak{p})$, we are going define a new one $(\mu', r') = \Phi(\mathfrak{p}')$ on the top right corner of (C.18). This enables us to go backward vertically in (C.18) and to translate the growth of the paths into that of rigged configurations starting from $\Phi(\emptyset) = \emptyset$. Here the growth means $w_1, w_1w_2, w_1w_2w_3, \dots$ in the notation of (C.2).

(i)' Set $\ell^{(d)} = \infty$. Do the following procedure for $a = d-1, d-2, \dots, k$ in this order. Find the longest color $a$ singular string $(a, l, p_l^{(a)})$ in $(\mu^{(a)}, r^{(a)})$ such that $l \le \ell^{(a+1)}$ and continue with $a-1$.<fnref n="21" /> If there is no such string, set $\ell^{(a)} = 0$ and continue with $a-1$.<fnref n="22" />

(ii)' Let $(d-1, \ell^{(d-1)}, *), (d-2, \ell^{(d-2)}, *), \dots, (k, \ell^{(k)}, *)$ be the so found singular, or 'length 0' strings, where * denotes the rigging equal to the respective vacancy.<fnref n="23" />

(iii)' The new rigged configuration $(\mu', r')$ such that $\delta(\mu', r') = (\mu, r)$ is obtained by replacing them with $(d-1, \ell^{(d-1)} + 1, \sharp), (d-2, \ell^{(d-2)} + 1, \sharp), \dots, (k, \ell^{(k)} + 1, \sharp)$ keeping all the other strings unchanged.<fnref n="24" /> Here riggings $\sharp$ are taken so that these strings become singular, namely, they are chosen to be respective vacancy in the new environment $\mathcal{B} \otimes B^{k,1}$. When $k=d$ in particular, we do not have any string to replace, hence $(\mu', r') = (\mu, r)$.

**Example C.1.** Let us show $\Phi(\mathfrak{p}_0) = S_0$, where $\mathfrak{p}_0$ is given by (C.6) and $S_0$ is by (C.11) or graphically (C.12). We utilize the detailed representation (C.13) and display how $\mathfrak{p}_0 = \Phi^{-1}(S_0)$ is constructed by applying $\delta$ successively. In each line, $rk(\mu, r)$ is calculated and it is newly inscribed in the box of the path in the next line which is just the removed one in $\mathcal{B}$ by rb. The leftmost data $\mathcal{B}$ corresponds to the blank boxes of the path which are yet to be determined. The image of $\Phi^{-1}$ of the rigged configuration on each line is the part of $\mathfrak{p}_0$ (C.6) that corresponds

<footnote n="19">¹⁹If there are more than one such strings, pick any one of it.</footnote>

<footnote n="20">²⁰The selected strings with length $\ell^{(c)} = 1$ (if any) are to be just removed.</footnote>

<footnote n="21">²¹If there are more than one such strings, pick any one of it.</footnote>

<footnote n="22">²²Once this happens, the rest becomes $\ell^{(a)} = \ell^{(a-1)} = \dots = \ell^{(k)} = 0$ by the definition.</footnote>

<footnote n="23">²³In case $\ell^{(a)} = 0$, * is undefined but it does not matter.</footnote>

<footnote n="24">²⁴When $\ell^{(a)} = 0$, replacing $(a, \ell^{(a)}, *)$ by $(a, \ell^{(a)} + 1, \sharp)$ is to create a color $a$ length 1 singular string.</footnote>
to the unfilled boxes of the path on the same line.

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

(C.20)

After this, we will only get vac$^{\otimes 3}$ and end up with $p_0$ in (C.6).

Let us finish with a remark on the situation of our cBBS. The states are elements of $B^{\otimes L}$ of the form vac$^{\otimes j-1} \otimes b_j \otimes \cdots \otimes b_l \otimes \text{vac}^{\otimes L-l}$ satisfying the boundary condition $1 \ll j \ll l \ll L$. Taking $j \gg 1$ assures the highest condition (C.3) automatically for any $b_j \otimes \cdots \otimes b_l$. So we can consider the image of the BBS states by $\Phi$. From the definition of $B$ in (2.21), the vacancy (C.10) becomes (3.3). Thus $l \ll L$ implies $\forall p_j^{(a)} \gg 1$, which allows us to increase the rigging of $\Phi(\text{state})$ without breaking the condition (C.8). As the above example indicates, if
to the unfilled boxes of the path on the same line.

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

(C.20)

After this, we will only get vac$^{\otimes 3}$ and end up with $p_0$ in (C.6).

Let us finish with a remark on the situation of our cBBS. The states are elements of $B^{\otimes L}$ of the form vac$^{\otimes j-1} \otimes b_j \otimes \cdots \otimes b_l \otimes \text{vac}^{\otimes L-l}$ satisfying the boundary condition $1 \ll j \ll l \ll L$. Taking $j \gg 1$ assures the highest condition (C.3) automatically for any $b_j \otimes \cdots \otimes b_l$. So we can consider the image of the BBS states by $\Phi$. From the definition of $B$ in (2.21), the vacancy (C.10) becomes (3.3). Thus $l \ll L$ implies $\forall p_j^{(a)} \gg 1$, which allows us to increase the rigging of $\Phi(\text{state})$ without breaking the condition (C.8). As the above example indicates, if
$\Phi(\mathfrak{p}) = (\mu, r) \in \text{RC}(B^{\otimes L}, \lambda)$, supplement of the vacuum tail does not influence the rigged configuration in the sense that $\Phi(\mathfrak{p} \otimes \text{vac}^{\otimes l}) = (\mu, r) \in \text{RC}(B^{\otimes L+l}, \lambda + ((n-1)l, (n-2)l, \dots, l, 0))$. Similarly, supplement of the vacuum in front just causes a uniform shift of the rigging $\Phi(\text{vac}^{\otimes l} \otimes \mathfrak{p}) = (\mu, r+l) \in \text{RC}(B^{\otimes L+l}, \lambda + ((n-1)l, (n-2)l, \dots, l, 0))$, where $r+l$ means that the rigging of every string gets increased by $l$.

## Acknowledgments

A.K. is supported by Grants-in-Aid for Scientific Research No. 16H03922, 18H01141 and 18K03452 from JSPS.

## References

[1] D. Takahashi and J. Satsuma, *A soliton cellular automaton*, J. Phys. Soc. Jpn. **59**, 3514 (1990), doi:10.1143/JPSJ.59.3514.

[2] R. Baxter, *Exactly solved models in statistical mechanics*, Dover Publications, Mineola, New York, ISBN 9780486462714 (2008).

[3] H. Bethe, *Zur Theorie der Metalle*, Z. Physik **71**, 205 (1931), doi:10.1007/BF01341708.

[4] M. Kashiwara, *On crystal bases of the Q-analogue of universal enveloping algebras*, Duke Math. J. **63**, 465 (1991), doi:10.1215/S0012-7094-91-06321-0.

[5] R. Inoue, A. Kuniba and T. Takagi, *Integrable structure of box-ball systems: crystal, Bethe ansatz, ultradiscretization and tropical geometry*, J. Phys. A: Math. Theor. **45**, 073001 (2012), doi:10.1088/1751-8113/45/7/073001.

[6] C. N. Yang and C. P. Yang, *Thermodynamics of a one-dimensional system of bosons with repulsive delta-function interaction*, J. Math. Phys. **10**, 1115 (1969), doi:10.1063/1.1664947.

[7] M. Takahashi, *Thermodynamics of one-dimensional solvable models*, Cambridge University Press, Cambridge, United Kingdom, ISBN 9780521551434 (1999), doi:10.1017/CBO9780511524332.

[8] B. Bertini, M. Collura, J. De Nardis and M. Fagotti, *Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents*, Phys. Rev. Lett. **117**, 207201 (2016), doi:10.1103/PhysRevLett.117.207201.

[9] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura, *Emergent hydrodynamics in integrable quantum systems out of equilibrium*, Phys. Rev. X **6**, 041065 (2016), doi:10.1103/PhysRevX.6.041065.

[10] B. Doyon, *Lecture notes on generalised hydrodynamics*, SciPost Phys. Lect. Notes, 18 (2020), doi:10.21468/SciPostPhysLectNotes.18.

[11] A. Kuniba, G. Misguich and V. Pasquier, *Generalized hydrodynamics in box-ball system*, J. Phys. A: Math. Theor. **53**, 404001 (2020), doi:10.1088/1751-8121/abadb9.

[12] W. Fulton, *Young tableaux: With applications to representation theory and geometry*, London Mathematical Society Student Texts, Cambridge University Press, Cambridge, ISBN 9780521567244 (1996).
$\Phi(\mathfrak{p}) = (\mu, r) \in \text{RC}(B^{\otimes L}, \lambda)$, supplement of the vacuum tail does not influence the rigged configuration in the sense that $\Phi(\mathfrak{p} \otimes \text{vac}^{\otimes l}) = (\mu, r) \in \text{RC}(B^{\otimes L+l}, \lambda + ((n-1)l, (n-2)l, \dots, l, 0))$. Similarly, supplement of the vacuum in front just causes a uniform shift of the rigging $\Phi(\text{vac}^{\otimes l} \otimes \mathfrak{p}) = (\mu, r+l) \in \text{RC}(B^{\otimes L+l}, \lambda + ((n-1)l, (n-2)l, \dots, l, 0))$, where $r+l$ means that the rigging of every string gets increased by $l$.

## Acknowledgments

A.K. is supported by Grants-in-Aid for Scientific Research No. 16H03922, 18H01141 and 18K03452 from JSPS.

## References

[1] D. Takahashi and J. Satsuma, *A soliton cellular automaton*, J. Phys. Soc. Jpn. **59**, 3514 (1990), doi:10.1143/JPSJ.59.3514.

[2] R. Baxter, *Exactly solved models in statistical mechanics*, Dover Publications, Mineola, New York, ISBN 9780486462714 (2008).

[3] H. Bethe, *Zur Theorie der Metalle*, Z. Physik **71**, 205 (1931), doi:10.1007/BF01341708.

[4] M. Kashiwara, *On crystal bases of the Q-analogue of universal enveloping algebras*, Duke Math. J. **63**, 465 (1991), doi:10.1215/S0012-7094-91-06321-0.

[5] R. Inoue, A. Kuniba and T. Takagi, *Integrable structure of box-ball systems: crystal, Bethe ansatz, ultradiscretization and tropical geometry*, J. Phys. A: Math. Theor. **45**, 073001 (2012), doi:10.1088/1751-8113/45/7/073001.

[6] C. N. Yang and C. P. Yang, *Thermodynamics of a one-dimensional system of bosons with repulsive delta-function interaction*, J. Math. Phys. **10**, 1115 (1969), doi:10.1063/1.1664947.

[7] M. Takahashi, *Thermodynamics of one-dimensional solvable models*, Cambridge University Press, Cambridge, United Kingdom, ISBN 9780521551434 (1999), doi:10.1017/CBO9780511524332.

[8] B. Bertini, M. Collura, J. De Nardis and M. Fagotti, *Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents*, Phys. Rev. Lett. **117**, 207201 (2016), doi:10.1103/PhysRevLett.117.207201.

[9] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura, *Emergent hydrodynamics in integrable quantum systems out of equilibrium*, Phys. Rev. X **6**, 041065 (2016), doi:10.1103/PhysRevX.6.041065.

[10] B. Doyon, *Lecture notes on generalised hydrodynamics*, SciPost Phys. Lect. Notes, 18 (2020), doi:10.21468/SciPostPhysLectNotes.18.

[11] A. Kuniba, G. Misguich and V. Pasquier, *Generalized hydrodynamics in box-ball system*, J. Phys. A: Math. Theor. **53**, 404001 (2020), doi:10.1088/1751-8121/abadb9.

[12] W. Fulton, *Young tableaux: With applications to representation theory and geometry*, London Mathematical Society Student Texts, Cambridge University Press, Cambridge, ISBN 9780521567244 (1996).
[13] I. G. Macdonald, *Symmetric functions and Hall polynomials*, Oxford Classic Texts in the Physical Sciences, Oxford University Press, Oxford, New York, ISBN 9780198739128 (1995).

[14] D. Takahashi, *On some soliton systems defined by using boxes and balls*, Proc. Inter. Symp. Nonlin. Theory Appl. **93**, 555 (1993).

[15] G. Hatayama, K. Hikami, R. Inoue, A. Kuniba, T. Takagi and T. Tokihiro, *The $A_M^{(1)}$ automata related to crystals of symmetric tensors*, J. Math. Phys. **42**, 274 (2001), doi:10.1063/1.1322077.

[16] K. Fukuda, Y. Yamada and M. Okado, *Energy functions in box ball systems*, Int. J. Mod. Phys. A **15**, 1379 (2000), doi:10.1142/S0217751X00000616.

[17] A. Kuniba, H. Lyu and M. Okado, *Randomized box-ball systems, limit shape of rigged configurations and thermodynamic Bethe ansatz*, Nucl. Phys. B **937**, 240 (2018), doi:10.1016/j.nuclphysb.2018.10.008.

[18] G. Misguich, N. Pavloff and V. Pasquier, *Domain wall problem in the quantum XXZ chain and semiclassical behavior close to the isotropic point*, SciPost Phys. **7**, 025 (2019), doi:10.21468/SciPostPhys.7.2.025.

[19] B. Doyon and H. Spohn, *Dynamics of hard rods with initial domain wall state*, J. Stat. Mech. 073210 (2017), doi:10.1088/1742-5468/aa7abf.

[20] D. Yamada, *Box ball system associated with antisymmetric tensor crystals*, J. Phys. A: Math. Gen. **37**, 9975 (2004), doi:10.1088/0305-4470/37/42/010.

[21] A. Kuniba, M. Okado and Y. Yamada, *Box-Ball system with reflecting end*, J. Nonlinear Math. Phys. **12**, 475 (2005), doi:10.2991/jnmp.2005.12.4.4.

[22] P. A. Ferrari and D. Gabrielli, *BBS invariant measures with independent soliton components*, Electron. J. Probab. **25**, 1 (2020), doi:10.1214/20-EJP475.

[23] P. A. Ferrari, C. Nguyen, L. T. Rolla and M. Wang, *Soliton decomposition of the Box-Ball system*, (2018), arXiv:1806.02798.

[24] D. A. Croydon and M. Sasada, *Invariant measures for the box-ball system based on stationary Markov chains and periodic Gibbs measures*, J. Math. Phys. **60**, 083301 (2019), doi:10.1063/1.5095622.

[25] D. A. Croydon and M. Sasada, *Generalized hydrodynamic limit for the Box-Ball system*, Commun. Math. Phys. **383**, 427 (2021), doi:10.1007/s00220-020-03914-x.

[26] J. Lewis, H. Lyu, P. Pylyavskyy and A. Sen, *Scaling limit of soliton lengths in a multicolor box-ball system*, (2019), arXiv:1911.04458.

[27] A. Kuniba and H. Lyu, *Large deviations and one-sided scaling limit of randomized multicolor Box-Ball system*, J. Stat. Phys. **178**, 38 (2020), doi:10.1007/s10955-019-02417-x.

[28] V. G. Kac, *Infinite-dimensional Lie algebras*, Cambridge University Press, Cambridge, United Kingdom, ISBN 9780521372152 (1990), doi:10.1017/CBO9780511626234.

[29] V. G. Drinfeld, *Quantum groups*, In Proceedings of the International Congress of Mathematicians, Amer. Math. Soc. **1-2**, 798 (1986).
[13] I. G. Macdonald, *Symmetric functions and Hall polynomials*, Oxford Classic Texts in the Physical Sciences, Oxford University Press, Oxford, New York, ISBN 9780198739128 (1995).

[14] D. Takahashi, *On some soliton systems defined by using boxes and balls*, Proc. Inter. Symp. Nonlin. Theory Appl. **93**, 555 (1993).

[15] G. Hatayama, K. Hikami, R. Inoue, A. Kuniba, T. Takagi and T. Tokihiro, *The $A_M^{(1)}$ automata related to crystals of symmetric tensors*, J. Math. Phys. **42**, 274 (2001), doi:10.1063/1.1322077.

[16] K. Fukuda, Y. Yamada and M. Okado, *Energy functions in box ball systems*, Int. J. Mod. Phys. A **15**, 1379 (2000), doi:10.1142/S0217751X00000616.

[17] A. Kuniba, H. Lyu and M. Okado, *Randomized box-ball systems, limit shape of rigged configurations and thermodynamic Bethe ansatz*, Nucl. Phys. B **937**, 240 (2018), doi:10.1016/j.nuclphysb.2018.10.008.

[18] G. Misguich, N. Pavloff and V. Pasquier, *Domain wall problem in the quantum XXZ chain and semiclassical behavior close to the isotropic point*, SciPost Phys. **7**, 025 (2019), doi:10.21468/SciPostPhys.7.2.025.

[19] B. Doyon and H. Spohn, *Dynamics of hard rods with initial domain wall state*, J. Stat. Mech. 073210 (2017), doi:10.1088/1742-5468/aa7abf.

[20] D. Yamada, *Box ball system associated with antisymmetric tensor crystals*, J. Phys. A: Math. Gen. **37**, 9975 (2004), doi:10.1088/0305-4470/37/42/010.

[21] A. Kuniba, M. Okado and Y. Yamada, *Box-Ball system with reflecting end*, J. Nonlinear Math. Phys. **12**, 475 (2005), doi:10.2991/jnmp.2005.12.4.4.

[22] P. A. Ferrari and D. Gabrielli, *BBS invariant measures with independent soliton components*, Electron. J. Probab. **25**, 1 (2020), doi:10.1214/20-EJP475.

[23] P. A. Ferrari, C. Nguyen, L. T. Rolla and M. Wang, *Soliton decomposition of the Box-Ball system*, (2018), arXiv:1806.02798.

[24] D. A. Croydon and M. Sasada, *Invariant measures for the box-ball system based on stationary Markov chains and periodic Gibbs measures*, J. Math. Phys. **60**, 083301 (2019), doi:10.1063/1.5095622.

[25] D. A. Croydon and M. Sasada, *Generalized hydrodynamic limit for the Box-Ball system*, Commun. Math. Phys. **383**, 427 (2021), doi:10.1007/s00220-020-03914-x.

[26] J. Lewis, H. Lyu, P. Pylyavskyy and A. Sen, *Scaling limit of soliton lengths in a multicolor box-ball system*, (2019), arXiv:1911.04458.

[27] A. Kuniba and H. Lyu, *Large deviations and one-sided scaling limit of randomized multicolor Box-Ball system*, J. Stat. Phys. **178**, 38 (2020), doi:10.1007/s10955-019-02417-x.

[28] V. G. Kac, *Infinite-dimensional Lie algebras*, Cambridge University Press, Cambridge, United Kingdom, ISBN 9780521372152 (1990), doi:10.1017/CBO9780511626234.

[29] V. G. Drinfeld, *Quantum groups*, In Proceedings of the International Congress of Mathematicians, Amer. Math. Soc. **1-2**, 798 (1986).
[30] M. Jimbo, A *q*-difference analogue of $U(g)$ and the Yang-Baxter equation, Lett. Math. Phys. **10**, 63 (1985), doi:10.1007/BF00704588.

[31] A. N. Kirillov and N. Yu. Reshetikhin, *Representations of Yangians and multiplicities of occurrence of the irreducible components of the tensor product of representations of simple Lie algebras*, J. Math. Sci. **52**, 3156 (1990), doi:10.1007/BF02342935.

[32] S.-J. Kang, M. Kashiwara, K. C. Misra, T. Miwa, T. Nakashima and A. Nakayashiki, *Perfect crystals of quantum affine Lie algebras*, Duke Math. J. **68**, 499 (1992), doi:10.1215/S0012-7094-92-06821-9.

[33] S.-J. Kang, M. Kashiwara, K. C. Misra, T. Miwa, T. Nakashima and A. Nakayashiki, *Affine crystals and vertex models*, Int. J. Mod. Phys. A **07**, 449 (1992), doi:10.1142/S0217751X92003896.

[34] M. Shimozono, *Affine type a crystal structure on tensor products of rectangles, Demazure characters, and nilpotent varieties*, J. Algebraic Comb. **15**, 151 (2002), doi:10.1023/A:1013894920862.

[35] M. Okado, *X = M conjecture*, Math. Soc. Japan Memoirs, ISBN 9784931469372 (2007), doi:10.2969/msjmemoirs/01701C030.

[36] D. Takahashi and J. Matsukidaira, *Box and ball system with a carrier and ultradiscrete modified KdV equation*, J. Phys. A: Math. Gen. **30**, L733 (1997), doi:10.1088/0305-4470/30/21/005.

[37] A. Kuniba, T. Nakanishi and J. Suzuki, *T-systems and Y-systems in integrable systems*, J. Phys. A: Math. Theor. **44**, 103001 (2011), doi:10.1088/1751-8113/44/10/103001.

[38] A. Kuniba, M. Okado, R. Sakamoto, T. Takagi and Y. Yamada, *Crystal interpretation of Kerov-Kirillov-Reshetikhin bijection*, Nucl. Phys. B **740**, 299 (2006), doi:10.1016/j.nuclphysb.2006.02.005.

[39] R. Sakamoto, *Kirillov-Schilling-Shimozono bijection as energy functions of crystals*, Int. Math. Res. Not. 579 (2009), doi:10.1093/imrn/rnn140.

[40] G. Hatayama, A. Kuniba, M. Okado, T. Takagi and Y. Yamada, *Remarks on fermionic formula*, American Mathematical Society, Providence, Rhode Island, ISBN 9780821811993 (1999), doi:10.1090/conm/248.

[41] V. E. Zakharov, *Kinetic equation for solitons*, Sov. JETP **33**, 538 (1971), http://www.jetp.ac.ru/cgi-bin/e/index/e/33/3/p538?a=list.

[42] G. A. El, *The thermodynamic limit of the Whitham equations*, Phys. Lett. A **311**, 374 (2003), doi:10.1016/S0375-9601(03)00515-2.

[43] B. Doyon, T. Yoshimura and J.-S. Caux, *Soliton gases and generalized hydrodynamics*, Phys. Rev. Lett. **120**, 045301 (2018), doi:10.1103/PhysRevLett.120.045301.

[44] A. Kuniba and T. Takagi, *Bethe Ansatz, inverse scattering transform and tropical Riemann theta function in a periodic soliton cellular automaton for $A_n^{(1)}$, SIGMA* (2010), doi:10.3842/SIGMA.2010.013.

[45] J. De Nardis, D. Bernard and B. Doyon, *Hydrodynamic diffusion in integrable systems*, Phys. Rev. Lett. **121**, 160603 (2018), doi:10.1103/PhysRevLett.121.160603.
[30] M. Jimbo, A *q*-difference analogue of $U(g)$ and the Yang-Baxter equation, Lett. Math. Phys. **10**, 63 (1985), doi:10.1007/BF00704588.

[31] A. N. Kirillov and N. Yu. Reshetikhin, *Representations of Yangians and multiplicities of occurrence of the irreducible components of the tensor product of representations of simple Lie algebras*, J. Math. Sci. **52**, 3156 (1990), doi:10.1007/BF02342935.

[32] S.-J. Kang, M. Kashiwara, K. C. Misra, T. Miwa, T. Nakashima and A. Nakayashiki, *Perfect crystals of quantum affine Lie algebras*, Duke Math. J. **68**, 499 (1992), doi:10.1215/S0012-7094-92-06821-9.

[33] S.-J. Kang, M. Kashiwara, K. C. Misra, T. Miwa, T. Nakashima and A. Nakayashiki, *Affine crystals and vertex models*, Int. J. Mod. Phys. A **07**, 449 (1992), doi:10.1142/S0217751X92003896.

[34] M. Shimozono, *Affine type a crystal structure on tensor products of rectangles, Demazure characters, and nilpotent varieties*, J. Algebraic Comb. **15**, 151 (2002), doi:10.1023/A:1013894920862.

[35] M. Okado, *X = M conjecture*, Math. Soc. Japan Memoirs, ISBN 9784931469372 (2007), doi:10.2969/msjmemoirs/01701C030.

[36] D. Takahashi and J. Matsukidaira, *Box and ball system with a carrier and ultradiscrete modified KdV equation*, J. Phys. A: Math. Gen. **30**, L733 (1997), doi:10.1088/0305-4470/30/21/005.

[37] A. Kuniba, T. Nakanishi and J. Suzuki, *T-systems and Y-systems in integrable systems*, J. Phys. A: Math. Theor. **44**, 103001 (2011), doi:10.1088/1751-8113/44/10/103001.

[38] A. Kuniba, M. Okado, R. Sakamoto, T. Takagi and Y. Yamada, *Crystal interpretation of Kerov-Kirillov-Reshetikhin bijection*, Nucl. Phys. B **740**, 299 (2006), doi:10.1016/j.nuclphysb.2006.02.005.

[39] R. Sakamoto, *Kirillov-Schilling-Shimozono bijection as energy functions of crystals*, Int. Math. Res. Not. 579 (2009), doi:10.1093/imrn/rnn140.

[40] G. Hatayama, A. Kuniba, M. Okado, T. Takagi and Y. Yamada, *Remarks on fermionic formula*, American Mathematical Society, Providence, Rhode Island, ISBN 9780821811993 (1999), doi:10.1090/conm/248.

[41] V. E. Zakharov, *Kinetic equation for solitons*, Sov. JETP **33**, 538 (1971), http://www.jetp.ac.ru/cgi-bin/e/index/e/33/3/p538?a=list.

[42] G. A. El, *The thermodynamic limit of the Whitham equations*, Phys. Lett. A **311**, 374 (2003), doi:10.1016/S0375-9601(03)00515-2.

[43] B. Doyon, T. Yoshimura and J.-S. Caux, *Soliton gases and generalized hydrodynamics*, Phys. Rev. Lett. **120**, 045301 (2018), doi:10.1103/PhysRevLett.120.045301.

[44] A. Kuniba and T. Takagi, *Bethe Ansatz, inverse scattering transform and tropical Riemann theta function in a periodic soliton cellular automaton for $A_n^{(1)}$, SIGMA* (2010), doi:10.3842/SIGMA.2010.013.

[45] J. De Nardis, D. Bernard and B. Doyon, *Hydrodynamic diffusion in integrable systems*, Phys. Rev. Lett. **121**, 160603 (2018), doi:10.1103/PhysRevLett.121.160603.
[46] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur, *Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems*, Phys. Rev. B **98**, 220303 (2018), doi:10.1103/PhysRevB.98.220303.

[47] M. Medenjak, K. Klobas and T. Prosen, *Diffusion in deterministic interacting lattice systems*, Phys. Rev. Lett. **119**, 110603 (2017), doi:10.1103/PhysRevLett.119.110603.

[48] K. Klobas, M. Medenjak, T. Prosen and M. Vanicat, *Time-dependent matrix product ansatz for interacting reversible dynamics*, Commun. Math. Phys. **371**, 651 (2019), doi:10.1007/s00220-019-03494-5.

[49] M. Kollar, F. A. Wolf and M. Eckstein, *Generalized Gibbs ensemble prediction of prethermalization plateaus and their relation to nonthermal steady states in integrable systems*, Phys. Rev. B **84**, 054304 (2011), doi:10.1103/PhysRevB.84.054304.

[50] A. N. Kirillov, A. Schilling and M. Shimozono, *A bijection between Littlewood-Richardson tableaux and rigged configurations*, Sel. Math. **8**, 67 (2002), doi:10.1007/s00029-002-8102-6.

[51] S. V. Kerov, A. N. Kirillov and N. Yu. Reshetikhin, *Combinatorics, Bethe Ansatz, and representations of the symmetric group*, J. Math. Sci. **41**, 916 (1988), doi:10.1007/BF01247087.

[52] A. N. Kirillov and N. Yu. Reshetikhin, *The Bethe Ansatz and the combinatorics of Young tableaux*, J. Math. Sci. **41**, 925 (1988), doi:10.1007/BF01247088.

[53] G. Hatayama, A. Kuniba, M. Okado, T. Takagi and Z. Tsuboi, *Paths, crystals and fermionic formulae*, in MathPhys Odyssey 2001, Birkhäuser Boston, Boston, MA, ISBN 9781461266051 (2002), doi:10.1007/978-1-4612-0087-1_9.

[54] P. Di Francesco and R. Kedem, *Proof of the combinatorial Kirillov-Reshetikhin conjecture*, Int. Math. Res. Not. 006 (2008), doi:10.1093/imrn/rnn006.
[46] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur, *Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems*, Phys. Rev. B **98**, 220303 (2018), doi:10.1103/PhysRevB.98.220303.

[47] M. Medenjak, K. Klobas and T. Prosen, *Diffusion in deterministic interacting lattice systems*, Phys. Rev. Lett. **119**, 110603 (2017), doi:10.1103/PhysRevLett.119.110603.

[48] K. Klobas, M. Medenjak, T. Prosen and M. Vanicat, *Time-dependent matrix product ansatz for interacting reversible dynamics*, Commun. Math. Phys. **371**, 651 (2019), doi:10.1007/s00220-019-03494-5.

[49] M. Kollar, F. A. Wolf and M. Eckstein, *Generalized Gibbs ensemble prediction of prethermalization plateaus and their relation to nonthermal steady states in integrable systems*, Phys. Rev. B **84**, 054304 (2011), doi:10.1103/PhysRevB.84.054304.

[50] A. N. Kirillov, A. Schilling and M. Shimozono, *A bijection between Littlewood-Richardson tableaux and rigged configurations*, Sel. Math. **8**, 67 (2002), doi:10.1007/s00029-002-8102-6.

[51] S. V. Kerov, A. N. Kirillov and N. Yu. Reshetikhin, *Combinatorics, Bethe Ansatz, and representations of the symmetric group*, J. Math. Sci. **41**, 916 (1988), doi:10.1007/BF01247087.

[52] A. N. Kirillov and N. Yu. Reshetikhin, *The Bethe Ansatz and the combinatorics of Young tableaux*, J. Math. Sci. **41**, 925 (1988), doi:10.1007/BF01247088.

[53] G. Hatayama, A. Kuniba, M. Okado, T. Takagi and Z. Tsuboi, *Paths, crystals and fermionic formulae*, in MathPhys Odyssey 2001, Birkhäuser Boston, Boston, MA, ISBN 9781461266051 (2002), doi:10.1007/978-1-4612-0087-1_9.

[54] P. Di Francesco and R. Kedem, *Proof of the combinatorial Kirillov-Reshetikhin conjecture*, Int. Math. Res. Not. 006 (2008), doi:10.1093/imrn/rnn006.