# Exact Solution of Bipartite Fluctuations in One-Dimensional Fermions

Kazuya Fujimoto  Department of Physics, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo 152-8551, Japan  Tomohiro Sasamoto  Department of
Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8551, Japan

###### Abstract

Emergence of hydrodynamics in quantum many-body systems has recently garnered
growing interest. The recent experiment of ultracold atoms [J. F. Wienand et
al., arXiv:2306.11457] studied emergent hydrodynamics in hard-core bosons
using a bipartite fluctuation, which quantifies how the particle number
fluctuates in a subsystem. In this Letter, we theoretically study the variance
of a bipartite fluctuation in one-dimensional noninteracting fermionic
dynamics starting from an alternating state, deriving the exact solution of
the variance and its asymptotic linear growth law for the long-time dynamics.
To compare the theoretical prediction with the experiment, we generalize our
exact solution by incorporating the incompleteness of the initial alternating
state, deriving the general linear growth law analytically. We find that it
shows quantitative agreement with the experimentally observed variance growth
without any fitting parameters.

Introduction.— Relaxation to an equilibrium state is ubiquitous in quantum
many-body dynamics, bringing up fundamental and intriguing questions, e.g.,
how an isolated quantum system relaxes to a thermal equilibrium state
Polkovnikov  _et al._ (2011); Eisert  _et al._ (2015); Nandkishore and Huse
(2015); Torres-Herrera  _et al._ (2015); Luca D’Alessio and Rigol (2016); Mori
_et al._ (2018); Abanin  _et al._ (2019); Gong and Hamazaki (2022). Over
decades, such quantum thermalization has been intensively studied from broad
perspectives, such as the eigenstate thermalization hypothesis Deutsch (1991);
Srednicki (1994), integrability Kinoshita  _et al._ (2006); Gring  _et al._
(2012); Cazalilla (2006); Rigol  _et al._ (2006, 2007); Calabrese  _et al._
(2011); Fagotti and Essler (2013); Langen  _et al._ (2015); Ilievski  _et al._
(2015); Tang  _et al._ (2018), many-body localization Basko  _et al._ (2006);
Žnidarič  _et al._ (2008); Pal and Huse (2010); Bardarson  _et al._ (2012);
Serbyn  _et al._ (2013); Huse  _et al._ (2014); Schreiber  _et al._ (2015);
Potter  _et al._ (2015); Vasseur and Moore (2016); yoon Choi  _et al._ (2016);
Lüschen  _et al._ (2017); Kohlert  _et al._ (2019), and quantum scars Bernien
_et al._ (2017); Shiraishi and Mori (2017); Turner  _et al._ (2018a, b); Lin
and Motrunich (2019); Shibata  _et al._ (2020). To date, many experiments
involving ultracold atoms and trapped ions have observed various phenomena
related to quantum thermalization Kinoshita  _et al._ (2006); Gring  _et al._
(2012); Trotzky  _et al._ (2012); Langen  _et al._ (2013); Schreiber  _et al._
(2015); Langen  _et al._ (2015); Clos  _et al._ (2016); Kaufman  _et al._
(2016); Neill  _et al._ (2016); Lüschen  _et al._ (2017); Bernien  _et al._
(2017); Tang  _et al._ (2018); Kohlert  _et al._ (2019). For example, S.
Trotzky et al. realized an isolated bosonic system on a one-dimensional
lattice, observing the relaxation dynamics starting from an alternating state
where the bosons occupy every other site Trotzky  _et al._ (2012).

Recently, beyond the conventional quantum thermalization, hydrodynamic
description based on local equilibrium has rapidly developed in quantum many-
body systems. This situation is exemplified by the recent observation of
electron fluids in strongly interacting systems Moll  _et al._ (2016);
Bandurin  _et al._ (2016); Gooth  _et al._ (2018); Krishna Kumar  _et al._
(2017); Sulpizio  _et al._ (2019); Aharon-Steinberg  _et al._ (2022) and the
success of generalized hydrodynamics being valid for integrable quantum models
Castro-Alvaredo  _et al._ (2016); Bertini  _et al._ (2016); Bulchandani  _et
al._ (2017, 2018); Doyon and Yoshimura (2017); Doyon  _et al._ (2017, 2018);
Collura  _et al._ (2018); De Nardis  _et al._ (2018); Schemmer  _et al._
(2019); Gopalakrishnan and Vasseur (2019); Doyon (2020); Alba  _et al._
(2021); Malvania  _et al._ (2021); Bouchoule and Dubail (2022); Essler (2023).
When considering such hydrodynamic description, it is of essence to scrutinize
how a system approaches a local equilibrium state. One of the most useful
quantities for diagnosing local equilibrium is a bipartite fluctuation that
quantifies how a physical quantity in a subsystem fluctuates. In a recent
experiment of ultracold atoms for hard-core bosons Wienand  _et al._ (2023),
J. F. Wienand et al. utilized the bipartite fluctuation of the particle number
in quench dynamics starting from the alternating state, studying the local
equilibrium and the emergence of fluctuating hydrodynamics in the quantum
many-body system.

Despite the recent strong interest in emergent hydrodynamics in a quantum
regime, exact and analytical results of the bipartite fluctuation with the
initial alternating state have been elusive even for the noninteracting
fermions Groha  _et al._ (2018); com , although there are several numerical
works Fujimoto  _et al._ (2020); Jin  _et al._ (2020); Fujimoto  _et al._
(2021, 2022); Cecile  _et al._ (2024); Aditya and Roy (2024); Bhakuni and Lev
(2023). Hence, its fundamental and analytical understanding is highly
desirable now, and it is of great importance to explain the recently observed
bipartite fluctuation Wienand  _et al._ (2023) quantitatively via analytical
solutions.

Figure 1:  Schematic illustration of the main result. (a) One-dimensional
fermions addressed in this work. The blue circles represent fermions occupying
every other site at the initial time. This initial state is referred to as an
alternating state. Under unitary time evolution, the fermions hop to
neighboring sites. The quantity of our interest is the variance of the
particle number in a subsystem $A$ for a time-evolved state. (b) Schematic for
the variance growth. The abscissa and ordinate, respectively, denote time $t$
and a subtracted variance, which is defined by the variance from which one
subtracts its initial value. The two solid lines represent our analytical
results, namely $2t/\pi$ and $2(n_{\rm even}-n_{\rm odd})^{2}t/\pi$, obtained
by considering the variances growing from the complete and incomplete
alternating states. Here, $n_{\rm even}$ and $n_{\rm odd}$ are, respectively,
averaged initial particle numbers at even and odd sites, quantifying the
incompleteness of the initial alternating state.

In this Letter, we consider one-dimensional fermions starting from the
alternating state, theoretically studying time evolution for the variance of
the particle numbers in a subsystem as schematically shown in Fig. 1 (a) and
compare our theoretical prediction with the recent experimental data of Ref.
Wienand  _et al._ (2023). First, in the noninteracting fermions, we exactly
compute infinite series involving the $n$th-order Bessel function $J_{n}(x)$
of the first kind, obtaining the exact and simple expression of the variance.
We apply asymptotic analysis to the exact result, analytically deriving the
linear growth law of the variance and its sub-leading correction for the long-
time dynamics. Second, we compare our analytical result with the recent
experiment of Ref. Wienand  _et al._ (2023). For this purpose, we focus on the
incompleteness of the initial alternating state realized in the experiment. We
incorporate the incompleteness into our exact solution, deriving the general
linear growth law of the variance. This general law can quantitatively explain
the variance growth observed in Ref. Wienand  _et al._ (2023) without any
fitting parameters. The impact of the incompleteness is schematically
displayed in Fig. 1 (b). Finally, we numerically investigate interaction
effects on the variance growth using the time-evolving decimation (TEBD)
method Vidal (2003, 2004); Schollwöck (2011); Paeckel  _et al._ (2019),
finding that the presence of the interaction substantially alters the behavior
of the variance.

Before ending the introduction, it is worth noting that a bipartite
fluctuation has been studied in terms of an integrated current and full
counting statistics. In the former case, most of previous works consider
dynamics starting from a domain-wall initial state, investigating a bipartite
fluctuation for a conserved quantity by connecting it to an integrated current
via a continuity equation Antal  _et al._ (2008); Eisler and Rácz (2013);
Ljubotina  _et al._ ; Moriya  _et al._ (2019); Gamayun  _et al._ (2020); Jin
_et al._ (2021); Gopalakrishnan and Vasseur (2023); Wei  _et al._ (2022);
Krajnik  _et al._ (2022a, b, 2024a, 2024b). One of the notable results is the
logarithmic growth of the variance in noninteracting fermions Antal  _et al._
(2008). In the latter case, a bipartite fluctuation has been employed to study
static and dynamic properties of entanglement entropy Klich and Levitov
(2009); Song  _et al._ (2010, 2012); Rachel  _et al._ (2012); Parez  _et al._
(2021); Oshima and Fuji (2023); Bertini  _et al._ (2023a) and fluctuations of
physical quantities Schönhammer (2007); Collura  _et al._ (2017); Stéphan and
Pollmann (2017); Najafi and Rajabpour (2017); Arzamasovs and Gangardt (2019);
Calabrese  _et al._ (2020); Ares  _et al._ (2021); Smith  _et al._ (2021);
McCulloch  _et al._ (2023); Hercé  _et al._ (2023); Bertini  _et al._ (2023b).

Setup.— Let us consider fermions on a one-dimensional lattice labeled by
$\Lambda\coloneqq\\{-L,-L+1,...,L\\}$ with a positive even number $L$. We
denote the fermionic creation and annihilation operators at a site
$j\in\Lambda$ by $\hat{a}_{j}^{\dagger}$ and $\hat{a}_{j}$. Then, the
Hamiltonian is given by

$\displaystyle\hat{H}\coloneqq-\sum_{j=-L}^{L-1}\left(\hat{a}_{j+1}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{j+1}\right)+U\sum_{j=-L}^{L-1}\hat{n}_{j+1}\hat{n}_{j}$
(1)

with the particle-number operator
$\hat{n}_{j}\coloneqq\hat{a}_{j}^{\dagger}\hat{a}_{j}$ and the interaction
parameter $U$. Here, the boundary condition is basically a periodic boundary
condition ($\hat{a}_{L}=\hat{a}_{-L}$), but we use an open boundary condition
($\hat{a}_{L}=0$) in several cases. In the latter case, we will explicitly
specify the boundary condition. The initial state is the alternating state
defined by $\ket{\psi_{\rm
alt}}\coloneqq\prod_{j=-L/2}^{L/2-1}\hat{a}_{2j+1}^{\dagger}\ket{0}$. We
denote a quantum state at time $t$ by $\ket{\psi(t)}$ and assume that the
state obeys the Schrödinger equation $\displaystyle{\rm
i}d\ket{\psi(t)}/dt=\hat{H}\ket{\psi(t)}$. Here, the Dirac constant $\hbar$ is
set to be unity.

