# Current correlations, Drude weights and large deviations in a box-ball
system

Atsuo Kuniba1, Grégoire Misguich2 and Vincent Pasquier2 1 Institute of
Physics, University of Tokyo, Komaba, Tokyo 153-8902, Japan  2 Université
Paris-Saclay, CNRS, CEA, Institut de physique théorique, 91191, Gif-sur-
Yvette, France<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We explore several aspects of the current fluctuations and correlations in the
box-ball system (BBS), an integrable cellular automaton in one space
dimension. The state we consider is an ensemble of microscopic configurations
where the box occupancies are independent random variables (i.i.d. state),
with a given mean ball density. We compute several quantities exactly in such
homogeneous stationary state: the mean value and the variance of the number of
balls $N_{t}$ crossing the origin during time $t$, and the scaled cumulants
generating function associated to $N_{t}$. We also compute two spatially
integrated current-current correlations. The first one, involving the long-
time limit of the current-current correlations, is the so-called Drude weight
and is obtained with thermodynamic Bethe Ansatz (TBA). The second one,
involving equal time current-current correlations is calculated using a
transfer matrix approach. A family of generalized currents, associated to the
conserved charges and to the different time evolutions of the models are
constructed. The long-time limits of their correlations generalize the Drude
weight and the second cumulant of $N_{t}$ and are found to obey nontrivial
symmetry relations. They are computed using TBA and the results are found to
be in good agreement with microscopic simulations of the model. TBA is also
used to compute explicitly the whole family of flux Jacobian matrices.
Finally, some of these results are extended to a (non-i.i.d.) two-temperatures
generalized Gibbs state (with one parameter coupled to the total number of
balls, and another one coupled to the total number of solitons).

,

###### Contents

  1. 1 Introduction
  2. 2 Current correlations and Drude weights
    1. 2.1 Mean current
    2. 2.2 Second cumulant
      1. 2.2.1 Correlations and sum rules.
      2. 2.2.2 Direct calculation.
      3. 2.2.3 Second cumulant using TBA.
    3. 2.3 Drude weights
    4. 2.4 Generalized current correlations and their symmetries
    5. 2.5 Numerical calculation of the Drude weights
    6. 2.6 Flux Jacobian and related matrices
      1. 2.6.1 Normal modes
      2. 2.6.2 Family of flux Jacobians
      3. 2.6.3 Covariance and Drude matrices
    7. 2.7 Equal-time current-current correlations
  3. 3 Large deviation function
    1. 3.1 General method
    2. 3.2 Scaled cumulants generating function
    3. 3.3 $l=\infty$
  4. 4 Conclusions
  5. 5 Acknowledgements
  6. A Correlation sum rules
  7. B Transfer matrix calculation for the current fluctuations
  8. C Drude weight and SCGF in a two-temperature GGE
  9. D Microscopic definition of $\hat{\eta}^{(l)}_{i}(x)$

##  1 Introduction

The box ball system (BBS) is an integrable cellular automaton introduced in
1990 by Takahashi and Satsuma [1]. In this model some “balls” occupy the sites
(“boxes”) of a one-dimensional lattice and propagate according to some simple
deterministic rules, as explained below. We start from an initial
configuration of balls divided into $L$ boxes, with at most one ball per box.
The configuration at the next time step is obtained by letting a “carrier”
travel through the system from left to right. Doing so, each time the carrier
passes over an occupied box and if it has not reached its maximum load $l$, it
loads the ball and leaves the box empty. Each time the carrier carries at
least one ball and passes over an empty box it unloads a ball in the box. An
example is given in Fig. 1.111The dynamics can also be defined on a periodic
system by initializing the carrier load to some suitable value (Proposition
5.1 of [2]).

Figure 1: Top line: ball configuration. Bottom line: configuration at the next
time step. We assumed here a carrier capacity $l\geq 3$.

Despite its apparent simplicity this model possesses stable solitons with
nontrivial scattering under collisions. It has an infinity of conserved
quantities, and its very rich mathematical and integrable structures have
attracted a lot of interest. For example, the above combinatorial rule for the
time evolution has its origin in the quantum $R$ matrix at $q=0$ as explained
in D (see [2] for a review).

An interesting family of problems arises when considering a statistical
ensemble of random microscopic ball configurations [3, 4, 5, 6, 7, 8, 9, 10,
11]. In the simplest case the occupancies of the boxes can be taken to be
independent and identically distributed (i.i.d.) random variables, and
parameterized with a single parameter (the mean ball density). Such a
statistical ensemble corresponds to an homogeneous and stationary state, and
many quantities, like the densities of the various types of solitons and their
associated mean velocities, can be computed exactly using the thermodynamic
Bethe Ansatz (TBA) [5, 11, 12].

It has also been shown that a hydrodynamic approach can accurately capture the
long-time and large distance evolution of some inhomogeneous states. Due to
the extensive number of conserved quantities in the model, the appropriate
framework is the so-called generalized hydrodynamics (GHD). Introduced a few
years ago [13, 14, 15], the GHD is based on the assumption that the system is
locally in an homogeneous state with maximum entropy (generalized Gibbs
ensemble (GGE) [16]) and the approach amounts to constructing and solving the
set of continuity equations involving the local currents and local densities
associated to the conserved quantities of an integrable system. It has been
applied successfully to many integrable quantum and classical systems. On the
classical side, we can mention for instance the application of GHD to the
hard-rods model [17, 18], to the Toda model [19], to the sinh-Gordon model
[20], to the BBS [10, 12] and to the (higher rank) complete BBS [21].

In the context of BBS, the evolution triggered by an initial domain-wall state
with two different ball densities in the left half and the right half has been
studied in details using GHD [12, 21]. In such a setup the BBS develops a
series of plateaux in the variable $\zeta=x/t$ (space coordinate divided by
time). The ball density and soliton content inside each plateau as well as the
location $\zeta_{0}$, $\zeta_{1}$, $\cdots$ of the steps between consecutive
plateaux could be determined exactly using GHD. In this domain-wall problem it
has also been possible to go beyond the simplest hydrodynamic description by
investigating some fluctuations effects. Due to the fact that the velocity of
a given soliton is affected by the fluctuations of the densities of the other
species of soliton, these velocities fluctuate and the step between two
consecutive plateaux is not perfectly sharp but broadened. This broadening has
a diffusive scaling and the width of the step number $k$ behaves as $\delta
x\simeq t^{1/2}\Sigma_{k}$ with constants $\Sigma_{k}$ that have been
calculated analytically. All these results for the domain wall problem have
been also checked accurately using numerical simulations.

In the present work we are interested in the probability distribution of the
number $N_{t}$ of balls passing through the origin during a time $t$. We may
take an initial domain wall state where the left half is an i.i.d. state with
ball density $p_{\rm left}>0$, and the right half initially empty ($p_{\rm
right}=0$). In such a case the number of transferred balls is also equal the
number of balls in the right half at time $t$. But the ball propagation takes
place only in the right direction and there is no “blocking effect” from
existing balls in the BBS dynamics. So, $N_{t}$ does in fact not depend on the
initial state in the right half of the system. The probability distribution of
$N_{t}$ is therefore the same i) for the domain wall setup with ($p_{\rm
left}=p$ and $p_{\rm right}=0$) and ii) in a uniform i.i.d. state with ball
density $p_{\rm left}=p_{\rm right}=p$.

In the first step (Sec. 2.2) we are interested in the second cumulant $c_{2}$
of $N_{t}$. The result is obtained by two different methods, by a direct
calculation and with the TBA. Next (Sec. 2.3) we compute the Drude weight,
which is defined here in terms of the long-time limit of a spatially
integrated current-current correlation. This correlation is formally similar
to the one defining the second cumulant of $N_{t}$ and the Drude weight is
obtained using a TBA calculation which parallels that of the second cumulant.

A particularity of the model is to have a commuting family of temporal
evolutions which are in duality with the conserved quantities.222Although the
commuting time evolutions should in principle be a common feature of
integrable systems, their implementation can be quite involved in general.
This is one reason why the BBS is a precious model, since it allows concrete
descriptions of the different time evolution, as well as their simulations. In
Sec. 2.4 the currents associated with the conserved charges under the
different temporal evolutions are introduced. These depend on one index
associated to a conserved charge, and another one associated to a time
evolution. The generalized current operators turn out to be symmetric under
the exchange of the two indices, highlighting this duality. We compute the
long-time limit of the correlations associated to these generalized currents.
They generalize $c_{2}$ and $D$, and are shown to enjoy nontrivial symmetry
relations among them. In Sec. 2.5 we report on numerical calculations of the
Drude weights and the generalized current correlations using microscopic
simulations of the BBS, and a good agreement is found with the TBA results. In
Sec. 2.6 we compute the so-called flux Jacobian, which is a matrix playing an
important role in generalized hydrodynamics. The Sec. 2.7 discusses another
quantity defined in terms of spatially integrated current-current correlation,
the variance of the total current. Contrary to the Drude weight the
correlations are now taken at equal time, and this variance is obtained using
a transfer matrix approach.

In the second part of the paper (Sec. 3) we compute the scaled cumulants
generating function (SCGF) associated to $N_{t}$. It characterizes all the
cumulants of the probability distribution [22]. The Legendre transform of the
SCGF, called large deviation rate function, describes not only the typical
fluctuations around its mean value, but also the rare events. The large
deviation rate function is compared with the results of numerical simulations.
Finally, several results (second cumulant, Drude weight and SCGF) are
generalized to a more complex GGE, with two temperatures (C), one coupled to
the total number of balls and the other one coupled to the total number of
solitons.

##  2 Current correlations and Drude weights

We consider an i.i.d. homogeneous stationary state (characterized by some ball
density $p=z/(z+1)<1/2$) and we study the probability distribution of the
number $N_{t}$ of balls crossing the origin between time 0 and time $t$, in
the long time limit. The capacity of the carrier which induces the dynamics is
denoted by $l$.

###  2.1 Mean current

The mean value of $N_{t}$ grows linearly with time at a rate given by the mean
ball current:

$\left<N_{t}\right>\sim t\;j^{(l)}.$ (2.1)

The mean ball current $j^{(l)}$ is a simple function of the soliton currents
$j_{k}^{(l)}$

$j^{(l)}(z)=\sum_{k=1}^{\infty}k\;j_{k}^{(l)}.$ (2.2)

The factor $k$ above reflects the fact that each type-$k$ solitons carries $k$
balls. The mean soliton currents are products of the soliton velocities
$v_{k}^{(l)}$ by the soliton densities $\rho_{k}$:

$j_{k}^{(l)}=\rho_{k}v_{k}^{(l)}.$ (2.3)

In integrable systems such densities can be obtained using the TBA [23]. This
approach, applied to the BBS [5, 12], leads to a simple expression for the
mean soliton densities in terms of the ball fugacity $z=p/(1-p)$:

$\rho_{k}=\frac{z^{k}(1-z)^{3}(1+z^{k+1})}{(1+z)(1-z^{k})(1-z^{k+1})(1-z^{k+2})}.$
(2.4)

As for the velocities, they obey a set of equations which reflect the fact
that the mean velocity of each soliton species is affected by the collisions
with the solitons of the other species. The essential ideas were proposed in
[24, 25, 26], to describe soliton gases in the context of the Korteweg-de
Vries and nonlinear Schrödinger equations. In the case of BBS the equation for
the effective velocities reads [27] (see also [10, 12])

$v^{(l)}_{i}=\kappa^{(l)}_{i}+\sum_{k=1}^{\infty}M_{i,k}(v^{(l)}_{i}-v^{(l)}_{k})\rho_{k}\qquad(\forall
i\geq 1),$ (2.5)

where the bare soliton velocities are

$\kappa^{(l)}_{i}=\min(i,l)$ (2.6)

and the matrix $M$ (with indices from $[1,\infty]\times[1,\infty]$) encoding
the shifts experienced by solitons during collisions is

$M_{i,k}=2\min(i,k).$ (2.7)

For i.i.d. states these equations could be solved explicitly [12] and the
following expression for the velocities were obtained:

$v^{(l)}_{k}=\frac{1+z^{l+1}}{1-z^{l+1}}v_{\min(k,l)},\qquad
v_{k}=\frac{1+z}{1-z}k-\frac{2z(1+z)(1-z^{k})}{(1-z)^{2}(1+z^{k+1})}.$ (2.8)

Remark: $v_{k}^{(\infty)}=v_{k}$. Combining the results above the mean ball
current finally reduces to

$j^{(l)}(z)=\frac{z}{1-z}-\left(l+1\right){\frac{{z}^{l+1}}{1-{z}^{l+1}}}.$
(2.9)

###  2.2 Second cumulant

Next we are interested in the second cumulant of $N_{t}$, which turns out to
be linear in $t$:

$\left<\delta
N_{t}^{2}\right>=\left<N_{t}^{2}\right>-\left<N_{t}\right>^{2}\sim t\,c_{2}.$
(2.10)

This cumulant can be expressed as a time-integrated current-current
correlation:

$\displaystyle N_{t}$ $\displaystyle=\int_{0}^{t}j(0,s)ds,$ (2.11)
$\displaystyle\left<\delta N_{t}^{2}\right>$
$\displaystyle=\int_{0}^{t}\int_{0}^{t}dsds^{\prime}\left<j(0,s)j(0,s^{\prime})\right>^{c}\sim
t\,c_{2},$ (2.12)

where $j(x,s)$ is the current flowing from the site $x-1$ to the site $x$ at
time $s$. The superscript $(l)$ specifying the dynamics $T_{l}$ will often be
omitted in what follows. A continuous time notation is used here for clarity,
but the BBS is a discrete time model and $\int_{0}^{t}ds$ is in fact
equivalent to $\sum_{s=0}^{t}$. In the carrier picture for BBS, $j(x,s)$ is
nothing but the load of the carrier on this link. A general method to obtain
such cumulants in integrable systems is described in [28, 29]. It applies when
all the cumulants scale linearly in time. In what follows we apply this method
to the case of the BBS, where it simplifies considerably.

####  2.2.1 Correlations and sum rules.

We start from the general sum rule (A.12) derived in A and specialize it to
$f(x)=|x|$:

$\displaystyle\sum_{x=-\infty}^{\infty}|x|\left<\left(n(x,t)-n(x,0)\right)\left(n(0,t)-n(0,0)\right)\right>^{c}$
$\displaystyle=-2\int_{0}^{t}\int_{0}^{t}dsds^{\prime}\left<j(0,s)j(0,s^{\prime})\right>^{c}.$
(2.13)

Using the translation invariance (in space and time) of the correlators
appearing in the l.h.s we get

$\displaystyle\left<\delta N_{t}^{2}\right>=$
$\displaystyle\frac{1}{2}\sum_{x=-\infty}^{\infty}|x|\left(\left<n(x,t)n(0,0)\right>^{c}+\left<n(-x,t)n(0,0)\right>^{c}-2\left<n(x,0)n(0,0)\right>^{c}\right)$
$\displaystyle=\sum_{x=-\infty}^{\infty}|x|\left(\left<n(x,t)n(0,0)\right>^{c}-\left<n(x,0)n(0,0)\right>^{c}\right).$
(2.14)

The equation above is essentially equivalent to (2.23) of [30] (see also
(A.33) of [29]). The second cumulant is then expressed as

$c_{2}=\lim_{t\to\infty}\frac{1}{t}\sum_{x=-\infty}^{\infty}|x|\left(\left<n(x,t)n(0,0)\right>^{c}-\left<n(x,0)n(0,0)\right>^{c}\right).$
(2.15)

The connected equal time correlation $\left<n(x,0)n(0,0)\right>^{c}$ vanishes
for $x\neq 0$ in the i.i.d. state, so that

$c_{2}=\lim_{t\to\infty}\frac{1}{t}\sum_{x=-\infty}^{\infty}|x|\langle
n(x,t)n(0,0)\rangle^{c}.$ (2.16)

Since the propagation only takes place in the right direction in the BBS,
causality implies that the connected correlation $\langle
n(x,t)n(0,0)\rangle^{c}$ vanishes if $x<0$ (for $t>0$). So we have:

$\sum_{x=-\infty}^{\infty}|x|\langle
n(x,t)n(0,0)\rangle^{c}=\sum_{x=-\infty}^{\infty}x\langle
n(x,t)n(0,0)\rangle^{c}.$ (2.17)

We differentiate with respect to time, and use the local conservation law
$\frac{d}{dt}n(x,t)=j(x,t)-j(x+1,t)$:

$\displaystyle\frac{d}{dt}\sum_{x=-\infty}^{\infty}|x|\langle
n(x,t)n(0,0)\rangle^{c}$
$\displaystyle=\sum_{x=-\infty}^{\infty}x\langle\left(j(x,t)-j(x+1,t)\right)n(0,0)\rangle^{c}$
(2.18) $\displaystyle=\sum_{x=-\infty}^{\infty}\langle
j(x,t)n(0,0)\rangle^{c}=\sum_{x=-\infty}^{\infty}\langle
j(0,0)n(x,-t)\rangle^{c}.$ (2.19)