In this Letter, we shall study the fluctuation of the particle number in a
subsystem. The operator for the bipartite particle-number is defined by
$\hat{N}_{M}\coloneqq\sum_{j=0}^{M-1}\hat{a}_{j}^{\dagger}\hat{a}_{j}$ with a
positive integer $M$. Then, a generating function for the $n$th moment of
$\hat{N}_{M}$ is given by
$G_{M}(\lambda,t)\coloneqq\braket{e^{\lambda\hat{N}_{M}}}_{t}$, where we
introduce the notation $\braket{\bullet}_{t}\coloneqq{\rm
Tr}[\hat{\rho}(t)\bullet]$ with the density matrix
$\hat{\rho}(t)\coloneqq\ket{\psi(t)}\bra{\psi(t)}$. We can compute the moments
by differentiating $G_{M}(\lambda,t)$ with respect to $\lambda$.

Exact solution for the variance of the noninteracting fermions.— We study the
time evolution for the variance of $\hat{N}_{M}$ in the noninteracting
fermions ($U=0$) under the thermodynamic limit ($L\rightarrow\infty$). As
described in Sec. I of the Supplemental Material SM , the generating function
$G_{M}(\lambda,t)$ becomes Eisler and Rácz (2013); Schönhammer (2007); Parez
_et al._ (2021)

$\displaystyle G_{M}(\lambda,t)={\rm
det}\left[\delta_{j,k}+(e^{\lambda}-1)\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}\right]_{j,k=0}^{M-1}.$
(2)

Here, the two-point correlator $\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}$
in the thermodynamic limit ($L\rightarrow\infty$) is given by Flesch  _et al._
(2008)

$\displaystyle\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}=\frac{\delta_{j,k}}{2}-\frac{{\rm
i}^{j+k}}{2}J_{k-j}(4t).$ (3)

We here derive the exact and simple expression for the variance of
$\hat{N}_{M}$ under the limit $M\rightarrow\infty$. The first step is to
express the variance via the Bessel function $J_{n}(x)$ of the first kind.
Differentiating Eq. (2) with respect to $\lambda$, we obtain the variance
$\sigma_{M}(t)^{2}\coloneqq\braket{\hat{N}_{M}^{2}}_{t}-\braket{\hat{N}_{M}}_{t}^{2}=\partial^{2}G_{M}(\lambda,t)/\partial\lambda^{2}|_{\lambda=0}-(\partial
G_{M}(\lambda,t)/\partial\lambda|_{\lambda=0})^{2}$ as

$\displaystyle\sigma_{M}(t)^{2}=\frac{M}{4}\left(1-J_{0}(4t)^{2}-2\sum_{k=1}^{M-1}J_{k}(4t)^{2}\right)+\frac{1}{2}\sum_{k=1}^{M-1}kJ_{k}(4t)^{2}.$
(4)

The detailed derivation of Eq. (4) is given in Sec. II of the Supplemental
Material SM . The next task is to take the limit $M\rightarrow\infty$ in Eq.
(4). As proved in Sec. III of the Supplemental Material SM , we have
$\lim_{M\rightarrow\infty}M\left(1-J_{0}(4t)^{2}-2\sum_{k=1}^{M-1}J_{k}(4t)^{2}\right)/4=0$
for $t>0$ and
$\lim_{M\rightarrow\infty}\sum_{k=1}^{M-1}kJ_{k}(4t)^{2}/2=4t^{2}\left(J_{0}(4t)^{2}+J_{1}(4t)^{2}\right)-tJ_{0}(4t)J_{1}(4t)$.
We utilize these formulae and Eq. (4), deriving the following exact and simple
expression of the variance
$\sigma(t)^{2}\coloneqq\lim_{M\rightarrow\infty}\sigma_{M}(t)^{2}$,

$\displaystyle\sigma(t)^{2}=4t^{2}\left(J_{0}(4t)^{2}+J_{1}(4t)^{2}\right)-tJ_{0}(4t)J_{1}(4t)$
(5)

for $t>0$. Note that Eq. (5) is still valid at $t=0$ since we can show
$\sigma_{M}(0)^{2}=0$ and $\sigma(0)^{2}=0$ from Eqs. (4) and (5),
respectively. Thus, we can relax the constraint $t>0$ for Eq. (5) to $t\geq
0$.

Figure 2:  Numerical verification of Eq. (6). The time evolution of the
variance denoted by the blue square symbols is numerically obtained using Eqs.
(2) and (3) with $M=200$. In the upper and lower panels, we compare Eq. (6)
with the numerical result.

We apply asymptotic analysis to Eq. (5), deriving the asymptotic forms of
$\sigma(t)^{2}$ both for the short-time ($t\ll 1$) and long-time ($t\gg 1$)
dynamics. Using the asymptotic forms of the Bessel functions of the first kind
mat (2015), we obtain

$\displaystyle\sigma(t)^{2}\simeq\begin{dcases}2t^{2}&(t\ll 1)\\\
\frac{2}{\pi}t-\frac{1}{64\pi t}\left(2\sin(8t)+1\right)&(t\gg
1).\end{dcases}$ (6)

The detailed derivation of Eq. (6) is given in Sec. IV of the Supplemental
Material SM . The essential consequence of Eq. (6) is that the variance is
proportional to time $t$ for $t\gg 1$. This linear growth law is entirely
different from the variance growth in noninteracting fermionic dynamics
starting from the domain-wall state $\ket{\psi_{\rm
DW}}\coloneqq\prod_{j\in\mathbb{Z}_{<0}}\hat{a}_{j}^{\dagger}\ket{0}$, where
the variance was shown to obey the logarithmic growth Antal  _et al._ (2008);
Moriya  _et al._ (2019). Note that the linear growth law was numerically
confirmed in Ref. Fujimoto  _et al._ (2022) and Eq. (6) gives its analytical
explanation.

Finally, we numerically verify Eq. (6). Figure 2 shows the time evolution of
the variance obtained by Eqs. (2) and (3) numerically. One can see the
excellent agreement with Eq. (6). Note that the linear growth of the variance
appears in $t\gtrsim 0.4$, though it is originally derived under $t\gg 1$.
This behavior is favorable for observing the linear growth experimentally.

Comparison with the experimental result of Ref. Wienand  _et al._ (2023).— We
consider whether our theoretical prediction of the variance for the
noninteracting fermions can explain the experimental result reported in Ref.
Wienand  _et al._ (2023). For this purpose, we must note that the experiment
does not realize the complete alternating initial state $\ket{\psi_{\rm
alt}}$. To make this point clear, let us denote averaged particle numbers at
the even and odd sites by $n_{\rm even}$ and $n_{\rm odd}$. The observed
imbalance parameter $I\coloneqq(n_{\rm odd}-n_{\rm even})/(n_{\rm even}+n_{\rm
odd})$ is about $0.91$ Wienand  _et al._ (2023), which means the existence of
deviation from the complete alternating state since $\ket{\psi_{\rm alt}}$ has
$I=1$. This observation strongly suggests that we need to incorporate the
incompleteness of the initial alternating state into the analytical results of
Eqs. (5) and (6). In what follows, we describe our analytical solution with
the incomplete alternating state, comparing it with the experimental data.

Before calculating the variance, we first specify an initial state for the
incomplete alternating state by considering how the alternating state is
experimentally prepared. The experiment prepares the initial state by ramping
up a strongly tilted double-well potential in one direction (see Sec. I. B. of
the supplementary material of Ref. Wienand  _et al._ (2023)). Thus, we expect
that the experimental initial state is well described by a product state,
assuming that the initial density matrix is approximated by

$\displaystyle\hat{\rho}_{\rm
alt}\coloneqq\dfrac{1}{\Xi}\exp\left(-\beta\hat{H}_{\rm
alt}+\beta\mu\hat{N}_{\rm tot}\right)$ (7)

with the inverse temperature $\beta$, the chemical potential $\mu$, and the
partition function $\Xi\coloneqq{\rm Tr}\exp(-\beta\hat{H}_{\rm
alt}+\beta\mu\hat{N}_{\rm tot})$. The operators $\hat{H}_{\rm alt}$ and
$\hat{N}_{\rm tot}$ are defined by $\hat{H}_{\rm
alt}\coloneqq\sum_{j=-L/2}^{L/2-1}(\hat{a}_{2j}^{\dagger}\hat{a}_{2j}-\hat{a}_{2j+1}^{\dagger}\hat{a}_{2j+1})$
and $\hat{N}_{\rm
tot}\coloneqq\sum_{j=-L}^{L-1}\hat{a}_{j}^{\dagger}\hat{a}_{j}$. The averaged
particle numbers at the even and odd sites become $n_{\rm
even}=1/[e^{\beta(1-\mu)}+1]$ and $n_{\rm odd}=1/[e^{\beta(-1-\mu)}+1]$. The
parameters $\beta$ and $\mu$ are determined by the filling factor
$\nu\coloneqq(n_{\rm even}+n_{\rm odd})/2$ and the imbalance parameter $I$,
both of which are observable in the experiment.

We analytically derive the exact and asymptotic forms of the variance for the
dynamics starting from the incomplete alternating state of Eq. (7). As
detailed in Sec. V of the Supplemental Material SM , the two-point correlator
becomes

$\displaystyle\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}=\frac{1}{2}\left(n_{\rm
even}+n_{\rm odd}\right)\delta_{j,k}+\frac{1}{2}\left(n_{\rm even}-n_{\rm
odd}\right){\rm i}^{j+k}J_{k-j}(4t)$ (8)

in the thermodynamic limit ($L\rightarrow\infty$). Here the density matrix is
given by $\hat{\rho}(t)=e^{-{\rm i}\hat{H}t}\hat{\rho}_{\rm alt}e^{{\rm
i}\hat{H}t}$ with Eq. (1) and $U=0$. Following the derivation given in Sec. VI
of the Supplemental Material SM , we obtain

$\displaystyle\sigma_{\rm
sub}(t)^{2}\coloneqq\lim_{M\rightarrow\infty}\left(\sigma_{M}(t)^{2}-\delta\sigma_{M}(t)^{2}\right)$
$\displaystyle=\left(n_{\rm even}-n_{\rm
odd}\right)^{2}\Bigl{[}4t^{2}\left(J_{0}(4t)^{2}+J_{1}(4t)^{2}\right)-tJ_{0}(4t)J_{1}(4t)\Bigl{]}$
(9)

with the function $\delta\sigma_{M}(t)^{2}\coloneqq M[n_{\rm even}(1-n_{\rm
even})+n_{\rm odd}(1-n_{\rm odd})]/2+J_{0}(4t)[n_{\rm even}(1-n_{\rm
even})-n_{\rm odd}(1-n_{\rm odd})][\sum_{m=0}^{M-1}(-1)^{m}]/2$. The time-
dependent term of $\delta\sigma_{M}(t)^{2}$ is much smaller than the time-
independent term when $M$ is large. Hence, we have
$\delta\sigma_{M}(t)^{2}\simeq\sigma_{M}(0)^{2}$ for $M\gg 1$, and thus
$\sigma_{\rm sub}(t)^{2}$ can be interpreted as the variance from which one
subtracts its initial value. Applying the asymptotic analysis used in Eq. (6)
to Eq. (9), we derive

$\displaystyle\sigma_{\rm sub}(t)^{2}\simeq\begin{dcases}2\left(n_{\rm
even}-n_{\rm odd}\right)^{2}t^{2}&(t\ll 1)\\\ \dfrac{2}{\pi}\left(n_{\rm
even}-n_{\rm odd}\right)^{2}t&(t\gg 1).\end{dcases}$ (10)

This result elucidates that the incompleteness of the alternating state
substantially affects the coefficient of the variance, but the exponents of
time are robust against it.

In addition to Eq. (10), we shall compare the numerical calculations in a
finite system with the experiment. Our numerical calculation is implemented
for $2L=40$, $\nu=0.44$, and $I=0.91$ with the initial state of Eq. (7). These
concrete values are taken from Ref. Wienand  _et al._ (2023) (see Sec. II. C.
of the supplementary material of Ref. Wienand  _et al._ (2023)). We adopt the
Hamiltonian defined by

$\displaystyle\hat{H}^{\prime}\coloneqq-\sum_{j=-L}^{L-1}\left(\hat{a}_{j+1}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{j+1}\right)+\sum_{j=-L}^{L-1}V_{j}\hat{n}_{j}$
(11)

with the open boundary condition ($\hat{a}_{L}=0$). Following the numerical
simulation in Ref. Wienand  _et al._ (2023), we here add a random potential
$V_{j}$, which is sampled from a uniform distribution with the range
$[-\Delta,\Delta]$. Here, $\Delta\geq 0$ denotes the strength of the
randomness. Under this setup, we investigate the variance
$\sigma^{\prime}_{M}(t)^{2}\coloneqq\overline{\braket{\hat{N}_{M}^{\prime
2}}}-\left(~{}\overline{\braket{\hat{N}_{M}^{\prime}}}~{}\right)^{2}$ with
$\hat{N}_{M}^{\prime}\coloneqq\sum_{m=-M/2}^{M/2-1}\hat{a}^{\dagger}_{m}\hat{a}_{m}$.
Here, the overline $\overline{\bullet}$ denotes the ensemble average over the
random potentials, and we use 10000 samples for this.

Figure 3:  Comparison of our theoretical prediction with the experimental
data. The parameters are $2L=40$, $\nu=0.44$, $I=0.91$, and $M=8$, which are
taken from Ref. Wienand  _et al._ (2023). The time evolution of the variance,
denoted by circle and pentagon markers, is numerically obtained using Eq. (11)
with $\Delta=0$ and $1$. For the comparison with the experiment, we plot the
subtracted variance $2\sigma^{\prime}_{M}(t)^{2}-2\sigma^{\prime}_{M}(0)^{2}$
for the theoretical results and $\sigma_{\rm exp}(t)^{2}-\sigma_{\rm
exp}(t_{0})^{2}$ for the experimental result because the experiment studies
two-ladder lattice systems and thus the observed variance is twice as large as
ours. Here, the experimental data $\sigma_{\rm exp}(t)^{2}$ is taken from Fig.
S8 of Ref. Wienand  _et al._ (2023), and $t_{0}$ is about $0.0006$. The dotted
and dashed lines correspond to our analytical expressions of Eqs. (6) and (10)
for $t\gg 1$, respectively. The plus marker denotes the numerical data of the
6th-order perturbative calculation.

Figure 3 compares our theoretical results with the experimental data of Ref.
Wienand  _et al._ (2023). Our numerical result with $\Delta=0$ can capture the
variance growth in the early stage of the dynamics, but the deviation from the
experimental data grows in time. On the other hand, when we turn on the random
potential with $\Delta=1$, the numerical result well reproduces the
experimental data. The saturation of the variance for $t\gtrsim 2$ is caused
by the finite size effect and the disorder potential. These findings were
reported in the numerical simulations of Ref. Wienand  _et al._ (2023). We
also implement the numerical perturbative calculation with $\Delta$, finding
that the 6th-order perturbative result describes the time evolution before the
saturation (see Sec. VII of the Supplemental Material for the detailed
numerical method). These results show that the disorder effect is irrelevant
in $t\lesssim 1$. For an arbitrary strength $\Delta$, we can identify such a
time scale at which the disorder begins to affect the dynamics, by computing
the variances for $\Delta=0$ and $\Delta\neq 0$.

The dotted and dashed lines in Fig. 3 show the analytical results for the
linear growth of Eqs. (6) and (10). We find that the latter including the
incompleteness of the initial alternating state exhibits the reasonable
agreement with the experimental linear growth, while the former not including
the incompleteness systematically deviates from the experimental data. We here
emphasize that our analytical result (10) quantitatively describes the
experimental data without any fitting parameters under the assumption of (7).
Note that one can see that the disorder effect emerges before the finite $M$
effect by comparing the numerical data for $\Delta=0$ and 1 with Eq. (10)
obtained in the limit ($L,M\rightarrow\infty$), though the localization length
is larger than $M=8$ Wienand  _et al._ (2023). In general, one can identify
which the finite $M$ effect or the disorder effect emerges earlier, in the
same manner with the help of Eq. (10).

Figure 4:  Numerical results for the variance
$\sigma_{M}^{\prime\prime}(t)^{2}$ with $M=40$. The number of the lattice is
$2L=120$. We numerically solve the Schrödinger equation using the TEBD method
Vidal (2003, 2004); Schollwöck (2011); Paeckel  _et al._ (2019), calculating
the variances for $U=0,-1,-2$, and $-3$. The dashed line denotes the leading
term of Eq. (6) for $t\gg 1$, and the dotted and dot-dashed lines are curves
for $t^{2/3}$ and $t^{1/5}$, respectively. We insert these curves with the
exponents $2/3$ and $1/5$ for reference.

Numerical study for the variance growth of the interacting fermions.— We
numerically study the interaction effect on the time evolution of the
variance. Our numerical method is the time-evolving decimation method Vidal
(2003, 2004); Schollwöck (2011); Paeckel  _et al._ (2019) with Eq. (1) and the
open boundary condition ($\hat{a}_{L}=0$). We set $L$ to be $60$, and compute
the variance for $U=0,-1,-2$, and $-3$. In the language of the XXZ chain,
$U=-2$ corresponds to the critical point at which the model is identical to
the XXX chain Franchini (2017). To weaken the boundary effect, we use the
variance
$\sigma^{\prime\prime}_{M}(t)^{2}\coloneqq\langle(\hat{N}_{M}^{\prime})^{2}\rangle_{t}-\langle\hat{N}_{M}^{\prime}\rangle_{t}^{2}$
with $M=40$.

Figure 4 displays the time evolution of the variance
$\sigma^{\prime\prime}_{M}(t)^{2}$. In the noninteracting fermions ($U=0$),
our numerical result for the finite system well reproduces the leading term of
Eq. (6) for the infinite system. This demonstrates that the boundary effect is
negligible. In the interacting cases ($U\neq 0$), the time evolution exhibits
different growth behaviors, and we cannot find the clear ballistic property,
particularly for the $U=-2$ and $-3$ cases.

Note that Cecile et al. recently reported similar numerical results for the
variance in the XXZ chain Cecile  _et al._ (2024). For example, they confirmed
the signature of the anomalous power law growth
($\sigma^{\prime\prime}_{M}(t)^{2}\propto t^{2/3}$) in the dynamics starting
from the N$\rm\acute{e}$el state identical to the alternating state. The
similar signature was also discussed in Ref. Fujimoto  _et al._ (2020). We
here display Fig. 4 to show that our analytical result for the noninteracting
fermions does not work in the interacting cases.

Conclusions and prospects.— We theoretically studied the time evolution for
the variance of the particle numbers in the subsystem of the fermions on the
one-dimensional lattice, quantitatively explaining the recent experimental
result of Ref. Wienand  _et al._ (2023). The initial state used in this work
was the alternating state. In the noninteracting fermions, we analytically
derived the exact and simple expression of the variance, elucidating that the
variance linearly grows in the long-time dynamics. Developing the exact
analysis, we incorporated the incompleteness of the initial alternating state
into our analytical expression of the variance, finding that the improved
expression shows the reasonable agreement with the experimental result of Ref.
Wienand  _et al._ (2023) without any fitting parameters. Finally, our
numerical results based on the TEBD method Vidal (2003, 2004); Schollwöck
(2011); Paeckel  _et al._ (2019) find that the presence of the interaction
qualitatively alters the variance growth.

As a prospect, it is intriguing to study the variance growth in other models
accompanying quantum phase transitions, such as the transverse Ising model and
the Su-Schrieffer-Heeger model. These models are mapped into noninteracting
fermions, and one may analytically study the variance growth from the
perspective of the universal nature of the phase transitions. Furthermore,
studying the variance of other physical quantities, such as energy and spins,
is fascinating. In another direction, it is worth studying the variance growth
in open quantum systems. Depending on kinds of interactions with environments,
an open quantum system approaches a nonequilibrium stationary state that is
completely different from the equilibrium state addressed in this work. Thus,
uncovering features of the variance growth in such a case is fundamentally
interesting.

###### Acknowledgements.

KF is grateful to Ryusuke Hamazaki and Yuki Kawaguchi for fruitful discussions
and comments on the manuscript, Masaya Kunimi for helpful discussions on TEBD
with a conserved quantity, and Hiroki Moriya for helpful discussions. KF and
ST are grateful to Monika Aidelsburger, Immanuel Bloch, Alexander Impertro,
Simon Karch, Christian Schweizer, and Julian F. Wienand for sharing the
experimental data of Ref. Wienand  _et al._ (2023) used in Fig. 3. The work of
KF has been supported by JSPS KAKENHI Grant No. JP23K13029. The work of TS has
been supported by JSPS KAKENHI Grants No. JP21H04432, No. JP22H01143.