We recognize the conserved charge $Q=\sum_{x=-\infty}^{\infty}n(x,-t)$ in the
r.h.s., which means that (2.18)-(2.19) is independent of time. It can be
evaluated at $t=0$, or, instead, averaged over time. This yields two
equivalent formulations of $c_{2}$ in terms of integrated current-density
correlations:

$\displaystyle c_{2}$ $\displaystyle=\sum_{x=-\infty}^{\infty}\langle
j(x,0)n(0,0)\rangle^{c}$ (2.20)
$\displaystyle=\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}ds\sum_{x=-\infty}^{\infty}\langle
j(x,s)n(0,0)\rangle^{c}.$ (2.21)

####  2.2.2 Direct calculation.

In a GGE where $\beta$ is the inverse temperature associated to the total
number of balls (and $z=\exp(-\beta)$), (2.20) can be written as

$c_{2}=-\frac{\partial j}{\partial\beta}=z\frac{\partial j}{\partial z}.$
(2.22)

Using (2.9) we directly obtain

$c_{2}=\frac{z}{(1-z)^{2}}-(l+1)^{2}\frac{z^{l+1}}{(1-z^{l+1})^{2}}.$ (2.23)

The above formula has been compared with numerical simulation of the BBS. The
mean value and the second cumulant turn out to be in good agreement with the
theoretical values (Tab. 1). A generalization of this result to a two-
temperature GGE is given in (C.13).

$p$ | $l$ | $t$ |  $c_{2}$ (2.23) |  $\left<\delta N_{t}^{2}\right>/t$ numerics
---|---|---|---|---
0.2 | 4 | 6000 | 0.41998 | 0.419
0.4 | 4 | 2000 | 1.63352 | 1.63
0.3 | 10 | 2000 | 1.30166 | 1.29
Table 1: Variance of the number of transferred balls for different values of
the density $p$ and capacity $l$. Comparison between (2.23) and numerical
results (last column).

####  2.2.3 Second cumulant using TBA.

We will now compute $c_{2}$ by a different approach. The total current $J$ and
the total number of balls $Q$ can be decomposed into a mean value plus a
fluctuating part: $J=Lj+\delta J$, $Q=Lp+\delta Q$. With these definitions
(2.21) gives

$c_{2}=L^{-1}\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}ds\left<\delta J(s)\delta
Q(0)\right>.$ (2.24)

The pseudoenergy
$\epsilon_{i}=-\ln\left(\frac{2E_{i}-E_{i+1}-E_{i-1}}{L-2E_{i}}\right)$ can be
defined for each microscopic configuration, from the energies $E_{i}$,
$E_{i-1}$ and $E_{i+1}$. See D or [eq. (2.7), [12]] for the explanation of
$E_{i}$.111 The pseudoenergies can be used to write down the (fermionic) free
energy $F=-\sum_{k}\ln(1+e^{-\epsilon_{k}})$ [(3.15) in [12]] and the mode
occupancies are
$n_{k}=dF/d(\epsilon_{k})=(1+e^{\epsilon_{k}})^{-1}=(1+1/y_{k})^{-1}$ with
$y_{k}=e^{-\epsilon_{k}}$.  These energies are conserved under the time
evolution, and so are the pseudoenergies. The pseudoenergies can however
fluctuate from configuration to configuration in a GGE. If we define
$\delta\epsilon_{i}=L^{1/2}(\epsilon_{i}-\bar{\epsilon}_{i})$, with
$\bar{\epsilon}_{i}=\left<\epsilon_{i}\right>$, these fluctuations have
diagonal correlations [12]

$\left<\delta\epsilon_{i}\delta\epsilon_{j}\right>=\delta_{ij}(1+e^{\bar{\epsilon}_{i}})/\sigma_{i},$
(2.25)

where

$\sigma_{i}=1-\sum_{k}M_{ik}\rho_{k}$ (2.26)

is the hole density and $M$ was given in (2.7).222The values of $\sigma_{i}$
in the i.i.d. state are given by (C.4) with $a=z$ (or (3.25) in [12]). Note
also that the expectation values of the pseudoenergies are given by
$\bar{\epsilon}_{i}=-\ln\left(\rho_{i}/\sigma_{i}\right)$. The relation (2.25)
was checked numerically (see Tab. 2).

i | j | $\delta_{ij}(1+e^{\bar{\epsilon}_{i}})/\sigma_{i}$ |  $\left<\delta\epsilon_{i}\delta\epsilon_{j}\right>$ numerics
---|---|---|---
1 | 1 | 8.01282 | 8.0117
1 | 2 | 0 | 0.00027
1 | 3 | 0 | -0.00079
2 | 2 | 27.2183 | 27.223
2 | 3 | 0 | -0.0053
3 | 3 | 65.5367 | 65.588
4 | 4 | 131.789 | 132.07
5 | 5 | 238.454 | 239.37
Table 2:  Pseudoenergy correlations for ball density $p=0.4$: comparison
between the TBA result (r.h.s of (2.25)), and numerical simulations (system
size $L=8\,10^{4}$ and average over $N_{\rm samples}=10^{8}$ configurations).
The numerical results correspond to
$L\left<\ln\left(\frac{2E_{i}-E_{i+1}-E_{i-1}}{L-2E_{i}}\right)\ln\left(\frac{2E_{j}-E_{j+1}-E_{j-1}}{L-2E_{j}}\right)\right>^{c}$.
As for the expectation values $\langle\epsilon_{i}\rangle$, the simulation
results agree with $-\ln\left(\rho_{i}/\sigma_{i}\right)$ with very good
accuracy (data not shown).

In turn, the fluctuations of the other quantities can be related to the
fluctuations $\delta\epsilon_{j}$. In a large system $\delta J/L$ and $\delta
Q/L$ are typically small, of order $\mathcal{O}(L^{-1/2})$, and we may
linearize the relation between the ball density fluctuation and the
pseudoenergies,

$\delta
Q/\sqrt{L}=\sum_{i=1}^{\infty}\delta\epsilon_{i}\frac{\partial\rho_{b}}{\partial\bar{\epsilon}_{i}},$
(2.27)

as well as the relation between the ball current density fluctuation and the
pseudo energies

$\delta J/\sqrt{L}=\sum_{i=1}^{\infty}\delta\epsilon_{i}\frac{\partial
j_{b}}{\partial\bar{\epsilon}_{i}}+\cdots.$ (2.28)

In the equations above the ball density was denoted by $\rho_{b}$ and the ball
current by $j_{b}$. In (2.28) we have decomposed a current fluctuation, which
is not conserved in time, in terms of fluctuations $\delta\epsilon_{i}$ of the
pseudoenergies. Since each $\epsilon_{i}$ is a function of the $\\{E_{k}\\}$,
the pseudoenergies and their fluctuations $\delta\epsilon_{i}$ are
configuration-dependent but independent of time. So, by decomposing $\delta J$
over the $\delta\epsilon_{i}$ we have dropped all the time dependence in
$\delta J$ and only the conserved component of the current has been kept
(hence the dots in (2.28)). Following the idea of hydrodynamics projection
onto the space of conserved quantities [15] we are here restricting ourselves
to the part of the current which is constant in time and can be expressed in
terms of the conserved energies. It is of course legitimate to do so in order
to compute $c_{2}$, thanks to the long-time limit in (2.24). Note that, in
contrast, nothing was dropped in (2.27) since the total number of balls (and
thus also $\delta Q$) is a conserved quantity.

Replacing (2.27), (2.28) and (2.25) in (2.24) yields

$c_{2}=\sum_{i=1}^{\infty}\frac{\partial\rho_{b}}{\partial\bar{\epsilon}_{i}}\,\frac{\partial
j_{b}}{\partial\bar{\epsilon}_{i}}\,\frac{1+e^{\bar{\epsilon}_{i}}}{\sigma_{i}}.$
(2.29)

What remains to be done is to compute the derivatives
$\frac{\partial\rho_{b}}{\partial\bar{\epsilon}_{i}}$ and $\frac{\partial
j_{b}}{\partial\bar{\epsilon}_{i}}$.

Once the pseudoenergies associated to the GGE are known we can introduce the
“dressing” operation, which is standard in the framework of TBA. A set of
quantities $o_{i}$ labelled by some index $i\in[1,\infty]$ representing a
soliton size can be grouped into a column vector $o=(o_{1},o_{2},\cdots)$. We
then define the “dressing matrix” $G$ and construct a new dressed vector
$o^{\rm dr}$ [12]:

$\displaystyle o^{\rm dr}$ $\displaystyle=Go,$ (2.30) $\displaystyle G$
$\displaystyle=\left(1+M\hat{y}\right)^{-1},$ (2.31)

where $\hat{y}$ is the diagonal matrix with elements
$y_{j}=\exp(-\bar{\epsilon}_{j})$.333The diagonal entries of the matrix $M$ do
not enter the velocity equation (2.5) and there is therefore some freedom to
redefine its diagonal part. This leads to some freedom in the definition of
the dressing operation and the present choice (also used in [12]) is not the
same as in [15]. With (2.31) we have $o^{\rm dr}=o-M{\hat{y}}o^{\rm dr}$. In
[15] the dressing operation (denoted with a prime) is instead defined by
$o^{\rm dr^{\prime}}=o+T{\hat{n}}o^{\rm dr^{\prime}}$ where $T=1-M$ and
$\hat{n}$ is the diagonal matrix with entries
$n_{i}=\rho_{i}/(\sigma_{i}+\rho_{i})$. The two dressing definitions lead to
the same effective soliton speeds $v_{i}^{(l)}=\kappa_{i}^{(l){\rm
dr}}/(1^{\rm dr})_{i}=\kappa_{i}^{(l){\rm dr^{\prime}}}/(1^{\rm
dr^{\prime}})_{i}$ and can both be used to obtain physical quantities.

The basic relations that will be used in the following calculation are

$\displaystyle
y_{i}=e^{-\bar{\epsilon}_{i}}=\frac{\rho_{i}}{\sigma_{i}},\qquad\sigma_{i}v^{(l)}_{i}=(G\kappa^{(l)})_{i}\;\bigl{(}=(\kappa^{(l)\mathrm{dr}})_{i}\bigr{)}.$
(2.32)

From $G+GM\hat{y}=\mathrm{Id}$, we find that $GM$ is symmetric:

$\displaystyle
GM=\hat{y}^{-1}-G\hat{y}^{-1}=\hat{y}^{-1}-(\hat{y}+\hat{y}M\hat{y})^{-1}=(GM)^{t}.$
(2.33)

Its $(i,j)$ element is concretely expressed as

$\displaystyle(GM)_{ij}$
$\displaystyle=2\sum_{k}G_{ik}\min(k,j)=2(G\kappa^{(j)})_{i}=2\sigma_{i}v^{(j)}_{i},$
(2.34)

where the invariance of the last expression under the interchange
$i\leftrightarrow j$ can also be confirmed at the level of explicit formulae,
see (2.72). Using $\frac{\partial
y_{n}}{\partial\bar{\epsilon}_{i}}=-\delta_{in}y_{i}$, we get

$\displaystyle\frac{\partial G_{jk}}{\partial\bar{\epsilon}_{i}}$
$\displaystyle=-\sum_{m,n}G_{jm}\frac{\partial(M\hat{y})_{mn}}{\partial\bar{\epsilon}_{i}}G_{nk}=\sum_{m}G_{jm}M_{mi}y_{i}G_{ik}$
$\displaystyle=(GM)_{ji}y_{i}G_{ik}=(GM)_{ij}y_{i}G_{ik}=2\sigma_{i}v^{(j)}_{i}y_{i}G_{ik}=2\rho_{i}v^{(j)}_{i}G_{ik}.$
(2.35)

Let us introduce $\eta^{(l)}_{j}$, which includes, due to $v^{(l=1)}_{k}=1$,
the density and the current of balls as special cases:

$\displaystyle\eta^{(l)}_{j}=\sum_{k}\min(j,k)\rho_{k}v^{(l)}_{k},\qquad\eta^{(1)}_{\infty}=\rho_{b},\qquad\eta^{(l)}_{\infty}=j_{b}.$
(2.36)

It is expressed in terms of $G$ as

$\displaystyle\eta^{(l)}_{j}$
$\displaystyle=\sum_{k}\min(j,k)y_{k}\sigma_{k}v^{(l)}_{k}=\frac{1}{2}\sum_{k}(M\hat{y})_{jk}(G\kappa^{(l)})_{k}$
$\displaystyle=\frac{1}{2}\sum_{k}(-\mathrm{Id}+G^{-1})_{jk}(G\kappa^{(l)})_{k}=-\frac{1}{2}\sum_{k}G_{jk}\kappa^{(l)}_{k}+\frac{1}{2}\kappa^{(l)}_{j}=\eta^{(j)}_{l},$
(2.37)

where the last equality is due to the symmetry (2.33). The symmetry of $\eta$
under the exchange of the upper and lower indices will play an important role
in Sec. 2.4. By means of (2.2.3), the $\bar{\epsilon}_{i}$ derivative is
calculated as

$\displaystyle\frac{\partial\eta^{(l)}_{j}}{\partial\bar{\epsilon}_{i}}=-\frac{1}{2}\sum_{k}\frac{\partial
G_{jk}}{\partial\bar{\epsilon}_{i}}\kappa^{(l)}_{k}=-\sum_{k}\rho_{i}v^{(j)}_{i}G_{ik}\kappa^{(l)}_{k}=-\rho_{i}\sigma_{i}v^{(j)}_{i}v^{(l)}_{i}.$
(2.38)

From this with $j=\infty$ and (2.36) it follows that

$\displaystyle\frac{\partial\rho_{b}}{\partial\bar{\epsilon}_{i}}$
$\displaystyle=-\rho_{i}\sigma_{i}v^{(\infty)}_{i},$ (2.39)
$\displaystyle\frac{\partial j_{b}}{\partial\bar{\epsilon}_{i}}$
$\displaystyle=-\rho_{i}\sigma_{i}v^{(\infty)}_{i}v^{(l)}_{i}.$ (2.40)

The results (2.39) and (2.40) can be inserted into (2.29):

$c_{2}=\sum_{i\geq
1}\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})(v^{(\infty)}_{i})^{2}v^{(l)}_{i}.$
(2.41)

This expression can be checked to agree with (2.23).

###  2.3 Drude weights

The Drude weight is an important quantity in the field of transport. This
linear response coefficient characterizes the increase of the current in
presence of a force which couples to the charges. The recent years have seen
very important progress in the understanding of this quantity in one-
dimensional integrable systems [31, 32, 33, 34, 35, 36]. $D$ can be defined in
terms of a time average current-current correlation function :444Using the sum
rule (A.12) with the choice $f(x)=x^{2}$, (2.43) can also be written
$D=\lim_{t\to\infty}\frac{1}{t^{2}}\sum_{x}x^{2}\left<n(x,t)n(0,0)\right>^{c},$
(2.42) which is another useful definition of $D$.

$D=\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}ds\sum_{x}\left<j(x,s)j(0,0)\right>^{c}.$
(2.43)

Noting the close similarity with (2.21), we repeat the hydrodynamic projection
of the current fluctuations onto the pseudoenergy fluctuations, which lead to
(2.29). It yields to

$D=\sum_{i=1}^{\infty}\left(\frac{\partial
j_{b}}{\partial\bar{\epsilon}_{i}}\right)^{2}\frac{1+e^{\bar{\epsilon}_{i}}}{\sigma_{i}}$
(2.44)

and the derivatives $\frac{\partial j_{b}}{\partial\bar{\epsilon}_{i}}$ are
again given by (2.40). We thus obtain a TBA expression for the Drude weight

$D=\sum_{i\geq
1}\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})(v^{(\infty)}_{i}v^{(l)}_{i})^{2}.$
(2.45)

The expression for the density correlation
$R=\sum_{x}\left<n(x,t)n(0,0)\right>^{c}$ has a form similar to that of
$c_{2}$ and $D$. Thanks to the conservation of the total number of balls $R$
is independent of time and in an i.i.d. state it is trivially equal to
$p(1-p)$. Following the reasoning leading to (2.41) and (2.45) we get

$p(1-p)=\sum_{i\geq
1}\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})(v^{(\infty)}_{i})^{2}.$ (2.46)

For more details about the above formula and their expression in matrix forms
we refer the readers to Sec. 2.6.3. The structure of these formulae is very
similar to the results for $c_{2}$ and $D$ in the Lieb-Liniger model, obtained
by Doyon and Spohn (eqs. (1.2) and (1.3) of [33]). In the latter work a
central idea is to implement the long-time limit in the definition of $D$ via
the hydrodynamic projection (see (2.18) of [33]). This approach amounts to
projecting the observable of interest (for $D$, the current) via a suitable
scalar product into the time-invariant subspace, spanned by the conserved
charges. As mentioned in the paragraph before (2.29), in the approach
presented here the projection is implemented by decomposing the current
fluctuation over the conserved variables $\epsilon_{i}$. The relation (2.25)
shows that the pseudoenergies form a convenient orthogonal basis in the
subspace of conserved quantities.

Aided by some computer algebra system it was possible to obtain explicit
formulae for the Drude weight as a function of $z$ for a few values of the
carrier capacity $l$, denoted by $D^{(l)}$ :

$\displaystyle D^{(2)}$
$\displaystyle=\frac{z(1-z)(1-z^{2})^{3}(1+11z+11z^{3}+z^{4})}{(1-z^{3})^{3}(1-z^{4})},$
(2.47) $\displaystyle D^{(3)}$
$\displaystyle=\frac{z(1-z)(1-z^{2})^{2}(1-z^{3})}{(1-z^{4})^{3}(1-z^{6})}(1+11z+44z^{2}+29z^{3}+30z^{4}+29z^{5}+44z^{6}+11z^{7}+z^{8}),$
(2.48) $\displaystyle D^{(4)}$
$\displaystyle=\frac{z(1-z^{2})^{4}(1-z^{3})}{(1-z^{5})^{3}(1-z^{6})(1-z^{8})}(1+10z+35z^{2}+117z^{3}+68z^{4}+254z^{5}+95z^{6}+357z^{7}$
$\displaystyle\qquad\qquad+126z^{8}+357z^{9}+95z^{10}+254z^{11}+68z^{12}+117z^{13}+35z^{14}+10z^{15}+z^{16}),$
(2.49) $\displaystyle D^{(5)}$
$\displaystyle=\frac{z(1-z)(1-z^{2})^{2}(1-z^{4})(1-z^{5})}{(1-z^{6})^{3}(1-z^{8})(1-z^{10})}(1+11z+44z^{2}+140z^{3}+355z^{4}+406z^{5}+480z^{6}$
$\displaystyle+443z^{7}+633z^{8}+714z^{9}+896z^{10}+714z^{11}+633z^{12}+443z^{13}+480z^{14}+406z^{15}$
$\displaystyle+355z^{16}+140z^{17}+44z^{18}+11z^{19}+z^{20}).$ (2.50)

The functions above as well as $D^{(\infty)}$ are plotted as a function of the
ball density $p=z/(1+z)$ in Fig. 2. When the density is slightly below 0.5 we
note a rapid increase of the Drude weight with $l$. Since $l$ acts as a cutoff
on the effective speed of the solitons of size $k\geq l$, this indicates that
the large solitons have a dominant contribution to the Drude weight.

Figure 2: Drude weight $D^{(l)}$ for $l=2,3,4,5$ and $l=\infty$, plotted as a
function of the ball density $p$. For $l\geq 2$ the Drude weight is linear in
$p$ at low density: $D^{(l)}=p+11p^{2}+\mathcal{O}(p^{3})$.

Going further one can obtain an expression for $D^{(l)}$ for generic $l$ with
finite number of terms. To this end we first set
$W_{i}=\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})v^{2}_{i}$. Then we have

$\begin{split}\sum_{i<l}W_{i}&=\frac{zA_{l}}{(1+z)^{2}(1-z^{l})^{2}(1-z^{l+1})^{2}},\\\
A_{l}&=1+3z^{2l}+8z^{2l+1}+3z^{2l+2}+z^{4l+2}-(l+1)^{2}(z^{l+2}+z^{3l})\\\
&+(l^{2}-4)(z^{l+1}+z^{3l+1})+(l-1)(l+3)(z^{l}+z^{3l+2})-l^{2}(z^{l-1}+z^{3l+3}).\end{split}$
(2.51)

Note that $A_{0}=A_{1}=0$ indeed holds. It leads to

$\displaystyle D^{(1)}=\sum_{i\geq 1}W_{i}(v^{(1)}_{i})^{2}=\sum_{i\geq
1}W_{i}=\frac{zA_{\infty}}{(1+z)^{2}}=\frac{z}{(1+z)^{2}}.$ (2.52)

By means of this, the infinite sum (2.45) is reduced to a finite one as

$\displaystyle D^{(l)}$
$\displaystyle=\sum_{i<l}W_{i}(v^{(l)}_{i})^{2}+(v^{(l)}_{l})^{2}\sum_{i\geq
l}W_{i}$ (2.53)
$\displaystyle=\sum_{i<l}W_{i}(v^{(l)}_{i})^{2}+(v^{(l)}_{l})^{2}\Bigl{(}\frac{z}{(1+z)^{2}}-\sum_{i<l}W_{i}\Bigr{)},$
(2.54)

where the second term is known by (2.51). Or equivalently,

$\displaystyle
D^{(l)}=\left(\frac{1+z^{l+1}}{1-z^{l+1}}\right)^{2}\left(\sum_{i<l}W_{i}(v_{i}^{2}-v_{l}^{2})+v_{l}^{2}\frac{z}{(1+z)^{2}}\right).$
(2.55)

We note that in the half filled limit the Drude weight tends to a simple value

$\displaystyle\lim_{z\rightarrow 1}D^{(l)}=\left(\frac{l(l+2)}{6}\right)^{2}.$
(2.56)

This is consistent with the observation that $D^{(\infty)}$ diverges when
$z\to 1$ (see Fig. 2). To conclude this section we mention that some
generalization of these results to a two-temperature GGE is discussed in C.

###  2.4 Generalized current correlations and their symmetries

Thus far we have considered $c_{2}$ (2.24), $D$ (2.43) and $R$ (2.46). They
are all associated with the number balls, which is one special case
$E_{\infty}$ of the conserved quantities $E_{j}$ with $j=1,2,\ldots$ [12].
Here we discuss a natural generalization corresponding to the whole family
$\\{E_{j}\\}$.

Let us denote by $\hat{\eta}^{(l)}_{j}(x)$ the microscopic operator measuring
the current associated to the $j^{\rm th}$ energy, under the time evolution
$T_{l}$, and at position $x$. To define this operator microscopically in a
given state $s$ of the BBS we need to consider two successive time evolutions
$s\rightarrow T_{l}(s)\rightarrow T_{i}T_{l}(s)$, as illustrated in Fig. 3.
Consider the carrier inducing the first $T_{l}$ and let $u(x)=(u_{0},u_{1})$
be its state at position $x$, where $u_{0}$ and $u_{1}$ are the numbers of
empty space and balls in it, respectively. Similarly, let
$u^{\prime}(x)=(u^{\prime}_{0},u^{\prime}_{1})$ be the state of another
carrier for the second time evolution $T_{i}$ at position $x$. By the
definition of the carrier capacity $u_{0}+u_{1}=l$ and
$u^{\prime}_{0}+u^{\prime}_{1}=i$. Now ${\hat{\eta}}^{(l)}_{i}(x)$ of the
state $s$ is defined by ${\hat{\eta}}^{(l)}_{i}(x)=\min(u^{\prime}_{0},u_{1})$
(red integers in Fig. 3). By extending the argument in Sec. 2.2 in [12], it
can be shown that these generalized currents are in fact symmetric in the two
capacities $l$ and $i$: $\hat{\eta}^{(l)}_{i}(x)=\hat{\eta}^{(i)}_{l}(x)$, and
$\hat{\eta}^{(l)}_{\infty}(x)=j^{(l)}(x)$ and ${\hat{\eta}}^{(1)}_{k}(x)=$
local term for the energy $E_{k}$. We refer the reader to D for more details
about the construction of these currents.

00011100110$T_{3}$21001232123211000011001$\hat{\eta}^{(3)}_{2}(x)$100012210121$T_{2}$12210001210101110000110

00011100110$T_{2}$10001221012110000111001$\hat{\eta}^{(2)}_{3}(x)$100012210121$T_{3}$23210012321201110000110

Figure 3: Top: The generalized current $\eta^{(3)}_{2}(x)$ is shown in red
letters. According to (D.7) and (D.3), it is evaluated as
$\min(2-\alpha_{1},\beta_{1})$, where $\alpha_{1}$ and $\beta_{1}$ are numbers
just under and above it, respectively. Bottom: Similarly, the generalized
current $\eta^{(2)}_{3}(x)$ is defined using $\min(3-\alpha_{1},\beta_{1})$.
One can see on this example that the two currents densities are in fact equal:
$\eta^{(3)}_{2}(x)=\eta^{(2)}_{3}(x)$.

The mean values of these energy currents are expressed in terms of
$\eta^{(l)}_{j}$ (2.36) as

$\displaystyle\eta^{(1)}_{j}=\left<\hat{\eta}^{(1)}_{j}\right>=\sum_{k\geq
1}\min(j,k)\rho_{k},\quad\eta^{(l)}_{j}=\left<\hat{\eta}^{(l)}_{j}\right>=\sum_{k\geq
1}\min(j,k)\rho_{k}v^{(l)}_{k},$ (2.57)

where $v^{(l=1)}_{k}=1$ is used. The mean ball current (2.2) corresponds to
$\eta^{(l)}_{\infty}$, and the ball density is $\eta^{(1)}_{\infty}$.
Comparing (2.57) and (2.36), we see that the index $l$ allows interpolating
between the density and the current, utilizing the variety of time evolutions
in the system.

Now consider the time-averaged correlation

$\displaystyle{\mathcal{C}}^{m,l,n}_{i,j}=\lim_{t_{n}\rightarrow\infty}\frac{1}{t_{n}}\int_{0}^{t_{n}}ds\sum_{x}\left<\hat{\eta}^{(m)}_{i}(x,s)\hat{\eta}^{(l)}_{j}(0,0)\right>^{c},$
(2.58)

which is also expected to coincide with the limit

$\displaystyle{\mathcal{C}}^{m,l,n}_{i,j}=\lim_{t_{n}\rightarrow\infty}\sum_{x}\left<\hat{\eta}^{(m)}_{i}(x,t_{n})\hat{\eta}^{(l)}_{j}(0,0)\right>^{c}.$
(2.59)

In the definitions above we have denoted the time variable by $t_{n}$ to
emphasize that the time evolution from time 0 to time $t_{n}$ is computed with
$T_{n}$. We however conjecture that, for sufficiently large $n$ this quantity
becomes independent of $n$, and we set

$\displaystyle
C^{l,m}_{i,j}=\lim_{n\rightarrow\infty}{\mathcal{C}}^{m,l,n}_{i,j}.$ (2.60)

This property is easy to check if at least one of the indices $m,l,i,j$ is
equal to one. Using the symmetry between the upper and lower indices in
$\hat{\eta}^{(l)}_{i}$, the index which is equal to one can be moved to an
upper position. The correlation (2.59) can then be formulated using, say,
$\sum_{x}\hat{\eta}^{(1)}_{i}(x)$, which is the conserved energy $E_{i}$. The
time evolution therefore drops and ${\mathcal{C}}^{m,l,n}_{i,j}$ appears to be
independent of $n$ in such a case. As we will see in Sec. 2.5, in more general
cases the numerical simulations support the fact that
${\mathcal{C}}^{m,l,n}_{i,j}$ is independent of $n$ for $n\geq\min(m,l)$.

The family of generalized correlations $C^{l,m}_{i,j}$ includes three
quantities we have defined previously:

$\displaystyle c_{2}=C^{1,l}_{\infty,\infty},\quad
D=C^{l,l}_{\infty,\infty},\quad R=C^{1,1}_{\infty,\infty}.$ (2.61)

As seen here, the superscripts $l,m$ in (2.60) are restricted to $1$ (for
density) or $l$ (for current) in the original setting. However, keeping them
general reveals an interesting symmetry of the problem as we see below.
Admitting the validity of the previous TBA argument leads to

$\displaystyle C^{l,m}_{i,j}$
$\displaystyle=\sum_{k=1}^{\infty}\frac{\partial\eta^{(l)}_{i}}{\partial\bar{}\epsilon_{k}}\frac{\partial\eta^{(m)}_{j}}{\partial\bar{\epsilon}_{k}}\frac{1+e^{\bar{\epsilon}_{k}}}{\sigma_{k}}$
(2.62) $\displaystyle=\sum_{k\geq
1}\rho_{k}\sigma_{k}(\rho_{k}+\sigma_{k})v^{(i)}_{k}v^{(j)}_{k}v^{(l)}_{k}v^{(m)}_{k}.$
(2.63)

It tells that $C^{l,m}_{i,j}$ is completely symmetric in the four indices.

The symmetry relating $C^{l,m}_{i,j}$ to $C^{i,m}_{l,j}$ and $C^{l,j}_{i,m}$
derives from the equality $\hat{\eta}_{i}^{(l)}=\hat{\eta}_{l}^{(i)}$. For
example, the ball density correlation $R=C^{1,1}_{\infty,\infty}$ coincides
with $C^{\infty,\infty}_{1,1}$ which is the soliton current correlation under
the time evolution $T_{\infty}$. In this particular case the symmetry can be
checked directly at the microscopic level.555To this end one compares the
$1-$energy density (or soliton density) $e_{1}(x,t)$ at times $t$ and $t+1$
(according to the $T_{\infty}$ evolution) to deduce the associated current
$j_{1}(x,t)$ satisfying the discrete continuity equation
$e_{1}(x,t+1)-e_{1}(x,t)+j_{1}(x+1/2,t)-j_{1}(x-1/2,t)=0$. Doing so one
realizes that the soliton current is identical to the ball density, up to a
half lattice spacing shift: $j_{1}(x,t)=n(x-1/2,t)$. This explains the
equality between their correlations and provides an explicit check of the
general symmetry mentioned above.

The symmetry relating $C^{l,m}_{i,j}$ to $C^{m,l}_{j,i}$ follows from the fact
that, in the correlation
$\left<\hat{\eta}^{(m)}_{i}(x,s)\hat{\eta}^{(l)}_{j}(0,0)\right>^{c}$ one can
replace $(x,s)$ by $(-x,-s)$. The possibility to exchange, say, the indices
$l$ and $m$ in $C^{l,m}_{i,j}$ is apparent in the TBA result (2.63) but it
does not follow from a simple microscopic symmetry. It will be interesting to
seek such an enhanced symmetry among the transport characteristics also in
other integrable models admitting the GHD approach.

We note that the family of characteristics
$C^{m,l}_{\infty,\infty}\,(m=1,2,\ldots,l)$ interpolating $c_{2}$ and $D$ in
(2.61) are all expressible as a finite sum similar to (2.54):

$\displaystyle
C^{m,l}_{\infty,\infty}=\sum_{k<m}W_{k}v^{(l)}_{k}(v^{(m)}_{k}-v^{(m)}_{m})+c_{2}v^{(m)}_{m}.$
(2.64)

When $i,j,l,m$ are all finite, $C^{l,m}_{i,j}$ is evaluated as a rational
function by means of (C.15) with $a=z$ as

$\displaystyle
C^{l,m}_{i,j}=\sum_{k<r}\rho_{k}\sigma_{k}(\rho_{k}+\sigma_{k})v^{(i)}_{k}v^{(j)}_{k}v^{(l)}_{k}v^{(m)}_{k}+\frac{z^{r}(1-z)^{4}(1+z^{2r+1})}{(1+z)^{3}(1-z^{r})^{2}(1-z^{r+1})^{2}}v^{(i)}_{i}v^{(j)}_{j}v^{(l)}_{l}v^{(m)}_{m},$
(2.65)

where $r=\max(i,j,l,m)$.

###  2.5 Numerical calculation of the Drude weights

In this section we present some direct numerical calculations of the
correlation ${\mathcal{C}}^{m,l,n}_{i,j}$. The results are displayed in Figs.
4, 5 and 6.

Fig. 4 illustrates the dependence of ${\mathcal{C}}^{m,l,n}_{i,j}$ on the
dynamical parameter $n$ in the case $m=l$ and $i=j=99$ (almost $\infty$).
Since $T_{1}$ is a simple translation we have
$\hat{\eta}^{(l)}_{i}(x,t_{1})=\hat{\eta}^{(l)}_{j}(x-t_{1},0)$ and (2.59)
reduces to an equal-time correlation when $n=1$ (no long-time limit). This
equal-time current-current correlation is computed in Sec. 2.7 and the
numerical data for $n=1$ agree with the analytical result. For general values
of $n$ we can compare the numerical results for
${\mathcal{C}}^{l,l,n}_{\infty,\infty}$ with the Drude weight
${\mathcal{C}}^{l,l}_{\infty,\infty}$, and the data plotted in Fig. 4 suggest
that we have ${\mathcal{C}}^{l,l,n}_{\infty,\infty}=C^{l,l}_{\infty,\infty}$
for $n\geq l$. We are however unaware of a method to obtain
${\mathcal{C}}^{l,l,n}_{\infty,\infty}$ in the intermediate regime $1<n<l$.

Fig. 5 illustrates the dependence on the dynamical parameter $n$ in a case
where $m\neq l$ (and still $i=j=99$). For $n=1$ the correlation
${\mathcal{C}}^{m,l,1}_{\infty,\infty}$ is again an equal-time correlation
between two generalized currents. For larger values of the dynamical parameter
$n$ the data displayed in Fig. 5 suggest that we have
${\mathcal{C}}^{m,l,n}_{\infty,\infty}=C^{m,l}_{\infty,\infty}$ for
$n\geq\min(m,l)$, even though we have no proof.

Finally, Fig. 6 illustrates the convergence of current-current correlations as
a function of time and for different time evolutions. This figure shows that
it converges to ${\mathcal{C}}^{l,l,n}_{\infty,\infty}$ after only a
relatively short relaxation time. In this particular example ($l=5$ and ball
density $p=0.3$) the relaxation time appears to decrease when $n$ increases.
As for the long-time limit of the correlator, it coincides with the TBA
expression (2.63) only for $n\geq l=5$, in agreement with the conjectured
proposed above.