## References

  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Eisert  _et al._ (2015) J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nature Physics 11, 124 (2015).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Torres-Herrera  _et al._ (2015) E. J. Torres-Herrera, D. Kollmar, and L. F. Santos, Relaxation and thermalization of isolated many-body quantum systems, Physica Scripta 2015, 014018 (2015).
  * Luca D’Alessio and Rigol (2016) A. P. Luca D’Alessio, Yariv Kafri and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Mori  _et al._ (2018) T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Thermalization and prethermalization in isolated quantum systems: a theoretical overview, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 112001 (2018).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Gong and Hamazaki (2022) Z. Gong and R. Hamazaki, Bounds in nonequilibrium quantum dynamics, International Journal of Modern Physics B 36, 2230007 (2022).
  * Deutsch (1991) J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki (1994) M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Kinoshita  _et al._ (2006) T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum newton’s cradle, Nature 440, 900 (2006).
  * Gring  _et al._ (2012) M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, Relaxation and prethermalization in an isolated quantum system, Science 337, 1318 (2012).
  * Cazalilla (2006) M. A. Cazalilla, Effect of suddenly turning on interactions in the luttinger model, Phys. Rev. Lett. 97, 156403 (2006).
  * Rigol  _et al._ (2006) M. Rigol, A. Muramatsu, and M. Olshanii, Hard-core bosons on optical superlattices: Dynamics and relaxation in the superfluid and insulating regimes, Phys. Rev. A 74, 053616 (2006).
  * Rigol  _et al._ (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, 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 (2007).
  * Calabrese  _et al._ (2011) P. Calabrese, F. H. L. Essler, and M. Fagotti, Quantum quench in the transverse-field ising chain, Phys. Rev. Lett. 106, 227203 (2011).
  * Fagotti and Essler (2013) M. Fagotti and F. H. L. Essler, Reduced density matrix after a quantum quench, Phys. Rev. B 87, 245107 (2013).
  * Langen  _et al._ (2015) T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Experimental observation of a generalized gibbs ensemble, Science 348, 207 (2015).
  * Ilievski  _et al._ (2015) E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Complete generalized gibbs ensembles in an interacting theory, Phys. Rev. Lett. 115, 157201 (2015).
  * Tang  _et al._ (2018) Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Thermalization near integrability in a dipolar quantum newton’s cradle, Phys. Rev. X 8, 021030 (2018).
  * Basko  _et al._ (2006) D. Basko, I. Aleiner, and B. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006).
  * Žnidarič  _et al._ (2008) M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the heisenberg XXZ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Bardarson  _et al._ (2012) J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ (2013) M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013).
  * Huse  _et al._ (2014) D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * Schreiber  _et al._ (2015) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Potter  _et al._ (2015) A. C. Potter, R. Vasseur, and S. A. Parameswaran, Universal properties of many-body delocalization transitions, Phys. Rev. X 5, 031033 (2015).
  * Vasseur and Moore (2016) R. Vasseur and J. E. Moore, Nonequilibrium quantum dynamics and transport: from integrability to many-body localization, Journal of Statistical Mechanics: Theory and Experiment 2016, 064010 (2016).
  * yoon Choi  _et al._ (2016) J. yoon Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Lüschen  _et al._ (2017) H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017).
  * Kohlert  _et al._ (2019) T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. Das Sarma, I. Bloch, and M. Aidelsburger, Observation of many-body localization in a one-dimensional system with a single-particle mobility edge, Phys. Rev. Lett. 122, 170403 (2019).
  * Bernien  _et al._ (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Shiraishi and Mori (2017) N. Shiraishi and T. Mori, Systematic construction of counterexamples to the eigenstate thermalization hypothesis, Phys. Rev. Lett. 119, 030601 (2017).
  * Turner  _et al._ (2018a) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018a).
  * Turner  _et al._ (2018b) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Quantum scarred eigenstates in a rydberg atom chain: Entanglement, breakdown of thermalization, and stability to perturbations, Phys. Rev. B 98, 155134 (2018b).
  * Lin and Motrunich (2019) C.-J. Lin and O. I. Motrunich, Exact quantum many-body scar states in the rydberg-blockaded atom chain, Phys. Rev. Lett. 122, 173401 (2019).
  * Shibata  _et al._ (2020) N. Shibata, N. Yoshioka, and H. Katsura, Onsager’s scars in disordered spin chains, Phys. Rev. Lett. 124, 180604 (2020).
  * Trotzky  _et al._ (2012) S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional bose gas, Nature Physics 8, 325 (2012).
  * Langen  _et al._ (2013) T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J. Schmiedmayer, Local emergence of thermal correlations in an isolated quantum many-body system, Nature Physics 9, 640 (2013).
  * Clos  _et al._ (2016) G. Clos, D. Porras, U. Warring, and T. Schaetz, Time-resolved observation of thermalization in an isolated quantum system, Phys. Rev. Lett. 117, 170401 (2016).
  * Kaufman  _et al._ (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Neill  _et al._ (2016) C. Neill, P. Roushan, M. Fang, Y. Chen, M. Kolodrubetz, Z. Chen, A. Megrant, R. Barends, B. Campbell, B. Chiaro, A. Dunsworth, E. Jeffrey, J. Kelly, J. Mutus, P. J. J. O’Malley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, A. Polkovnikov, and J. M. Martinis, Ergodic dynamics and thermalization in an isolated quantum system, Nature Physics 12, 1037 (2016).
  * Moll  _et al._ (2016) P. J. W. Moll, P. Kushwaha, N. Nandi, B. Schmidt, and A. P. Mackenzie, Evidence for hydrodynamic electron flow in pdcoo¡sub¿2¡/sub¿, Science 351, 1061 (2016).
  * Bandurin  _et al._ (2016) D. A. Bandurin, I. Torre, R. K. Kumar, M. B. Shalom, A. Tomadin, A. Principi, G. H. Auton, E. Khestanova, K. S. Novoselov, I. V. Grigorieva, L. A. Ponomarenko, A. K. Geim, and M. Polini, Negative local resistance caused by viscous electron backflow in graphene, Science 351, 1055 (2016).
  * Gooth  _et al._ (2018) J. Gooth, F. Menges, N. Kumar, V. Sü$\beta$, C. Shekhar, Y. Sun, U. Drechsler, R. Zierold, C. Felser, and B. Gotsmann, Thermal and electrical signatures of a hydrodynamic electron fluid in tungsten diphosphide, Nature Communications 9, 4093 (2018).
  * Krishna Kumar  _et al._ (2017) R. Krishna Kumar, D. A. Bandurin, F. M. D. Pellegrino, Y. Cao, A. Principi, H. Guo, G. . H. Auton, M. Ben Shalom, L. A. Ponomarenko, G. Falkovich, K. Watanabe, T. Taniguchi, I. . V. Grigorieva, L. S. Levitov, M. Polini, and A. . K. Geim, Superballistic flow of viscous electron fluid through graphene constrictions, Nature Physics 13, 1182 (2017).
  * Sulpizio  _et al._ (2019) J. A. Sulpizio, L. Ella, A. Rozen, J. Birkbeck, D. J. Perello, D. Dutta, M. Ben-Shalom, T. Taniguchi, K. Watanabe, T. Holder, R. Queiroz, A. Principi, A. Stern, T. Scaffidi, A. K. Geim, and S. Ilani, Visualizing poiseuille flow of hydrodynamic electrons, Nature 576, 75 (2019).
  * Aharon-Steinberg  _et al._ (2022) A. Aharon-Steinberg, T. Völkl, A. Kaplan, A. K. Pariari, I. Roy, T. Holder, Y. Wolf, A. Y. Meltzer, Y. Myasoedov, M. E. Huber, B. Yan, G. Falkovich, L. S. Levitov, M. Hücker, and E. Zeldov, Direct observation of vortices in an electron fluid, Nature 607, 74 (2022).
  * Castro-Alvaredo  _et al._ (2016) O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ (2016) B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Bulchandani  _et al._ (2017) V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore, Solvable hydrodynamics of quantum integrable systems, Phys. Rev. Lett. 119, 220604 (2017).
  * Bulchandani  _et al._ (2018) V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore, Bethe-boltzmann hydrodynamics and spin transport in the XXZ chain, Phys. Rev. B 97, 045407 (2018).
  * Doyon and Yoshimura (2017) B. Doyon and T. Yoshimura, A note on generalized hydrodynamics: inhomogeneous fields and other concepts, SciPost Phys. 2, 014 (2017).
  * Doyon  _et al._ (2017) B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Large-scale description of interacting one-dimensional bose gases: Generalized hydrodynamics supersedes conventional hydrodynamics, Phys. Rev. Lett. 119, 195301 (2017).
  * Doyon  _et al._ (2018) B. Doyon, T. Yoshimura, and J.-S. Caux, Soliton gases and generalized hydrodynamics, Phys. Rev. Lett. 120, 045301 (2018).
  * Collura  _et al._ (2018) M. Collura, A. De Luca, and J. Viti, Analytic solution of the domain-wall nonequilibrium stationary state, Phys. Rev. B 97, 081111 (2018).
  * De Nardis  _et al._ (2018) J. De Nardis, D. Bernard, and B. Doyon, Hydrodynamic diffusion in integrable systems, Phys. Rev. Lett. 121, 160603 (2018).
  * Schemmer  _et al._ (2019) M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Generalized hydrodynamics on an atom chip, Phys. Rev. Lett. 122, 090601 (2019).
  * Gopalakrishnan and Vasseur (2019) S. Gopalakrishnan and R. Vasseur, Kinetic theory of spin diffusion and superdiffusion in XXZ spin chains, Phys. Rev. Lett. 122, 127202 (2019).
  * Doyon (2020) B. Doyon, Lecture notes on Generalised Hydrodynamics, SciPost Phys. Lect. Notes , 18 (2020).
  * Alba  _et al._ (2021) V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, Generalized-hydrodynamic approach to inhomogeneous quenches: correlations, entanglement and quantum effects, Journal of Statistical Mechanics: Theory and Experiment 2021, 114004 (2021).
  * Malvania  _et al._ (2021) N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1d bose gases, Science 373, 1129 (2021).
  * Bouchoule and Dubail (2022) I. Bouchoule and J. Dubail, Generalized hydrodynamics in the one-dimensional bose gas: theory and experiments, Journal of Statistical Mechanics: Theory and Experiment 2022, 014003 (2022).
  * Essler (2023) F. H. Essler, A short introduction to generalized hydrodynamics, Physica A: Statistical Mechanics and its Applications 631, 127572 (2023), lecture Notes of the 15th International Summer School of Fundamental Problems in Statistical Physics.
  * Wienand  _et al._ (2023) J. F. Wienand, S. Karch, A. Impertro, C. Schweizer, E. McCulloch, R. Vasseur, S. Gopalakrishnan, M. Aidelsburger, and I. Bloch, Emergence of fluctuating hydrodynamics in chaotic quantum systems (2023), arXiv:2306.11457 [cond-mat.quant-gas] .
  * Groha  _et al._ (2018) S. Groha, F. H. L. Essler, and P. Calabrese, Full counting statistics in the transverse field Ising chain, SciPost Phys. 4, 043 (2018).
  * (68) In Ref. Groha  _et al._ (2018), a probability function for the bipartite fluctuation of a transverse Ising chain was analytically studied, but this model is different from the experimental situation .
  * Fujimoto  _et al._ (2020) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Family-vicsek scaling of roughness growth in a strongly interacting bose gas, Phys. Rev. Lett. 124, 210604 (2020).
  * Jin  _et al._ (2020) T. Jin, A. Krajenbrink, and D. Bernard, From stochastic spin chains to quantum Kardar-Parisi-Zhang dynamics, Phys. Rev. Lett. 125, 040603 (2020).
  * Fujimoto  _et al._ (2021) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Dynamical scaling of surface roughness and entanglement entropy in disordered fermion models, Phys. Rev. Lett. 127, 090601 (2021).
  * Fujimoto  _et al._ (2022) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Impact of dissipation on universal fluctuation dynamics in open quantum systems, Phys. Rev. Lett. 129, 110403 (2022).
  * Cecile  _et al._ (2024) G. Cecile, J. De Nardis, and E. Ilievski, Squeezed ensembles and anomalous dynamic roughening in interacting integrable chains, Phys. Rev. Lett. 132, 130401 (2024).
  * Aditya and Roy (2024) S. Aditya and N. Roy, Family-vicsek dynamical scaling and kardar-parisi-zhang-like superdiffusive growth of surface roughness in a driven one-dimensional quasiperiodic model, Phys. Rev. B 109, 035164 (2024).
  * Bhakuni and Lev (2023) D. S. Bhakuni and Y. B. Lev, Dynamic scaling relation in quantum many-body systems (2023), arXiv:2309.03273 [cond-mat.dis-nn] .
  * Vidal (2003) G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Phys. Rev. Lett. 91, 147902 (2003).
  * Vidal (2004) G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev. Lett. 93, 040502 (2004).
  * Schollwöck (2011) U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019).
  * Antal  _et al._ (2008) T. Antal, P. L. Krapivsky, and A. Rákos, Logarithmic current fluctuations in nonequilibrium quantum spin chains, Phys. Rev. E 78, 061115 (2008).
  * Eisler and Rácz (2013) V. Eisler and Z. Rácz, Full counting statistics in a propagating quantum front and random matrix spectra, Phys. Rev. Lett. 110, 060602 (2013).
  * (82) M. Ljubotina, M. Žnidarič, and T. Prosen, .
  * Moriya  _et al._ (2019) H. Moriya, R. Nagao, and T. Sasamoto, Exact large deviation function of spin current for the one dimensional XX spin chain with domain wall initial condition, Journal of Statistical Mechanics: Theory and Experiment 2019, 063105 (2019).
  * Gamayun  _et al._ (2020) O. Gamayun, O. Lychkovskiy, and J.-S. Caux, Fredholm determinants, full counting statistics and Loschmidt echo for domain wall profiles in one-dimensional free fermionic chains, SciPost Phys. 8, 036 (2020).
  * Jin  _et al._ (2021) T. Jin, T. Gautié, A. Krajenbrink, P. Ruggiero, and T. Yoshimura, Interplay between transport and quantum coherences in free fermionic systems, Journal of Physics A: Mathematical and Theoretical 54, 404001 (2021).
  * Gopalakrishnan and Vasseur (2023) S. Gopalakrishnan and R. Vasseur, Anomalous transport from hot quasiparticles in interacting spin chains, Reports on Progress in Physics 86, 036502 (2023).
  * Wei  _et al._ (2022) D. Wei, A. Rubio-Abadal, B. Ye, F. Machado, J. Kemp, K. Srakaew, S. Hollerith, J. Rui, S. Gopalakrishnan, N. Y. Yao, I. Bloch, and J. Zeiher, Quantum gas microscopy of kardar-parisi-zhang superdiffusion, Science 376, 716 (2022).
  * Krajnik  _et al._ (2022a) Ž. Krajnik, E. Ilievski, and T. Prosen, Absence of normal fluctuations in an integrable magnet, Phys. Rev. Lett. 128, 090604 (2022a).
  * Krajnik  _et al._ (2022b) Ž. Krajnik, J. Schmidt, V. Pasquier, E. Ilievski, and T. Prosen, Exact anomalous current fluctuations in a deterministic interacting model, Phys. Rev. Lett. 128, 160601 (2022b).
  * Krajnik  _et al._ (2024a) Ž. Krajnik, J. Schmidt, E. Ilievski, and T. Prosen, Dynamical criticality of magnetization transfer in integrable spin chains, Phys. Rev. Lett. 132, 017101 (2024a).
  * Krajnik  _et al._ (2024b) Ž. Krajnik, J. Schmidt, V. Pasquier, T. Prosen, and E. Ilievski, Universal anomalous fluctuations in charged single-file systems, Phys. Rev. Res. 6, 013260 (2024b).
  * Klich and Levitov (2009) I. Klich and L. Levitov, Quantum noise as an entanglement meter, Phys. Rev. Lett. 102, 100502 (2009).
  * Song  _et al._ (2010) H. F. Song, S. Rachel, and K. Le Hur, General relation between entanglement and fluctuations in one dimension, Phys. Rev. B 82, 012405 (2010).
  * Song  _et al._ (2012) H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Bipartite fluctuations as a probe of many-body entanglement, Phys. Rev. B 85, 035409 (2012).
  * Rachel  _et al._ (2012) S. Rachel, N. Laflorencie, H. F. Song, and K. Le Hur, Detecting quantum critical points using bipartite fluctuations, Phys. Rev. Lett. 108, 116401 (2012).
  * Parez  _et al._ (2021) G. Parez, R. Bonsignori, and P. Calabrese, Quasiparticle dynamics of symmetry-resolved entanglement after a quench: Examples of conformal field theories and free fermions, Phys. Rev. B 103, L041104 (2021).
  * Oshima and Fuji (2023) H. Oshima and Y. Fuji, Charge fluctuation and charge-resolved entanglement in a monitored quantum circuit with U(1) symmetry, Phys. Rev. B 107, 014308 (2023).
  * Bertini  _et al._ (2023a) B. Bertini, P. Calabrese, M. Collura, K. Klobas, and C. Rylands, Nonequilibrium full counting statistics and symmetry-resolved entanglement from space-time duality, Phys. Rev. Lett. 131, 140401 (2023a).
  * Schönhammer (2007) K. Schönhammer, Full counting statistics for noninteracting fermions: Exact results and the levitov-lesovik formula, Phys. Rev. B 75, 205329 (2007).
  * Collura  _et al._ (2017) M. Collura, F. H. L. Essler, and S. Groha, Full counting statistics in the spin-1/2 heisenberg XXZ chain, Journal of Physics A: Mathematical and Theoretical 50, 414002 (2017).
  * Stéphan and Pollmann (2017) J.-M. Stéphan and F. Pollmann, Full counting statistics in the haldane-shastry chain, Phys. Rev. B 95, 035119 (2017).
  * Najafi and Rajabpour (2017) K. Najafi and M. A. Rajabpour, Full counting statistics of the subsystem energy for free fermions and quantum spin chains, Phys. Rev. B 96, 235109 (2017).
  * Arzamasovs and Gangardt (2019) M. Arzamasovs and D. M. Gangardt, Full counting statistics and large deviations in a thermal 1d bose gas, Phys. Rev. Lett. 122, 120401 (2019).
  * Calabrese  _et al._ (2020) P. Calabrese, M. Collura, G. D. Giulio, and S. Murciano, Full counting statistics in the gapped XXZ spin chain, Europhysics Letters 129, 60007 (2020).
  * Ares  _et al._ (2021) F. Ares, M. A. Rajabpour, and J. Viti, Exact full counting statistics for the staggered magnetization and the domain walls in the XY spin chain, Phys. Rev. E 103, 042107 (2021).
  * Smith  _et al._ (2021) N. R. Smith, P. L. Doussal, S. N. Majumdar, and G. Schehr, Full counting statistics for interacting trapped fermions, SciPost Phys. 11, 110 (2021).
  * McCulloch  _et al._ (2023) E. McCulloch, J. De Nardis, S. Gopalakrishnan, and R. Vasseur, Full counting statistics of charge in chaotic many-body quantum systems, Phys. Rev. Lett. 131, 210402 (2023).
  * Hercé  _et al._ (2023) G. Hercé, J.-P. Bureik, A. Ténart, A. Aspect, A. Dareau, and D. Clément, Full counting statistics of interacting lattice gases after an expansion: The role of condensate depletion in many-body coherence, Phys. Rev. Res. 5, L012037 (2023).
  * Bertini  _et al._ (2023b) B. Bertini, K. Klobas, M. Collura, P. Calabrese, and C. Rylands, Dynamics of charge fluctuations from asymmetric initial states (2023b), arXiv:2306.12404 [cond-mat.stat-mech] .
  * (110) See Supplemental Material for (I) Derivation of Eq. (2), (II) Derivation of Eq. (4), (III) Derivation of the two limiting formulae, (IV) Derivation of Eq. (6), (V) Derivation of Eq. (8), (VI) Derivation of Eq. (9), and (VII) Numerical perturbative calculation for the disorder potential .
  * Flesch  _et al._ (2008) A. Flesch, M. Cramer, I. P. McCulloch, U. Schollwöck, and J. Eisert, Probing local relaxation of cold atoms in optical superlattices, Phys. Rev. A 78, 033608 (2008).
  * mat (2015) 8 - special functions, in  _Table of Integrals, Series, and Products (Eighth Edition)_, edited by D. Zwillinger, V. Moll, I. Gradshteyn, and I. Ryzhik (Academic Press, Boston, 2015) eighth edition ed., pp. 867–1013.
  * Franchini (2017) F. Franchini,  _An introduction to integrable techniques for one-dimensional quantum systems_ , Vol. 940 (Springer, 2017).
  * Watson (1922) G. N. Watson,  _A treatise on the theory of Bessel functions_ , Vol. 2 (The University Press, 1922).
  * Ifantis and Siafarikas (1990) E. Ifantis and P. Siafarikas, Inequalities involving bessel and modified bessel functions, Journal of Mathematical Analysis and Applications 147, 214 (1990).
  * Anderson and Qiu (1997) G. D. Anderson and S.-L. Qiu, A monotoneity property of the gamma function, Proc. Amer. Math. Soc. 125, 3355 (1997).

## Supplemental Material for “Exact solution of bipartite fluctuations in one-
dimensional fermions”

Kazuya Fujimoto and Tomohiro Sasamoto

Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-
ku, Tokyo 152-8551, Japan

This Supplemental Material describes the following:

  * (I) Derivation of Eq. (2),

  * (II) Derivation of Eq. (4),

  * (III) Derivation of the two limiting formulae,

  * (IV) Derivation of Eq. (6),

  * (V) Derivation of Eq. (8),

  * (VI) Derivation of Eq. (9),

  * (VII) Numerical perturbative calculation for the disorder potential.

##  I Derivation of Eq. (2)

We derive the determinantal formula for the generating function, namely Eq.
(2) of the main text. The essential ingredient of the derivation is the fact
that the Wick theorem can be applicable since the quantum state in our setup
is the Gaussian state. The detailed calculation is given by

$\displaystyle G_{M}(\lambda,t)$
$\displaystyle=\left\langle\prod_{j=0}^{M-1}e^{\lambda\hat{a}_{j}^{\dagger}\hat{a}_{j}}\right\rangle_{t}$
$\displaystyle=\left\langle\prod_{j=0}^{M-1}\biggl{[}1+(e^{\lambda}-1)\hat{a}_{j}^{\dagger}\hat{a}_{j}\biggl{]}\right\rangle_{t}$
$\displaystyle=1+\sum_{n=1}^{M}(e^{\lambda}-1)^{n}\sum_{\begin{subarray}{c}j_{1}<j_{2}<...<j_{n}\\\
j_{k}\in\\{0,1,2,...,M-1\\}\end{subarray}}\left\langle\prod_{k=1}^{n}\hat{a}_{j_{k}}^{\dagger}\hat{a}_{j_{k}}\right\rangle_{t}$
$\displaystyle=1+\sum_{n=1}^{M}(e^{\lambda}-1)^{n}\sum_{\begin{subarray}{c}j_{1}<j_{2}<...<j_{n}\\\
j_{k}\in\\{0,1,2,...,M-1\\}\end{subarray}}{\rm
det}\left[\braket{\hat{a}_{j_{k}}^{\dagger}\hat{a}_{j_{l}}}_{t}\right]_{k,l=1}^{n}$
$\displaystyle={\rm
det}\left[\delta_{j,k}+(e^{\lambda}-1)\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}\right]_{j,k=0}^{M-1}.$
(S-1)

In the fourth line, we use the Wick theorem. This kind of determinantal form
was derived in several previous works Schönhammer (2007); Eisler and Rácz
(2013); Parez  _et al._ (2021).