If the long-time limit of a current-current correlation (that is, the Drude
weight) is subtracted, the time integral of the remaining part of the
correlation is an element of the Onsager matrix, which characterizes the
diffusive corrections to the hydrodynamics [15, 36]. We see from Fig. 6 that
such Onsager coefficient would be dominated by a few short-time value of the
correlator. The detailed study of the Onsager matrix in the BBS is left for
future studies.666The diffusive broadening of density steps in a domain-wall
problem was studied analytically and numerically for the BBS in [12, 21].

In the hydrodynamic projection scheme the Drude weights and its
generalizations (2.59) depend on the projections of the current operators onto
the space spanned by the conserved quantities. The energies $E_{i}$ are
conserved whatever the dynamical parameter $n$ and their contribution to the
long-time limit of the current-current correlations was included in (2.63),
via the pseudoenergy variables. However, the numerical observation that
${\mathcal{C}}^{m,l,n}_{i,j}$ can be different from (2.63) suggest that
$T_{n}$ may not always be sufficiently mixing inside each energy sector fixed
by the $\\{E_{i}\\}$ data, in particular for small $n$. This is obvious for
$n=1$ (no mixing at all), but as a less trivial example consider two solitons
of sizes $p$ and $p^{\prime}$ under the evolution $T_{n}$, with $n\leq
p<p^{\prime}$. Because the effective soliton velocities (2.8) are independent
of the soliton size when they are larger than $n$ the mean distance between
two big solitons will be almost constant in time (up to small fluctuations due
to the collisions with smaller and slower solitons). This may be a hint why
some minimal value of $n$ seems required for a current-current correlation to
converge to the TBA result (2.63).

Figure 4: Numerical results for the Drude weight $\mathcal{C}^{l,l,n}_{99,99}$
(2.59) for different values of the ball densities ($p=0.4$ and $p=0.2$), two
values of $l$ and different values of the parameter $n$ characterizing the
time evolution. The simulations were carried out on a periodic system of size
$L=60000$, with time $t_{n}=1000$ and $N_{\rm samples}=5\,10^{5}$ random
initial configurations (top and middle panel) and with $N_{\rm
samples}=8\,10^{6}$ in the bottom panel. The red horizontal lines represent
the TBA result (2.45) and the green lines represent (B.7).  Figure 5:
Numerical results for ${\mathcal{C}}^{3,6,n}_{99,99}$ (2.59) plotted as a
function of $n$. The data suggests that
${\mathcal{C}}^{m,l,n}_{\infty,\infty}$ coincides with the TBA result (2.63)
for $C^{m,l}_{\infty,\infty}$ when $n\geq\min(m,l)$. Simulations performed
with (periodic) system size $L=30000$, time $t_{n}=1000$ and $1.5\,10^{6}$
random initial conditions with ball density $p=0.3$. Figure 6: Correlation
$\sum_{x}\langle\eta^{(5)}_{99}(x,t_{n})\eta^{(5)}_{99}(0,0)\rangle^{c}$
plotted as a function of time for different time evolutions $T_{n}$. Since
$T_{1}$ is a simple translation the dynamics is trivial and the correlation is
independent of time (purple curve, $n=1$). For $n\geq 2$ the dynamics is
however nontrivial, and the correlation relaxes to a finite value at
sufficiently long times. For $n\geq 5$ this long-time limit coincides with the
TBA value (2.45) (horizontal red line). Simulations performed with (periodic)
system size $L=20000$, and $36\,10^{6}$ random initial conditions with ball
density $p=0.3$.

###  2.6 Flux Jacobian and related matrices

In this section we present concrete forms of the flux Jacobians corresponding
to the commuting family of time evolutions in our BBS. We will also describe
the relations with the matrices of 2nd cumulants and the Drude weights, as
well as their generalization $C^{l,m}_{i,j}$.

####  2.6.1 Normal modes

In BBS, the role of the conserved charge $q_{i}$ and the current $j_{i}$ in
the flux Jacobian [eq. (26) of [15]] are played by $\eta^{(1)}_{i}$ and
$\eta^{(l)}_{i}$ if the time evolution is taken as $T_{l}$. From [eq. (3.9),
[12]], $\eta^{(1)}_{i}$ is related to the hole density as
$\sigma_{i}=1-2\eta^{(1)}_{i}$. The associated conserved currents for the time
evolution $T_{l}$ are given by $\sigma_{i}v^{(l)}_{i}$. Therefore, we have the
equations of conservation:

$\displaystyle\partial_{t_{l}}\sigma_{i}+\partial_{x}(\sigma_{i}v^{(l)}_{i})=0.$
(2.66)

We use the fact that ${\sigma}_{i}$ and $v^{(l)}_{i}$ are functions of $y_{j}$
given by $\sigma_{i}=\sum_{k=1}^{\infty}G_{ik}$ and
$\sigma_{i}v^{(l)}_{i}=(G\kappa^{(l)})_{i}=(GM)_{il}/2$. Their partial
derivatives with respect to $y_{j}$ are calculated as

$\displaystyle{\partial\sigma_{i}\over\partial y_{j}}$
$\displaystyle=-\sum_{k}y^{-1}_{j}\frac{\partial
G_{ik}}{\partial\bar{\epsilon}_{j}}=-\sum_{k}y^{-1}_{j}2\sigma_{j}v^{(i)}_{j}y_{j}G_{jk}=-2\sigma_{j}^{2}v^{(i)}_{j}=-(GM)_{ij}\sigma_{j},$
(2.67) $\displaystyle{\partial(\sigma_{i}v^{(l)}_{i})\over\partial y_{j}}$
$\displaystyle=-\frac{1}{2}\sum_{k}y^{-1}_{j}\frac{\partial
G_{ik}}{\partial\bar{\epsilon}_{j}}M_{kl}=-\sum_{k}y^{-1}_{j}\sigma_{j}v^{(i)}_{j}y_{j}G_{jk}M_{kl}$
$\displaystyle=-2\sigma^{2}_{j}v^{(i)}_{j}v^{(l)}_{j}=-(GM)_{ij}\sigma_{j}v^{(l)}_{j},$
(2.68)

where (2.2.3) is used. Substituting them into (2.66), we deduce that $y_{j}$
are normal modes [eq. (5.20), [12]]:

$\displaystyle\partial_{t_{l}}y_{j}+v^{(l)}_{j}\partial_{x}y_{j}=0.$ (2.69)

####  2.6.2 Family of flux Jacobians

Let us introduce the infinite dimensional matrices

$\displaystyle V=(v^{(i)}_{j})_{i,j\geq 1},\quad
V^{(l)}_{d}=\text{diag}(v^{(l)}_{1},v^{(l)}_{2},\ldots),\quad{\hat{\sigma}}=\text{diag}(\sigma_{1},\sigma_{2},\ldots),$
(2.70)

where the effective velocity is related to the hole density [eq. (4.19), [12]]
as

$\displaystyle
v^{(l)}_{i}=\sum_{k=1}^{\min(l,i)}\frac{\sigma_{l}}{\sigma_{k-1}\sigma_{k}}$
(2.71)

with $\sigma_{0}=1$. We therefore have an explicit expression of the hole
currents

$\displaystyle j_{i}^{(l)\rm
hole}=\sigma_{i}v^{(l)}_{i}=\sum_{k=1}^{\min(l,i)}\frac{\sigma_{i}\sigma_{l}}{\sigma_{k-1}\sigma_{k}}$
(2.72)

in terms of the hole densities $\sigma_{i}$. In the GHD context such an
expression is an equation of state [([15]), eq. (22)]. The results (2.67) and
(2.68) are expressed in matrix forms as

$\displaystyle\frac{\partial\sigma}{\partial
y}=-GM{\hat{\sigma}},\qquad\frac{\partial(\sigma v^{(l)})}{\partial
y}=-GMV^{(l)}_{d}{\hat{\sigma}}.$ (2.73)

If we use $\sigma_{i}$ as variables to express the currents, (2.66) is written
as

$\displaystyle\partial_{t_{l}}\sigma_{i}+\sum_{j=1}^{\infty}A^{(l)j}_{i}\partial_{x}\sigma_{j}=0.$
(2.74)

The matrix $A^{(l)}=(A^{(l)j}_{i})_{i,j\geq 1}$ (row $i$, column $j$) is the
flux Jacobian, where

$\displaystyle A^{(l)j}_{i}$
$\displaystyle={\partial{(\sigma_{i}v_{i}^{(l)})}\over\partial\sigma_{j}}=\sum_{k}\frac{\partial(\sigma_{i}v^{(l)}_{i})}{\partial
y_{k}}\frac{\partial y_{k}}{\partial\sigma_{j}}.$ (2.75)

In terms of matrices in (2.73) this is equivalent to

$\displaystyle A^{(l)}$
$\displaystyle=GMV^{(l)}_{d}{\hat{\sigma}}(GM{\hat{\sigma}})^{-1}=GMV^{(l)}_{d}M^{-1}G^{-1}=MG^{t}V^{(l)}_{d}(G^{t})^{-1}M^{-1},$
(2.76)

where $GM=MG^{t}$ is used. Note from (2.33) and (2.34) that
$GM=2{\hat{\sigma}}V^{t}=2V{\hat{\sigma}}$. Thus the above result is also
expressed as

$\displaystyle
A^{(l)}=VV^{(l)}_{d}V^{-1},\quad\text{i.e.}\quad\sum_{k=1}^{\infty}A_{i}^{(l)k}v^{(k)}_{j}$
$\displaystyle=v^{(l)}_{j}v^{(i)}_{j}\qquad(\forall i,j,l\in{\mathbb{Z}}_{\geq
1}).$ (2.77)

This confirms the fact that the effective velocities $v^{(l)}_{j}\,(1\leq
j\leq l)$ are the eigenvalues of $A^{(l)}$ [eq. (30) of [15]].

We remark that our flux Jacobians form a commuting family:

$\displaystyle[A^{(l)},A^{(l^{\prime})}]=0\qquad(A^{(1)}=1).$ (2.78)

This comes from the fact that the matrix $V$ which diagonalizes $A^{(l)}$ does
not depend on $l$ (only the eigenvalues depend on $l$). This property reflects
the commutativity $[T_{l},T_{l^{\prime}}]=0$ of the time evolutions in BBS.

Substitution of (2.71) into the first expression of (2.75) yields

$\displaystyle
A_{i}^{(l)j}=\sum_{k=1}^{\min(l,i)}\frac{\partial}{\partial\sigma_{j}}\Bigl{(}\frac{\sigma_{i}\sigma_{l}}{\sigma_{k-1}\sigma_{k}}\Bigr{)}.$
(2.79)

A little inspection of this shows that $A^{(l)}$ has the block structure whose
top left block is of size $l\times l$, the top right one is zero and the
bottom right one is $v^{(l)}_{l}$ times the identity matrix of infinite size:

$\displaystyle A^{(l)}$ $\displaystyle=\begin{pmatrix}\mathscr{A}^{(l)}&0\\\
\mathscr{U}&v^{(l)}_{l}\mathrm{I}\end{pmatrix},$ (2.80)
$\displaystyle\mathscr{A}^{(l)}$
$\displaystyle=\Bigl{(}\frac{\partial(\sigma_{i}v^{(l)}_{i})}{\partial\sigma_{j}}\Bigr{)}_{1\leq
i,j\leq
l},\quad\mathscr{U}=\begin{pmatrix}\sigma_{l+1}u_{1}&\sigma_{l+1}u_{2}\,\ldots&\sigma_{l+1}u_{l}\\\
\sigma_{l+2}u_{1}&\sigma_{l+2}u_{2}\,\ldots&\sigma_{l+2}u_{l}\\\
&\cdots&\end{pmatrix},\quad u_{j}=\frac{\partial
v^{(l)}_{l}}{\partial\sigma_{j}},$ (2.81)

For instance $\mathscr{A}^{(3)}$ reads as

$\displaystyle\mathscr{A}^{(3)}=\left(\begin{array}[]{ccc}0&0&1\\\
-\frac{(\sigma_{2}+1)\sigma_{3}}{\sigma_{1}^{2}}&\frac{\sigma_{3}}{\sigma_{1}}&\frac{\sigma_{2}+1}{\sigma_{1}}\\\
-\frac{(\sigma_{2}+1)\sigma_{3}^{2}}{\sigma_{1}^{2}\sigma_{2}}&-\frac{\sigma_{3}(\sigma_{1}+\sigma_{3})}{\sigma_{1}\sigma_{2}^{2}}&\frac{\sigma_{1}+2\sigma_{2}\sigma_{3}+2\sigma_{3}}{\sigma_{1}\sigma_{2}}\\\
\end{array}\right).$ (2.85)

####  2.6.3 Covariance and Drude matrices

Let us introduce the matrices $C=(C_{i,j})_{i,j\geq 1}$,
$B^{(l)}=(B^{(l)}_{i,j})_{i,j\geq 1}$ and $D^{(l)}=(D^{(l)}_{i,j})_{i,j\geq
1}$ whose elements are special cases of $C^{l,m}_{i,j}$ in (2.63):

$\displaystyle C_{i,j}$
$\displaystyle=C^{1,1}_{i,j}=\sum_{p}\rho_{p}\sigma_{p}(\rho_{p}+\sigma_{p})v^{(i)}_{p}v^{(j)}_{p},$
(2.86) $\displaystyle B^{(l)}_{i,j}$
$\displaystyle=C^{1,l}_{i,j}=\sum_{p}\rho_{p}\sigma_{p}(\rho_{p}+\sigma_{p})v^{(l)}_{p}v^{(i)}_{p}v^{(j)}_{p},$
(2.87) $\displaystyle D^{(l)}_{i,j}$
$\displaystyle=C^{l,l}_{i,j}=\sum_{p}\rho_{p}\sigma_{p}(\rho_{p}+\sigma_{p})(v^{(l)}_{p})^{2}v^{(i)}_{p}v^{(j)}_{p}.$
(2.88)

The quantities $C^{(l)}_{i,j}$ and $D^{(l)}_{i,j}$ are static covariance and
Drude weights [eqs.(158), (163), of [15]].

By using (2.77) it is easy to see

$\displaystyle A^{(l)}C=CA^{(l)t},\quad B^{(l)}=A^{(l)}C,\quad
D^{(l)}=\bigl{(}A^{(l)}\bigr{)}^{2}C=A^{(l)}CA^{(l)t}.$ (2.89)

In fact, the first relation is a special case of
$\sum_{k}A_{i}^{(l)k}{\tilde{C}}_{k,j}=\sum_{k}{\tilde{C}}_{i,k}A_{j}^{(l)k}$
for ${\tilde{C}}_{i,j}=\sum_{p}w_{p}v^{(i)}_{p}v^{(j)}_{p}$ with arbitrary
parameters $w_{k}$. The equality of the coefficient of each $w_{p}$ follows
directly from (2.77).

Recall the matrix form of the flux Jacobian in (2.76), i.e.,

$\displaystyle A^{(l)}=M(1+\hat{y}M)^{-1}V^{(l)}_{d}(1+\hat{y}M)M^{-1}.$
(2.90)

Introduce further diagonal matrices

$\displaystyle\hat{\rho}=\mathrm{diag}(\rho_{1},\rho_{2},\ldots),\qquad{\hat{g}}=1+{\hat{y}}=\mathrm{diag}(1+y_{1},1+y_{2},\ldots).$
(2.91)

Then the formula

$\displaystyle C_{i,j}$
$\displaystyle=\sum_{p}\rho_{p}(1+y_{p})(\sigma_{p}v^{(i)}_{p})(\sigma_{p}v^{(j)}_{p})=\frac{1}{4}\sum_{p}(GM)_{pi}\rho_{p}(1+y_{p})(GM)_{pj}$
(2.92)

derived from (2.33) and (2.34) is interpreted as
$C=\frac{1}{4}MG^{t}{\hat{\rho}}{\hat{g}}GM$ due to $GM=MG^{t}$. Applying
(2.89) to this and (2.90) we find

$\displaystyle C$
$\displaystyle=\frac{1}{4}M(1+\hat{y}M)^{-1}\hat{\rho}{\hat{g}}(1+M\hat{y})^{-1}M,$
(2.93) $\displaystyle B^{(l)}$
$\displaystyle=\frac{1}{4}M(1+\hat{y}M)^{-1}V^{(l)}_{d}\hat{\rho}{\hat{g}}(1+M\hat{y})^{-1}M,$
(2.94) $\displaystyle D^{(l)}$
$\displaystyle=\frac{1}{4}M(1+\hat{y}M)^{-1}(V^{(l)}_{d})^{2}\hat{\rho}{\hat{g}}(1+M\hat{y})^{-1}M.$
(2.95)

The factor $\frac{1}{4}$ originates in the coefficient 2 in
$M_{ij}=2\min(i,j)$, which is a non-essential artifact. The results (2.90) and
(2.93) – (2.95) agree with [eqs.(165)–(168), [15]].777 The correspondence with
eqs. (127), (149), (157) etc. in [15] as follows:
$\rho_{p}(k)\leftrightarrow\rho_{k}$, $\rho_{s}(k)\leftrightarrow\sigma_{k}$,
$n(k)\leftrightarrow(1+1/y_{k})^{-1}$, $f(p)\leftrightarrow(1+y_{p})^{-1}$,
$E^{\prime}(p)\leftrightarrow\min(l,p)$ (for $T_{l}$ BBS dynamics),
$h_{i}(p)\leftrightarrow\min(i,p)$,
$(E^{\prime})^{dr}(p)\leftrightarrow\sigma_{p}v^{(l)}_{p}$,
$h^{dr}_{i}(p)\leftrightarrow\sigma_{p}v^{(i)}_{p}$,
$1^{dr}(p)\leftrightarrow\sigma_{p}$,
$\int\frac{dp}{2\pi}\leftrightarrow\sum_{p=1}^{\infty}$.  We remark a natural
generalization

$\displaystyle C^{(l,m)}:=(C^{l,m}_{i,j})_{i,j\geq
1}=A^{(l)}C\bigl{(}A^{(m)}\bigr{)}^{t}(=C^{(m,l)})$ (2.96)

covers $C=C^{(1,1)},B^{(l)}=C^{(l,1)}=C^{(1,l)}$ and $D^{(l)}=C^{(l,l)}$.

###  2.7 Equal-time current-current correlations

As a slight digression, we consider here the spatially integrated current

$J=\sum_{x}j(x,0)$ (2.97)

and its variance

$f=L^{-1}\left(\left<J^{2}\right>-\left<J\right>^{2}\right)=\sum_{x}\left<j(x,0)j(0,0)\right>^{c},$
(2.98)

where the factor $L^{-1}$ insures that $f$ is finite in the thermodynamic
limit. The r.h.s of (2.98) is analogous to the current-current correlation
(2.43) which defines the Drude weight, except for the important fact it is an
equal-time correlation.

It is possible to compute $f$ using a transfer matrix approach, as explained
in B. $f$ of course differs from the Drude weight, but the numerical values
turn out to be close (see Tab. 3). Both diverge when approaching the half-
filled limit $p\to 1/2$ at $l=\infty$, and both are equal to $p(1-p)$ for
$l=1$ (in which case the BBS dynamics reduces to translations). We are however
unaware of a simple way to obtain $f$ using TBA. The fact that $f\geq D$ is a
general property, which follows from the Cauchy-Scwharz inequality associated
to a suitable scalar product between observables (see (3.14) of [37]).

Ball density $p$ | $l$ |  $D$ (2.44) (2.45) |  $f$ (2.98) (B.7)
---|---|---|---
$p$ | 1 | $p(1-p)$ | $p(1-p)$
0.2 | 2 | 0.6383506 | 0.672498
0.2 | 5 | 1.744263 | 1.857532
0.2 | 10 | 1.901627 | 2.023662
0.4 | 2 | 1.606536 | 1.790640
0.4 | 5 | 22.570957 | 26.04813
0.4 | 10 | 117.6194 | 134.2497390
Table 3:  Drude weight $D$ for a few values of the ball density $p$ and
carrier capacity $l$. For comparison the total current variance $f$ is also
given. For $l=1$ the BBS dynamics reduces to translations. In such a case the
spatially integrated current is trivially independent of time and $f=D$
(compare (2.43) and (2.98)).

##  3 Large deviation function

We are interested in the fluctuations of a quantity which, on average, grows
linearly with time. In such a case the large deviation function is a useful
tool to characterize the fluctuations. In general the large deviation function
contains the information about all the cumulants of the fluctuating quantity
(full counting statistics).111See however [35, 38] for an example where is
this not the case. It also describes the rate at which large fluctuations
occur [22]. The large deviation function is in general a difficult quantity to
obtain, but some important progress has been made recently in the context of
integrable systems [28, 29], where a general method to obtain the large
deviation function associated to the transport of a conserved charge has been
constructed. In this section we apply these ideas to the case of the BBS,
where the calculations simplify considerably.

###  3.1 General method

We want to compute the large deviation associated to the number $N_{t}$ of
balls transferred from the left to the right during time $t$, in some simple
GGE stationary state.222We refer the reader to C for the large deviation
function associated to the joint distribution of the number of transferred
balls and the number of transferred solitons. By definition, if the principle
of large deviation is obeyed, we have in the large time limit:

$F(\lambda)=\lim_{t\to\infty}\frac{1}{t}\ln\left[\langle
e^{\lambda{N_{t}}}\rangle_{\beta}\right].$ (3.1)

Or, in a more compact way:

$e^{tF(\lambda)}\sim\langle e^{\lambda{N_{t}}}\rangle_{\beta},$ (3.2)

where ${N_{t}}=\sum_{s=0}^{t}{j(0,s)}$ is the number of balls passing through
the origin during the time interval $[0,t]$. In the above expression the
expectation value $\left<f(t)\right>_{\beta}$ of an operator $f$ at time $t$
means:

$\langle f(t)\rangle_{\beta}={{\sum_{c}f(c(t))e^{-\beta
Q(c)}}\over\sum_{c}e^{-\beta Q(c)}}$ (3.3)

where $c(t)$ is the ball configuration $c$ evolved up to time $t$, and $Q(c)$
is (conserved) total number of balls in $c$. The parameter $\beta$ describing
the i.i.d. state is related to the ball fugacity by $e^{-\beta}=z$. We focus
here on a single conserved quantity (the total number of balls), but the
approach can be generalized to a GGE with several $\beta^{i}$ coupled to
several conserved energies, and a function $F$ of several variables
$\lambda^{i}$. See C for a two-$\beta$ case. Expanding $F$ in powers of
$\lambda$ gives access to the scaled cumulants $c_{n}$:

$F(\lambda)=\sum_{n=1}^{\infty}\frac{\lambda^{n}}{n!}c_{n}\;\;{\rm
and}\;\;c_{n}=\lim_{t\to\infty}\frac{1}{t}\left<(N_{t})^{n}\right>^{c}_{\beta}.$
(3.4)

From the expression above it appears clearly that $F$ is well-defined only if
all the cumulants $\left<(N_{t})^{n}\right>^{c}_{\beta}$ have the same
$t$-linear scaling. It should be noted that there are some integrable models
where the above relation is not obeyed [39, 38].

Taking the derivative with respect to $\lambda$, we have:

${dF(\lambda)\over d\lambda}=\lim_{t\to\infty}{1\over t}{{\langle
N_{t}e^{\lambda{N_{t}}}\rangle_{\beta}}\over{\langle e^{\lambda
N_{t}}\rangle_{\beta}}}$ (3.5)

where we assume that the above large-time limit exists.

Let us denote by $Q^{L}_{t}$ and $Q^{R}_{t}$ the charge in the left and right
halves of the system at time $t$. We have $Q^{L}_{t}+Q^{R}_{t}=Q$, independent
of time. Thus,

$N_{t}=Q^{L}_{0}-Q^{L}_{t}=Q-Q_{0}^{R}-Q_{t}^{L}.$ (3.6)

The relation (3.5) can be written

${dF(\lambda)\over d\lambda}=\lim_{t\to\infty}{1\over
t}{{\langle{{N_{t}}}e^{-\lambda(Q_{0}^{R}+Q_{t}^{L})}\rangle_{\beta-\lambda}}\over\langle
e^{-\lambda(Q_{0}^{R}+Q_{t}^{L})}\rangle_{\beta-\lambda}},$ (3.7)

where the factor $\exp(\lambda Q)$ has been converted into a shift
$\beta\to\beta-\lambda$.

In the i.i.d. states we consider the equal-time connected correlation of two
local operators, $\left<O(x_{1},t)O(x_{2},t)\right>^{c}$ vanish if $x_{1}\neq
x_{2}$. But since the ball propagation only takes place in the $x>0$
direction, we also know, by causality, that a two-time correlation of the form
$\left<O(x_{1},t_{1})O(x_{2},t_{2})\right>^{c}$ vanishes if
$(x_{1}-x_{2})(t_{1}-t_{2})<0$. This means that the local terms appearing in
$Q_{0}^{R}+Q_{t}^{L}$ (acting at $x>0$ and at time $0$ or at $x<0$ and time
$t>0$), and the terms in $N_{t}$ ($x=0$ and $t>0$) are uncorrelated. The term
$e^{-\lambda(Q_{0}^{R}+Q_{t}^{L})}$ thus decouples from $N_{t}$ (Fig. 7) and
cancels between the numerator and denominator of (3.7). This leads to

${dF(\lambda)\over d\lambda}=\lim_{t\to\infty}{1\over
t}\langle{N_{t}}\rangle_{\beta-\lambda}.$ (3.8)

The state defined by $\langle\cdots\rangle_{\beta-\lambda}$ is an i.i.d.
stationary state, so the expectation value of the currents is independent of
time, and we get

${dF(\lambda)\over d\lambda}={\langle j(0)}\rangle_{\beta-\lambda}$ (3.9)

where $j(0)$ is the ball current at the origin (and time zero).

The equation (3.9) is closely related to the so-called “flow equation” [29,
28]. As explained in these two works, for a general integrable system the
derivative ${dF(\lambda)\over d\lambda}$ of the SCGF is the expectation value
$\left<j\right>_{\tilde{\beta}(\lambda)}$ of the current in a modified GGE
state parameterized with modified inverse temperatures
$\tilde{\beta}^{i}(\lambda)$. These inverse temperatures are determined by
integrating the flow equation $\frac{d\tilde{\beta}^{i}}{d\lambda}=-{\rm
sign}(A(\lambda))_{i^{*}}^{i}$ from the initial condition where
$\tilde{\beta}^{i}(0)=\beta^{i}$ are the parameters of the original GGE
(single parameter $e^{-\beta}=z$ in the i.i.d. case we consider). The index
$i$ labels the conserved charges and $i^{*}$ is the index of the particular
charge for which the SCGF is computed. We compute here the SCGE associated to
the number of balls, so $i^{\ast}=\infty$. The matrix
$A_{ik}(\lambda)={\partial\left<j_{i}\right>_{\beta(\lambda)}}/{\partial
q_{k}}$, is the flux Jacobian (see Sec. 2.6), it appears when linearizing the
Euler equations. $A(\lambda)$ is defined by the derivatives of the current
densities with respect to the charges densities, evaluated in the
$\lambda$-modified state. Its eigenvalues are the effective velocities $v_{i}$
of the normal hydrodynamical modes (the soliton velocities in the BBS case).
${\rm sign}(A)$ is defined as the matrix with the same eigenvectors as $A$ but
replacing each eigenvalue $v_{i}$ by ${\rm sign}(v_{i})=\pm 1$. Here comes a
drastic simplification in the BBS case: all the soliton velocities are
positive and, therefore, ${\rm sign}(A)={\rm Id}$ and the flow equation
becomes $\frac{d\tilde{\beta}^{i}}{d\lambda}=-\delta_{i^{*},i}$. Since in the
present case $i^{\ast}=\infty$ the index $i$ is thus also fixed to $i=\infty$.
It follows that the $\lambda-$modification of the state is a simple shift of
the parameter $\beta$ associated to the ball density. The $\lambda-$modified
state remains i.i.d. but with $\tilde{\beta}(\lambda)=\beta-\lambda$.

Figure 7: The density operators $n(x>0,t=0)$ appearing in $Q^{R}_{0}$ are
uncorrelated from the current operators $j(x=0,t>0)$ in $N_{t}$. In the same
way, the density operators $n(x<0,t>0)$ appearing in $Q^{L}_{t}$ are
uncorrelated from $N_{t}$ and from $Q^{R}_{0}$.

###  3.2 Scaled cumulants generating function

In the r.h.s of (3.9) we need to consider the mean ball current in a modified
i.i.d. state with parameter $\beta^{\prime}=\beta-\lambda$, and hence a
modified ball fugacity $z^{\prime}=e^{-\beta+\lambda}$.

${{dF}(\lambda)\over d\lambda}=j(e^{-\beta+\lambda}).$ (3.10)

Writing the integration we obtain

$F(\lambda)=\int_{0}^{\lambda}j(e^{-\beta+\lambda^{\prime}})d\lambda^{\prime}=\int_{z}^{ze^{\lambda}}\frac{dz^{\prime}}{z^{\prime}}j(z^{\prime}).$
(3.11)

Reintroducing explicitly the parameter $l$ of the dynamics (carrier capacity)
and using the explicit form of the ball current (2.9), the integration in
(3.11) can be carried out explicitly and gives

${F^{(l)}(\lambda)}=\ln\left({1-(ze^{\lambda})^{l+1}\over
1-ze^{\lambda}}\right)-\ln\left({1-z^{l+1}\over 1-z}\right).$ (3.12)

As a sanity check we can compute its second derivative at $\lambda=0$:

$c_{2}=\left.\frac{d^{2}F^{(l)}}{d\lambda^{2}}\right|_{\lambda=0}$ (3.13)

and we recover (2.23).

The C presents some generalization of the above results to a two-temperature
GGE, with one parameter coupled to the total number of balls and another one
to the total number of solitons. In particular a generalization of (3.12) is
given in (C.23).

It is useful to consider the Legendre transform of $F$ with respect to
$\lambda$, the so-called large-deviation rate function:

$G(j)=j\lambda(j)-F(\lambda(j)),$ (3.14)

where $\lambda(j)$ is a solution of

$F^{\prime}(\lambda(j))=j.$ (3.15)

The probability to observe a transferred charge $N=j\,t$ at time $t$ is then

$P_{t}(N)\sim\exp\left(-t\,G(N/t)\right).$ (3.16)

$G(j)$ is defined for $0<j<l$, since the maximum possible value of the ball
current is $l$. It is a convex function obeying
$G(0)=\ln(({1-z^{l+1})/(1-z)})$, $G(l)=G(0)-l\ln(z)$ reaching its maximum zero
at the mean current $j=F^{\prime}(0)$. The function $G$ is represented in Fig.
8, in the vicinity of its maximum, for one particular case ($p=z/(1+z)=0.3$
and $l=10$).

Figure 8:  Probability distribution of the number $N_{t}$ of balls having
crossed the origin between $t=0$ and $t=400$ (blue), and between $t=0$ and
$t=2000$ (red). Ball density $p=0.3$ and capacity $l=10$. The numerical data
at $t=400$ are obtained using $6\,10^{6}$ random initial conditions and the
data at $t=2000$ involve $2\,10^{6}$ random initial conditions. These results
are compared with the Gaussian obtained by doing a least-square fit of the
data (dashed lines). Such a fit has two free parameters: mean and variance.
The data at $t=400$ are also compared to the Gaussian
$\exp[-t(j-<j>)^{2}/(2c_{2})]$ (magenta line), without any adjustable
parameter, where $<j>=0.7490144$ is the exact mean value of the current (2.9)
and $c_{2}=1.3016578$ is the exact second cumulant (2.23). Finally, the data
are compared to the exact large deviation rate function (continuous lines)
$\exp(-tG(N/t))$. At $t=400$ it is possible to see that the numerical data
departs from a simple Gaussian and are instead consistent with the theoretical
prediction involving the large deviation rate function $G$ (3.14).

###  3.3 $l=\infty$

In the limit $l\to\infty$ the SCGF (3.12) simplifies to

${F^{(\infty)}(\lambda)}=\ln\left({1-z\over 1-ze^{\lambda}}\right).$ (3.17)

The interval where $F$ is defined is $\lambda\in[-\infty,\beta]$, which is
equivalent to $ze^{\lambda}<1$ (recall that $z=e^{-\beta}$). For
$\lambda\to-\infty$ the current $j=F^{\prime}(\lambda)$ tends to zero, and for
$\lambda\to\beta$ we have instead $j=F^{\prime}(\lambda)\to+\infty$. The full
range of physical values for the ball current is covered by $F^{\prime}$, as
it should. A generalization of (3.17) to a two-temperature GGE is given in
(C.26).

In that case the large-deviation rate function can be obtained explicitly:

$G^{(\infty)}(j)=-\ln(1-z)-j\ln(z)-(1+j)\ln(1+j)+j\ln(j).$ (3.18)

This function is plotted in Fig. 9 for a few different values of the ball
fugacity. The generalization to a two-temperature GGE is given in (C.28). When
$z\to 1^{-}$ $F^{(\infty)}(\lambda)$ becomes singular in $\lambda=0$ and the
associated cumulants diverge. For instance: $c_{2}=z/(1-z)^{2}$ at $l=\infty$.
This can be interpreted as a phase transition when the ball density approaches
$1/2$. Some other properties of this transition were discussed in [3].

Figure 9: Large-deviation rate function $G^{(l=\infty)}(j)$ plotted as a
function of the current $j$, for different values of the ball fugacity $z$.
The mean current (minimum of $G$) and the second cumulant $c_{2}$ diverge when
$z\to 1^{-}$.

To conclude this section we note that the large deviation principle for the
energies $E_{k}$ in i.i.d. states was studied for a multicolor BBS by
regarding the history of carriers going through the states as a Markov process
[11]. In that previous study the system size is the variable which plays a
role similar to the role played by time here.

##  4 Conclusions