##  II Derivation of Eq. (4)

We express the variance $\sigma_{M}(t)^{2}$ by the Bessel function $J_{n}(x)$
of the first kind. Using the property of the generating function, we can
obtain

$\displaystyle\sigma_{M}(t)^{2}$
$\displaystyle=\left.\dfrac{\partial^{2}}{\partial\lambda^{2}}G_{M}(\lambda,t)\right|_{\lambda=0}-\left(\left.\dfrac{\partial}{\partial\lambda}G_{M}(\lambda,t)\right|_{\lambda=0}\right)^{2}$
(S-2)
$\displaystyle=\sum_{k=0}^{M-1}\braket{\hat{a}_{k}^{\dagger}\hat{a}_{k}}_{t}-\sum_{j=0}^{M-1}\sum_{k=0}^{M-1}|\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}|^{2}.$
(S-3)

The first term on the right-hand side of Eq. (S-3) becomes

$\displaystyle\sum_{k=0}^{M-1}\braket{\hat{a}_{k}^{\dagger}\hat{a}_{k}}_{t}=\frac{M}{2}-\frac{J_{0}(4t)}{2}\sum_{k=0}^{M-1}(-1)^{k},$
(S-4)

where we use the following expression of the two-point correlator,

$\displaystyle\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}=\frac{\delta_{j,k}}{2}-\frac{{\rm
i}^{j+k}}{2}J_{k-j}(4t).$ (S-5)

Similarly, the second term on the right-hand side of Eq. (S-3) becomes

$\displaystyle\sum_{j=1}^{M-1}\sum_{k=1}^{M-1}|\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}|^{2}$
$\displaystyle=\frac{M}{4}-\frac{J_{0}(4t)}{2}\sum_{k=0}^{M-1}(-1)^{k}+\frac{1}{4}\sum_{j=0}^{M-1}\sum_{k=0}^{M-1}J_{j-k}(4t)^{2}.$
(S-6)

The double summation on the right-hand side of Eq. (S-6) becomes

$\displaystyle\sum_{j=0}^{M-1}\sum_{k=0}^{M-1}J_{j-k}(4t)^{2}$
$\displaystyle=\sum_{l=0}^{M-1}\sum_{m=0}^{M-l-1}J_{l}(4t)^{2}+\sum_{l=-M+1}^{-1}\sum_{m=-l}^{M-1}J_{l}(4t)^{2}$
(S-7)
$\displaystyle=MJ_{0}(4t)^{2}+2\sum_{l=1}^{M-1}\left(M-l\right)J_{l}(4t)^{2}$
(S-8)
$\displaystyle=M\left(J_{0}(4t)^{2}+2\sum_{l=1}^{M-1}J_{l}(4t)^{2}\right)-2\sum_{l=1}^{M-1}lJ_{l}(4t)^{2}.$
(S-9)

In Eq. (S-7), we introduce $l=j-k$ and divide the single double summation into
the two double summations. Figure S-1 schematically describes this procedure
in detail. We also use the formula $J_{-l}(4t)^{2}=J_{l}(4t)^{2}$ in Eq.
(S-8). Finally, putting Eq. (S-9) into Eq. (S-6), we get

$\displaystyle\sum_{j=1}^{M-1}\sum_{k=1}^{M-1}|\braket{\hat{a}_{j}^{\dagger}\hat{a}_{k}}_{t}|^{2}$
$\displaystyle=\frac{M}{4}-\frac{J_{0}(4t)}{2}\sum_{k=0}^{M-1}(-1)^{k}+\frac{M}{4}\left(J_{0}(4t)^{2}+2\sum_{l=1}^{M-1}J_{l}(4t)^{2}\right)-\dfrac{1}{2}\sum_{l=1}^{M-1}lJ_{l}(4t)^{2}.$
(S-10)

Substituting Eqs. (S-4) and (S-10) into (S-3), we derive Eq. (4) of the main
text, namely

$\displaystyle\sigma_{M}(t)^{2}=\frac{M}{4}\left(1-J_{0}(4t)^{2}-2\sum_{k=1}^{M-1}J_{k}(4t)^{2}\right)+\frac{1}{2}\sum_{k=1}^{M-1}kJ_{k}(4t)^{2}.$
(S-11)

Figure S-1:  Schematic illustration in deriving Eq. (S-7). The filled circles
denotes the points for the double summation
$\sum_{j=0}^{M-1}\sum_{k=0}^{M-1}\bullet$. In the blue region specified by the
set $\\{(k,j)\in\\{0,1,...,M-1\\}^{2}\mid j\geq k\\}$, the $j$-intercept
$l=j-k$ can take the values from $\\{0,1,...,M-1\\}$. On the other hand, the
$j$-intercept takes the values from $\\{-M+1,-M+2,...,-1\\}$ in the red region
specified by the set $\\{(k,j)\in\\{0,1,...,M-1\\}^{2}\mid j<k\\}$. Using the
$j$-intercept $l$, we can transform the double summation of the left-hand side
of Eq. (S-7) into the right-hand side of Eq. (S-7).

##  III Derivation of the two limiting formulae

In the derivation of Eq. (5), we use the following formulae,

$\displaystyle\lim_{M\rightarrow\infty}\frac{M}{4}\left(1-J_{0}(4t)^{2}-2\sum_{k=1}^{M-1}J_{k}(4t)^{2}\right)=0~{}~{}~{}(t>0),$
(S-12)
$\displaystyle\lim_{M\rightarrow\infty}\frac{1}{2}\sum_{k=1}^{M-1}kJ_{k}(4t)^{2}=4t^{2}\left(J_{0}(4t)^{2}+J_{1}(4t)^{2}\right)-tJ_{0}(4t)J_{1}(4t).$
(S-13)

In this section, we shall prove them.

###  III.1  Proof of Eq. (S-12)

For the proof of Eq. (S-12), we first derive the following inequality,

$\displaystyle\frac{M}{4}\left|\left(1-J_{0}(4t)^{2}-2\sum_{k=1}^{M-1}J_{k}(4t)^{2}\right)\right|<\dfrac{1}{2}M^{(\gamma-1)M+1}(M+1)^{(\gamma-1)(M+1)+1}(2t)^{2M}$
(S-14)

for $t>0$. Here $\gamma\simeq 0.577$ is the Euler-Mascheroni constant. The
proof of Eq. (S-14) is described in the following. We employ the recursion
relation $2dJ_{n}(x)/dx=J_{n-1}(x)-J_{n+1}(x)$ for the Bessel function of the
first kind, deriving

$\displaystyle\frac{d}{dx}\sum_{n=1}^{M-1}J_{n}(x)^{2}$
$\displaystyle=\sum_{n=1}^{M-1}\left(J_{n}(x)J_{n-1}(x)-J_{n}(x)J_{n+1}(x)\right)$
(S-15) $\displaystyle=J_{1}(x)J_{0}(x)-J_{M-1}(x)J_{M}(x)$ (S-16)
$\displaystyle=-\frac{1}{2}\frac{d}{dx}J_{0}(x)^{2}-J_{M-1}(x)J_{M}(x),$
(S-17)

where we use the formula $J_{1}(x)=-dJ_{0}(x)/dx$ in the last equality.
Integrating Eq. (S-17) from $0$ to $x$ leads to

$\displaystyle
J_{0}(x)^{2}+2\sum_{k=1}^{M-1}J_{k}(x)^{2}-1=-2\int_{0}^{x}J_{M-1}(y)J_{M}(y)dy,$
(S-18)

where we use $J_{0}(0)=1$ and $J_{l}(0)=0~{}(l>0)$. We next consider the
absolute value of Eq. (S-18) and then obtain

$\displaystyle\left|\left(1-J_{0}(x)^{2}-2\sum_{k=1}^{M-1}J_{k}(x)^{2}\right)\right|$
(S-19) $\displaystyle\leq
2\int_{0}^{x}\left|J_{M-1}(y)\right|\left|J_{M}(y)\right|dy$ (S-20)
$\displaystyle<2M^{(\gamma-1)M+1}(M+1)^{(\gamma-1)(M+1)+1}\int_{0}^{x}\left(\frac{y}{2}\right)^{2M-1}\exp\left(-\frac{y^{2}}{4M}-\frac{y^{2}}{4M+4}\right)dy$
(S-21)
$\displaystyle<2M^{(\gamma-1)M+1}(M+1)^{(\gamma-1)(M+1)+1}\int_{0}^{x}\left(\frac{y}{2}\right)^{2M-1}dy$
(S-22)
$\displaystyle=2M^{(\gamma-1)M}(M+1)^{(\gamma-1)(M+1)+1}\left(\frac{x}{2}\right)^{2M}.$
(S-23)

In Eq. (S-21), we use the following inequality,

$\displaystyle
J_{M}(x)<\left(\frac{x}{2}\right)^{M}(M+1)^{(\gamma-1)(M+1)+1}\exp\left(-\frac{x^{2}}{4M+4}\right)$
(S-24)

for $x>0$. This inequality can be derived using the two inequalities given by

$\displaystyle
J_{\nu}(x)<\left(\frac{x}{2}\right)^{\nu}\frac{1}{\Gamma(\nu+1)}\exp\left(-\frac{x^{2}}{4\nu+4}\right)$
(S-25)

for $\nu\geq 0$ and $x>0$, and

$\displaystyle x^{(1-\gamma)x-1}<\Gamma(x)<x^{x-1}$ (S-26)

for $x\in(1,\infty)$. Here, $\Gamma(x)$ denotes the Gamma function. The
inequalities of Eqs. (S-25) and (S-26) were given in Refs. Watson (1922);
Ifantis and Siafarikas (1990) and Anderson and Qiu (1997), respectively. Thus,
employing the inequality of Eq. (S-23) with $x=4t$, we can readily derive the
inequality of Eq. (S-14).

Finally, we take the limit $M\rightarrow\infty$ in Eq. (S-14). For this
purpose, we consider the logarithmic of the right-hand side of Eq. (S-14),
which reads

$\displaystyle-\log 2+\left[(\gamma-1)M+1\right]\log
M+\left[(\gamma-1)(M+1)+1\right]\log(M+1)+2M\log(2t).$ (S-27)

Thus, when $M$ is much larger than unity, the above for a fixed nonzero time
$t$ can be approximated to be $2(\gamma-1)M\log M$. This becomes $-\infty$ for
$M\rightarrow\infty$ because the Euler-Mascheroni constant $\gamma\simeq
0.577$ is smaller than unity. Hence, the left-hand side of Eq. (S-14) for
$M\rightarrow\infty$ can be bounded by zero. This proves Eq. (S-12).

###  III.2  Proof of Eq. (S-13)

We shall prove Eq. (S-13) by noting the recursion relation
$2kJ_{k}(x)/x=J_{k-1}(x)+J_{k+1}(x)$ for the Bessel function of the first
kind. This relation leads to

$\displaystyle\frac{2k}{x}J_{k}(x)J_{k}(y)=J_{k-1}(x)J_{k}(y)+J_{k+1}(x)J_{k}(y),$
(S-28)
$\displaystyle\frac{2k}{y}J_{k}(x)J_{k}(y)=J_{k-1}(y)J_{k}(x)+J_{k+1}(y)J_{k}(x).$
(S-29)