We have computed analytically several quantities related to the current and
density fluctuations in stationary i.i.d. states. We have obtained the SCGF
associated to the number $N_{t}$ of balls crossing the origin during time $t$,
from which all the cumulants can be extracted. The Legendre transform of this
function could be compared with the probability distribution of $N_{t}$
extracted from numerical simulations. This is one of the very few interacting
and deterministic models where the SCGF could be computed exactly ([38] is
another recent example).

The Drude weights – defined as the long-time limit of a spatially integrated
current-current correlation – could also be obtained using TBA combined with
hydrodynamical projection. Explicit analytical expressions for the Drude
weights could be compared successfully with numerical simulations.

The existence of a family of commuting time evolutions is an important
property of integrable systems, although implementing them in actual systems
can be cumbersome in practice. The BBS is a distinguished example of
integrable cellular automata where the all commuting time evolutions have a
simple and neat implementation. Exploiting these time evolutions a set of new
generalized currents correlations (or generalized Drude weights) was
constructed and shown to enjoy unexpected symmetry relations. Some of these
symmetry relations can be explained at the microscopic level, while others
only emerge in the long-time limit. We observed that the numerical results for
the current correlations at long times coincide with the TBA results only for
certain time evolution, namely $T_{n}$ with sufficiently large $n$. This
suggests that $T_{n}$ for small $n$ is insufficiently mixing and the observed
correlations in this regime still escape our understanding. Understanding
these results and reconciling them with the hydrodynamic projection ideas
would certainly deserve further study.

We also stress that a number of these results could be generalized to a larger
family of (non i.i.d.) GGE states with two temperatures, respectively coupled
to the number of balls and to the number of solitons.

It is quite remarkable that so many explicit formulae could be obtained for
nontrivial quantities related to long-distance and long-time limit of
correlations in such an out-of-equilibrium interacting problem. These could
prove to be useful to compare BBS with other models, either integrable or non-
integrable, and to shed some light about fundamental questions like the
emergence of hydrodynamics.

##  5 Acknowledgements

V. P. thanks E. Ilievski, Z. Krajnik, T. Prosen and J. Schmidt for numerous
discussions and collaboration on related subjects. We also acknowledge the DRF
of CEA for providing us with CPU time on the supercomputer TOPAZE at CCRT.

##  Appendix A Correlation sum rules

We derive a sum rule connecting density-density correlations to current-
current ones. The argument is directly inspired from [30]. We consider the
following quantity

$\sum_{x=-\infty}^{\infty}f(x)\left(n(x,t)-n(x,0)\right),$ (A.1)

where $n$ is the ball density operator, and $f$ some test function defined on
the lattice sites. Next we introduce the charge of the region $]-\infty,x]$:

$\displaystyle N^{\langle}(x,t)$ $\displaystyle=\sum_{y=-\infty}^{x}n(y,t)$
(A.2) $\displaystyle n(x,t)$
$\displaystyle=N^{\langle}(x,t)-N^{\langle}(x-1,t)$ (A.3)

(A.1) can be rewritten as

$\displaystyle\sum_{x=-\infty}^{\infty}f(x)\left(n(x,t)-n(x,0)\right)$
$\displaystyle=\sum_{x=-\infty}^{\infty}f(x)\left(N^{\langle}(x,t)-N^{\langle}(x-1,t)-N^{\langle}(x,0)+N^{\langle}(x-1,0)\right)$
(A.4)
$\displaystyle=\sum_{x=-\infty}^{\infty}\left(f(x)-f(x+1)\right)\left(N^{\langle}(x,t)-N^{\langle}(x,0)\right)$
(A.5)
$\displaystyle=-\sum_{x=-\infty}^{\infty}\left(f(x)-f(x+1)\right)\int_{0}^{t}j(x+1,s)ds$
(A.6)
$\displaystyle=\sum_{x=-\infty}^{\infty}\left(f(x)-f(x-1)\right)\int_{0}^{t}j(x,s)ds.$
(A.7)

The current can also be used to write

$n(0,t)-n(0,0)=\int_{0}^{t}\left(j(0,s)-j(1,s)\right)ds.$ (A.8)

We then multiply (A.7) by (A.8)

$\displaystyle\sum_{x=-\infty}^{\infty}f(x)\left(n(x,t)-n(x,0)\right)\left(n(0,t)-n(0,0)\right)$
$\displaystyle=\sum_{x=-\infty}^{\infty}\left(f(x)-f(x-1)\right)\int_{0}^{t}j(x,s)ds\int_{0}^{t}\left(j(0,s^{\prime})-j(1,s^{\prime})\right)ds^{\prime}$
(A.9)

and take the connected average:

$\displaystyle\sum_{x=-\infty}^{\infty}f(x)\left<\left(n(x,t)-n(x,0)\right)\left(n(0,t)-n(0,0)\right)\right>^{c}$
$\displaystyle=\sum_{x=-\infty}^{\infty}\int_{0}^{t}\int_{0}^{t}dsds^{\prime}\left(f(x)-f(x-1)\right)\left<j(x,s)\left(j(0,s^{\prime})-j(1,s^{\prime})\right)\right>^{c}$
(A.10)
$\displaystyle=\sum_{x=-\infty}^{\infty}\int_{0}^{t}\int_{0}^{t}dsds^{\prime}\left(f(x)-f(x-1)\right)\left<\left(j(x,s)-j(x-1,s)\right)j(0,s^{\prime})\right>^{c}$
(A.11)

where, in the last equality, we have used the translation invariance of the
current-current correlator. This can finally be rewritten

$\displaystyle\sum_{x=-\infty}^{\infty}f(x)\left<\left(n(x,t)-n(x,0)\right)\left(n(0,t)-n(0,0)\right)\right>^{c}$
$\displaystyle=\sum_{x=-\infty}^{\infty}\int_{0}^{t}\int_{0}^{t}dsds^{\prime}\left(2f(x)-f(x+1)-f(x-1)\right)\left<j(x,s)j(0,s^{\prime})\right>^{c}.$
(A.12)

##  Appendix B Transfer matrix calculation for the current fluctuations

This section presents a transfer matrix calculation of the fluctuations of the
current (2.98). It is an extension of the calculation presented in the
Appendix E of [12].

Here we consider the periodic BBS model on a chain of length $L$. It admits a
discrete set of commuting evolutions characterized by a row to row transfer
propagator $T_{l}$ labelled by a positive integer $l$ representing the
capacity of the carrier. The propagator $T_{l}$ takes the form of a vertex
transfer matrix (see Fig. 10) where the horizontal links can contain up to $l$
balls which are auxiliary variables: $0\leq n^{c}\leq l$. The vertical links
contain zero or one ball: $0\leq n^{b}\leq 1$ which are the BBS variables.
Time flows down, so that at each time step, the south vertical links are
occupied according to the north vertical link configurations. We can see the
evolution as the result of the passage of a carrier transporting up to $l$
balls from west to east and updating each vertical link successively passing
through the vertices. If the north vertical link is empty ($n^{b}_{N}=0$) and
the carrier has at least one ball ($n^{c}_{W}>0$), it leaves one of its balls
to the south vertical link during the passage ($n^{b}_{S}=1$,
$n^{c}_{E}=n^{c}_{W}-1$). If it does not carry balls ($n^{c}_{W}=0$), it
passes without changing either the vertical occupation or its load ($n_{S}=0$,
$n^{c}_{E}=0$). If the north vertical link is filled ($n^{b}_{N}=1$) and the
carrier carries strictly less than $l$ balls ($n^{c}_{W}<l$), it picks up a
ball and leaves the south vertical link empty ($n^{c}_{E}=n^{c}_{W}+1$,
$n^{b}_{S}=0$). If it carries $l$ balls ($n^{c}_{W}=l$), it passes without
changing either the vertical occupation or its load ($n^{b}_{S}=1$,
$n^{c}_{E}=l$). On a periodic chain, it can be shown (Proposition 5.1 of [2])
that the periodicity condition on the horizontal links uniquely determines the
load of the carrier.111To be precise: it is so when the density is not exactly
equal to 1/2.

$n^{c}_{W}=0$$n^{c}_{E}=n^{c}_{W}$$n^{b}_{S}=0$$n^{b}_{N}=0$

$n^{c}_{W}<l$$n^{c}_{E}=n^{c}_{W}+1$$n^{b}_{S}=0$$n^{b}_{N}=1$

$n^{c}_{W}>0$$n^{c}_{E}=n^{c}_{W}-1$$n^{b}_{S}=1$$n^{b}_{N}=0$

$n^{c}_{W}=l$$n^{c}_{E}=n^{c}_{W}$$n^{b}_{S}=1$$n^{b}_{N}=1$

Figure 10: Four examples of vertices.

As a result, the transfer propagator can be expressed as:

$(T_{l})^{n^{b}(t)}_{n^{b}(t+1)}=\left<n^{b}(t+1)\right|{\rm
Tr}\left(\prod_{k=1}^{L}\mathbb{L}_{k}\right)\left|n^{b}(t)\right>$ (B.1)

where $n_{b}(t)$ stands for the configuration $\\{n^{b}(x,t),1\leq x\leq L\\}$
and $\langle n^{b}_{S}|\mathbb{L}|n^{b}_{N}\rangle=\langle
n^{b}_{S}|\mathbb{L}|n^{b}_{N}\rangle_{n^{c}_{W},n^{c}_{E}}$ is a
$(l+1)\times(l+1)$ permutation matrix with nonzero matrix elements equal to
one whenever the corresponding vertex is allowed.

Due to charge conservation, $n_{N}+n_{W}=n_{S}+n_{E}$, the ball density
$n^{b}(x,t)$ coincides with the number of balls on the vertical links and the
ball current $j(x,t)=n^{c}(x,t)$ to the number of balls on the horizontal
links.

We can give a fugacity to the balls by inserting the operator
$z^{Q}=z^{\sum_{x=1}^{L}n^{b}(x,t)}$ which commutes with the propagator.
Denote $|0\rangle$ the state obtained by summing all the ball configurations
with weight one. The i.i.d. stationary state with ball density $p=z/(z+1)$ is
$z^{Q}|0\rangle$.

We can also give a fugacity to the total current at time $t$ by weighting the
horizontal links by $y^{J(t)}=y^{\sum_{x=1}^{L}j(x,t)}$. We denote by
$T_{l}(y)$ (or $\mathbb{L}(y)$ for a single site) the modified propagator
where we weight each vertex by $y^{n^{c}_{W}}$ instead of $1$, which has the
effect to weight each configuration by $y^{J(t)}$.

(B.2)

Consider the matrix element $Z=\langle 0|z^{Q}y^{J(t)}|0\rangle={\rm
Tr}(\mathbb{L}^{0}(y)+z\mathbb{L}^{1}(y))^{L}$ where $\mathbb{L}^{0}(y),\
\mathbb{L}^{1}(y)$ are the deformed vertex matrices with $n^{b}_{N}=0,1$
respectively ($n^{b}_{S}$ is redundant).

Denoting $\mathbb{L}^{z}(y)=\mathbb{L}^{0}(y)+z\mathbb{L}^{1}(y)$, we have:

$\displaystyle\langle 0|\mathbb{L}^{z}(y)$ $\displaystyle=\langle 0|+z\langle
1|$ $\displaystyle\langle n|\mathbb{L}^{z}(y)$ $\displaystyle=y^{n}\langle
n-1|+zy^{n}\langle n+1|\ {\rm if}\ n\neq 0,l$ $\displaystyle\langle
l|\mathbb{L}^{z}(y)$ $\displaystyle=y^{l}\langle l-1|+zy^{l}\langle l|.$ (B.3)

As an example, for $l=3$ we have

$\displaystyle\mathbb{L}^{z}(y)$ $\displaystyle=\begin{pmatrix}1&z&0&0\\\
y&0&zy&0\\\ 0&y^{2}&0&zy^{2}\\\ 0&0&y^{3}&zy^{3}\end{pmatrix}.$ (B.4)

The matrix element $\mathbb{L}^{z}(1)_{\alpha\beta}/(z+1)$ is the probability
to have $\beta$ balls on the east link knowing that there are $\alpha$ on the
west link. Therefore, the components of the left eigenvector of
$\mathbb{L}^{z}(1)/(z+1)$ with eigenvalue one, $\sum_{k=0}^{l}z^{k}\langle
k|$, are the unnormalized probabilities for the carrier to contain $k$ balls.
The current formula (2.9) is obtained as the average number of balls in the
carrier.

The second cumulant $c_{2}$ and $f$, the integrated current-current connected
correlation, are then given by second derivatives of $\ln(Z)$:

$c_{2}=\frac{1}{L}\left(z\frac{\partial}{\partial
z}\right)\left(y\frac{\partial}{\partial y}\right)\ln(Z)|_{y=1},$ (B.5)
$f=\langle J(0)j(0,0)\rangle^{c}=\frac{1}{L}\left(y\frac{\partial}{\partial
y}\right)^{2}\ln(Z)|_{y=1}.$ (B.6)

By analyzing the values of $f$ for various values of $l$, it has also been
possible to conjecture the following analytical expression for $f$:

$\begin{split}f=&\frac{z(1+6z+z^{2})(1-z^{3l+3})}{(1-z)^{4}(1-z^{l+1})^{3}}-\frac{24z^{l+5}(1-z^{l-3})}{(1-z)^{4}(1-z^{l+1})^{3}}-\frac{3z^{l+2}(g_{l+2}+7zg_{l}+8z^{2}g_{l-2})}{(1-z)^{3}(1-z^{l+1})^{3}}\\\
&-\frac{(l+1)^{2}z^{l+1}(g_{l+2}+2l(1+z)g_{l+1}+3zg_{l})}{(1-z)(1-z^{l+1})^{3}},\end{split}$
(B.7)

where $g_{j}=1+z^{j}$.

##  Appendix C Drude weight and SCGF in a two-temperature GGE

The BBS has the conserved quantities $E_{k}=L\sum_{j\geq 1}\min(j,k)\rho_{j}$
for $k=1,2,\ldots$ [12]. $E_{\infty}$ is the total number of balls with which
the main text is concerned. In this appendix we present a partial
generalization of the results to the two-temperature GGE corresponding to the
statistical weight $e^{-\beta_{1}E_{1}-\beta_{\infty}E_{\infty}}$. The
conserved quantity $E_{1}$ is the number of solitons. The inverse temperature
$\beta$ in the main text is denoted by $\beta_{\infty}$ here.

Following [12] (3.20) we parameterize the temperatures by $a,z$ as

$\displaystyle\mathrm{e}^{\frac{1}{2}\beta_{1}}=\frac{a^{\frac{1}{2}}-a^{-\frac{1}{2}}}{z^{\frac{1}{2}}-z^{-\frac{1}{2}}},\qquad\mathrm{e}^{-\beta_{\infty}}=z.$
(C.1)

The single temperature GGE($\beta_{\infty}$) corresponds to the limit
$a\rightarrow z$. Derivatives by the temperatures are expressed in terms of
$a$ and $z$ as in [eq. (3.22) and (3.23) in[12]]:

$\displaystyle\frac{\partial}{\partial\beta_{1}}$
$\displaystyle=\frac{\partial a}{\partial\beta_{1}}\frac{\partial}{\partial
a}+\frac{\partial z}{\partial\beta_{1}}\frac{\partial}{\partial
z}=-\frac{a(1-a)}{1+a}\frac{\partial}{\partial a},$ (C.2)
$\displaystyle\frac{\partial}{\partial\beta_{\infty}}$
$\displaystyle=\frac{\partial
a}{\partial\beta_{\infty}}\frac{\partial}{\partial a}+\frac{\partial
z}{\partial\beta_{\infty}}\frac{\partial}{\partial
z}=-\frac{a(1-a)(1+z)}{(1+a)(1-z)}\frac{\partial}{\partial
a}-z\frac{\partial}{\partial z}.$ (C.3)

Densities, effective velocity and currents have been obtained in [12] as
follows ((C.7) was not included therein):

$\displaystyle\sigma_{k}$
$\displaystyle=\frac{(1-a)(1+az^{k})}{(1+a)(1-az^{k})},\qquad\rho_{k}=\frac{az^{k-1}(1-a)(1-z)^{2}(1+az^{k})}{(1+a)(1-az^{k-1})(1-az^{k})(1-az^{k+1})},$
(C.4) $\displaystyle v^{(l)}_{k}$
$\displaystyle=\frac{1+az^{l}}{1-az^{l}}v_{\min(k,l)},\qquad
v_{k}=\frac{1+a}{1-a}k-\frac{2a(1+z)(1-z^{k})}{(1-a)(1-z)(1+az^{k})}.$ (C.5)
$\displaystyle\text{ball current}:\;j^{(l)}_{\infty}=\sum_{k\geq
1}k\rho_{k}v^{(l)}_{k}=\frac{a(1+z)(1-z^{l})}{(1+a)(1-z)(1-az^{l})}-\frac{laz^{l}}{1-az^{l}},$
(C.6) $\displaystyle\text{soliton current}:\;j^{(l)}_{1}=\sum_{k\geq
1}\rho_{k}v^{(l)}_{k}=\frac{a(1-z^{l})}{(1+a)(1-az^{l})}.$ (C.7)

The results (C.6) and (C.7) are the $j=\infty$ and $j=1$ cases of (2.57):

$\displaystyle\eta^{(l)}_{j}=\frac{a(1+z)(1-z^{\min(j,l)})(1+az^{\max(j,l)})}{(1+a)(1-z)(1-az^{j})(1-az^{l})}-\frac{\min(j,l)a(z^{j}+z^{l})}{(1-az^{j})(1-az^{l})}.$
(C.8)

Set $W_{i}=\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})v^{2}_{i}$ as in the main
text. Then the sum formula (2.51) admits the following generalization:

$\displaystyle\sum_{i<l}W_{i}$
$\displaystyle=\frac{a(1-a)(1+z)A_{l}}{(1+a)^{3}(1-z)(1-az^{l-1})^{2}(1-az^{l})^{2}},$
(C.9) $\displaystyle A_{l}$
$\displaystyle=1+2a(1+a^{2})z^{2l-1}+a^{4}z^{4l-2}+a^{2}z^{2l-2}(1+8z+z^{2})-az^{l-1}(2+a^{3}z^{2l-1})g_{l}(z)$
$\displaystyle-z^{l}(1+2a^{3}z^{2l-1})g_{l}(z^{-1})-a^{2}z^{l-1}h_{l}(z)-a^{2}z^{3l-1}h_{l}(z^{-1}),$
(C.10) $\displaystyle g_{l}(z)$
$\displaystyle=\frac{z+\bigl{(}1-l(1-z)\bigr{)}^{2}}{1+z},\quad
h_{l}(z)=(1+l)^{2}+l^{2}z+\frac{(1-2l+z)(3+z+2lz)}{1+z}.$ (C.11)

This is easily shown by induction on $l$ and $A_{0}=A_{1}=0$. The above
$A_{l}$ precisely reduces to (2.51) in the single temperature case $a=z$.

Define the second cumulant $c_{2}^{(l)}$ and the Drude weight $D^{(l)}$ by
(2.41) and (2.45) with $\rho_{i},\sigma_{i},v^{(l)}_{i}$ specified in
(C.4)–(C.5). From $v^{(1)}_{i}=1$ we have

$\displaystyle D^{(1)}$ $\displaystyle=\sum_{i\geq
1}W_{i}=\frac{a(1-a)(1+z)A_{\infty}}{(1+a)^{3}(1-z)}=\frac{a(1-a)(1+z)}{(1+a)^{3}(1-z)},$
(C.12) $\displaystyle c^{(l)}_{2}$
$\displaystyle=\frac{a(1-a)(1+z)^{2}(1-z^{l})(1+a^{2}z^{l})}{(1+a)^{3}(1-z)^{2}(1-az^{l})^{2}}+\frac{2az(1-z^{l})}{(1+a)(1-z)^{2}(1-az^{l})}$
$\displaystyle\quad-\frac{alz^{l}}{(1-az^{l})^{2}}\left(l+\frac{2(1-a)(1+z)}{(1+a)(1-z)}\right)$
(C.13) $\displaystyle=-\frac{\partial
j^{(l)}_{\infty}}{\partial\beta_{\infty}}.$ (C.14)

For $D^{(l)}$ with general $l$, formally the same formula as (2.54) with
$\frac{z}{(1+z)^{2}}$ replaced by (C.12) is valid. These results reduce to the
single temperature case at $a=z$. Another useful sum formula is

$\displaystyle\sum_{i\geq r}\rho_{i}\sigma_{i}(\rho_{i}+\sigma_{i})$
$\displaystyle=\frac{a(1-a)^{3}(1-z)z^{r-1}\left(1+a^{2}z^{2r-1}\right)}{(1+a)^{3}\left(1-az^{r-1}\right)^{2}\left(1-az^{r}\right)^{2}}.$
(C.15)

The ball density $p$ in the two-temperature GGE is known to be
$p=\frac{a}{1+a}$ in [[12] (3.23)]. Unlike (2.46), the result (C.12) does not
coincide with the $p(1-p)$ reflecting the fact that the two-temperature GGE
under consideration is not i.i.d.

It is natural to introduce the joint cumulant generating function

$\displaystyle F(\lambda,\mu)=\lim_{t\rightarrow\infty}\frac{1}{t}\ln\langle
e^{\lambda N_{\infty,t}+\mu N_{1,t}}\rangle_{\beta_{\infty},\beta_{1}},$
(C.16)

where the superscript $l$ is suppressed and
$N_{\infty,t}=\int_{0}^{t}j^{(l)}_{\infty}(0,\tau)d\tau$ and
$N_{1,t}=\int_{0}^{t}j^{(l)}_{1}(0,\tau)d\tau$. By the same argument as
before, we have

$\displaystyle\frac{\partial
F(\lambda,\mu)}{\partial\lambda}=j^{(l)}_{\infty}(\beta_{\infty}-\lambda,\beta_{1}-\mu),\qquad\frac{\partial
F(\lambda,\mu)}{\partial\mu}=j^{(l)}_{1}(\beta_{\infty}-\lambda,\beta_{1}-\mu).$
(C.17)

Let $(\zeta,\alpha)$ be the parameters $(z,a)$ corresponding to
$(\beta_{\infty}-\lambda,\beta_{1}-\mu)$ in the sense of (C.1). Then from
(C.6) and (C.7) the currents are expressed as

$\displaystyle j^{(l)}_{\infty}(\beta_{\infty}-\lambda,\beta_{1}-\mu)$
$\displaystyle=\frac{\alpha(1+\zeta)(1-\zeta^{l})}{(1+\alpha)(1-\zeta)(1-\alpha\zeta^{l})}-\frac{l\alpha\zeta^{l}}{1-\alpha\zeta^{l}},$
(C.18) $\displaystyle j^{(l)}_{1}(\beta_{\infty}-\lambda,\beta_{1}-\mu)$
$\displaystyle=\frac{\alpha(1-\zeta^{l})}{(1+\alpha)(1-\alpha\zeta^{l})}.$
(C.19)

Here $(\zeta,\alpha)$ are regarded as functions of $(\lambda,\mu)$ including
$(z,a)$ as parameters. Namely,

$\displaystyle\frac{\alpha^{\frac{1}{2}}-\alpha^{-\frac{1}{2}}}{\zeta^{\frac{1}{2}}-\zeta^{-\frac{1}{2}}}=\mathrm{e}^{\frac{1}{2}(\beta_{1}-\mu)}=\frac{a^{\frac{1}{2}}-a^{-\frac{1}{2}}}{z^{\frac{1}{2}}-z^{-\frac{1}{2}}}e^{-\frac{1}{2}\mu}\qquad\zeta=\mathrm{e}^{-\beta_{\infty}+\lambda}=ze^{\lambda}.$
(C.20)

They imply the derivative relations similar to (C.2) and (C.3):

$\displaystyle\frac{\partial}{\partial\lambda}=\frac{\alpha(1-\alpha)(1+\zeta)}{(1+\alpha)(1-\zeta)}\frac{\partial}{\partial\alpha}+\zeta\frac{\partial}{\partial\zeta},\qquad\frac{\partial}{\partial\mu}=\frac{\alpha(1-\alpha)}{1+\alpha}\frac{\partial}{\partial\alpha}.$
(C.21)

By using (C.18), (C.19) and (C.21), one can directly check the consistency of
(C.17):

$\displaystyle\frac{\partial
j^{(l)}_{\infty}(\beta_{\infty}-\lambda,\beta_{1}-\mu)}{\partial\mu}=\frac{\partial
j^{(l)}_{1}(\beta_{\infty}-\lambda,\beta_{1}-\mu)}{\partial\lambda}.$ (C.22)

The solution to (C.17) satisfying $F(0,0)=0$ is given by

$\displaystyle
F(\lambda,\mu)=\ln\left(\frac{1-\alpha\zeta^{l}}{1-\alpha}\right)-\ln\left(\frac{1-az^{l}}{1-a}\right),$
(C.23)

where $\alpha=\alpha(\lambda,\mu)$ and
$\zeta=\zeta(\lambda,\mu)(=ze^{\lambda})$ are specified as the solution to
(C.20) satisfying $\alpha>0$. Note that
$\alpha(\lambda,\mu=0)|_{a=z}=ze^{\lambda}$. Thus, the result (C.23) provides
a generalization of (3.12) reproducing the latter as
$F(\lambda,\mu=0)|_{a=z}=F^{(l)}(\lambda)$.

One can check
$\left.\frac{\partial^{2}F}{\partial\lambda^{2}}\right|_{\lambda=\mu=0}=c_{2}^{(l)}$.
The other second order scaled cumulants are given by

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\langle
N_{\infty,t}N_{1,t}\rangle^{c}_{\beta_{\infty},\beta_{1}}$
$\displaystyle=\left.\frac{\partial^{2}F}{\partial\lambda\partial\mu}\right|_{\lambda=\mu=0}=\frac{a(1-a)\bigl{(}(1+z)(1-z^{l})(1+a^{2}z^{l})-(1+a)^{2}(1-z)lz^{l}\bigr{)}}{(1+a)^{3}(1-z)(1-az^{l})^{2}},$
(C.24) $\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\langle
N_{1,t}^{2}\rangle^{c}_{\beta_{\infty},\beta_{1}}$
$\displaystyle=\left.\frac{\partial^{2}F}{\partial\mu^{2}}\right|_{\lambda=\mu=0}=\frac{a(1-a)(1-z^{l})(1+a^{2}z^{l})}{(1+a)^{3}(1-az^{l})^{2}}.$
(C.25)

We leave an interesting problem of studying the Hessian of $F$ in relation to
the convexity of $F$ for a future work.

When $l=\infty$, (C.23) in the regime $ze^{\lambda}<1$ simplifies to

$\displaystyle
F^{(\infty)}(\lambda,\mu)=\ln\left(\frac{1-a}{1-\alpha}\right).$ (C.26)

This is still a non-trivial function of $\lambda$ and $\mu$ via (C.20). From
(C.17) – (C.19) with $l=\infty$, the equation $\Bigl{(}\frac{\partial
F^{(\infty)}(\lambda,\mu)}{\partial\lambda},\frac{\partial
F^{(\infty)}(\lambda,\mu)}{\partial\mu}\Bigr{)}=(\mathscr{J}_{\infty},\mathscr{J}_{1})$
has the solution

$\displaystyle\lambda_{\ast}=\ln\left(\frac{\mathscr{J}_{\infty}-\mathscr{J}_{1}}{z(\mathscr{J}_{\infty}+\mathscr{J}_{1})}\right),\quad\mu_{\ast}=\ln\left(\frac{4(1-\mathscr{J}_{1})\mathscr{J}_{1}^{3}(a+a^{-1}-2)}{(\mathscr{J}_{\infty}^{2}-\mathscr{J}_{1}^{2})(1-2\mathscr{J}_{1})^{2}(z+z^{-1}-2)}\right),$
(C.27)

which is deduced from
$\alpha_{\ast}=\frac{\mathscr{J}_{1}}{1-\mathscr{J}_{1}}$,
$\zeta_{\ast}=\frac{\mathscr{J}_{\infty}-\mathscr{J}_{1}}{\mathscr{J}_{\infty}+\mathscr{J}_{1}}$.
Now the large deviation rate function
$G^{(\infty)}(\mathscr{J}_{\infty},\mathscr{J}_{1})=\mathscr{J}_{\infty}\lambda_{\ast}+\mathscr{J}_{1}\mu_{\ast}-F^{(\infty)}(\lambda_{\ast},\mu_{\ast})$
is given by

$\begin{split}G^{(\infty)}(\mathscr{J}_{\infty},\mathscr{J}_{1})&=\mathscr{J}_{\infty}\ln\left(\frac{\mathscr{J}_{\infty}-\mathscr{J}_{1}}{z(\mathscr{J}_{\infty}+\mathscr{J}_{1})}\right)+\mathscr{J}_{1}\ln\left(\frac{4(1-\mathscr{J}_{1})\mathscr{J}_{1}^{3}(a+a^{-1}-2)}{(\mathscr{J}_{\infty}^{2}-\mathscr{J}_{1}^{2})(1-2\mathscr{J}_{1})^{2}(z+z^{-1}-2)}\right)\\\
&-\ln\left(\frac{(1-a)(1-\mathscr{J}_{1})}{1-2\mathscr{J}_{1}}\right).\end{split}$
(C.28)

Since every soliton carries at least one ball, the ball current
$\mathscr{J}_{\infty}$ and the soliton current $\mathscr{J}_{1}$ are to be
considered in the domain $\mathscr{J}_{\infty}>\mathscr{J}_{1}$. Similarly,
from (C.7) with $l\rightarrow\infty$, one should suppose $\mathscr{J}_{1}<1$.
Note that $G^{(\infty)}(\mathscr{J}_{\infty},\mathscr{J}_{1})$ is finite at
$\mathscr{J}_{1}=\frac{1}{2}$.

The single temperature case in the main text corresponds to setting $a=z$,
$\mu=0$. In fact, (C.20) enforces $\alpha_{\ast}=\zeta_{\ast}$ leading to
$\mathscr{J}_{1}=\frac{\mathscr{J}_{\infty}}{1+2\mathscr{J}_{\infty}}$ and
$\mu_{\ast}=0$. Then (C.28) reduces $G^{(\infty)}(j=\mathscr{J}_{\infty})$ in
the main text.

##  Appendix D Microscopic definition of $\hat{\eta}^{(l)}_{i}(x)$

Here we outline the proof of the properties of the generalized current
$\hat{\eta}^{(l)}_{i}(x)$ claimed in the beginning of Section 2.4. We begin by
recalling the combinatorial $R$ in the crystal base theory, a theory of
quantum groups at $q=0$, for which readers are referred to [Section 2.2, [2]]
and the references therein.

For a positive integer $l$, define a set
$B_{l}=\\{\alpha=(\alpha_{0},\alpha_{1})\in({\mathbb{Z}}_{\geq
0})^{2}\mid\alpha_{0}+\alpha_{1}=l\\}$. Introduce an infinite set known as an
affine crystal
$\mathrm{Aff}(B_{l})=B_{l}\times\zeta^{{\mathbb{Z}}}=\\{\alpha\zeta^{n}\mid\alpha\in
B_{l},n\in{\mathbb{Z}}\\}$, where $\zeta$ is an indeterminate. Define a map
$R_{l,m}:\mathrm{Aff}(B_{l})\otimes\mathrm{Aff}(B_{m})\rightarrow\mathrm{Aff}(B_{m})\otimes\mathrm{Aff}(B_{l})$
by222Tensor product $\otimes$ in this appendix can just be regarded as product
of sets.

$\displaystyle R_{l,m}:$
$\displaystyle\;\alpha\zeta^{a}\otimes\beta\zeta^{b}\mapsto\tilde{\beta}\zeta^{b+H(\alpha\otimes\beta)}\otimes\tilde{\alpha}\zeta^{a-H(\alpha\otimes\beta)}$
(D.1)

for any $a,b\in{\mathbb{Z}}$. Here for $\alpha=(\alpha_{0},\alpha_{1})\in
B_{l}$, $\beta=(\beta_{0},\beta_{1})\in B_{m}$, the image
$\tilde{\alpha}=(\tilde{\alpha}_{0},\tilde{\alpha}_{1})\in B_{l}$,
$\tilde{\beta}=(\tilde{\beta}_{0},\tilde{\beta}_{1})\in B_{m}$ and
$H(\alpha\otimes\beta)\in{\mathbb{Z}}$ are specified by

$\displaystyle\tilde{\alpha}_{i}=\alpha_{i}+\min(\alpha_{i+1},\beta_{i})-\min(\alpha_{i},\beta_{i+1}),\;\;\tilde{\beta}_{i}=\beta_{i}-\min(\alpha_{i+1},\beta_{i})+\min(\alpha_{i},\beta_{i+1}),$
(D.2) $\displaystyle H(\alpha\otimes\beta)=\min(\alpha_{0},\beta_{1})$ (D.3)

with all the indices in ${\mathbb{Z}}_{2}$. The map $R_{l,m}$ is called an
(affine) combinatorial R. It is the quantum $R$ matrix for $(\text{spin $l/2$
rep})\otimes(\text{spin $m/2$ rep})$ of $U_{q}(\widehat{sl}_{2})$ at $q=0$
with respect to the crystal base, which retains all the combinatorial essence.
The indeterminate $\zeta$ is a remnant of the spectral parameter of the $R$
matrices. The simpler version forgetting it (the one formally corresponding to
$\zeta=1$) is called (classical) combinatorial $R$. In what follows the both
versions will simply be denoted by $R$. The relation (D.1) is customarily
depicted as

$\alpha\zeta^{a}$$\beta\zeta^{b}$$\tilde{\alpha}\zeta^{a-H}\qquad(H=H(\alpha\otimes\beta))$.$\tilde{\beta}\zeta^{b+H}$

It satisfies the inversion and the Yang-Baxter relations

$\displaystyle R_{m,l}R_{l,m}=\mathrm{Id},$ (D.4) $\displaystyle(1\otimes
R_{k,l})(R_{k,m}\otimes 1)(1\otimes R_{l,m})=(R_{l,m}\otimes 1)(1\otimes
R_{k,m})(R_{k,l}\otimes 1),$ (D.5)