Subtracting them and taking the summation for the index $k$ from zero to a
positive integer $N$, we obtain

$\displaystyle\left(\frac{2}{x}-\frac{2}{y}\right)\sum_{k=0}^{N}kJ_{k}(x)J_{k}(y)$
$\displaystyle=\sum_{k=1}^{N}(J_{k+1}(x)J_{k}(y)-J_{k+1}(y)J_{k}(x)+J_{k-1}(x)J_{k}(y)-J_{k-1}(y)J_{k}(x))$
(S-30)
$\displaystyle=J_{N+1}(x)J_{N}(y)-J_{N+1}(y)J_{N}(x)+J_{0}(x)J_{1}(y)-J_{0}(y)J_{1}(x),$
(S-31)

where we use the property of the telescoping series to derive the last
expression. Taking the limit $N\rightarrow\infty$, we get

$\displaystyle\sum_{k=0}^{\infty}kJ_{k}(x)J_{k}(y)=\frac{xy}{2}\frac{J_{0}(x)J_{1}(y)-J_{0}(y)J_{1}(x)}{y-x}.$
(S-32)

Finally, we take the limit $y\rightarrow x$ in Eq. (S-32) via the L’Hôpital’s
rule, obtaining

$\displaystyle\sum_{k=0}^{\infty}kJ_{k}(x)^{2}=\frac{x^{2}}{2}J_{0}(x)^{2}+\frac{x^{2}}{2}J_{1}(x)^{2}-\frac{x}{2}J_{0}(x)J_{1}(x).$
(S-33)

Here, we use the following formulae,

$\displaystyle\frac{d}{dx}J_{k+1}(x)=J_{k}(x)-\frac{k+1}{x}J_{k+1}(x),$ (S-34)
$\displaystyle\frac{d}{dx}J_{k-1}(x)=-J_{k}(x)+\frac{k-1}{x}J_{k-1}(x).$
(S-35)

Finally, we put $x=4t$ into Eq. (S-33), completing the proof of Eq. (S-13).

##  IV  Derivation of Eq. (6)

We derive the asymptotic expression of the variance given by Eq. (6) of the
main text for $t\ll 1$ and $t\gg 1$. For the short-time dynamics ($t\ll 1$),
the Bessel functions $J_{0}(4t)$ and $J_{1}(4t)$ of the first kind is
approximated by

$\displaystyle J_{0}(4t)\simeq 1,$ (S-36) $\displaystyle J_{1}(4t)\simeq 2t,$
(S-37)

which can be readily derived using the infinite series of $J_{n}(x)$.
Substituting Eqs. (S-36) and (S-37) into Eq. (5) of the main text, we obtain,
up to the leading order,

$\displaystyle\sigma(t)^{2}\simeq 2t^{2}.$ (S-38)

For the long-time dynamics ($t\gg 1$), we use the asymptotic formula mat
(2015) given by

$\displaystyle
J_{n}(4t)=\frac{\displaystyle\cos\left(4t-\frac{\pi}{2}n-\frac{\pi}{4}\right)}{\displaystyle\sqrt{2\pi
t}}\left(1-\frac{\Gamma(n+5/2)}{128\Gamma(n-3/2)}t^{-2}\right)-\frac{\displaystyle\sin\left(4t-\frac{\pi}{2}n-\frac{\pi}{4}\right)}{\displaystyle\sqrt{2\pi}t^{3/2}}\frac{\Gamma(n+3/2)}{8\Gamma(n-1/2)}+\mathcal{O}(t^{-7/2}).$
(S-39)

This expression leads to

$\displaystyle J_{0}(4t)^{2}$ $\displaystyle\simeq\dfrac{1}{2\pi
t}\cos\left(4t-\dfrac{\pi}{4}\right)^{2}+\dfrac{1}{64\pi
t^{2}}\cos\left(8t\right)+\dfrac{1}{2048\pi
t^{3}}\sin\left(4t-\dfrac{\pi}{4}\right)^{2}-\dfrac{9}{2048\pi
t^{3}}\cos\left(4t-\dfrac{\pi}{4}\right)^{2},$ (S-40) $\displaystyle
J_{1}(4t)^{2}$ $\displaystyle\simeq\dfrac{1}{2\pi
t}\sin\left(4t-\dfrac{\pi}{4}\right)^{2}-\dfrac{3}{64\pi
t^{2}}\cos\left(8t\right)+\dfrac{9}{2048\pi
t^{3}}\cos\left(4t-\dfrac{\pi}{4}\right)^{2}+\dfrac{15}{2048\pi
t^{3}}\sin\left(4t-\dfrac{\pi}{4}\right)^{2},$ (S-41) $\displaystyle
J_{0}(4t)J_{1}(4t)$ $\displaystyle\simeq-\dfrac{1}{8\pi
t}\cos\left(8t\right)+\dfrac{3}{64\pi
t^{2}}\cos\left(4t-\dfrac{\pi}{4}\right)^{2}+\dfrac{1}{64\pi
t^{2}}\sin\left(4t-\dfrac{\pi}{4}\right)^{2}.$ (S-42)

Employing these expressions and Eq. (5), we derive

$\displaystyle\sigma(t)^{2}\simeq\dfrac{2}{\pi}t-\dfrac{1}{64\pi
t}\left(2\sin(8t)+1\right).$ (S-43)

##  V  Derivation of Eq. (8)

We shall derive the explicit expression of the two-point correlator with the
incomplete alternating initial state $\hat{\rho}_{\rm alt}$ defined by Eq. (7)
of the main text. We first define the two-point correlator $C_{m,n}(t)$ as

$\displaystyle C_{m,n}(t)\coloneqq{\rm Tr}\left[e^{-{\rm
i}\hat{H}t}\hat{\rho}_{\rm alt}e^{{\rm
i}\hat{H}t}\hat{a}^{\dagger}_{m}\hat{a}_{n}\right].$ (S-44)

The straightforward calculation leads to the equation of motion for
$C_{m,n}(t)$,

$\displaystyle{\rm
i}\dfrac{d}{dt}C_{m,n}(t)=C_{m+1,n}(t)+C_{m-1,n}(t)-C_{m,n+1}(t)-C_{m,n-1}(t).$
(S-45)

To solve the differential equation, we expand the correlator via the discrete
Fourier transformation, getting

$\displaystyle
C_{m,n}(t)=\dfrac{1}{4L^{2}}\sum_{\alpha=-L}^{L-1}\sum_{\beta=-L}^{L-1}D_{\alpha,\beta}(t)e^{{\rm
i}\pi(n\beta-m\alpha)/L},$ (S-46)

where $D_{\alpha,\beta}(t)$ is the coefficient of the expansion. Substituting
Eq. (S-46) into Eq. (S-45), we readily derive

$\displaystyle{\rm
i}\dfrac{d}{dt}D_{\alpha,\beta}(t)=\mathcal{E}_{\alpha,\beta}D_{\alpha,\beta}(t)$
(S-47)

with $\mathcal{E}_{\alpha,\beta}\coloneqq
2\cos(\pi\alpha/L)-2\cos(\pi\beta/L)$. The initial condition
$D_{\alpha,\beta}(0)$ is calculated as

$\displaystyle D_{\alpha,\beta}(0)$
$\displaystyle=\sum_{m=-L}^{L-1}\sum_{n=-L}^{L-1}C_{m,n}(0)e^{{\rm
i}\pi(-n\beta+m\alpha)/L}$ (S-48) $\displaystyle=n_{\rm
even}\sum_{m=-L/2}^{L/2-1}e^{{\rm i}2\pi m(\alpha-\beta)/L}+n_{\rm
odd}\sum_{m=-L/2}^{L/2-1}e^{{\rm i}\pi(2m+1)(\alpha-\beta)/L}$ (S-49)
$\displaystyle=Ln_{\rm
even}(\delta_{\alpha,\beta}+\delta_{\alpha,\beta+L}+\delta_{\alpha,\beta-L})+Ln_{\rm
odd}(\delta_{\alpha,\beta}-\delta_{\alpha,\beta+L}-\delta_{\alpha,\beta-L}),$
(S-50)

where we use $-2L+1\leq\alpha\pm\beta\leq 2L-1$ since the values of $\alpha$
and $\beta$ are restricted to the 1st Brillouin zone. As to the initial
condition $C_{m,n}(0)$, by the definition of $\hat{\rho}_{\rm alt}$, we can
obtain

$\displaystyle C_{m,n}(0)=\begin{dcases}n_{\rm even}&(n=m\wedge n~{}{\rm
is~{}even})\\\ n_{\rm odd}&(n=m\wedge n~{}{\rm is~{}odd})\\\ 0&({\rm
otherwise}),\end{dcases}$ (S-51)

where $n_{\rm even}$ and $n_{\rm odd}$ are averaged particle numbers at the
even and odd sites for the initial density matrix $\hat{\rho}_{\rm alt}$,
respectively. Solving the differential equation of Eq. (S-47), we get

$\displaystyle D_{\alpha,\beta}(t)$ $\displaystyle=D_{\alpha,\beta}(0)e^{-{\rm
i}\mathcal{E}_{\alpha,\beta}t}$ (S-52) $\displaystyle=Ln_{\rm
T}\delta_{\alpha,\beta}+Ln_{\rm
D}\left(\delta_{\alpha,\beta+L}+\delta_{\alpha,\beta-L}\right)e^{{\rm
i}4t\cos(\pi\beta/L)},$ (S-53)

where we define $n_{\rm T}\coloneqq n_{\rm even}+n_{\rm odd}$ and $n_{\rm
D}\coloneqq n_{\rm even}-n_{\rm odd}$. We substitute Eq. (S-53) into Eq.
(S-46), getting

$\displaystyle C_{m,n}(t)$
$\displaystyle=\dfrac{1}{4L^{2}}\sum_{\alpha=-L}^{L-1}\sum_{\beta=-L}^{L-1}D_{\alpha,\beta}(t)e^{{\rm
i}\pi(n\beta-m\alpha)/L}$ (S-54)
$\displaystyle=\dfrac{1}{4L}\sum_{\alpha=-L}^{L-1}\sum_{\beta=-L}^{L-1}\left(n_{\rm
T}\delta_{\alpha,\beta}+n_{\rm
D}\left(\delta_{\alpha,\beta+L}+\delta_{\alpha,\beta-L}\right)e^{{\rm
i}4t\cos(\pi\beta/L)}\right)e^{{\rm i}\pi(n\beta-m\alpha)/L}$ (S-55)
$\displaystyle=\dfrac{n_{\rm T}}{4L}\sum_{\alpha=-L}^{L-1}e^{{\rm
i}\pi(n-m)\alpha/L}+\dfrac{n_{\rm D}(-1)^{n}}{4L}\sum_{\alpha=0}^{L-1}e^{-{\rm
i}4t\cos(\pi\alpha/L)}e^{{\rm i}\pi(n-m)\alpha/L}+\dfrac{n_{\rm
D}(-1)^{n}}{4L}\sum_{\alpha=-L}^{-1}e^{-{\rm i}4t\cos(\pi\alpha/L)}e^{{\rm
i}\pi(n-m)\alpha/L}$ (S-56) $\displaystyle=\dfrac{n_{\rm
T}}{4L}\sum_{\alpha=-L}^{L-1}e^{{\rm i}\pi(n-m)\alpha/L}+\dfrac{n_{\rm
D}(-1)^{n}}{4L}\sum_{\alpha=-L}^{L-1}e^{-{\rm i}4t\cos(\pi\alpha/L)}e^{{\rm
i}\pi(n-m)\alpha/L}$ (S-57) $\displaystyle=\dfrac{n_{\rm
T}}{2}\delta_{m,n}+\dfrac{n_{\rm D}(-1)^{n}}{4L}\sum_{\alpha=-L}^{L-1}e^{-{\rm
i}4t\cos(\pi\alpha/L)}e^{{\rm i}\pi(n-m)\alpha/L},$ (S-58)