which are equalities of the maps
$\mathrm{Aff}(B_{l})\otimes\mathrm{Aff}(B_{m})\rightarrow\mathrm{Aff}(B_{l})\otimes\mathrm{Aff}(B_{m})$
and
$\mathrm{Aff}(B_{k})\otimes\mathrm{Aff}(B_{l})\otimes\mathrm{Aff}(B_{m})\rightarrow\mathrm{Aff}(B_{m})\otimes\mathrm{Aff}(B_{l})\otimes\mathrm{Aff}(B_{k})$
for any $k,l,m\in{\mathbb{Z}}_{\geq 1}$. Note that these relations include the
equality of the powers of $\zeta$ for each tensor component. For example, the
inversion relation tells

$\displaystyle H(\alpha\otimes\beta)=H(\tilde{\beta}\otimes\tilde{\alpha}).$
(D.6)

Now we consider BBS. An element $\alpha\in B_{l}$ can be interpreted as a
capacity $l$ carrier containing $\alpha_{1}$ balls. When $l=1$, it may also be
regarded as a local BBS state containing $\alpha_{1}(=0,1)$ ball. The BBS on
the length $L$ periodic lattice is a dynamical system on $B_{1}^{\otimes L}$.
In what follows, a BBS state $s=(s_{1},\ldots,s_{L})\in\\{0,1\\}^{L}$ is
identified with $(1-s_{1},s_{1})\otimes\cdots\otimes(1-s_{L},s_{L})\in
B^{\otimes L}_{1}$, which will also be denoted by $s_{1}\otimes\cdots\otimes
s_{L}$.

The time evolution $s^{\prime}=T_{l}(s)$ and the associated $l^{\text{th}}$
energy $E_{l}(s)$ are defined by the composition of the combinatorial $R$ as
follows:

$s_{1}\quad
s_{2}$$s^{\prime}_{1}\quad\,s^{\prime}_{2}$$u\zeta^{0}$.........$s_{L}$$s^{\prime}_{L}$$u\zeta^{-E_{l}(s)}$

Here the carrier $u\in B_{l}$ is determined uniquely from
$s=s_{1}\otimes\cdots\otimes s_{L}$ by the periodic boundary condition,
namely, by requiring that it comes back to $u$ itself after penetrating $s$
provided that the ball density is not exactly $1/2$.

Let $u^{\prime}\in B_{i}$ be another carrier for the time evolution
$T_{l}(s)\rightarrow T_{i}T_{l}(s)$. Set $R:u^{\prime}\zeta^{0}\otimes
u\zeta^{0}\mapsto{\tilde{u}}\zeta^{h}\otimes{\tilde{u}^{\prime}}\zeta^{-h}$
where $h=H(u^{\prime}\otimes u)$. From the Yang-Baxter relation we have

$\cdots$$\cdots$$\scriptstyle u$$\scriptstyle u^{\prime}$$\scriptstyle
s$$\scriptstyle T_{l}(s)$$\scriptstyle T_{i}T_{l}(s)$$\scriptstyle
u\zeta^{-E_{l}(s)}$$\scriptstyle
u^{\prime}\zeta^{-E_{i}(T_{l}(s))}$$\scriptstyle\tilde{u}\zeta^{-E_{l}(s)+h}$$\scriptstyle{\tilde{u}^{\prime}}\zeta^{-E_{i}(T_{l}(s))-h}$$=$$\scriptstyle
u^{\prime}$$\scriptstyle
u$$\scriptstyle{\tilde{u}}\zeta^{h}$$\scriptstyle{\tilde{u}^{\prime}}\zeta^{-h}$$\cdots$$\cdots$$\scriptstyle
s$$\scriptstyle T_{i}(s)$$\scriptstyle
T_{l}T_{i}(s)$$\scriptstyle{\tilde{u}^{\prime}}\zeta^{-E_{i}(s)-h}$$\scriptstyle{\tilde{u}}\zeta^{-E_{l}(T_{i}(s))+h}$

Comparing the two sides we get the commutativity $T_{l}T_{i}(s)=T_{i}T_{l}(s)$
and the energy conservation $E_{i}(s)=E_{i}(T_{l}(s))$,
$E_{l}(s)=E_{l}(T_{i}(s))$ as is well known.

Now consider the intermediate stage where the carriers have only gone through
the first $x-1$ local states $s_{<x}=s_{1}\otimes\cdots\otimes s_{x-1}$. To
systematize the notation we write $u,u^{\prime},h$ in the above diagram as
$u(0),u^{\prime}(0),h(0)$. Then the corresponding diagram looks as

$\cdots$$\cdots$$\scriptstyle B_{l}\,\ni\,u(0)$$\scriptstyle
B_{i}\,\ni\,u^{\prime}(0)$$\scriptstyle s_{<x}$$\scriptstyle
T_{l}(s_{<x})$$\scriptstyle T_{i}T_{l}(s_{<x})$$\scriptstyle
u(x)\zeta^{-E_{l}(s_{<x})}$$\scriptstyle
u^{\prime}(x)\zeta^{-E_{i}(T_{l}(s_{<x}))}$$\scriptstyle\tilde{u}(x)\zeta^{-E_{l}(s_{<x})+h(x)}$$\scriptstyle{\tilde{u}^{\prime}}(x)\zeta^{-E_{i}(T_{l}(s_{<x}))-h(x)}$$=$$\scriptstyle
B_{l}\,\ni\,u(0)$$\scriptstyle
B_{i}\,\ni\,u^{\prime}(0)$$\scriptstyle{\tilde{u}}(0)\zeta^{h(0)}$$\scriptstyle{\tilde{u}^{\prime}}(0)\zeta^{-h(0)}$$\cdots$$\cdots$$\scriptstyle
s_{<x}$$\scriptstyle T_{i}(s_{<x})$$\scriptstyle
T_{l}T_{i}(s_{<x})$$\scriptstyle{\tilde{u}^{\prime}}(x)\zeta^{-E_{i}(s_{<x})-h(0)}$$\scriptstyle{\tilde{u}}(x)\zeta^{-E_{l}(T_{i}(s_{<x}))+h(0)}$

where $h(x)=H(u^{\prime}(x)\otimes u(x))$. The carriers $u(x)$ and
$u^{\prime}(x)$ (resp. ${\tilde{u}}(x)$ and ${\tilde{u}^{\prime}}(x)$) at this
position do not yet have to return to the initial ones $u(0)$ and
$u^{\prime}(0)$ (resp. ${\tilde{u}}(0)$ and ${\tilde{u}^{\prime}}(0)$). Define
the local observables

$\displaystyle\hat{\eta}^{(l)}_{i}(x)=H(u^{\prime}(x)\otimes
u(x)),\qquad\hat{\eta}^{(i)}_{l}(x)=H({\tilde{u}}(x)\otimes{\tilde{u}^{\prime}}(x)).$
(D.7)

Then the following properties are satisfied.

(i) $\displaystyle\hat{\eta}^{(l)}_{i}(x)=\hat{\eta}^{(i)}_{l}(x),$ (D.8) (ii)
$\displaystyle\begin{cases}E_{l}(T_{i}(s_{<x}))-E_{l}(s_{<x})+{\hat{\eta}}^{(i)}_{l}(x)-{\hat{\eta}}^{(i)}_{l}(0)=0,\\\
E_{i}(T_{l}(s_{<x}))-E_{i}(s_{<x})+{\hat{\eta}}^{(l)}_{i}(x)-{\hat{\eta}}^{(l)}_{i}(0)=0.\end{cases}$
(D.9)

In fact (i) is a consequence of (D.6). As for (ii) both relations follow
simultaneously by comparing the powers of $\zeta$ in the above diagram using
(i). The upper relation, for instance, is nothing but the space
$[0,x]$-integrated equation of continuity for $E_{l}$ with respect to the time
evolution $T_{i}$, where $\hat{\eta}^{(i)}_{l}(x)$ plays the role of local
current at $x$.

As an example, $\hat{\eta}^{(3)}_{2}(x)$ associated with $T_{2}T_{3}$ are
given in red letters in the top panel of Fig. 3. On the other hand,
$\hat{\eta}^{(2)}_{3}(x)$ associated with $T_{3}T_{2}$ are given similarly in
the bottom panel of Fig. 3. One can observe the equality
$\hat{\eta}^{(3)}_{2}(x)=\hat{\eta}^{(2)}_{3}(x)$ everywhere. In this example,
$\hat{\eta}^{(2)}_{3}(x)=\hat{\eta}^{(2)}_{\infty}(x)$ holds and it coincides
with the carriers for $T_{2}$, which is indeed the ball current.

## Bibliography

## References

  * [1] Takahashi D and Satsuma J 1990 A Soliton Cellular Automaton J. Phys. Soc. Jpn. 59 3514–3519 URL https://journals.jps.jp/doi/10.1143/JPSJ.59.3514
  * [2] Inoue R, Kuniba A and Takagi T 2012 Integrable structure of box–ball systems: crystal, Bethe ansatz, ultradiscretization and tropical geometry J. Phys. A: Math. Theor. 45 073001 URL https://doi.org/10.1088%2F1751-8113%2F45%2F7%2F073001
  * [3] Levine L, Lyu H and Pike J 2020 Double jump phase transition in a soliton cellular automaton arXiv:1706.05621 URL http://arxiv.org/abs/1706.05621
  * [4] Croydon D A, Kato T, Sasada M and Tsujimoto S 2018 Dynamics of the box-ball system with random initial conditions via Pitman’s transformation arXiv:1806.02147 URL http://arxiv.org/abs/1806.02147
  * [5] Kuniba A, Lyu H and Okado M 2018 Randomized box–ball systems, limit shape of rigged configurations and thermodynamic Bethe ansatz Nucl. Phys. B 937 240–271 URL http://www.sciencedirect.com/science/article/pii/S055032131830289X
  * [6] Ferrari P A and Gabrielli D 2020 BBS invariant measures with independent soliton components Electron. J. Probab. 25 URL https://projecteuclid.org/euclid.ejp/1594432886
  * [7] Ferrari P A, Nguyen C, Rolla L T and Wang M 2020 Soliton decomposition of the Box-Ball System arXiv:1806.02798 URL http://arxiv.org/abs/1806.02798
  * [8] Lewis J, Lyu H, Pylyavskyy P and Sen A 2020 Scaling limit of soliton lengths in a multicolor box-ball system arXiv:1911.04458 URL http://arxiv.org/abs/1911.04458
  * [9] Croydon D A and Sasada M 2019 Invariant measures for the box-ball system based on stationary Markov chains and periodic Gibbs measures J. Math. Phys. 60 083301 URL https://aip.scitation.org/doi/10.1063/1.5095622
  * [10] Croydon D A and Sasada M 2021 Generalized Hydrodynamic Limit for the Box–Ball System Commun. Math. Phys. 383 427–463 URL https://doi.org/10.1007/s00220-020-03914-x
  * [11] Kuniba A and Lyu H 2020 Large Deviations and One-Sided Scaling Limit of Randomized Multicolor Box-Ball System J. Stat. Phys. 178 38–74 URL https://doi.org/10.1007/s10955-019-02417-x
  * [12] Kuniba A, Misguich G and Pasquier V 2020 Generalized hydrodynamics in box-ball system J. Phys. A: Math. Theor. 53 404001 URL https://doi.org/10.1088%2F1751-8121%2Fabadb9
  * [13] Castro-Alvaredo O A, Doyon B and Yoshimura T 2016 Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium Phys. Rev. X 6 041065 URL https://link.aps.org/doi/10.1103/PhysRevX.6.041065
  * [14] Bertini B, Collura M, De Nardis J and Fagotti M 2016 Transport in Out-of-Equilibrium XXZ Chains: Exact Profiles of Charges and Currents Phys. Rev. Lett. 117 207201 URL https://link.aps.org/doi/10.1103/PhysRevLett.117.207201
  * [15] Doyon B 2020 Lecture notes on Generalised Hydrodynamics SciPost Physics Lecture Notes 018 URL https://scipost.org/10.21468/SciPostPhysLectNotes.18
  * [16] Rigol M, Dunjko V, Yurovsky V and Olshanii M 2007 Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons Phys. Rev. Lett. 98 050405 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.050405
  * [17] Doyon B and Spohn H 2017 Dynamics of hard rods with initial domain wall state J. Stat. Mech. 2017 073210 URL https://doi.org/10.1088%2F1742-5468%2Faa7abf
  * [18] Doyon B, Yoshimura T and Caux J S 2018 Soliton Gases and Generalized Hydrodynamics Phys. Rev. Lett. 120 045301 URL https://link.aps.org/doi/10.1103/PhysRevLett.120.045301
  * [19] Doyon B 2019 Generalized hydrodynamics of the classical Toda system J. Math. Phys. 60 073302 URL https://aip.scitation.org/doi/10.1063/1.5096892
  * [20] Bastianello A, Doyon B, Watts G and Yoshimura T 2018 Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model SciPost Physics 4 045 URL https://scipost.org/10.21468/SciPostPhys.4.6.045
  * [21] Kuniba A, Misguich G and Pasquier V 2021 Generalized hydrodynamics in complete box-ball system for $U_{q}(\widehat{sl_{n}})$ SciPost Physics 10 095 URL https://www.scipost.org/SciPostPhys.10.4.095
  * [22] Touchette H 2009 The large deviation approach to statistical mechanics Physics Reports 478 1–69 URL https://www.sciencedirect.com/science/article/pii/S0370157309001410
  * [23] Yang C N and Yang C P 1969 Thermodynamics of a One-Dimensional System of Bosons with Repulsive Delta-Function Interaction J. Math. Phys. 10 1115–1122 URL https://aip.scitation.org/doi/10.1063/1.1664947
  * [24] Zakharov V E 1971 Kinetic Equation for Solitons JETP 33 538 URL http://www.jetp.ac.ru/cgi-bin/e/index/e/33/3/p538?a=list
  * [25] El G A 2003 The thermodynamic limit of the Whitham equations Phys. Lett. A 311 374–383 URL http://www.sciencedirect.com/science/article/pii/S0375960103005152
  * [26] El G A and Kamchatnov A M 2005 Kinetic Equation for a Dense Soliton Gas Phys. Rev. Lett. 95 204101 URL https://link.aps.org/doi/10.1103/PhysRevLett.95.204101
  * [27] Ferrari P A, Nguyen C, Rolla L T and Wang M 2021 Soliton decomposition of the Box-Ball System arXiv:1806.02798 URL http://arxiv.org/abs/1806.02798
  * [28] Myers J, Bhaseen J, Harris R J and Doyon B 2020 Transport fluctuations in integrable models out of equilibrium SciPost Physics 8 007 URL https://scipost.org/SciPostPhys.8.1.007
  * [29] Doyon B and Myers J 2020 Fluctuations in Ballistic Transport from Euler Hydrodynamics Ann. Henri Poincaré 21 255–302 URL https://doi.org/10.1007/s00023-019-00860-w
  * [30] Mendl C B and Spohn H 2015 Current fluctuations for anharmonic chains in thermal equilibrium J. Stat. Mech. Theor. Exp. 2015 P03007 URL https://doi.org/10.1088/1742-5468/2015/03/p03007
  * [31] Ilievski E and De Nardis J 2017 Microscopic Origin of Ideal Conductivity in Integrable Quantum Models Phys. Rev. Lett. 119 020602 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.020602
  * [32] Ilievski E and De Nardis J 2017 Ballistic transport in the one-dimensional Hubbard model: The hydrodynamic approach Phys. Rev. B 96 081118 URL https://link.aps.org/doi/10.1103/PhysRevB.96.081118
  * [33] Doyon B and Spohn H 2017 Drude Weight for the Lieb-Liniger Bose Gas SciPost Physics 3 039 URL https://www.scipost.org/10.21468/SciPostPhys.3.6.039
  * [34] Bulchandani V B, Vasseur R, Karrasch C and Moore J E 2018 Bethe-Boltzmann hydrodynamics and spin transport in the XXZ chain Phys. Rev. B 97 045407 URL https://link.aps.org/doi/10.1103/PhysRevB.97.045407
  * [35] Krajnik Z, Ilievski E, Prosen T and Pasquier V 2021 Anisotropic Landau-Lifshitz model in discrete space-time SciPost Physics 11 051 URL https://scipost.org/SciPostPhys.11.3.051
  * [36] Nardis J D, Doyon B, Medenjak M and Panfil M 2022 Correlation functions and transport coefficients in generalised hydrodynamics J. Stat. Mech. 2022 014002 URL https://doi.org/10.1088/1742-5468/ac3658
  * [37] Spohn H 2018 Interacting and noninteracting integrable systems J. Math. Phys. 59 091402 URL https://aip.scitation.org/doi/10.1063/1.5018624
  * [38] Krajnik Z, Schmidt J, Pasquier V, Ilievski E and Prosen T 2022 Exact anomalous current fluctuations in a deterministic interacting model arXiv:2201.05126 (to appear in Phys. Rev. Lett.) URL http://arxiv.org/abs/2201.05126
  * [39] Krajnik Z, Ilievski E and Prosen T 2022 Absence of Normal Fluctuations in an Integrable Magnet Phys. Rev. Lett. 128 090604 URL https://link.aps.org/doi/10.1103/PhysRevLett.128.090604