where we again use $-2L+1\leq\alpha\pm\beta\leq 2L-1$. Finally, we take the
thermodynamic limit ($L\rightarrow\infty$), obtaining

$\displaystyle C_{m,n}(t)$ $\displaystyle=\dfrac{n_{\rm
T}}{2}\delta_{m,n}+\dfrac{n_{\rm
D}(-1)^{n}}{4\pi}\int_{-\pi}^{\pi}d\theta~{}e^{{\rm i}(n-m)\theta-{\rm
i}4t\cos\theta}$ (S-59) $\displaystyle=\dfrac{n_{\rm
T}}{2}\delta_{m,n}+\dfrac{n_{\rm
D}(-1)^{m}}{4\pi}\int_{0}^{2\pi}d\theta~{}e^{{\rm i}(n-m)\theta+{\rm
i}4t\cos\theta}$ (S-60) $\displaystyle=\dfrac{n_{\rm
T}}{2}\delta_{m,n}+\dfrac{n_{\rm D}}{2}{\rm i}^{n+m}J_{n-m}(4t).$ (S-61)

In the last line, we use the integral formula for the Bessel function of the
first kind given by

$\displaystyle J_{n}(x)=\dfrac{1}{2\pi{\rm
i}^{n}}\int_{0}^{2\pi}d\theta~{}e^{{\rm i}n\theta+{\rm i}x\cos\theta}.$ (S-62)

##  VI Derivation of Eq. (9)

We shall derive Eq. (9) of the main text using Eq. (8) of the main text. Note
that the initial state $\hat{\rho}_{\rm alt}$ is the Gaussian state and thus
we can use the Wick theorem. As a result, the variance is given by

$\displaystyle\sigma_{M}(t)^{2}=\sum_{m=0}^{M-1}C_{m,m}(t)-\sum_{m=0}^{M-1}\sum_{n=0}^{M-1}|C_{m,n}(t)|^{2}.$
(S-63)

Following the almost same procedure as Sec. II, we obtain

$\displaystyle\sum_{m=0}^{M-1}C_{m,m}(t)$ $\displaystyle=\dfrac{M}{2}(n_{\rm
even}+n_{\rm odd})+\dfrac{1}{2}(n_{\rm even}-n_{\rm
odd})J_{0}(4t)\sum_{m=0}^{M-1}(-1)^{m},$ (S-64)
$\displaystyle\sum_{m=0}^{M-1}\sum_{n=0}^{M-1}|C_{m,n}(t)|^{2}$
$\displaystyle=\dfrac{M}{4}(n_{\rm even}+n_{\rm odd})^{2}+\dfrac{1}{2}(n_{\rm
even}^{2}-n_{\rm odd}^{2})J_{0}(4t)\sum_{m=0}^{M-1}(-1)^{m}$ (S-65)
$\displaystyle+\dfrac{1}{4}(n_{\rm even}-n_{\rm
odd})^{2}\sum_{m=0}^{M-1}\sum_{n=0}^{M-1}J_{n-m}^{2}(4t)$ (S-66)
$\displaystyle=\dfrac{M}{4}(n_{\rm even}+n_{\rm odd})^{2}+\dfrac{1}{2}(n_{\rm
even}^{2}-n_{\rm odd}^{2})J_{0}(4t)\sum_{m=0}^{M-1}(-1)^{m}$ (S-67)
$\displaystyle+\dfrac{M}{4}(n_{\rm even}-n_{\rm
odd})^{2}\left(J_{0}(4t)^{2}+2\sum_{m=1}^{M-1}J_{m}(4t)^{2}\right)-\dfrac{1}{2}(n_{\rm
even}-n_{\rm odd})^{2}\sum_{m=1}^{M-1}mJ_{m}(4t)^{2}.$ (S-68)

Thus, the variance becomes

$\displaystyle\sigma_{M}(t)^{2}$ $\displaystyle=\dfrac{M}{4}(n_{\rm
even}+n_{\rm odd})(2-n_{\rm even}-n_{\rm odd})+\dfrac{1}{2}(n_{\rm
even}-n_{\rm odd}-n_{\rm even}^{2}+n_{\rm
odd}^{2})J_{0}(4t)\sum_{m=0}^{M-1}(-1)^{m}$ (S-69)
$\displaystyle-\dfrac{M}{4}(n_{\rm even}-n_{\rm
odd})^{2}\left(J_{0}(4t)^{2}+2\sum_{m=1}^{M-1}J_{m}(4t)^{2}\right)+\dfrac{1}{2}(n_{\rm
even}-n_{\rm odd})^{2}\sum_{m=1}^{M-1}mJ_{m}(4t)^{2}.$ (S-70)

Next, we define $\delta\sigma_{M}(t)^{2}$ as

$\displaystyle\delta\sigma_{M}(t)^{2}\coloneqq\dfrac{M}{2}\bigl{(}n_{\rm
even}(1-n_{\rm even})+n_{\rm odd}(1-n_{\rm
odd})\bigl{)}+\dfrac{1}{2}J_{0}(4t)\bigl{(}n_{\rm even}(1-n_{\rm even})-n_{\rm
odd}(1-n_{\rm odd})\bigl{)}\sum_{m=0}^{M-1}(-1)^{m}.$ (S-71)

As a result, we can obtain

$\displaystyle\sigma_{M}(t)^{2}-\delta\sigma_{M}(t)^{2}=(n_{\rm even}-n_{\rm
odd})^{2}\left[\dfrac{M}{4}\left(1-J_{0}(4t)^{2}-2\sum_{m=1}^{M-1}J_{m}(4t)^{2}\right)+\dfrac{1}{2}\sum_{m=1}^{M-1}mJ_{m}(4t)^{2}\right].$
(S-72)

Finally, applying Eqs. (S-12) and (S-13) into the above, we derive Eq. (9) of
the main text, namely

$\displaystyle\lim_{M\rightarrow\infty}\bigl{(}\sigma_{M}(t)^{2}-\delta\sigma_{M}(t)^{2}\bigl{)}=(n_{\rm
even}-n_{\rm
odd})^{2}\bigl{[}4t^{2}\left(J_{0}(4t)^{2}+J_{1}(4t)^{2}\right)-tJ_{0}(4t)J_{1}(4t)\bigl{]}.$
(S-73)

##  VII Numerical perturbative calculation for the disorder potential

We shall explain how to implement the numerical perturbative calculation for
the disorder potential. In the main text, we define the random potential
$V_{j}$, which is independently sampled from the uniform distribution with the
range $[-\Delta,\Delta]$. Here, $\Delta\geq 0$ denotes the strength of the
randomness. This definition is inconvenient to proceed perturbative
calculation for the disorder potential, so thus we introduce a new random
variable $v_{n}$, which is independently sampled from the uniform distribution
with the range $[-1,1]$. Then, the Hamiltonian $\hat{H}^{\prime}$ of Eq. (11)
in the main text is written as

$\displaystyle\hat{H}^{\prime}=-\sum_{j=-L}^{L-1}\left(\hat{a}_{j+1}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{j+1}\right)+\Delta\sum_{j=-L}^{L-1}v_{j}\hat{n}_{j}.$
(S-74)

Using this notation, we can derive the equation of motion for the two-point
correlator $C_{m,n}(t)$

$\displaystyle{\rm
i}\dfrac{d}{dt}C_{m,n}(t)=C_{m+1,n}(t)+C_{m-1,n}(t)-C_{m,n+1}(t)-C_{m,n-1}(t)+\Delta(v_{n}-v_{m})C_{m,n}(t).$
(S-75)

In what follows, we regard $\Delta$ as a small parameter and construct the
perturbative equation of motion for the two-point correlator.

We expand the two-point correlator with the parameter $\Delta$ as

$\displaystyle C_{m,n}(t)=C_{m,n}^{(0)}(t)+\Delta
C_{m,n}^{(1)}(t)+\Delta^{2}C_{m,n}^{(2)}(t)+...~{}~{}.$ (S-76)

Substituting Eq. (S-76) into Eq. (S-75), we obtain the $\alpha$th-order
equation of motion as

$\displaystyle{\rm
i}\dfrac{d}{dt}C_{m,n}^{(\alpha)}(t)=C_{m+1,n}^{(\alpha)}(t)+C_{m-1,n}^{(\alpha)}(t)-C_{m,n+1}^{(\alpha)}(t)-C_{m,n-1}^{(\alpha)}(t)+(v_{n}-v_{m})C_{m,n}^{(\alpha-1)}(t),$
(S-77)

where we define $C_{m,n}^{-1}(t)\coloneqq 0$. When you calculate the variance
up to the $\alpha$th-order, we numerically solve Eq. (S-77) for
$C_{m,n}^{(0)}(t),C_{m,n}^{(1)}(t),...,C_{m,n}^{(\alpha)}(t)$, by the 4th-
order Runge-Kutta method, obtaining the $\alpha$th-order two-point correlator
as

$\displaystyle\bar{C}_{m,n}(t)\coloneqq\sum_{\beta=0}^{\alpha}\Delta^{\beta}C_{m,n}^{(\beta)}(t).$
(S-78)

Then, the variance $\sigma_{M}^{\prime}(t)^{2}$ defined in the main text can
be computed using

$\displaystyle\braket{\hat{N}_{M}^{\prime}}=\sum_{m=-M/2}^{M/2-1}\bar{C}_{m,m}(t),$
(S-79) $\displaystyle\braket{\hat{N}_{M}^{\prime
2}}=\braket{\hat{N}_{M}^{\prime}}+\braket{\hat{N}_{M}^{\prime}}^{2}-\sum_{m=-M/2}^{M/2-1}\sum_{n=-M/2}^{M/2-1}|\bar{C}_{m,n}(t)|^{2}.$
(S-80)

Figure S-2 displays the numerical results of the variance for the 2nd-, 4th-,
and 6th-order perturbative calculations with $\Delta=1$. One can see that the
agreement of the perturbative calculations with the exact one becomes better
as the order $\alpha$ increases.

Figure S-2:  Numerical results of the variance for the perturbative
calculations with $\Delta=1$. The parameters are the same as those used in
Fig. 3 of the main text. The pentagon, square, and triangle markers denote the
data of the variance $\sigma^{\prime}_{M}(t)^{2}$ for the $\alpha=2,4,$ and
$6$th-order perturbative calculations, respectively, and the circle marker
does the exact numerical data. Here, the exact result is numerically obtained
without using the perturbation calculation.