# Lieb-Robinson bound and almost-linear light-cone in interacting boson
systems

Tomotaka Kuwahara1 and Keiji Saito2 1 Mathematical Science Team, RIKEN Center
for Advanced Intelligence Project (AIP),1-4-1 Nihonbashi, Chuo-ku, Tokyo
103-0027, Japan  2Department of Physics, Keio University, Yokohama 223-8522,
Japan

###### Abstract

In this work, we investigate how quickly local perturbations propagate in
interacting boson systems with Bose-Hubbard-type Hamiltonians. In general,
these systems have unbounded local energies, and arbitrarily fast information
propagation may occur. We focus on a specific but experimentally natural
situation in which the number of bosons at any one site in the unperturbed
initial state is approximately limited. We rigorously prove the existence of
an almost-linear information-propagation light-cone, thus establishing a
Lieb–Robinson bound: the wave-front grows at most as $t\log^{2}(t)$. We prove
the clustering theorem for gapped ground states and study the time complexity
of classically simulating one-dimensional quench dynamics, a topic of great
practical interest.

Introduction.— In non-relativistic quantum many-body systems, the speed limit
of information propagation is characterized by the Lieb–Robinson bound [1, 2,
3], an effective light-cone outside which the amount of transferred
information rapidly decays with distance. In standard spin models such as the
transverse Ising model, the light-cone is linear over time and characterized
by the Lieb–Robinson velocity, which depends only on the system details. As a
fundamental restriction applied to generic many-body quantum systems, the
Lieb–Robinson bound has been utilized to establish the clustering theorem on
bi-partite correlations in ground states [4, 5, 6] and an efficient
classical/quantum algorithm to simulate quantum many-body dynamics [7, 8, 9].
It has featured in many fields of quantum many-body physics including
condensed matter theory [10, 11, 12, 13, 14, 15, 16], statistical mechanics
[17, 18, 19, 20, 21, 22, 23], high-energy physics [24, 25, 26, 27, 28, 29,
30], and quantum information [31, 32, 33, 34, 35].

The Lieb–Robinson bound and the existence of a linear light-cone are well-
understood under the following two conditions [36, 5, 6, 37, 3]: i) the
interaction is short-range, and ii) the Hamiltonian is locally bounded. If
either of these conditions is broken, as often happens in real-world quantum
systems, the shape of the linear light-cone becomes quite complicated. When
there are long-range interactions, breaking the first condition, a
comprehensive characterization of the shape of the light-cone has been
achieved [38, 39, 28, 40, 41, 42, 29]. However, it remains challenging to
clarify the Lieb–Robinson bound when the second condition breaks down.

Quantum boson systems are representative examples of the breakdown of this
second condition with locally unbounded Hamiltonians. The difficulty lies in
the fact that the standard approach for the Lieb–Robinson bound necessarily
results in a Lieb–Robinson velocity proportional to the norm of the local
energy. When $N$ bosons clump at a single location, the on-site energy can be
as large as ${\rm poly}(N)$, leading to an infinite Lieb–Robinson velocity as
$N\to\infty$. Even though it is quite unlikely that many bosons will clump
together in realistic experiments, the theoretical possibility of such
situations must be taken into account. If harmonic and anharmonic systems [43,
44, 45, 46, 47] and spin boson models [48, 49, 50] are considered, the
Lieb–Robinson bound with the linear light-cone has been established. However,
we have no hope of unconditionally proving the existence of a Lieb–Robinson
bound without restricting the form of Hamiltonians or initial states. (In Ref.
[51], Eisert and Gross provided 1D quantum boson systems with nearest-neighbor
interactions, inducing an exponential speed of information propagation.)

Recent experiments have focused on interacting bosonic systems of the Bose-
Hubbard type [52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66],
which typically appear in cold atom setups. Since the earliest experiments on
the Lieb–Robinson bound [55, 56], there have been many attempts to clarify
information propagation in these models rigorously. However, with a few
exceptions [67, 68], establishing the Lieb–Robinson bound in Bose-Hubbard-type
models remains an open problem. A previous rigorous study [67] showed that
initially concentrated bosons in the vacuum spread at a finite speed. In Ref.
[68], the Lieb–Robinson velocity was qualitatively improved from
$\mathcal{O}(N)$ to $\mathcal{O}(\sqrt{N})$ (still infinitely large in the
limit of $N\to\infty$), where $N$ is the total number of bosons. On the other
hand, numerical calculations and theoretical case studies indicate that a
linear light-cone should be observed in practical settings such as quench
dynamics [69, 70, 71, 72, 73, 74, 75, 76]. The most natural condition is to
require a finite number of bosons at any one site in the initial state, for
example, a Mott state. However, this condition can break down over time, and a
large bias in the boson distribution may cause an unexpected acceleration of
information propagation [67]. Until now, no theoretical tools have been
developed to overcome this obstacle.

In this work, we establish the Lieb–Robinson bound with an almost linear
light-cone when a local perturbation is added to quantum states that are
initially time-independent and have low boson density [see the condition (6)].
Our Lieb–Robinson bound characterizes a wave-front that propagates as
$t\log^{2}(t)$ with time. As a practical application, we derive the clustering
theorem for non-critical ground states by extending the technique in [4, 5,
6]. In addition, we extend our theory to analyze the time complexity of
computing quantum dynamics by quenching the Hamiltonian parameter, a topic of
major research interest [77, 69, 78, 79, 71, 72, 73, 75, 80, 81, 82, 70, 83,
84, 85, 72, 86, 87, 88, 89, 90, 74, 75, 91, 92, 93, 76]. We rigorously
establish the time complexity of $e^{t\log^{3}(t)}$ to simulate local quench
dynamics for one-dimensional Bose-Hubbard-type Hamiltonians.

Setup and main result.— We consider a quantum system on a finite-dimensional
lattice (graph), where bosons interact with each other. An unbounded number of
bosons can sit on each of the sites, and the local Hilbert dimension is thus
infinitely large. We denote by $\Lambda$ the set of all sites on the lattice.
For an arbitrary partial set $X\subseteq\Lambda$, we denote the cardinality
(the number of sites contained in $X$) by $|X|$. For arbitrary subsets
$X,Y\subseteq\Lambda$, we define $d_{X,Y}$ as the shortest path-length on the
graph that connects $X$ and $Y$. For a subset $X\subseteq\Lambda$, we define
the extended subset $X[r]$ by length $r$ as

$\displaystyle X[r]:=\\{i\in\Lambda|d_{X,i}\leq r\\},$ (1)

where $X[0]=X$ and $r$ is an arbitrary positive number (i.e.,
$r\in\mathbb{R}^{+}$).

We define $b_{i}$ and $b_{i}^{\dagger}$ as the annihilation and creation
operators of the boson, respectively. We also define
$\hat{n}_{i}:=b_{i}^{\dagger}b_{i}$ as the number operator of bosons on site
$i$. We consider a Hamiltonian of the form

$\displaystyle H:=\sum_{\langle i,j\rangle}J_{i,j}(b_{i}b_{j}^{\dagger}+{\rm
h.c.})+\sum_{Z\subset\Lambda:|Z|\leq k}v_{Z},$ (2)

where $|J_{i,j}|\leq\bar{J}$ and $\sum_{\langle i,j\rangle}$ denotes summation
over all pairs of adjacent sites $\\{i,j\\}$ on the lattice. Here, $v_{Z}$
consists of finite-range boson-boson interactions on subset $Z$. We now assume
that $v_{Z}$ is given as a function of the number operators
$\\{\hat{n}_{i}\\}_{i\in Z}$. The simplest example is the Bose-Hubbard model:

$\displaystyle H=\sum_{\langle i,j\rangle}J(b_{i}b_{j}^{\dagger}+{\rm
h.c.})+\frac{U}{2}\sum_{i\in\Lambda}\hat{n}_{i}(\hat{n}_{i}-1)-\mu\sum_{i\in\Lambda}\hat{n}_{i},$

where $U$ and $\mu$ are $\mathcal{O}(1)$ constants. For an arbitrary operator
$O$, the time-evolution due to another operator $A$ is

$\displaystyle O(A,t):=e^{iAt}Oe^{-iAt}.$ (3)

(We abbreviate $O(H,t)$ as $O(t)$ for simplicity.)

Let $\rho_{0}$ be a time-independent quantum state, i.e., $[\rho_{0},H]=0$. We
consider propagation of a local perturbation to $\rho_{0}$ such as $\rho\to
O_{i_{0}}\rho_{0}O_{i_{0}}^{\dagger},$ where $i_{0}\in\Lambda$ and $O_{i_{0}}$
can take the form of a projection onto site $i_{0}$. We are interested in how
fast this perturbation propagates. Mathematically, after the time evolution,
$\rho(t)$ is given by $O_{i_{0}}(t)\rho_{0}O_{i_{0}}(t)^{\dagger}$. Thus, we
must estimate the approximation error of

$\displaystyle O_{i_{0}}(t)\rho_{0}\approx O^{(t)}_{i_{0}[R]}\rho_{0},$ (4)

where $O^{(t)}_{i_{0}[R]}$ is an appropriate operator supported on subset
$i_{0}[R]$ [see the notation (1)]. Our main result concerns the approximation
error for finite $R$ (see Sec. S.II. in Supplementary materials [94] for the
formal expression).

Following Ref. [41], we define the shape of the light cone in the following
sense. We say that the Hamiltonian dynamics $e^{-iHt}$ have an effective light
cone with velocity $v_{t,\delta}$ if the following inequality holds for an
arbitrary error $\delta\in\mathbb{R}$ and $t$:

$\displaystyle\|[O_{i_{0}}(t)-O^{(t)}_{i_{0}[R]}]\|\leq\delta\left\|O_{i_{0}}\right\|\quad\textrm{for}\quad
R\geq v_{t,\delta}|t|.$ (5)

When $v_{t,\delta}$ converges to a finite value for $t\to\infty$ (i.e.,
$v_{\infty,\delta}={\rm const.}$), we say that the effective light cone is
linear. From the definition, the amount of information propagation is smaller
than $\delta$ outside the region separated by the distance $v_{t,\delta}|t|$.

Main Theorem. Let us assume that the number of boson creations by $O_{i_{0}}$
is finitely bounded. Then, for an arbitrary time-independent quantum state
$\rho_{0}$ satisfying the low-boson-density condition

$\displaystyle\max_{i\in\Lambda}{\rm
tr}(e^{c_{0}(\hat{n}_{i}-\bar{q})}\rho_{0})\leq 1\quad c_{0}\leq 1,$ (6)

we can approximate $O_{i_{0}}(t)\rho_{0}$ by another operator
$O^{(t)}_{i_{0}[R]}$ supported on $i_{0}[R]$ with the following approximation
error:

$\displaystyle\left\|\left(O_{i_{0}}(t)-O^{(t)}_{i_{0}[R]}\right)\rho_{0}\right\|$
$\displaystyle\leq\|O_{i_{0}}\|\exp\left(c_{0}\bar{q}-C_{1}\frac{R}{t\log(R)}+C_{2}\log(R)\right),$
(7)

where $t\geq 1$, and $C_{1}$ and $C_{2}$ are constants of $\mathcal{O}(1)$
that are independent of $\bar{q}$ and only depend on the details of the
system. For a general operator $O_{X_{0}}$, we can obtain a similar inequality
by slightly changing (Lieb-Robinson bound and almost-linear light-cone in
interacting boson systems).

Condition (6) ensures that the probability for many bosons to be concentrated
on one site is exponentially small in the initial state $\rho_{0}$. We notice
that the condition can break down as time increases. By applying the
inequality (Lieb-Robinson bound and almost-linear light-cone in interacting
boson systems) to (5), we obtain
$v_{t,\delta}\propto\log^{2}(t)[\log(1/\delta)+c_{0}\bar{q}]$. Hence,
information propagation is restricted in the region that is separated from
$i_{0}$ by at most $\mathcal{O}(\bar{q})t\log^{2}(t)$. Therefore, we can
ensure that the acceleration of information propagation observed in Ref. [51]
cannot occur in our model, because the speed of information becomes at most
polylogarithmically large with time, i.e., $\leq\log^{2}(t)$.

Clustering Theorem.— As an immediate application of the main theorem, we
consider the exponential decay of bi-partite correlations in gapped ground
states, i.e., the clustering theorem. Here, we denote the non-degenerate
ground state by $|E_{0}\rangle$ and the spectral gap by $\Delta E$. We prove
an upper-bound on the correlation function ${\rm Cor}(O_{X},O_{Y}):=\langle
E_{0}|O_{X}O_{Y}|E_{0}\rangle-\langle E_{0}|O_{X}|E_{0}\rangle\langle
E_{0}|O_{Y}|E_{0}\rangle$, where $O_{X}$ and $O_{Y}$ are operators supported
on $X$ and $Y$. For simplicity, we let $\bar{q}=\mathcal{O}(1)$. Then, the
following inequality holds if $|E_{0}\rangle$ satisfies condition (6) (see
Sec. S.III. in Supplementary materials [94]):

$\displaystyle{\rm Cor}(O_{X},O_{Y})\leq
C_{3}\|O_{X}\|\cdot\|O_{Y}\|\exp\left(-\sqrt{\frac{C_{3}^{\prime}\Delta
E}{\log(R)}R}\right),$ (8)

where $C_{3}$, $C_{3}^{\prime}$ and $C_{3}^{\prime\prime}$ are
$\mathcal{O}(1)$ constants. From the inequality, the bi-partite correlations
decay beyond $R\approx\tilde{\mathcal{O}}(1/\Delta E)$. This sub-exponential
decay, which is weaker than the exponential decay described in Ref. [4, 5, 6],
is a consequence of the asymptotic form of $e^{-\mathcal{O}(R/(t\log R))}$ in
our Lieb–Robinson bound (Lieb-Robinson bound and almost-linear light-cone in
interacting boson systems).

Application to quench dynamics.— We next consider the application of our
results to quench dynamics, the most popular setup in the study of non-
equilibrium quantum systems. Here, a system is initially prepared in a steady
state $\rho_{0}$ (e.g., the ground state), and then evolves unitarily in time
under the sudden change of the Hamiltonian $H\to H^{\prime}$. We consider the
case where the Hamiltonian $H^{\prime}$ is given by $H^{\prime}=H+h_{X_{0}}$,
where we assume $H^{\prime}$ still has the form of Eq. (2). In addition, the
interaction $h_{X_{0}}$ includes only polynomials of finite degree in
$\\{\hat{n}_{i}\\}_{i\in X_{0}}$, such as $\hat{n}_{i}^{2}$ and
$\hat{n}_{i}^{2}\hat{n}_{j}^{3}$, etc.

Our purpose is to find an appropriate unitary operator $U_{i_{0}[R]}$
supported on $i_{0}[R]$ that gives $\rho_{0}(H^{\prime},t)\approx
U_{i_{0}[R]}\rho_{0}U^{\dagger}_{i_{0}[R]}$. We can prove the following
theorem (see Sec. S.IX. in the Supplementary materials [94] for details):

Quench theorem. For initial state $\rho_{0}$ with the conditions
$[\rho_{0},H]=0$ and (6), we have

$\displaystyle\left\|\rho_{0}(H^{\prime},t)-U_{i_{0}[R]}\rho_{0}U^{\dagger}_{i_{0}[R]}\right\|_{1}$
$\displaystyle\leq\exp\left(c_{0}\bar{q}-C^{\prime}_{1}\frac{(R-r_{0})}{t\log(R)}+C^{\prime}_{2}\log(R)\right),$
(9)

where we define $r_{0}$ such that $X_{0}\in i_{0}[r_{0}]$ for an appropriate
$i_{0}\in\Lambda$, and $C^{\prime}_{1}$ and $C^{\prime}_{2}$ are constants of
$\mathcal{O}(1)$ that are independent of $\bar{q}$ and only depend on the
details of the system. Moreover, the computational cost of constructing the
unitary operator $U_{i_{0}[R]}$ is at most
$\exp\left[\mathcal{O}\left(R^{D}\log(R)\right)\right].$

This theorem immediately gives the following corollary on the time complexity
of preparing $U_{i_{0}[R]}$:
Corollary. The computational cost of calculating the quench dynamics on 1D
chains up to an error $\epsilon$ is at most

$\displaystyle\exp\left[t\log^{3}(t)+t\log(1/\epsilon)\log\log^{2}(1/\epsilon)\right],$
(10)

where we assume $r_{0}=\mathcal{O}(1)$ and $\bar{q}=\mathcal{O}(1)$. When the
error $\epsilon$ is fixed, we have a time complexity of $e^{t\log^{3}(t)}$.
This is the first rigorous result on the efficiency of the classical
simulation of interacting boson systems.

Proof of the main theorem.— For the proof, we connect the Lieb–Robinson bounds
for small time evolutions step by step, based on previous analyses of the
Lieb–Robinson bound in long-range interacting systems [95, 29]. The great
merit of this approach is that we have to derive the Lieb–Robinson bound only
for short-time evolution. We decompose the total time $t$ into $m_{t}$ pieces
and define $\Delta t:=t/m_{t}$ with $m_{t}=\mathcal{O}(t)$. Note that we can
make $\Delta t$ arbitrarily small by making $m_{t}$ sufficiently large. For a
fixed $R$, we define the subset $X_{m}$ as follows:

$\displaystyle X_{m}:=i_{0}[m\Delta r],\quad\Delta r=\left\lfloor
R/m_{t}\right\rfloor,$

where $X_{m}=X_{0}[m\Delta r]$ and $X_{m_{t}}\subseteq i_{0}[R]$.

We connect the step-by-step approximations of the short-time evolution to
reach the final approximation. Under the assumption of the time invariance of
$\rho_{0}$ (i.e., $\rho_{0}(t)=\rho$), we can derive the following inequality
[29]:

$\displaystyle\left\|\left[O_{i_{0}}(m_{t}\Delta
t)-O_{X_{m_{t}}}^{(m_{t})}\right]\rho_{0}\right\|_{1}$
$\displaystyle\leq\sum_{m=1}^{m_{t}}\left\|\left[O_{X_{m-1}}^{(m-1)}(\Delta
t)-O_{X_{m}}^{(m)}\right]\rho_{0}\right\|_{1},$ (11)

where $O_{X_{0}}^{(0)}=O_{X_{0}}$, and $O_{X_{m}}^{(m)}$ is recursively
defined by approximating $O_{X_{m}}^{(m)}(\Delta t)$. When $\rho_{0}$ depends
on the time, a severe modification is required in the inequality (Lieb-
Robinson bound and almost-linear light-cone in interacting boson systems) (see
Sec. IV. B in Supplementary materials [94]). In order to reduce (Lieb-Robinson
bound and almost-linear light-cone in interacting boson systems) to the main
inequality (Lieb-Robinson bound and almost-linear light-cone in interacting
boson systems), we need to obtain

$\displaystyle O_{X_{m-1}}^{(m-1)}(\Delta t)\rho_{0}\approx
U_{X_{m}}^{(m)\dagger}O_{X_{m-1}}^{(m-1)}U_{X_{m}}^{(m)}\rho_{0}=O_{X_{m}}^{(m)}\rho_{0},$
(12)

by using an appropriate unitary operator supported on $X_{m}$.

Therefore, our primary task is to estimate the approximation error of (12),
which gives the Lieb–Robinson bound for the short time $\Delta t$. We can
prove that, for a general operator $O_{X}$ with $X\subseteq i[r]$
($i\in\Lambda$), there exists a unitary operator $U_{X[\ell]}$ supported on
$X[\ell]$ such that

$\displaystyle\left\|\left(O_{X}(t)-U_{X[\ell]}^{\dagger}O_{X}U_{X[\ell]}\right)\rho_{0}\right\|_{1}$
$\displaystyle\leq\|O_{X}\|e^{c_{0}\bar{q}-\ell/\log(r)+C_{0}\log(r)},$ (13)

for $t\leq\Delta t_{0}$ (see Subtheorem 1 in Supplementary materials [94]),
where $C_{0}$ and $\Delta t_{0}$ are $\mathcal{O}(1)$ constants. We here
choose the time width $\Delta t$ such that $\Delta t\leq\Delta t_{0}$. By
using the inequality (Lieb-Robinson bound and almost-linear light-cone in
interacting boson systems) with $\ell=\Delta r$ and $t=\Delta t$, we can
reduce the inequality (Lieb-Robinson bound and almost-linear light-cone in
interacting boson systems) to the desired form (Lieb-Robinson bound and
almost-linear light-cone in interacting boson systems) by choosing $C_{1}$ and
$C_{2}$ appropriately. This completes the proof of the main theorem. $\square$

Figure 1: Boson density after time evolution. In Ref. [67], all bosons were
initially concentrated in a region $X$ on the vacuum and were shown to spread
beyond it with a finite velocity. However, if there is initially a finite
number of bosons outside $X$, the upper bound of the boson number increases
exponentially with $t$. This spoils the approach of Ref. [67] in general
setups for long-time evolution. Our approach only considers the short time
$t=\mathcal{O}(1)$, when the exponential increase $e^{\mathcal{O}(t)}$ is
still $\mathcal{O}(1)$. We then ensure that the boson number distribution for
$\hat{n}_{i}$ exponentially decays if the site $i$ of interest is sufficiently
separated from the region $X$.

Short-time Lieb–Robinson bound.— We have seen that the bosonic Lieb–Robinson
bound can be immediately derived if we can prove the inequality (Lieb-Robinson
bound and almost-linear light-cone in interacting boson systems), which
includes all the difficulties in our proof. We will now provide a sketch of
the proof; a fuller and more formal presentation can be found in the
Supplementary materials [94] (Secs. S.V., S.VI., S.VII. and S.VIII.).

We first consider the boson density after short-time evolution (see Sec. S.VI.
in Supplementary materials [94]). For this purpose, we need to estimate

$\displaystyle{\rm
tr}\left[\hat{n}_{i}^{s}\tilde{\rho}(t)\right],\quad\tilde{\rho}(t)=e^{-iHt}O_{X}\rho_{0}O_{X}^{\dagger}e^{iHt},$
(14)

with $s\in\mathbb{N}$. This quantity characterizes the influence of the
perturbation $O_{X}$ on the boson density after time evolution. In the state
$\tilde{\rho}(0)$, the boson number $\hat{n}_{i}$ ($i\notin X$) is
exponentially suppressed because of condition (6), while the bosons may be
highly concentrated in the region $X$. Time evolution will cause these
concentrated bosons to spread outside $X$ (see Fig. 1).

In order to characterize the dynamics of the bosons, we utilize the method in
Ref. [67]. We can prove that

$\displaystyle\frac{{\rm
tr}\left[\hat{n}_{i}^{s}\tilde{\rho}(t)\right]}{\|O_{X}\|^{2}}\leq
c^{\prime}_{1}e^{c_{0}\bar{q}}|X|^{3}(c_{1}s|X|)^{s}e^{-d_{i,X}}+c^{\prime\prime}_{1}e^{c_{0}\bar{q}}(c_{1}s)^{s},$
(15)

where $c_{1}$, $c_{1}^{\prime}$, and $c_{1}^{\prime\prime}$ depend on the time
as $e^{\mathcal{O}(t)}$. The above upper bound induces an exponential increase
of the boson density with time; hence we cannot use it for arbitrarily large
$t$. However, the key point of our proof-method is that we only need to treat
the short-time evolution, where the coefficients $c_{1}$, $c_{1}^{\prime}$,
and $c_{1}^{\prime\prime}$ are $\mathcal{O}(1)$ constants. By using Markov’s
inequality, we can ensure that the probability distribution of the boson
number $\hat{n}_{i}$ obeys

$\displaystyle P_{i,\geq z_{0}}^{(t)}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\|O_{X}\|^{2}\left(\frac{\tilde{c}_{1}d_{i,X}}{z_{0}}\right)^{\tilde{c}^{\prime}_{1}d_{i,X}/\log(r)}$
(16)

under the condition $d_{i,X}\gtrsim\log(r)$, where $P_{i,\geq z_{0}}^{(t)}$ is
the probability that $z_{0}$ or more bosons are observed at the site $i$.
(Recall that by definition $X\subseteq i[r]$.) Finally, we remark that it is
essential to the proof that the Hamiltonian be the form (2); if the
Hamiltonian includes interactions such as
$\hat{n}_{i}\hat{n}_{j}b_{i}b_{j}^{\dagger}$, the inequality (15) may break
down even for small $t$.

Figure 2:  Schematic picture of the region where the boson number is
truncated. In the region $L_{1}$ ($=X[\ell_{0}]$), we cannot restrict the
boson number distribution. On the other hand, as long as $\ell_{0}$ is
sufficiently large, the boson number outside $L_{1}$ can be truncated up to a
finite value $q$. We perform the boson number truncation in the shaded region
$\tilde{L}=L_{2}\setminus L_{1}$ with $L_{2}=X[2\ell_{0}]$. By using an
effective Hamiltonian $\tilde{H}[\tilde{L},q]$ as in Eq. (18), we can
approximate the exact dynamics by choosing the appropriate $q$ [i.e.,
$\eta\ell_{0}$; see the inequality (Lieb-Robinson bound and almost-linear
light-cone in interacting boson systems)].

In the second technique, we construct an effective Hamiltonian that has
bounded local energy in a specific region and approximates the exact dynamics
(see Sec. S.VII. in Supplementary materials [94]). The inequality (16) implies
that the boson number $\hat{n}_{i}$ is strongly suppressed when the site $i$
is sufficiently separated from the region $X$. Hence, we expect that, in the
original Hamiltonian $H$, the maximum boson number at one site can be
truncated during short-time evolution. We first define two regions
$L_{1}:=X[\ell_{0}]$ and $L_{2}:=X[2\ell_{0}]$, where the length $\ell_{0}$ is
appropriately chosen. We then consider the boson truncation in the region
$\tilde{L}$ which is defined as (see Fig. 2)

$\displaystyle\tilde{L}:=L_{2}\setminus L_{1}.$ (17)

We now define $\bar{\Pi}_{\tilde{L},q}$ as the projection onto the eigenspace
such that the boson number $\hat{n}_{i}$ ($\forall i\in\tilde{L}$) is
truncated up to $q$, i.e., $\|\hat{n}_{i}\bar{\Pi}_{\tilde{L},q}\|\leq q$. We
then approximate the time-evolution operator $e^{-iHt}$ by using an effective
Hamiltonian $\tilde{H}[\tilde{L},q]$, defined by

$\displaystyle\tilde{H}[\tilde{L},q]:=\bar{\Pi}_{\tilde{L},q}H\bar{\Pi}_{\tilde{L},q},$
(18)

with a bounded local energy in the region $\tilde{L}$. In general, the time
evolution $O_{X}(t)$ cannot be approximated by
$O_{X}(\tilde{H}[\tilde{L},q],t)$ at all, where we have used the notation (3).
However, we are only interested in the norm difference between
$O_{X}(t)\rho_{0}$ and $O_{X}(\tilde{H}[\tilde{L},q],t)\rho_{0}$. We can prove

$\displaystyle\left\|\left[O_{X}(t)-O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\right]\rho_{0}\right\|_{1}$
$\displaystyle\leq\frac{\|O_{X}\|}{2}e^{c_{0}\bar{q}}e^{-2\ell_{0}/\log(r)}$
(19)

for $q=\eta\ell_{0}$ and $\ell_{0}\geq C_{0}\log^{2}(r)$, where $\eta$ and
$C_{0}$ are $\mathcal{O}(1)$ constants which are independent of $\bar{q}$, and
$r$ has been defined by $X\subseteq i[r]$. From this upper bound, we can see
that the error exponentially decreases with the number of the boson
truncation. Thus, by using the Hamiltonian
$\tilde{H}[\tilde{L},\eta\ell_{0}]$, the biggest obstacle, namely the
unboundedness of the interaction norms, has been removed, at least in the
region $\tilde{L}$. However, outside this region, the norm is still unbounded.
We thus need to consider how to derive the Lieb–Robinson bound for
$e^{-i\tilde{H}[\tilde{L},\eta\ell_{0}]t}$ only from the finiteness of the
Hamiltonian norm in the region $\tilde{L}$.

Our final task is to approximate the time evolution
$O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)$ by
$U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}$, where $U_{L_{2}}$ is an appropriate
unitary operator supported on the subset $L_{2}$ (see Sec. S.VIII. in
Supplementary materials [94]). By a careful calculation based on the standard
approach to deriving the Lieb–Robinson bound, we can show that the
approximation error obeys

$\displaystyle\left\|O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right\|\leq\frac{\|O_{X}\|}{2}e^{-2\ell_{0}/\log(r)},$
(20)

assuming $t\leq\Delta t_{0}$ with $\Delta t_{0}$ an $\mathcal{O}(1)$ constant.
Therefore, under the conditions $\ell_{0}\geq C_{0}\log^{2}(r)$ and
$t\leq\Delta t_{0}$, we have the inequalities (Lieb-Robinson bound and almost-
linear light-cone in interacting boson systems) and (20), which together yield
the desired inequality (Lieb-Robinson bound and almost-linear light-cone in
interacting boson systems) since $e^{c_{0}\bar{q}}\geq 1$.

Conclusion.— In this work, we have established the Lieb–Robinson bound (Lieb-
Robinson bound and almost-linear light-cone in interacting boson systems) with
an almost-linear light cone $R\propto t\log^{2}(t)$ for arbitrary initial
steady states under the condition (6). Our bound leads to the clustering
theorem (8) for gapped ground states and the efficient simulation of the
quench dynamics as in (10). Our result gives the first rigorous
characterization of the light cone of interacting boson systems under
experimentally realistic conditions.

Nevertheless, this Lieb–Robinson bound might be further improved. First, the
asymptotic form $e^{-R/(t\log R)}$ in (Lieb-Robinson bound and almost-linear
light-cone in interacting boson systems) could be changed to $e^{-R+vt}$,
which would induce a strictly linear light cone for information propagation.
Second, there remains the challenge to clarify the class of quantum states
that rigorously satisfy the assumption (6). Third, regarding the time
independence of $\rho_{0}$, we conjecture that an information wave-front of at
least polynomial form (i.e., $R\propto t^{\zeta}$, $\zeta\geq 1$) can be
derived when $\rho_{0}$ is time-dependent. Although our current techniques
cannot immediately accommodate these improvements, we hope to develop a better
Lieb–Robinson bound for interacting bosons in the future.

Note added. For the readers information, we would like to refer to a
subsequent study by Yin and Lucas [96], which proves the linear light cone for
interacting boson systems in another specific setup.

###### Acknowledgements.

The work of T. K. was supported by the RIKEN Center for AIP and JSPS KAKENHI
(Grant No. 18K13475). TK gives thanks to God for his wisdom. K.S. was
supported by JSPS Grants-in-Aid for Scientific Research (JP16H02211 and
JP19H05603).

## References

  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, Communications in Mathematical Physics 28, 251 (1972).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-Robinson Bounds and the Generation of Correlations and Topological Quantum Order, Phys. Rev. Lett. 97, 050401 (2006).
  * Nachtergaele  _et al._ [2019] B. Nachtergaele, R. Sims, and A. Young, Quasi-locality bounds for quantum lattice systems. I. Lieb-Robinson bounds, quasi-local maps, and spectral flow automorphisms, Journal of Mathematical Physics 60, 061101 (2019).
  * Hastings [2004a] M. B. Hastings, Locality in Quantum and Markov Dynamics on Lattices and Networks, Phys. Rev. Lett. 93, 140402 (2004a).
  * Hastings and Koma [2006] M. B. Hastings and T. Koma, Spectral Gap and Exponential Decay of Correlations, Communications in Mathematical Physics 265, 781 (2006).
  * Nachtergaele and Sims [2006] B. Nachtergaele and R. Sims, Lieb-Robinson Bounds and the Exponential Clustering Theorem, Communications in Mathematical Physics 265, 119 (2006).
  * Osborne [2006] T. J. Osborne, Efficient Approximation of the Dynamics of One-Dimensional Quantum Spin Systems, Phys. Rev. Lett. 97, 157202 (2006).
  * Haah  _et al._ [2018] J. Haah, M. Hastings, R. Kothari, and G. H. Low, Quantum Algorithm for Simulating Real Time Evolution of Lattice Hamiltonians, in  _2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS)_ (2018) pp. 350–360.
  * Tran  _et al._ [2019] M. C. Tran, A. Y. Guo, Y. Su, J. R. Garrison, Z. Eldredge, M. Foss-Feig, A. M. Childs, and A. V. Gorshkov, Locality and Digital Quantum Simulation of Power-Law Interactions, Phys. Rev. X 9, 031006 (2019).
  * Hastings [2004b] M. B. Hastings, Lieb-Schultz-Mattis in higher dimensions, Phys. Rev. B 69, 104431 (2004b).
  * Nachtergaele and Sims [2007] B. Nachtergaele and R. Sims, A Multi-Dimensional Lieb-Schultz-Mattis Theorem, Communications in Mathematical Physics 276, 437 (2007).
  * Hastings and Wen [2005] M. B. Hastings and X.-G. Wen, Quasiadiabatic continuation of quantum states: The stability of topological ground-state degeneracy and emergent gauge invariance, Phys. Rev. B 72, 045141 (2005).
  * Hastings [2007a] M. B. Hastings, An area law for one-dimensional quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007a), arXiv:0705.2024 .
  * Bravyi  _et al._ [2010] S. Bravyi, M. B. Hastings, and S. Michalakis, Topological quantum order: Stability under local perturbations, Journal of Mathematical Physics 51, 093512 (2010).
  * Haegeman  _et al._ [2013] J. Haegeman, S. Michalakis, B. Nachtergaele, T. J. Osborne, N. Schuch, and F. Verstraete, Elementary Excitations in Gapped Quantum Spin Systems, Phys. Rev. Lett. 111, 080401 (2013).
  * Hastings and Michalakis [2014] M. B. Hastings and S. Michalakis, Quantization of Hall Conductance for Interacting Electrons on a Torus, Communications in Mathematical Physics 334, 433 (2014).
  * Müller  _et al._ [2015] M. P. Müller, E. Adlam, L. Masanes, and N. Wiebe, Thermalization and Canonical Typicality in Translation-Invariant Quantum Lattice Systems, Communications in Mathematical Physics 340, 499 (2015).
  * Damanik  _et al._ [2014] D. Damanik, M. Lemm, M. Lukic, and W. Yessen, New anomalous lieb-robinson bounds in quasiperiodic $xy$ chains, Phys. Rev. Lett. 113, 127202 (2014).
  * Iyoda  _et al._ [2017] E. Iyoda, K. Kaneko, and T. Sagawa, Fluctuation Theorem for Many-Body Pure Quantum States, Phys. Rev. Lett. 119, 100601 (2017).
  * Kuwahara  _et al._ [2016] T. Kuwahara, T. Mori, and K. Saito, Floquet-Magnus theory and generic transient dynamics in periodically driven many-body quantum systems, Annals of Physics 367, 96 (2016).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. m. c. Huveneers, Effective Hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Phys. Rev. B 95, 014112 (2017).
  * Bachmann  _et al._ [2017] S. Bachmann, W. De Roeck, and M. Fraas, Adiabatic Theorem for Quantum Spin Systems, Phys. Rev. Lett. 119, 060201 (2017).
  * Kuwahara  _et al._ [2021] T. Kuwahara, A. M. Alhambra, and A. Anshu, Improved Thermal Area Law and Quasilinear Time Algorithm for Quantum Gibbs States, Phys. Rev. X 11, 011047 (2021).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 2016, 106 (2016).
  * Roberts and Swingle [2016] D. A. Roberts and B. Swingle, Lieb-Robinson Bound and the Butterfly Effect in Quantum Field Theories, Phys. Rev. Lett. 117, 091602 (2016).
  * Joshi  _et al._ [2020] M. K. Joshi, A. Elben, B. Vermersch, T. Brydges, C. Maier, P. Zoller, R. Blatt, and C. F. Roos, Quantum Information Scrambling in a Trapped-Ion Quantum Simulator with Tunable Range Interactions, Phys. Rev. Lett. 124, 240505 (2020).
  * Zhou  _et al._ [2020] T. Zhou, S. Xu, X. Chen, A. Guo, and B. Swingle, Operator Lévy Flight: Light Cones in Chaotic Long-Range Interacting Systems, Phys. Rev. Lett. 124, 180601 (2020).
  * Chen and Lucas [2019] C.-F. Chen and A. Lucas, Finite Speed of Quantum Scrambling with Long Range Interactions, Phys. Rev. Lett. 123, 250605 (2019).
  * Kuwahara and Saito [2021] T. Kuwahara and K. Saito, Absence of Fast Scrambling in Thermodynamically Stable Long-Range Interacting Systems, Phys. Rev. Lett. 126, 030604 (2021).
  * Chen [2021] C.-F. Chen, Concentration of otoc and lieb-robinson velocity in random hamiltonians (2021), arXiv:2103.09186 [quant-ph] .
  * Hastings [2007b] M. B. Hastings, Quantum belief propagation: An algorithm for thermal quantum systems, Phys. Rev. B 76, 201102 (2007b).
  * Eisert and Osborne [2006] J. Eisert and T. J. Osborne, General Entanglement Scaling Laws from Time Evolution, Phys. Rev. Lett. 97, 150404 (2006).
  * Ge  _et al._ [2016] Y. Ge, A. Molnár, and J. I. Cirac, Rapid Adiabatic Preparation of Injective Projected Entangled Pair States and Gibbs States, Phys. Rev. Lett. 116, 080503 (2016).
  * Deshpande  _et al._ [2018] A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig, and A. V. Gorshkov, Dynamical Phase Transitions in Sampling Complexity, Phys. Rev. Lett. 121, 030501 (2018).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of Trotter Error with Commutator Scaling, Phys. Rev. X 11, 011020 (2021).
  * Nachtergaele  _et al._ [2006] B. Nachtergaele, Y. Ogata, and R. Sims, Propagation of Correlations in Quantum Lattice Systems, Journal of Statistical Physics 124, 1 (2006).
  * Nachtergaele and Sims [2010] B. Nachtergaele and R. Sims, Lieb-Robinson bounds in quantum many-body physics, Contemp. Math 529, 141 (2010), arXiv:1004.2086 .
  * Eisert  _et al._ [2013] J. Eisert, M. van den Worm, S. R. Manmana, and M. Kastner, Breakdown of Quasilocality in Long-Range Quantum Lattice Models, Phys. Rev. Lett. 111, 260401 (2013).
  * Foss-Feig  _et al._ [2015] M. Foss-Feig, Z.-X. Gong, C. W. Clark, and A. V. Gorshkov, Nearly Linear Light Cones in Long-Range Interacting Quantum Systems, Phys. Rev. Lett. 114, 157201 (2015).
  * Tran  _et al._ [2020] M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo, A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov, and A. Lucas, Hierarchy of Linear Light Cones with Long-Range Interactions, Phys. Rev. X 10, 031009 (2020).
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions, Phys. Rev. X 10, 031010 (2020).
  * Tran  _et al._ [2021] M. C. Tran, A. Y. Guo, A. Deshpande, A. Lucas, and A. V. Gorshkov, Optimal State Transfer and Entanglement Generation in Power-Law Interacting Systems, Phys. Rev. X 11, 031016 (2021).
  * Cramer  _et al._ [2008a] M. Cramer, A. Serafini, and J. Eisert, Locality of dynamics in general harmonic quantum systems (2008a), arXiv:0803.0890 [quant-ph] .
  * Nachtergaele  _et al._ [2009] B. Nachtergaele, H. Raz, B. Schlein, and R. Sims, Lieb-Robinson Bounds for Harmonic and Anharmonic Lattice Systems, Communications in Mathematical Physics 286, 1073 (2009).
  * Raz and Sims [2009] H. Raz and R. Sims, Estimating the Lieb-Robinson Velocity for Classical Anharmonic Lattice Systems, Journal of Statistical Physics 137, 79 (2009).
  * Nachtergaele  _et al._ [2010] B. Nachtergaele, B. Schlein, R. Sims, S. Starr, and Z. Valentin, On the Existence of the Dynamics for Anharmonic Quantum Oscillator Systems, Reviews in Mathematical Physics 22, 207 (2010).
  * Nachtergaele and Sims [2014] B. Nachtergaele and R. Sims, On the dynamics of lattice systems with unbounded on-site terms in the Hamiltonian (2014), arXiv:1410.8174 [math-ph] .
  * Jünemann  _et al._ [2013] J. Jünemann, A. Cadarso, D. Pérez-García, A. Bermudez, and J. J. García-Ripoll, Lieb-Robinson Bounds for Spin-Boson Lattice Models and Trapped Ions, Phys. Rev. Lett. 111, 230404 (2013).
  * Woods  _et al._ [2015] M. P. Woods, M. Cramer, and M. B. Plenio, Simulating Bosonic Baths with Error Bars, Phys. Rev. Lett. 115, 130401 (2015).
  * Woods and Plenio [2016] M. P. Woods and M. B. Plenio, Dynamical error bounds for continuum discretisation via Gauss quadrature rules—A Lieb-Robinson bound approach, Journal of Mathematical Physics 57, 022105 (2016).
  * Eisert and Gross [2009] J. Eisert and D. Gross, Supersonic Quantum Communication, Phys. Rev. Lett. 102, 240501 (2009).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Sherson  _et al._ [2010] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Single-atom-resolved fluorescence imaging of an atomic Mott insulator, Nature 467, 68 (2010).
  * Bakr  _et al._ [2010] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and M. Greiner, Probing the Superfluid–to–Mott Insulator Transition at the Single-Atom Level, Science 329, 547 (2010).
  * Cheneau  _et al._ [2012] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, Light-cone-like spreading of correlations in a quantum many-body system, Nature 481, 484 (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).
  * Braun  _et al._ [2015] S. Braun, M. Friesdorf, S. S. Hodgman, M. Schreiber, J. P. Ronzheimer, A. Riera, M. del Rey, I. Bloch, J. Eisert, and U. Schneider, Emergence of coherence and the dynamics of quantum phase transitions, Proceedings of the National Academy of Sciences 112, 3641 (2015).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Choi  _et al._ [2016] J.-y. 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).
  * Meinert  _et al._ [2016] F. Meinert, M. J. Mark, K. Lauber, A. J. Daley, and H.-C. Nägerl, Floquet Engineering of Correlated Tunneling in the Bose-Hubbard Model with Ultracold Atoms, Phys. Rev. Lett. 116, 205301 (2016).
  * Baier  _et al._ [2016] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, Extended Bose-Hubbard models with ultracold magnetic atoms, Science 352, 201 (2016).
  * Ye  _et al._ [2019] Y. Ye, Z.-Y. Ge, Y. Wu, S. Wang, M. Gong, Y.-R. Zhang, Q. Zhu, R. Yang, S. Li, F. Liang, J. Lin, Y. Xu, C. Guo, L. Sun, C. Cheng, N. Ma, Z. Y. Meng, H. Deng, H. Rong, C.-Y. Lu, C.-Z. Peng, H. Fan, X. Zhu, and J.-W. Pan, Propagation and Localization of Collective Excitations on a 24-Qubit Superconducting Processor, Phys. Rev. Lett. 123, 050502 (2019).
  * Yan  _et al._ [2019] Z. Yan, Y.-R. Zhang, M. Gong, Y. Wu, Y. Zheng, S. Li, C. Wang, F. Liang, J. Lin, Y. Xu, C. Guo, L. Sun, C.-Z. Peng, K. Xia, H. Deng, H. Rong, J. Q. You, F. Nori, H. Fan, X. Zhu, and J.-W. Pan, Strongly correlated quantum walks with a 12-qubit superconducting processor, Science 364, 753 (2019).
  * Rubio-Abadal  _et al._ [2020] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Floquet Prethermalization in a Bose-Hubbard System, Phys. Rev. X 10, 021044 (2020).
  * Yang  _et al._ [2020] B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache, J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan, Observation of gauge invariance in a 71-site Bose–Hubbard quantum simulator, Nature 587, 392 (2020).
  * Takasu  _et al._ [2020] Y. Takasu, T. Yagami, H. Asaka, Y. Fukushima, K. Nagao, S. Goto, I. Danshita, and Y. Takahashi, Energy redistribution and spatiotemporal evolution of correlations after a sudden quench of the Bose-Hubbard model, Science Advances 6, 10.1126/sciadv.aba9255 (2020).
  * Schuch  _et al._ [2011] N. Schuch, S. K. Harrison, T. J. Osborne, and J. Eisert, Information propagation for interacting-particle systems, Phys. Rev. A 84, 032309 (2011).
  * Wang and Hazzard [2020] Z. Wang and K. R. Hazzard, Tightening the Lieb-Robinson Bound in Locally Interacting Systems, PRX Quantum 1, 010303 (2020).
  * Läuchli and Kollath [2008] A. M. Läuchli and C. Kollath, Spreading of correlations and entanglement after a quench in the one-dimensional Bose–Hubbard model, Journal of Statistical Mechanics: Theory and Experiment 2008, P05018 (2008).
  * Barmettler  _et al._ [2012] P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath, Propagation front of correlations in an interacting Bose gas, Phys. Rev. A 85, 053625 (2012).
  * Natu and Mueller [2013] S. S. Natu and E. J. Mueller, Dynamics of correlations in shallow optical lattices, Phys. Rev. A 87, 063616 (2013).
  * Carleo  _et al._ [2014] G. Carleo, F. Becca, L. Sanchez-Palencia, S. Sorella, and M. Fabrizio, Light-cone effect and supersonic correlations in one- and two-dimensional bosonic superfluids, Phys. Rev. A 89, 031602 (2014).
  * Lo Gullo and Dell’Anna [2015] N. Lo Gullo and L. Dell’Anna, Spreading of correlations and Loschmidt echo after quantum quenches of a Bose gas in the Aubry-André potential, Phys. Rev. A 92, 063619 (2015).
  * Bernier  _et al._ [2018] J.-S. Bernier, R. Tan, L. Bonnes, C. Guo, D. Poletti, and C. Kollath, Light-Cone and Diffusive Propagation of Correlations in a Many-Body Dissipative System, Phys. Rev. Lett. 120, 020401 (2018).
  * Fitzpatrick and Kennett [2018] M. R. C. Fitzpatrick and M. P. Kennett, Light-cone-like spreading of single-particle correlations in the Bose-Hubbard model after a quantum quench in the strong-coupling regime, Phys. Rev. A 98, 053618 (2018).
  * Mokhtari-Jazi  _et al._ [2021] A. Mokhtari-Jazi, M. R. C. Fitzpatrick, and M. P. Kennett, Phase and group velocities for correlation spreading in the Mott phase of the Bose-Hubbard model in dimensions greater than one, Phys. Rev. A 103, 023334 (2021).
  * Kollath  _et al._ [2007] C. Kollath, A. M. Läuchli, and E. Altman, Quench Dynamics and Nonequilibrium Phase Diagram of the Bose-Hubbard Model, Phys. Rev. Lett. 98, 180601 (2007).
  * Cramer  _et al._ [2008b] M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Exact Relaxation in a Class of Nonequilibrium Quantum Lattice Systems, Phys. Rev. Lett. 100, 030602 (2008b).
  * Cramer  _et al._ [2008c] M. Cramer, A. Flesch, I. P. McCulloch, U. Schollwöck, and J. Eisert, Exploring Local Quantum Many-Body Relaxation by Atoms in Optical Superlattices, Phys. Rev. Lett. 101, 063001 (2008c).
  * Roux [2009] G. Roux, Quenches in quantum many-body systems: One-dimensional Bose-Hubbard model reexamined, Phys. Rev. A 79, 021608 (2009).
  * Navez and Schützhold [2010] P. Navez and R. Schützhold, Emergence of coherence in the Mott-insulator–superfluid quench of the Bose-Hubbard model, Phys. Rev. A 82, 063603 (2010).
  * Enss and Sirker [2012] T. Enss and J. Sirker, Light cone renormalization and quantum quenches in one-dimensional Hubbard models, New Journal of Physics 14, 023008 (2012).
  * Bernier  _et al._ [2012] J.-S. Bernier, D. Poletti, P. Barmettler, G. Roux, and C. Kollath, Slow quench dynamics of Mott-insulating regions in a trapped Bose gas, Phys. Rev. A 85, 033641 (2012).
  * Sorg  _et al._ [2014] S. Sorg, L. Vidmar, L. Pollet, and F. Heidrich-Meisner, Relaxation and thermalization in the one-dimensional Bose-Hubbard model: A case study for the interaction quantum quench from the atomic limit, Phys. Rev. A 90, 033606 (2014).
  * Bernier  _et al._ [2014] J.-S. Bernier, R. Citro, C. Kollath, and E. Orignac, Correlation Dynamics During a Slow Interaction Quench in a One-Dimensional Bose Gas, Phys. Rev. Lett. 112, 065301 (2014).
  * Geiger  _et al._ [2014] R. Geiger, T. Langen, I. E. Mazets, and J. Schmiedmayer, Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional Bose gas, New Journal of Physics 16, 053034 (2014).
  * Krutitsky  _et al._ [2014] K. V. Krutitsky, P. Navez, F. Queisser, and R. Schützhold, Propagation of quantum correlations after a quench in the Mott-insulator regime of the Bose-Hubbard model, EPJ Quantum Technology 1, 12 (2014).
  * Andraschko and Sirker [2015] F. Andraschko and J. Sirker, Propagation of a single-hole defect in the one-dimensional Bose-Hubbard model, Phys. Rev. B 91, 235132 (2015).
  * Shen  _et al._ [2017] H. Shen, P. Zhang, R. Fan, and H. Zhai, Out-of-time-order correlation at a quantum phase transition, Phys. Rev. B 96, 054503 (2017).
  * Liu  _et al._ [2018] F. Liu, J. R. Garrison, D.-L. Deng, Z.-X. Gong, and A. V. Gorshkov, Asymmetric Particle Transport and Light-Cone Dynamics Induced by Anyonic Statistics, Phys. Rev. Lett. 121, 250404 (2018).
  * Pietraszewicz  _et al._ [2019] J. Pietraszewicz, M. Stobińska, and P. Deuar, Correlation evolution in dilute bose-einstein condensates after quantum quenches, Phys. Rev. A 99, 023620 (2019).
  * Despres  _et al._ [2019] J. Despres, L. Villa, and L. Sanchez-Palencia, Twofold correlation spreading in a strongly correlated lattice Bose gas, Scientific Reports 9, 4135 (2019).
  * Villa  _et al._ [2020] L. Villa, J. Despres, S. J. Thomson, and L. Sanchez-Palencia, Local quench spectroscopy of many-body quantum systems, Phys. Rev. A 102, 033337 (2020).
  * [94] See Supplemental Material for the details of the rigorous proof of the main statements, which includes Refs. [97, 98]. [97] [1st reference in Supplemental Material not already in paper], [98] [8th reference in Supplemental Material not already in paper] .
  * Kuwahara [2016] T. Kuwahara, Exponential bound on information spreading induced by quantum many-body dynamics with long-range interactions, New Journal of Physics 18, 053034 (2016).
  * Yin and Lucas [2021] C. Yin and A. Lucas, Finite speed of quantum information in models of interacting bosons at finite density (2021), arXiv:2106.09726 [quant-ph] .
  * Arad  _et al._ [2016] I. Arad, T. Kuwahara, and Z. Landau, Connecting global and local energy distributions in quantum spin models on a lattice, Journal of Statistical Mechanics: Theory and Experiment 2016, 033301 (2016).
  * Cramer and Eisert [2006] M. Cramer and J. Eisert, Correlations, spectral gap and entanglement in harmonic quantum systems on generic lattices, New Journal of Physics 8, 71 (2006).

Supplementary Material for “Lieb-Robinson bound and almost-linear light-cone
in interacting boson systems”
Tomotaka Kuwahara1, Keiji Saito2

1Mathematical Science Team, RIKEN Center for Advanced Intelligence Project
(AIP),1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan
2Department of Physics, Keio University, Yokohama 223-8522, Japan

###### Contents

  1. S.I Set up
    1. S.I.1 Boson operators
    2. S.I.2 Bose-Hubbard type Hamiltonian
    3. S.I.3 Initial condition for the boson density
  2. S.II Main result: Lieb–Robinson bound for interacting bosons
  3. S.III Clustering theorem for gapped ground states
    1. S.III.1 Proof of Corollary 2
  4. S.IV Proof outline of the main theorem (Fig. 4)
    1. S.IV.1 Lieb–Robinson bound for short-time evolution
    2. S.IV.2 Proof of Theorem 1 based on Subtheorem 1
  5. S.V Proof of Subtheorem 1
    1. S.V.1 Boson density after time evolution (Proposition 3)
      1. S.V.1.1 Proof of Corollary 4
    2. S.V.2 Effective Hamiltonian (Proposition 5)
    3. S.V.3 Lieb–Robinson bound for the effective Hamiltonian (Proposition 6)
    4. S.V.4 Completing the proof
  6. S.VI Proof of Proposition 3: boson density after time evolution
    1. S.VI.1 Restatement
    2. S.VI.2 Proof
      1. S.VI.2.1 Proof of Lemma 7
      2. S.VI.2.2 Proof of Lemma 8
      3. S.VI.2.3 Proof of Lemma 9
      4. S.VI.2.4 Proof of Lemma 10
      5. S.VI.2.5 Proof of Lemma 11
  7. S.VII Proof of Proposition 5: effective Hamiltonian
    1. S.VII.1 Restatement
    2. S.VII.2 Proof
      1. S.VII.2.1 Proof of Lemma 12
      2. S.VII.2.2 Proof of Lemma 13
  8. S.VIII Proof of Proposition 6: Lieb–Robinson bound for the effective Hamiltonian
    1. S.VIII.1 Restatement
    2. S.VIII.2 Proof
      1. S.VIII.2.1 Proof of Lemma 14
      2. S.VIII.2.2 Proof of Lemma 15
  9. S.IX Quench
    1. S.IX.1 Setup of the quench dynamics
    2. S.IX.2 Main theorem
    3. S.IX.3 Proof of Theorem 2
    4. S.IX.4 Proof of Proposition 16
      1. S.IX.4.1 Proof of Lemma 17

##  S.I Set up

We first describe the setup, which is described in less detail in the main
text.

We consider a quantum system on a $D$-dimensional lattice (graph), where
bosons interact with each other. An unbounded number of bosons occupy each
site, and hence the local Hilbert dimension is infinitely large. We define
$d_{G}$ as the maximum degree of the lattice (graph). We denote the set of
total sites by $\Lambda$. For an arbitrary partial set $X\subseteq\Lambda$, we
denote the cardinality (i.e., the number of sites contained in $X$) by $|X|$.

For arbitrary subsets $X,Y\subseteq\Lambda$, we define $d_{X,Y}$ as the
shortest path length on the graph that connects $X$ and $Y$; that is, if
$X\cap Y\neq\emptyset$, $d_{X,Y}=0$. When $X$ contains only one element (i.e.,
$X=\\{i\\}$), we abbreviate $d_{\\{i\\},Y}$ as $d_{i,Y}$ for simplicity. We
also denote the complementary subset and surface subset of $X$ by $X^{\rm
c}:=\Lambda\setminus X$ and $\partial X:=\\{i\in X|d_{i,X^{\rm c}}=1\\}$,
respectively. We also define ${\rm diam}(X)$ as follows:

$\displaystyle{\rm diam}(X):=1+\max_{i,j\in X}(d_{i,j}).$ (S.2)

For a subset $X\subseteq\Lambda$, we define the extended subset $X[r]$ as

$\displaystyle X[r]:=\\{i\in\Lambda|d_{X,i}\leq r\\},$ (S.3)

where $X[0]=X$, and $r$ is an arbitrary positive number (i.e.,
$r\in\mathbb{R}^{+}$). From the definition, $i[r]$ is given by a ball region
with radius $r$ centered at the site $i$. We introduce a geometric parameter
$\gamma$, which is determined only by the lattice structure. We define
$\gamma\geq 1$ as a constant of $\mathcal{O}(1)$ that satisfies the following
inequalities:

$\displaystyle|i[r]|\leq\gamma r^{D}\quad(r\geq 1).$ (S.4)

We define the constant $\lambda_{0}$ as follows:

$\displaystyle\lambda_{0}:=\max_{i\in\Lambda}\sum_{j\in\Lambda}e^{-d_{i,j}}.$
(S.5)

The parameter $\lambda_{0}$ depends on $\gamma$ and $D$ because, when the
parameter $\gamma$ is used, it is upper-bounded as

$\displaystyle\sum_{j\in\Lambda}e^{-d_{i,j}}$
$\displaystyle=1+\sum_{x=1}^{\infty}\sum_{j:d_{i,j}=x}e^{-x}$
$\displaystyle\leq 1+\gamma\sum_{x=1}^{\infty}x^{D}e^{-x}\leq
1+\gamma\int_{0}^{\infty}x^{D}e^{-x+1}dx=1+e\gamma D!,$ (S.6)

where we use the inequality $\\#\\{j:d_{i,j}=x\\}\leq\\#\\{j:d_{i,j}\leq
x\\}=|i[x]|\leq\gamma x^{D}$.

###  S.I.1 Boson operators

We define $b_{i}$ and $b_{i}^{\dagger}$ as the boson annihilation and creation
operators, respectively. We also define $\hat{n}_{i}$ as the boson number
operator on site $i$, $\hat{n}_{i}:=b_{i}^{\dagger}b_{i}$. We denote the boson
number on a subset $X$ by $\hat{n}_{X}$ as follows:

$\displaystyle\hat{n}_{X}=\sum_{i\in X}\hat{n}_{i}.$ (S.7)

For an arbitrary subset $X\subseteq\Lambda$, we define $\Pi_{X,q}$ as the
projection onto the eigenspace of $\hat{n}_{X}$ with eigenvalue $q$:

$\displaystyle\hat{n}_{X}\Pi_{X,q}=q\Pi_{X,q}.$ (S.8)

When $X$ includes only one site (i.e., $X=\\{i\\}$), we denote
$\Pi_{\\{i\\},q}$ by $\Pi_{i,q}$ for simplicity. We also define $\Pi_{X,\geq
q}$ as $\sum_{q^{\prime}=q}^{\infty}\Pi_{X,q^{\prime}}$.

###  S.I.2 Bose-Hubbard type Hamiltonian

We consider a Hamiltonian of the form

$\displaystyle H=H_{0}+V,$ $\displaystyle H_{0}:=\sum_{\langle
i,j\rangle}J_{i,j}(b_{i}b_{j}^{\dagger}+{\rm h.c.}),\quad
V:=\sum_{Z\subset\Lambda:|Z|\leq k}v_{Z},$
$\displaystyle|J_{i,j}|\leq\bar{J},$ (S.9)

where $\sum_{\langle i,j\rangle}$ represents the summation of all the pairs of
adjacent sites $\\{i,j\\}$ on the lattice, and $v_{Z}$ represents boson–boson
interactions in subset $Z$. We assume that $v_{Z}$ is now given by a function
of the number operators $\\{\hat{n}_{i}\\}_{i\in Z}$. For example, for
$Z=\\{i,j\\}$, $v_{Z}$ includes interactions such as $\hat{n}_{i}^{4}$,
$\hat{n}_{i}^{2}e^{\hat{n}_{j}}$, $e^{\hat{n}_{i}^{2}\hat{n}_{j}^{3}}$, and so
on. The simplest example is the Bose–Hubbard model:

$\displaystyle H=\sum_{\langle i,j\rangle}J(b_{i}b_{j}^{\dagger}+{\rm
h.c.})+\frac{U}{2}\sum_{i\in\Lambda}\hat{n}_{i}(\hat{n}_{i}-1)-\mu\sum_{i\in\Lambda}\hat{n}_{i},$
(S.10)

where $U$ and $\mu$ are $\mathcal{O}(1)$ constants.

For an arbitrary subset $X\subseteq\Lambda$, we define the subset Hamiltonians
$H_{0,X}$, $V_{X}$, and $H_{X}$ as follows:

$\displaystyle H_{X}=H_{0,X}+V_{X},$ $\displaystyle H_{0,X}:=\sum_{i,j\in
X}[J_{i,j}(b_{i}b_{j}^{\dagger}+{\rm h.c.})],\quad V_{X}:=\sum_{Z\subseteq
X:|Z|\leq k}v_{Z}.$ (S.11)

Note that they are supported on the subset $X$.

For an arbitrary operator $O$, we denote the time evolution by an operator $A$
as $O(A,t)$ as follows:

$\displaystyle O(A,t):=e^{iAt}Oe^{-iAt}.$ (S.12)

In particular, when $A=H$, we often abbreviate $O(H,t)$ as $O(t)$ for
simplicity.

We define $t_{0}$ as the unit of time, which is an $\mathcal{O}(1)$ constant;
for example, we can choose $t_{0}=1$.

###  S.I.3 Initial condition for the boson density

We here define the condition of low boson density as follows.

######  Assumption 1 (Low boson density).

For a quantum state $\sigma$, we say that the state $\sigma$ satisfies the
low-boson-density condition if there exist $\mathcal{O}(1)$ constants $c_{0}$
and $\bar{q}$ such that

$\displaystyle\max_{i\in\Lambda}{\rm
tr}(e^{c_{0}(\hat{n}_{i}-\bar{q})}\sigma)\leq 1\quad(c_{0}\leq 1).$ (S.13)

This condition ensures that the probability that many bosons are concentrated
on one site is exponentially small. Indeed, the probability that $\sigma$ has
more than $q$ ($q\in\mathbb{N}$) bosons on a site $i$ is upper-bounded by

$\displaystyle{\rm tr}(\Pi_{i,\geq q}\rho_{0})$ $\displaystyle={\rm
tr}\left(\Pi_{i,\geq
q}e^{-c_{0}(\hat{n}_{i}-\bar{q})/2}e^{c_{0}(\hat{n}_{i}-\bar{q})/2}\rho_{0}e^{c_{0}(\hat{n}_{i}-\bar{q})/2}e^{-c_{0}(\hat{n}_{i}-\bar{q})/2}\Pi_{i,\geq
q}\right)$ $\displaystyle\leq\|\Pi_{i,\geq
q}e^{-c_{0}(\hat{n}_{i}-\bar{q})/2}\|^{2}\cdot{\rm
tr}(e^{c_{0}(\hat{n}_{i}-\bar{q})}\rho_{0})\leq e^{-c_{0}(q-\bar{q})}.$ (S.14)

In the first inequality, we use the inequality ${\rm
tr}(O^{\dagger}AO)\leq\|O\|^{2}{\rm tr}(A)$ for an arbitrary positive operator
$A$. Therefore, the probability decays exponentially beyond $q\approx\bar{q}$.

Remark. This condition is expected to hold in real experimental setups,
although it would be difficult to prove rigorously in general. Quantum Gibbs
states are among the candidates that satisfy the condition; as a trivial
example, our theory can be applied to the infinite-temperature state. If we
specify a setup (e.g., that there exists a repulsive force between bosons), we
could prove the condition (S.13) for low-energy states by employing the
techniques in [97].

##  S.II Main result: Lieb–Robinson bound for interacting bosons

Let $\rho_{0}$ be a time-independent quantum state, $[\rho_{0},H]=0$. We then
consider the propagation of a perturbation to $\rho_{0}$ as

$\displaystyle\rho=O_{X_{0}}\rho_{0}O_{X_{0}}^{\dagger},\quad X\subseteq
i_{0}[r_{0}].$ (S.15)

We are now interested in how fast this perturbation propagates. After time
evolution, $\rho(t)$ is given mathematically by
$O_{X_{0}}(t)\rho_{0}O_{X_{0}}(t)^{\dagger}$, and hence we need to estimate
the approximation error of

$\displaystyle O_{X_{0}}(t)\rho_{0}\approx O^{(t)}_{i_{0}[R]}\rho_{0},$ (S.16)

where $O^{(t)}_{i_{0}[R]}$ is an appropriate operator supported on the subset
$i_{0}[R]$. We aim to upper-bound the error as a function of $R$ (Fig. 3).

In deriving the Lieb–Robinson bound, we need to assume that boson creation by
$O_{X_{0}}$ is not infinitely large; thus, we adopt the following condition:

$\displaystyle\|\Pi_{X_{0},q}O_{X_{0}}\Pi_{X_{0},q^{\prime}}\|=0\quad\textrm{for}\quad
q^{\prime}>q+q_{0}$ (S.17)

with $q_{0}=\mathcal{O}(1)$, where the projection $\Pi_{X_{0},q}$ has been
defined by Eq. (S.8). The above condition implies that the number of bosons
created by $O_{X_{0}}$ is less than or equal to $q_{0}$. We notice that the
above condition also implies

$\displaystyle\|\Pi_{X,q}O_{X_{0}}\Pi_{X,q^{\prime}}\|=0\quad\textrm{for}\quad
q^{\prime}>q+q_{0}$ (S.18)

for an arbitrary subset $X\supseteq X_{0}$*1*1*1If we consider a subset
$X_{1}\subset X_{0}$, it may break down. For example, for $X_{0}=X_{1}\sqcup
X_{2}$, let us consider an eigenstate $|q_{X_{1}},q_{X_{2}}\rangle$, where
$\hat{n}_{X_{1}}|q_{X_{1}},q_{X_{2}}\rangle=q_{X_{1}}|q_{X_{1}},q_{X_{2}}\rangle$,
$\hat{n}_{X_{2}}|q_{X_{1}},q_{X_{2}}\rangle=q_{X_{2}}|q_{X_{1}},q_{X_{2}}\rangle$,
and $q_{X_{1}}+q_{X_{2}}=q$. If we have
$O_{X_{0}}|q_{X_{1}},q_{X_{2}}\rangle=|q_{X_{1}}+q_{1},q_{X_{2}}-q_{2}\rangle$
with $q_{1}-q_{2}=q_{0}$, the operator $O_{X_{0}}$ satisfies Eq. (S.17).
However, for $q_{2}>0$, the equation
$\displaystyle\|\Pi_{X_{1},q}O_{X_{0}}\Pi_{X_{1},q^{\prime}}\|=0\quad\textrm{for}\quad
q^{\prime}>q+q_{0}$ (S.19) does not hold..

Our main result gives the Lieb–Robinson bound for an arbitrary quantum state
$\rho_{0}$ satisfying the low-boson-density condition (1).

######  Theorem 1.

Let $O_{X_{0}}$ be an arbitrary operator supported on a subset $X_{0}\subseteq
i_{0}[r_{0}]$ ($i_{0}\in\Lambda$). We assume that the number of bosons created
by $O_{X_{0}}$ is finite, as in Eq. (S.17). Then, for an arbitrary steady
quantum state $\rho_{0}$ satisfying the assumption (1), the operator
$O_{X_{0}}(t)\rho_{0}$ is approximated by using another operator
$O^{(t)}_{i_{0}[R]}$ supported on $i_{0}[R]$ with the following approximation
error:

$\displaystyle\left\|\left(O_{X_{0}}(t)-O^{(t)}_{i_{0}[R]}\right)\rho_{0}\right\|_{1}\leq\|O_{X_{0}}\|\exp\left(c_{0}\bar{q}-C_{1}\frac{(R-r_{0})}{t\log(R)}+C_{2}\log(R)\right)\quad(t\geq
1),$ (S.20)

where $C_{1}$ and $C_{2}$ are constants of $\mathcal{O}(1)$ which are
independent of $\bar{q}$ and depend only on the details of the system.

Figure 3: Schematic illustration of Theorem 1.

Remark. From the above result, we can identify the shape of the effective
light cone. We here assume $t\geq e$ and $X_{0}=i_{0}$ ($r_{0}=0$) for
simplicity. As shown in the main text, the velocity $v_{t,\delta}$ to
characterize the effective light cone has been defined by the following
inequality:

$\displaystyle\|[O_{i_{0}}(t)-O^{(t)}_{i_{0}[R]}]\|\leq\delta\left\|O_{i_{0}}\right\|\quad\textrm{for}\quad
R\geq v_{t,\delta}|t|.$ (S.21)

In order to achieve the above inequality for a particular choices of $t$ and
$\delta$, we need to choose $R$ such that

$\displaystyle\exp\left(c_{0}\bar{q}-C_{1}\frac{R}{t\log(R)}+C_{2}\log(R)\right)\leq\delta$
$\displaystyle\longrightarrow C_{1}\frac{R}{t\log(R)}-C_{2}\log(R)\geq
c_{0}\bar{q}+\log(1/\delta)\longrightarrow\frac{R}{t\log(R)}-\tilde{C}_{2}\log(R)\geq
Q,$ (S.22)

where we define $\tilde{C}_{2}:=C_{2}/C_{1}$ and
$Q=[c_{0}\bar{q}+\log(1/\delta)]/C_{1}$. We now choose $R=\Gamma t\log^{2}(t)$
by using a parameter $\Gamma$ ($\geq 1$). Then, from $\Gamma t\leq R\leq\Gamma
t^{2}$ for $t\geq e$, we obtain the condition of

$\displaystyle\Gamma\log(t)-2\tilde{C}_{2}\log(t)-\tilde{C}_{2}\log(\Gamma)\geq
Q.$ (S.23)

Let $\Gamma_{c}$ be a constant depending on $\tilde{C}_{2}$ such that
$(1/2)\Gamma\geq\tilde{C}_{2}\log(\Gamma)$ for $\Gamma\geq\Gamma_{c}$. Then,
for $\Gamma\geq\Gamma_{c}$, the condition (S.25) reduces to

$\displaystyle\Gamma(\log(t)-1/2)\geq
2\tilde{C}_{2}\log(t)+Q\longrightarrow\Gamma\geq\frac{2\tilde{C}_{2}\log(t)+Q}{\log(t)-1/2}\longrightarrow\Gamma\geq
4\tilde{C}_{2}+2Q,$ (S.24)

where the last inequality results from $t\geq e$. Therefore, by choosing

$\displaystyle
R=t\log^{2}(t)\max(\Gamma_{c},4\tilde{C}_{2}+2Q)=t\log^{2}(t)\max\left(\Gamma_{c},4C_{2}/C_{1}+2[c_{0}\bar{q}+\log(1/\delta)]/C_{1}\right),$
(S.25)

the inequality (S.21) is satisfied. We thus conclude that
$v_{t,\delta}\propto\log^{2}(t)[c_{0}\bar{q}+\log(1/\delta)]$.

##  S.III Clustering theorem for gapped ground states

As a direct application of Theorem 1, we can discuss the clustering property
of the bipartite correlations in gapped ground states. Let
$\\{|E_{j}\rangle\\}_{j\geq 0}$ be eigenstates of the Hamiltonian $H$. We set
$E_{0}=0$ and define $\Delta E:=E_{1}$ as the spectral gap between the ground
and first excited energies. We can prove the following corollary.

######  Corollary 2.

Let $O_{X}$ and $O_{Y}$ be arbitrary operators satisfying the condition (S.17)
for a fixed $q_{0}$. Also, we here assume that $\bar{q}=\mathcal{O}(1)$. Then,
for the ground states $|E_{0}\rangle$ with low boson density according to Eq.
(1), the bipartite correlation between $O_{X}$ and $O_{Y}$ satisfies the
inequality

$\displaystyle{\rm Cor}(O_{X},O_{Y}):=\langle
E_{0}|O_{X}O_{Y}|E_{0}\rangle-\langle E_{0}|O_{X}|E_{0}\rangle\langle
E_{0}|O_{Y}|E_{0}\rangle\leq
C_{3}\cdot\|O_{X}\|\cdot\|O_{Y}\|\cdot\exp\left(-\sqrt{\frac{C_{3}^{\prime}\Delta
E}{\log(R)}R}\right)$ (S.26)

for $R\geq{\rm const.}\times(1/\Delta E)\log^{3}(1/\Delta E)$. Here, $C_{3}$
and $C_{3}^{\prime}$ are $\mathcal{O}(1)$ constants.

Remark. According to this corollary, the correlation decays sub-exponentially
with $R$ that is larger than $\tilde{\mathcal{O}}(1/\Delta E)$. This corollary
is weaker than the standard clustering theorem [4, 5, 6], which yields
exponential decay of the bipartite correlations as $e^{-\mathcal{O}(\Delta
ER)}$. The primary reason is that the asymptotic form of the Lieb–Robinson
bound in Theorem 1 is given by $e^{-\mathcal{O}(R/(t\log R))}$ instead of
$e^{-\mathcal{O}(R-vt)}$. This point is reflected in the choice of $T$ in Eqs.
(S.32), (S.34), and (S.35) below. For example, if we could improve the present
upper bound to $e^{-\mathcal{O}(R/(\log R)-vt)}$, we would be able to obtain
nearly exponential decay of the bipartite correlations as
$e^{-\mathcal{O}(\Delta ER/(\log R))}$.

###  S.III.1 Proof of Corollary 2

The proof relies on the method in Refs. [4, 5, 6]. We start with the equation

$\displaystyle\langle E_{0}|[O_{X}(t),O_{Y}]|E_{0}\rangle$
$\displaystyle=\sum_{s\geq 0}\left(e^{-iE_{s}t}\langle
E_{0}|O_{X}|E_{s}\rangle\langle E_{s}|O_{Y}|E_{0}\rangle-e^{iE_{s}t}\langle
E_{0}|O_{Y}|E_{s}\rangle\langle E_{s}|O_{X}|E_{0}\rangle\right)$
$\displaystyle=\sum_{s\geq 1}\left(e^{-iE_{s}t}\langle
E_{0}|O_{X}|E_{s}\rangle\langle E_{s}|O_{Y}|E_{0}\rangle-e^{iE_{s}t}\langle
E_{0}|O_{Y}|E_{s}\rangle\langle E_{s}|O_{X}|E_{0}\rangle\right).$ (S.27)

Using the function $K(t)$, where

$\displaystyle K(t)=\frac{i}{2\pi}\lim_{\epsilon\to+0}\frac{e^{-\alpha
t^{2}}}{t+i\epsilon},$ (S.28)

we have, from Ref. [5],

$\displaystyle\int_{-\infty}^{\infty}e^{-i\omega t}K(t)dt=\begin{cases}{\rm
const}\cdot e^{-\omega^{2}/(4\alpha)}&\quad\textrm{for}\quad\omega>0,\\\
1+{\rm const}\cdot
e^{-\omega^{2}/(4\alpha)}&\quad\textrm{for}\quad\omega<0.\end{cases}$

Here, $\alpha$ is a free parameter. Because $K(t)$ decays as $e^{-\alpha
t^{2}}$, we can obtain

$\displaystyle\int_{-T}^{T}e^{-i\omega t}K(t)dt$
$\displaystyle=\begin{cases}{\rm
const}\cdot\left(e^{-\omega^{2}/(4\alpha)}+c_{\alpha,T}e^{-\alpha
T^{2}}\right)&\quad\textrm{for}\quad\omega>0,\\\ 1+{\rm
const}\cdot\left(e^{-\omega^{2}/(4\alpha)}+c_{\alpha,T}e^{-\alpha
T^{2}}\right)&\quad\textrm{for}\quad\omega<0,\end{cases}$
$\displaystyle=\begin{cases}{\rm const}\cdot(e^{-T\omega^{2}/(2\Delta
E)}+c_{\Delta E,T}e^{-T\Delta E/2})&\quad\textrm{for}\quad\omega>0,\\\ 1+{\rm
const}\cdot(e^{-T\omega^{2}/(2\Delta E)}+c_{\Delta E,T}e^{-T\Delta
E/2})&\quad\textrm{for}\quad\omega<0,\end{cases}$ (S.29)

where $c_{\alpha,T}$ and $c_{\Delta E,T}$ are appropriate constants, and we
choose $\alpha=\Delta E/(2T)$. We note that $|c_{\Delta E,T}|\leq{\rm
const.}/\sqrt{\alpha}={\rm const.}\sqrt{2T/\Delta E}$.

By applying Eq. (S.29) to Eq. (S.27), we obtain

$\displaystyle\int_{-T}^{T}K(t)\langle E_{0}|[O_{X}(t),O_{Y}]|E_{0}\rangle dt$
$\displaystyle=\int_{-T}^{T}K(t)\sum_{s\geq 1}\left(e^{-iE_{s}t}\langle
E_{0}|O_{X}|E_{s}\rangle\langle E_{s}|O_{Y}|E_{0}\rangle-e^{iE_{s}t}\langle
E_{0}|O_{Y}|E_{s}\rangle\langle E_{s}|O_{X}|E_{0}\rangle\right)dt$
$\displaystyle=\sum_{s\geq 1}\left[K_{E_{s},\Delta E}\langle
E_{0}|O_{X}|E_{s}\rangle\langle E_{s}|O_{Y}|E_{0}\rangle+(1+K_{E_{s},\Delta
E})\langle E_{0}|O_{Y}|E_{s}\rangle\langle E_{s}|O_{X}|E_{0}\rangle\right]$
$\displaystyle=(\langle E_{0}|O_{X}Q_{T}O_{Y}|E_{0}\rangle+{\rm
c.c.})+\sum_{s\geq 1}\langle E_{0}|O_{Y}|E_{s}\rangle\langle
E_{s}|O_{X}|E_{0}\rangle,$ (S.30)

where we set $K_{\omega,\Delta E}:={\rm const}\cdot(e^{-T\omega^{2}/(2\Delta
E)}+c_{\Delta E,T}e^{-T\Delta E/2})$ and define $Q_{T}$ as

$\displaystyle Q_{T}:=\sum_{s\geq 1}K_{E_{s},\Delta E}|E_{s}\rangle\langle
E_{s}|.$ (S.31)

By combining Eqs. (S.29) and (S.30), we obtain

$\displaystyle\left|{\rm Cor}(O_{X},O_{Y})\right|$
$\displaystyle\leq\int_{-T}^{T}K(t)|\langle
E_{0}|[O_{X}(t),O_{Y}]|E_{0}\rangle|dt+{\rm const}\cdot c_{\Delta
E,T}\|O_{X}\|\cdot\|O_{Y}\|\cdot e^{-T\Delta E/2},$ (S.32)

where we use $\sum_{s\geq 1}\langle E_{0}|O_{Y}|E_{s}\rangle\langle
E_{s}|O_{X}|E_{0}\rangle={\rm Cor}(O_{X},O_{Y})$ and $\|Q_{T}\|\leq
K_{E_{1},\Delta E}={\rm const}\cdot c_{\Delta E,T}e^{-T\Delta E/2}$ from
$\Delta E=E_{1}$. From Theorem 1, we have

$\displaystyle|\langle
E_{0}|[O_{X}(t),O_{Y}]|E_{0}\rangle|\leq\|[O_{X}(t),O_{Y}]\|\leq\|O_{X}\|\cdot\|O_{Y}\|\cdot\exp\left(-{\rm
const}\cdot\frac{R}{t\log(R)}\right)$ (S.33)

for $R\geq{\rm const}\cdot t\log^{2}(t)$. Hence, for $R\geq{\rm const}\cdot
T\log^{2}(T)$, we obtain

$\displaystyle\int_{-T}^{T}K(t)|\langle
E_{0}|[O_{X}(t),O_{Y}]|E_{0}\rangle|dt\leq{\rm
const}\cdot\|O_{X}\|\cdot\|O_{Y}\|\cdot T\exp\left(-{\rm
const}\cdot\frac{R}{T\log(R)}\right).$ (S.34)

We thus choose $T$ as

$\displaystyle T=\sqrt{\frac{{\rm const}\cdot R}{\Delta E\log(R)}},$ (S.35)

and we obtain the main inequality (S.26). Finally, the condition $R\geq{\rm
const}\cdot T\log^{2}(T)$ reduces to

$\displaystyle R\log R\geq{\rm const}\cdot\frac{1}{\Delta
E}[\log^{4}(R)+\log^{4}(1/\Delta E)],$ (S.36)

which is satisfied if $R\geq{\rm const}\cdot(1/\Delta E)\log^{3}(1/\Delta E)$.
This completes the proof. $\square$

##  S.IV Proof outline of the main theorem (Fig. 4)

###  S.IV.1 Lieb–Robinson bound for short-time evolution

The key ingredient in our proof is the Lieb–Robinson bound for short-time
evolution. We consider a quantum state $\tilde{\rho}$, which is defined using
an operator $O_{X}$ supported on a subset $X\subset\Lambda$:

$\displaystyle\tilde{\rho}:=O_{X}\rho_{0}O_{X}^{\dagger},\quad X\in
i_{0}[r]\quad(i_{0}\in\Lambda,\ r\geq 3),$ (S.37)

where $O_{X}$ is given in the form

$\displaystyle
O_{X}=U_{X}^{\dagger}O_{X_{0}}U_{X},\quad[U_{X},\hat{n}_{X}]=0,$ (S.38)

where $U_{X}$ is a unitary operator that commutes with $\hat{n}_{X}$. Because
$O_{X_{0}}$ satisfies the condition (S.17), and $U_{X}$ does not change the
total number of bosons on $X$, we obtain

$\displaystyle\|\Pi_{X,q}O_{X}\Pi_{X,q^{\prime}}\|=0\quad\textrm{for}\quad
q^{\prime}>q+q_{0}.$ (S.39)

We are now interested in the approximation

$\displaystyle
O_{X}(t)\rho_{0}\approx(U_{X[\ell]}^{\dagger}O_{X}U_{X[\ell]})\rho_{0}$ (S.40)

for a sufficiently small $t=\mathcal{O}(1)$, where $U_{X[\ell]}$ is an
appropriate unitary operator defined on the subset $X[\ell]$. We can prove the
following statement about the approximation.

######  Subtheorem 1.

Let $O_{X}$ be an arbitrary operator as defined in Eq. (S.38), which is
supported on the subset $X\subseteq i_{0}[r]$ ($i_{0}\in\Lambda$). Then, for a
length $\ell$ that satisfies

$\displaystyle\ell\geq C_{0}\log^{2}(r),$ (S.41)

with $C_{0}=\mathcal{O}(1)$ which does not depend on $\bar{q}$, we can find a
unitary operator $U_{X[\ell]}$ that does not depend on the form of $O_{X}$
such that

$\displaystyle[U_{X[\ell]},\hat{n}_{X[\ell]}]=0$ (S.42)

and

$\displaystyle\left\|\left(O_{X}(t)-U_{X[\ell]}^{\dagger}O_{X}U_{X[\ell]}\right)\rho_{0}\right\|_{1}\leq\|O_{X}\|e^{c_{0}\bar{q}-\ell/\log(r)}$
(S.43)

for $t\leq\Delta t_{0}$, where $\Delta t_{0}=\mathcal{O}(1)$ and $C_{0}$ are
appropriately chosen. We can combine condition (S.41) with Ineq. (S.43) as
follows:

$\displaystyle\left\|\left(O_{X}(t)-U_{X[\ell]}^{\dagger}O_{X}U_{X[\ell]}\right)\rho_{0}\right\|_{1}\leq\|O_{X}\|e^{c_{0}\bar{q}-\ell/\log(r)+C_{0}\log(r)},$
(S.44)

where the inequality holds trivially for $\ell\leq C_{0}\log^{2}(r)$.

Remark. In this subtheorem, we do not need to assume the time-independence of
the state $\rho_{0}$. Hence, only the assumption 1 is used. When $\rho_{0}$
does not satisfy $[\rho_{0},H]=0$, Ineq. (S.43) is replaced by

$\displaystyle\left\|\left(O_{X}(t)-U_{X[\ell]}^{\dagger}O_{X}U_{X[\ell]}\right)\rho_{0}(-t)\right\|_{1}\leq\|O_{X}\|e^{c_{0}\bar{q}-\ell/\log(r)}.$
(S.45)

By contrast, when we prove the main theorem, we need to connect the short-time
evolutions, and to perform this procedure, the time-independence of $\rho_{0}$
is required [see Eqs. (S.IV.2) and (S.53)].

It is necessary to generalize Eq. (S.45) when applying the subtheorem to
analyze the quench dynamics (see Sec. S.IX).

###  S.IV.2 Proof of Theorem 1 based on Subtheorem 1

For the convenience of readers, we present the approach to the proof in the
main text again with additional explanations. We use the connection of unitary
time evolutions addressed in Refs. [95, 29], which assumes the time-
independence of the initial state (i.e., $\rho_{0}(t)=\rho_{0}$).

We introduce $\Delta t\leq\Delta t_{0}$ and the following decompositions of
the time $t$ to $t/\Delta t$ pieces:

$\displaystyle t:=m_{t}\Delta t\quad{\rm with}\quad\Delta t\geq\min(t,\Delta
t_{0}/2).$ (S.46)

For fixed $R$, we define the subset $X_{m}$ as follows:

$\displaystyle X_{m}:=i_{0}[r_{0}+m\Delta r]=X_{0}[m\Delta r],\quad\Delta
r=\left\lfloor\frac{R-r_{0}}{m_{t}}\right\rfloor.$ (S.47)

Note that $X_{m_{t}}\subseteq i_{0}[R]$.

Then, assuming that $\rho_{0}$ is time-invariant, we can derive the following
inequality [29]:

$\displaystyle\left\|\left[O_{X_{0}}(m_{t}\Delta
t)-O_{X_{m_{t}}}^{(m_{t})}\right]\rho_{0}\right\|_{1}$
$\displaystyle\leq\sum_{m=1}^{m_{t}}\left\|\left[O_{X_{m-1}}^{(m-1)}(\Delta
t)-O_{X_{m}}^{(m)}\right]\rho_{0}\right\|_{1},$ (S.48)

where $O_{X_{0}}^{(0)}=O_{X_{0}}$, and $O_{X_{m}}^{(m)}$ is recursively
defined by approximating $O_{X_{m-1}}^{(m-1)}(\Delta t)$. To see why the time-
invariance of $\rho_{0}$ is essential, let us look at the derivation of Ineq.
(S.48) for $m_{t}=2$. For $m=1$, we define

$\displaystyle
O_{X_{1}}^{(1)}:=U_{X_{1}}^{(1)\dagger}O_{X_{0}}U_{X_{1}}^{(1)},$ (S.49)

where we choose the unitary operator $U_{X_{1}}$ by following Subtheorem 1.
Note that $O_{X_{1}}^{(1)}$ is now supported on the subset $X_{1}$. For $m=2$,
we consider the approximation $O_{X_{1}}^{(1)}(\Delta t)$ by

$\displaystyle
O_{X_{2}}^{(2)}:=U_{X_{2}}^{(2)\dagger}O_{X_{1}}^{(1)}U_{X_{2}}^{(2)},$ (S.50)

where the unitary operator $U_{X_{2}}^{(2)}$ is chosen according to Subtheorem
1. We then connect the two approximations:

$\displaystyle O_{X_{0}}(\Delta t)\xrightarrow{\rm
approximate}O_{X_{1}}^{(1)},\quad O_{X_{1}}^{(1)}(\Delta t)\xrightarrow{\rm
approximate}O_{X_{2}}^{(2)}.$ (S.51)

To obtain the approximation error, we need to consider

$\displaystyle\left\|\left[O_{X_{0}}(2\Delta
t)-O_{X_{2}}^{(2)}\right]\rho_{0}\right\|_{1}$
$\displaystyle\leq\left\|\left[O_{X_{0}}(2\Delta t)-O_{X_{1}}^{(1)}(\Delta
t)+O_{X_{1}}^{(1)}(\Delta t)-O_{X_{2}}^{(2)}\right]\rho_{0}\right\|_{1}$
$\displaystyle\leq\left\|\left[O_{X_{0}}(2\Delta t)-O_{X_{1}}^{(1)}(\Delta
t)\right]\rho_{0}\right\|_{1}+\left\|\left[O_{X_{1}}^{(1)}(\Delta
t)-O_{X_{2}}^{(2)}\right]\rho_{0}\right\|_{1}.$ (S.52)

The second term is upper-bounded according to Subtheorem 1. The first term is
given by

$\displaystyle\left\|\left[O_{X_{0}}(2\Delta t)-O_{X_{1}}^{(1)}(\Delta
t)\right]\rho_{0}\right\|_{1}=\left\|\left[O_{X_{0}}(\Delta
t)-O_{X_{1}}^{(1)}\right]\rho_{0}(-\Delta t)\right\|_{1}.$ (S.53)

If $\rho_{0}(-\Delta t)=\rho_{0}$, we can upper-bound the above quantity using
Subtheorem 1. However, when $\rho_{0}$ is time-dependent, the condition (1)
for low boson density may not be satisfied for $\rho_{0}(-\Delta t)$.
Therefore, to prove the main theorem for generic $\rho_{0}$, we need to prove
the low-boson-density condition for $\rho_{0}(-m\Delta t)$ ($m\leq m_{t}-1$).

We return to Ineq. (S.48). According to Subtheorem 1, we can find
$O_{X_{m}}^{(m)}$ such that

$\displaystyle
O_{X_{m}}^{(m)}:=U_{X_{m}}^{(m)\dagger}O_{X_{m-1}}^{(m-1)}U_{X_{m}}^{(m)}$
(S.54)

for each $m=1,2,\ldots,m_{t}$. From Ineq. (S.44), the unitary operator
$U_{X_{m}}^{(m)}$ satisfies

$\displaystyle\left\|\left(U_{X_{m}}^{(m)\dagger}O_{X_{m-1}}^{(m-1)}U_{X_{m}}^{(m)}-O_{X_{m-1}}^{(m-1)}(\Delta
t)\right)\rho_{0}\right\|_{1}\leq\|O_{X_{0}}\|e^{c_{0}\bar{q}-\Delta
r/\log(R)+C_{0}\log(R)},$ (S.55)

where we use $X_{m}\subseteq i_{0}[R]$ for all $m$ in applying Subtheorem 1.
We thus obtain

$\displaystyle\left\|\left[O_{X_{0}}(m_{t}\Delta
t)-O_{X_{m_{t}}}^{(m_{t})}\right]\rho_{0}\right\|_{1}$ $\displaystyle\leq
m_{t}\|O_{X_{0}}\|e^{c_{0}\bar{q}-\Delta r/\log(R)+C_{0}\log(R)}$
$\displaystyle\leq\|O_{X_{0}}\|\exp\left(c_{0}\bar{q}-\frac{\Delta
t(R-r_{0})}{t\log(R)}+\frac{1}{\log(R)}+(C_{0}+1)\log(R)\right),$ (S.56)

for $R\geq 2$, where in the second inequality we use $\Delta r\geq 1$ [or
$m_{t}\leq(R-r_{0})$]*2*2*2Otherwise, the upper bound is worse than the
trivial inequality, i.e., $\|[O_{X_{0}}(m_{t}\Delta
t)-O_{X_{m_{t}}}^{(m_{t})}]\rho_{0}\|_{1}\leq\|O_{X_{0}}\|$.. Thus, because of
$1/\log(R)\leq 3\log(R)$ for $R\geq 2$, by choosing $C_{1}=\Delta t$ and
$C_{2}=C_{0}+4$, the inequality (S.IV.2) reduces to Ineq. (S.20). This
completes the proof of Theorem 1. $\square$

In the following sections, we give the full proof of Subtheorem 1. Note that
the following proof repeats some of the explanations in the main text. We show
the outline of the proof in Fig. 4.

Figure 4: Schematic illustration of the proof of Theorem 1.

##  S.V Proof of Subtheorem 1

Here, we derive the inequality (S.45) in which the time-independence of
$\rho_{0}$ is not satisfied. Throughout the proof, we denote $\|O_{X_{0}}\|$
by $\zeta_{0}$. Because of the definition (S.38) of $O_{X}$, we have

$\displaystyle\|O_{X}\|=\zeta_{0}.$ (S.57)

We first show that the following simple analysis does not work in proving the
subtheorem. If short-time evolution is considered, the following simple Taylor
expansion is expected to work:

$\displaystyle\|[O_{X}(t),u_{i}]\rho_{0}\|_{1}=\|[O_{X},u_{i}(-t)]\rho_{0}\|_{1}=\left\|\sum_{s=0}^{\infty}\frac{(-t)^{s}}{s!}[{\rm
ad}_{H}^{s}(u_{i}),O_{X}]\rho_{0}\right\|_{1}\leq\sum_{s=0}^{\bar{s}}\frac{t^{s}}{s!}\left\|[{\rm
ad}_{H}^{s}(u_{i}),O_{X}]\rho_{0}\right\|_{1},$ (S.58)

where $u_{i}$ is an arbitrary unitary operator acting on a site $i\in\Lambda$.
Because the Hamiltonian is short-range, we have $[{\rm
ad}_{H}^{s}(u_{i}),O_{X}]=0$ for $s\leq\bar{s}$, where
$\bar{s}=\mathcal{O}(d_{i,X})$. Hence, we have

$\displaystyle\|[O_{X}(t),u_{i}]\rho_{0}\|_{1}\leq\sum_{s=\bar{s}+1}^{\infty}\frac{t^{s}}{s!}\left\|[{\rm
ad}_{H}^{s}(u_{i}),O_{X}]\rho_{0}\right\|_{1}.$ (S.59)

In standard spin models with bounded local energy, we have

$\displaystyle\|{\rm ad}_{H}^{s}(u_{i})\|\lesssim s!,$ (S.60)

and hence we can ensure the exponential convergence of the expansion (S.58)
for $t=\mathcal{O}(1)$. Unfortunately, this simple estimation cannot be used
to obtain Subtheorem 1. When we formally describe the Hamiltonian (S.9) as

$\displaystyle H=\sum_{|Z|\leq k}h_{Z},$ (S.61)

we need to consider the norm of

$\displaystyle h_{Z_{1}}h_{Z_{2}}\cdots
h_{Z_{s_{1}}}u_{i}h_{Z_{s_{1}+1}}\cdots h_{Z_{s}}O_{X}\rho_{0},$ (S.62)

where $h_{Z}$ consists of the boson hoppings $b_{i}b_{j}^{\dagger}$ and the
boson–boson interactions $v_{Z}$ in Eq. (S.9). Because the boson number
$\hat{n}_{i}$ on an arbitrary site is strongly suppressed in $\rho_{0}$ from
Assumption 1, we have $\|h_{Z}\rho_{0}\|_{1}=\mathcal{O}(1)$. By contrast, in
the state $O_{X}\rho_{0}O_{X}^{\dagger}$, all the bosons in the region $X$ can
be concentrated on one site, which may give
$\|\hat{n}_{i}O_{X}\rho_{0}\|_{1}\propto|X|$ for $i\in X$. We thus obtain
$\|h_{Z}O_{X}\rho_{0}\|_{1}\gtrsim|X|^{\nu}$ for $Z\cap X\neq\emptyset$, where
$\nu$ is a positive integer depending on the form of $v_{Z}$. Consequently,
the norm $\left\|[{\rm ad}_{H}^{s}(u_{i}),O_{X}]\rho_{0}\right\|_{1}$ has a
rather weaker upper bound:

$\displaystyle\left\|[{\rm
ad}_{H}^{s}(u_{i}),O_{X}]\rho_{0}\right\|_{1}\lesssim s!|X|^{\nu(s-\bar{s})}.$
(S.63)

By combining the above inequality with Eq. (S.59), we can ensure the
convergence of the expansion only for $t=\mathcal{O}(1/|X|^{\nu})$, which is
too weak to prove Theorem 1.

In the following, we take a different route to prove Subtheorem 1 by the three
steps in Secs. S.V.1, S.V.2, and S.V.3.

###  S.V.1 Boson density after time evolution (Proposition 3)

We first consider the boson number distribution after a short-time evolution.
To this end, we need to estimate

$\displaystyle{\rm tr}\left[\Pi_{i,\geq
z}\tilde{\rho}(t)\right],\quad\tilde{\rho}(t)=e^{-iHt}O_{X}\rho_{0}O_{X}^{\dagger}e^{iHt}=O_{X}(-t)\rho_{0}(t)O_{X}(-t)^{\dagger},$
(S.64)

where $\Pi_{i,z}$ has been defined by Eq. (S.8). In the state $\tilde{\rho}$,
the boson number $\hat{n}_{i}$ for $i\in X$ can be as large as
$\mathcal{O}(|X|)$, whereas the boson number $\hat{n}_{i}$ ($i\in X^{\rm c}$)
is exponentially suppressed, as shown in Eq. (S.I.3). During the time
evolution, the bosons concentrated on $X$ spread outside of $X$. We expect
that after a sufficiently small time, the exponential decay of the boson
number still holds for sites that are sufficiently separated from $X$.

The first proposition ensures that the boson density is not seriously affected
by the time evolution for $d_{i,X}\gg 1$ if the time $t$ is of
$\mathcal{O}(1)$ (see Sec. S.VI for the proof).

######  Proposition 3.

We first define the operator $M^{(s)}_{i}(t)$ as

$\displaystyle M^{(s)}_{i}(t):={\rm
tr}\left[\hat{n}_{i}^{s}\tilde{\rho}(t)\right].$ (S.65)

Then, for $t\leq t_{0}$, the following upper bound for $M^{(s)}_{i}(t)$ holds:

$\displaystyle M_{i}^{(s)}(t)$ $\displaystyle\leq
c^{\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}|X|^{3}(c_{1}s|X|)^{s}e^{-d_{i,X}}+c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}(c_{1}s)^{s},$
(S.66)

where $c_{1}$, $c^{\prime}_{1}$, and $c^{\prime\prime}_{1}$ are defined as

$\displaystyle c_{1}:=e^{8\bar{J}d_{G}t_{0}}/c_{0},$ $\displaystyle
c^{\prime}_{1}:=320c_{0}^{-3}e^{4\bar{J}d_{G}t_{0}+c_{0}(1+q_{0}/|X|)},$
$\displaystyle
c^{\prime\prime}_{1}:=80\lambda_{0}c_{0}^{-1}e^{4\bar{J}d_{G}t_{0}+c_{0}}.$
(S.67)

Notice that they are $\mathcal{O}(1)$ constants if $t_{0}=\mathcal{O}(1)$.

Remark. In the proof, we fully use the methods in Ref. [67], which treats the
case of $s=1$. The above upper bound increases exponentially with $t$, and
hence we cannot use it to upper-bound the boson density for general $t$. The
key point of the proof in Sec. S.IV.2 is that we need to treat only the short-
time evolution in this subtheorem. We afterward connect the short-time
evolutions step by step, as in Ineq. (S.48).

In this proposition, the form of the Hamiltonian, i.e., Eq. (S.9), is
essential; if the Hamiltonian includes interactions such as
$\hat{n}_{i}\hat{n}_{j}b_{i}b_{j}^{\dagger}$, the above proposition breaks
down even for small $t$.

By using Proposition 3, we can immediately derive the following corollary of
the boson number distribution.

######  Corollary 4.

Let us define the boson number distribution on a site $i\in\Lambda$ as
follows:

$\displaystyle P_{i,\geq z_{0}}^{(t)}:=\sum_{z\geq z_{0}}{\rm
tr}\left[\Pi_{i,z}\tilde{\rho}(t)\right],$ (S.68)

where $\Pi_{i,z}$ has been defined by Eq. (S.8). Then, for arbitrary
$i\in\Lambda$ such that

$\displaystyle\frac{\gamma^{3}c^{\prime}_{1}}{c^{\prime\prime}_{1}}r^{3D}\leq
e^{d_{i,X}/2}\to d_{i,X}\geq
2\log\left(\frac{\gamma^{3}c^{\prime}_{1}}{c^{\prime\prime}_{1}}\right)+6D\log(r)=\mathcal{O}(\log(r)),$
(S.69)

we obtain

$\displaystyle P_{i,\geq z_{0}}^{(t)}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}\left(\frac{\tilde{c}_{1}d_{i,X}}{z_{0}}\right)^{\tilde{c}^{\prime}_{1}d_{i,X}/\log(r)},$
(S.70)

where we have defined $r$ ($\geq 3$) by $X\subseteq i_{0}[r]$ for an
appropriate choice of $i_{0}\in\Lambda$ in Eq. (S.37), and $\tilde{c}_{1}$ and
$\tilde{c}^{\prime}_{1}$ are constants of $\mathcal{O}(1)$.

####  S.V.1.1 Proof of Corollary 4

Because $X\subseteq i_{0}[r]$, we have from Ineq. (S.4)

$\displaystyle|X|\leq\gamma r^{D}.$ (S.71)

We then choose $s$ such that

$\displaystyle
c^{\prime}_{1}|X|^{3}(c_{1}s|X|)^{s}e^{-d_{i,X}}\leq\gamma^{3}c^{\prime}_{1}r^{3D}(\gamma
c_{1}sr^{D})^{s}e^{-d_{i,X}}\leq c^{\prime\prime}_{1}(c_{1}s)^{s},$ (S.72)

which yields

$\displaystyle\frac{\gamma^{3}c^{\prime}_{1}}{c^{\prime\prime}_{1}}r^{3D}(\gamma
r^{D})^{s}\leq e^{d_{i,X}/2}(\gamma r^{D})^{s}\leq e^{d_{i,X}}$
$\displaystyle\longrightarrow s\leq\frac{d_{i,X}}{2\log(\gamma r^{D})},$
(S.73)

where we use the condition (S.69). We choose $s$ as

$\displaystyle s=\left\lfloor\frac{d_{i,X}}{2\log(\gamma
r^{D})}\right\rfloor\in\left[\tilde{c}^{\prime}_{1}\frac{d_{i,X}}{\log(r)},\tilde{c}_{1}^{\prime\prime}\frac{d_{i,X}}{\log(r)}\right],$
(S.74)

where $\tilde{c}^{\prime}_{1}$ and $\tilde{c}^{\prime\prime}_{1}$ are
constants which depend only on $\gamma$ and $D$. By using the inequality
(S.72), we reduce the Ineq. (S.66) to

$\displaystyle M_{i}^{(s)}(t)$ $\displaystyle\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}(c_{1}s)^{s},$ (S.75)

and hence Markov’s inequality yields

$\displaystyle P_{i,\geq z_{0}}^{(t)}\leq\frac{M_{i}^{(s)}(t)}{z_{0}^{s}}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}\left(\frac{c_{1}\tilde{c}^{\prime\prime}_{1}d_{i,X}}{z_{0}\log(r)}\right)^{\tilde{c}^{\prime}_{1}d_{i,X}/\log(r)}.$
(S.76)

By using $\log(r)\geq 1$ from $r\geq 3$ and defining
$\tilde{c}_{1}:=c_{1}\tilde{c}^{\prime\prime}_{1}$, we prove the main
inequality (S.70). This completes the proof. $\square$

###  S.V.2 Effective Hamiltonian (Proposition 5)

(a) Subsets $L_{1}$ and $L_{2}$ (b) Subset $\tilde{L}$

Figure 5:  Schematic illustrations of sets $L_{1}$, $L_{2}$, and $\tilde{L}$.
In the region $L_{1}$, we cannot ensure that the boson number distribution is
exponentially localized. By contrast, if $\ell_{0}$ is sufficiently large [see
the condition (S.83)], it is ensured that in the region
$\tilde{L}=L_{2}\setminus L_{1}$, the boson number can be truncated up to a
finite value. We define the effective Hamiltonian
$\tilde{H}[\tilde{L},\eta\ell_{0}]$ using Eqs. (S.79) and (S.V.2), which well
approximates the exact dynamics, as described by Ineq. (S.84).

Proposition 3 and Corollary 4 imply that during the time evolution, the boson
number $\hat{n}_{i}$ is strongly suppressed as long as site $i$ is separated
from the region $X$ by a sufficient distance. Hence, we expect that in the
Hamiltonian $H$, the maximum boson number on one site can be truncated. In
this subsection, we rigorously justify this procedure.

As shown in Fig. 5, we first define $L_{1}$ and $L_{2}$ as

$\displaystyle L_{1}=X[\ell_{0}],\quad L_{2}=X[2\ell_{0}].$ (S.77)

We also define $\tilde{L}$ as

$\displaystyle\tilde{L}:=L_{2}\setminus L_{1}.$ (S.78)

We now define $\bar{\Pi}_{L,q}$ ($L\subseteq\Lambda$) as the projection onto
the eigenspace such that for an arbitrary $i\in L$, the boson number
$\hat{n}_{i}$ is truncated up to $q$:

$\displaystyle\bar{\Pi}_{L,q}:=\prod_{i\in L}\Pi_{i,\leq q}.$ (S.79)

Note that $\|\hat{n}_{j}\bar{\Pi}_{L,q}\|\leq q$ for $\forall j\in L$. During
the time evolution of $\tilde{\rho}$ in Eq. (S.37), the boson number is
exponentially suppressed as long as site $i$ is sufficiently separated from
$X$ (see Lemma 12 below). We then aim to approximate the time evolution
$e^{-iHt}$ by using another Hamiltonian $\tilde{H}[L,q]$ as

$\displaystyle\tilde{H}[L,q]:=\tilde{H}_{0}[L,q]+\tilde{V}[L,q],$
$\displaystyle\tilde{H}_{0}[L,q]:=\bar{\Pi}_{L,q}H_{0}\bar{\Pi}_{L,q},\quad
V[L,q]:=\bar{\Pi}_{L,q}V\bar{\Pi}_{L,q},$ (S.80)

where the subset $L$ can be arbitrarily chosen in the definitions. In this
effective Hamiltonian, the boson number is truncated up to $q$ in the region
$L$.

In the following, we choose $L=\tilde{L}$ in Eq. (S.V.2) and consider
$\tilde{H}[\tilde{L},\eta\ell_{0}]$ as the effective Hamiltonian, where
$\tilde{L}$ is given by Eq. (S.78), and $\eta$ is chosen appropriately (see
Proposition 5 below). We usually cannot say that the time evolution by $H$ is
approximated by $e^{-i\tilde{H}[\tilde{L},\eta\ell_{0}]t}$, that is,

$\displaystyle\left\|O_{X}(t)-O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\right\|\approx\|O_{X}\|.$
(S.81)

However, when the above operator acts on the state $\rho_{0}(-t)$, the
difference can be small. Hence, we need to estimate the difference between

$\displaystyle O_{X}(t)\rho_{0}(-t)\quad{\rm and}\quad
O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\rho_{0}(-t)\quad(t\leq t_{0}).$
(S.82)

We can prove the following proposition (see Sec. S.VII for the proof).

######  Proposition 5.

Let us choose $\ell_{0}$ so that it satisfies

$\displaystyle\ell_{0}\geq c_{2}\log^{2}(r).$ (S.83)

Recall that $r$ was defined in Eq. (S.37). Then, there exists a constant
$\eta$ that gives

$\displaystyle\left\|\left[O_{X}(t)-O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\right]\rho_{0}(-t)\right\|_{1}\leq\frac{1}{2}\zeta_{0}e^{c_{0}\bar{q}}e^{-2\ell_{0}/\log(r)}$
(S.84)

for $t\leq t_{0}$. Here, $c_{2}$ and $\eta$ are $\mathcal{O}(1)$ constants
which do not depend on $\bar{q}$.

From this proposition, we can see that the error decreases exponentially with
the number of the boson truncations. In the Hamiltonian
$\tilde{H}[\tilde{L},\eta\ell_{0}]$, the greatest obstacle, namely, the
unboundedness of the interaction norms, has been removed, at least in the
region $\tilde{L}$. However, outside the region $\tilde{L}$, the norm is still
unbounded. In the following subsection, we consider how to derive the
Lieb–Robinson bound for $e^{-i\tilde{H}[\tilde{L},\eta\ell_{0}]t}$ only from
the finiteness of the norm in the region $\tilde{L}$.

###  S.V.3 Lieb–Robinson bound for the effective Hamiltonian (Proposition 6)

We are now interested in the Lieb–Robinson bound for the effective Hamiltonian
$\tilde{H}[\tilde{L},\eta\ell_{0}]$ defined in Proposition 5. In this section,
we adopt an additional condition for $\ell_{0}$, as follows:

$\displaystyle\ell_{0}\geq 8k,$ (S.85)

where $k$ represents the maximum interaction length in $V$ [see Eq. (S.9)].

We would like to calculate the norm

$\displaystyle\left\|O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right\|\quad(t\leq
t_{0}),$ (S.86)

where $U_{L_{2}}$ is an appropriate unitary operator supported on the subset
$L_{2}$. To this end, we first decompose
$e^{i\tilde{H}[\tilde{L},\eta\ell_{0}]t}$ as follows:

$\displaystyle
e^{-i\tilde{H}[\tilde{L},\eta\ell_{0}]t}=e^{-i\tilde{V}[\tilde{L},\eta\ell_{0}]t}\
\mathcal{T}\exp\left(-i\int_{0}^{t}e^{i\tilde{V}[\tilde{L},\eta\ell_{0}]x}\tilde{H}_{0}[\tilde{L},\eta\ell_{0}]e^{-i\tilde{V}[\tilde{L},\eta\ell_{0}]x}dx\right),$
(S.87)

where $\mathcal{T}$ is the time-ordering operator, and we use the definition
(S.V.2). Because $\tilde{V}[\tilde{L},\eta\ell_{0}]$ consists of operators
that commute with each other, the time-evolved operator

$\displaystyle
O_{X}(\tilde{V}[\tilde{L},\eta\ell_{0}],\eta\ell_{0}],t)=e^{i\tilde{V}_{X[k]}[\tilde{L},\eta\ell_{0}]t}O_{X}e^{-i\tilde{V}_{X[k]}[\tilde{L},\eta\ell_{0}]t}$
(S.88)

is supported on the subset $X[k]$, where
$\tilde{V}_{X[k]}[\tilde{L},\eta\ell_{0}]$ picks up all the interaction terms
$v_{Z}[\tilde{L},\eta\ell_{0}]$ such that $Z\subset X[k]$ [see also Eqs. (S.9)
and (S.11)]. Because of the condition (S.85), we have
$L_{1}=X[\ell_{0}]\supset X[k]$, and we write
$O_{X}(\tilde{V}[\tilde{L},\eta\ell_{0}],\eta\ell_{0}],t)$ as

$\displaystyle
O_{X}(\tilde{V}[\tilde{L},\eta\ell_{0}],\eta\ell_{0}],t)=\tilde{O}_{L_{1}}\quad(\|\tilde{O}_{L_{1}}\|=\|O_{X}\|=\zeta_{0}).$
(S.89)

In addition, for arbitrary $0\leq\tau\leq t$,
$e^{-i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}\tilde{H}_{0}[\tilde{L},\eta\ell_{0}]e^{i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}$
is formally described by

$\displaystyle\tilde{H}_{\tau}:=e^{i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}\tilde{H}_{0}[\tilde{L},\eta\ell_{0}]e^{-i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}$
$\displaystyle=\sum_{\langle
i,j\rangle}J_{i,j}e^{i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}\bar{\Pi}_{\tilde{L},\eta\ell_{0}}(b_{i}^{\dagger}b_{j}+{\rm
h.c.})\bar{\Pi}_{\tilde{L},\eta\ell_{0}}e^{-i\tilde{V}[\tilde{L},\eta\ell_{0}]\tau}$
$\displaystyle=:\sum_{Z\subset\Lambda:{\rm diam}(Z)\leq
2k}\tilde{h}_{Z,\tau},$ (S.90)

where we use the fact that $\tilde{V}[\tilde{L},\eta\ell_{0}]$ consists of
interaction terms with a maximum interaction length of $k$. For an arbitrary
time-dependent operator $A_{\tau}$, we adopt the following notations:

$\displaystyle U_{A_{\tau},x\to
t}=\mathcal{T}e^{-i\int_{x}^{t}A_{\tau}d\tau},$ $\displaystyle O(A_{\tau},x\to
t)=U_{A_{\tau},x\to t}^{\dagger}OU_{A_{\tau},x\to t}.$ (S.91)

Using these notations, we obtain

$\displaystyle
O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)=\tilde{O}_{L_{1}}(\tilde{H}_{\tau},0\to
t),$ (S.92)

where $\tilde{O}_{L_{1}}$ has been defined in Eq. (S.89). In the following, we
approximate $\tilde{O}_{L_{1}}(\tilde{H}_{\tau},0\to t)$ using the subset
Hamiltonian $\tilde{H}_{L,\tau}$:

$\displaystyle\tilde{O}_{L_{1}}(\tilde{H}_{\tau},0\to
t)\approx\tilde{O}_{L_{1}}(\tilde{H}_{L,\tau},0\to
t),\quad\tilde{H}_{L,\tau}=\sum_{Z\subset L}\tilde{h}_{Z,\tau},$ (S.93)

where we choose subset $L$ appropriately as described later. In the following
proposition, we estimate the approximation error (see Sec. S.VIII for the
proof).

######  Proposition 6.

When $L=L_{2}^{\prime}:=X[2\ell_{0}-2k]$ ($\subset L_{2}$) is chosen, the
approximation error in Eq. (S.93) is bounded from above by

$\displaystyle\left\|\tilde{O}_{L_{1}}(\tilde{H}_{\tau},0\to
t)-\tilde{O}_{L_{1}}(\tilde{H}_{L^{\prime}_{2},\tau},0\to t)\right\|\leq
2e^{3}\zeta_{0}c_{3}t|\partial L_{2}^{\prime}|\ell_{0}e^{-\ell_{0}/(2k)}$
(S.94)

under the condition

$\displaystyle t\leq\frac{1}{ec_{3}^{\prime}},$ (S.95)

where $c_{3}:=4\bar{J}\eta\gamma(2k)^{D}d_{G}$, and
$c_{3}^{\prime}:=16ekc_{3}\gamma(2k)^{D}$.

From the above proposition, by choosing $U_{L_{2}}$ as

$\displaystyle
U_{L_{2}}=e^{-i\tilde{V}_{X[k]}[\tilde{L},\eta\ell_{0}]t}U_{\tilde{H}_{L^{\prime}_{2},\tau},0\to
t},$ (S.96)

we find that

$\displaystyle[U_{L_{2}},\hat{n}_{L_{2}}]=0$ (S.97)

because of $[\tilde{V}_{X[k]}[\tilde{L},\eta\ell_{0}],\hat{n}_{L_{2}}]=0$ and
$[\tilde{H}_{L^{\prime}_{2},\tau},\hat{n}_{L_{2}}]=0$. We upper-bound the norm
(S.86) as

$\displaystyle\left\|O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right\|\leq
2e^{3}\zeta_{0}c_{3}t|\partial L_{2}^{\prime}|\ell_{0}e^{-\ell_{0}/(2k)}\leq
2e^{3}\zeta_{0}c_{3}\gamma t(2\ell_{0}+r)^{D}\ell_{0}e^{-\ell_{0}/(2k)},$
(S.98)

where the last inequality is derived from $\partial L_{2}^{\prime}\subset
L_{2}=X[2\ell_{0}]\subseteq i_{0}[2\ell_{0}+r]$.

###  S.V.4 Completing the proof

We now have all the ingredients to prove Subtheorem 1. First, we set the
parameter $\Delta t_{0}$ in the statement to $\Delta
t_{0}=1/(ec_{3}^{\prime})$, which is an $\mathcal{O}(1)$ constant from the
definition of $c_{3}^{\prime}$ in Proposition 6. By choosing $\ell_{0}$ such
that it satisfies the conditions (S.83) and (S.85), we obtain

$\displaystyle\left\|\left[O_{X}(t)-O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\right]\rho_{0}(-t)\right\|_{1}\leq\frac{1}{2}e^{c_{0}\bar{q}}\zeta_{0}e^{-2\ell_{0}/\log(r)}$
(S.99)

and

$\displaystyle\left\|O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right\|\leq
2e^{3}\zeta_{0}c_{3}\gamma t(2\ell_{0}+r)^{D}\ell_{0}e^{-\ell_{0}/(2k)},$
(S.100)

from Propositions 5 and 6, respectively. By combining them, we obtain

$\displaystyle\left\|\left[O_{X}(t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right]\rho_{0}(-t)\right\|_{1}\leq\frac{1}{2}\zeta_{0}e^{c_{0}\bar{q}}e^{-2\ell_{0}/\log(r)}+2e^{3}\zeta_{0}c_{3}\gamma
t(2\ell_{0}+r)^{D}\ell_{0}e^{-\ell_{0}/(2k)}.$ (S.101)

For the second term, there exists a constant $\delta c_{2}=\mathcal{O}(1)$
such that for $\ell_{0}\geq(c_{2}+\delta c_{2})\log^{2}(r)$,

$\displaystyle 2e^{3}\zeta_{0}c_{3}\gamma
t(2\ell_{0}+r)^{D}\ell_{0}e^{-\ell_{0}/(2k)}\leq\frac{1}{2}\zeta_{0}e^{-2\ell_{0}/\log(r)}\leq\frac{1}{2}e^{c_{0}\bar{q}}\zeta_{0}e^{-2\ell_{0}/\log(r)},$
(S.102)

which reduces Ineq. (S.101) to

$\displaystyle\left\|\left[O_{X}(t)-U_{L_{2}}^{\dagger}O_{X}U_{L_{2}}\right]\rho_{0}(-t)\right\|_{1}\leq
e^{c_{0}\bar{q}}\zeta_{0}e^{-2\ell_{0}/\log(r)}.$ (S.103)

Note that $\delta c_{2}$ does not depend on $\bar{q}$. By letting
$2\ell_{0}=\ell$ in this inequality (i.e., $L_{2}=X[2\ell_{0}]=X[\ell]$), we
obtain the main inequality (S.43). All the conditions for $\ell_{0}=\ell/2$
can be written in the form of (S.41) by choosing $C_{0}=\mathcal{O}(1)$
appropriately.

This completes the proof of Subtheorem 1. $\square$

##  S.VI Proof of Proposition 3: boson density after time evolution

###  S.VI.1 Restatement

Proposition 3. We first define the operator $M^{(s)}_{i}(t)$ as

$\displaystyle M^{(s)}_{i}(t):={\rm
tr}\left[\hat{n}_{i}^{s}\tilde{\rho}(t)\right].$ (S.104)

Then, for $t\leq t_{0}$, the following upper bound for $M^{(s)}_{i}(t)$ holds:

$\displaystyle M_{i}^{(s)}(t)$ $\displaystyle\leq
c^{\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}|X|^{3}(c_{1}s|X|)^{s}e^{-d_{i,X}}+c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}(c_{1}s)^{s},$
(S.105)

where $c_{1}$, $c^{\prime}_{1}$, and $c^{\prime\prime}_{1}$ are defined as

$\displaystyle c_{1}:=e^{8\bar{J}d_{G}t_{0}}/c_{0},$ $\displaystyle
c^{\prime}_{1}:=320c_{0}^{-3}e^{4\bar{J}d_{G}t_{0}+c_{0}(1+q_{0}/|X|)},$
$\displaystyle
c^{\prime\prime}_{1}:=80\lambda_{0}c_{0}^{-1}e^{4\bar{J}d_{G}t_{0}+c_{0}}.$
(S.106)

Note that they are $\mathcal{O}(1)$ constants if $t_{0}=\mathcal{O}(1)$.

###  S.VI.2 Proof

For the proof, we start from the differential equation for $M^{(s)}_{i}(t)$:

$\displaystyle\frac{d}{dt}M^{(s)}_{i}(t)$ $\displaystyle=-i{\rm
tr}\left(\hat{n}_{i}^{s}[H,\tilde{\rho}(t)]\right)=i{\rm
tr}\left([H,\hat{n}_{i}^{s}]\tilde{\rho}(t)\right).$ (S.107)

The form of the Hamiltonian (S.9) gives

$\displaystyle[H,\hat{n}_{i}^{s}]=[H_{0},\hat{n}_{i}^{s}]=\sum_{\langle
j,j^{\prime}\rangle}J_{j,j^{\prime}}[b_{j}b_{j^{\prime}}^{\dagger}+{\rm
h.c.},\hat{n}_{i}^{s}].$ (S.108)

Note that $[V,\hat{n}_{i}^{s}]=0$ because $[\hat{n}_{i},\hat{n}_{j}]=0$ for
$\forall j\in\Lambda$.

For the convenience of readers, we first consider the case of $s=1$, where
$M^{(1)}_{i}(t)$ gives the expectation value of $\hat{n}_{i}$. This case was
considered in Ref. [67]. We prove the following lemma.

######  Lemma 7.

For arbitrary site $i\in\Lambda$, the first-order moment $M^{(1)}_{i}(t)$ is
upper-bounded by

$\displaystyle M^{(1)}_{i}(t)$ $\displaystyle\leq
10\tilde{N}(d_{i,X})e^{3\bar{J}d_{G}t},\quad\tilde{N}(d_{i,X}):=N_{X}e^{-d_{i,X}}+n_{0}\lambda_{0},$
(S.109)

where we define $N_{X}$ and $n_{0}$ as $N_{X}:={\rm tr}(\rho_{0}\hat{n}_{X})$
and $n_{0}:=\max_{i\in\Lambda}{\rm tr}(\rho_{0}\hat{n_{i}})$, respectively.

####  S.VI.2.1 Proof of Lemma 7

Using the relation $[b_{i},\hat{n}_{i}]=b_{i}$ (or
$[b^{\dagger}_{i},\hat{n}_{i}]=-b_{i}^{\dagger}$), we obtain

$\displaystyle[H_{0},\hat{n}_{i}]=\sum_{j:d_{i,j}=1}J_{i,j}(b_{i}b_{j}^{\dagger}-{\rm
h.c.}),$ (S.110)

which can be used to reduce Eq. (S.107) with $s=1$ to

$\displaystyle\frac{d}{dt}M^{(1)}_{i}(t)$
$\displaystyle=-\sum_{j:d_{i,j}=1}2J_{i,j}{\rm Im}{\rm
tr}\left[\left(b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\right)\right].$ (S.111)

The Cauchy–Schwarz inequality gives an upper bound of

$\displaystyle\left|{\rm
tr}\left[\left(b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\right)\right]\right|\leq\sqrt{{\rm
tr}\left[\left(b_{i}^{\dagger}b_{i}\tilde{\rho}(t)\right)\right]\cdot{\rm
tr}\left[\left(b_{j}^{\dagger}b_{j}\tilde{\rho}(t)\right)\right]}\leq\frac{M^{(1)}_{i}(t)+M^{(1)}_{j}(t)}{2},$
(S.112)

and hence from $|J_{i,j}|\leq\bar{J}$

$\displaystyle\left|\frac{d}{dt}M^{(1)}_{i}(t)\right|$
$\displaystyle\leq\bar{J}\sum_{j:d_{i,j}=1}[M^{(1)}_{i}(t)+M^{(1)}_{j}(t)].$
(S.113)

We then give the upper bound of
$\vec{M}^{(1)}(t)=\\{M^{(1)}_{i}(t)\\}_{i=1}^{n}$ as

$\displaystyle\vec{M}^{(1)}(t)\leq
e^{\bar{J}d_{G}t}e^{\bar{J}\mathcal{M}t}\vec{M}^{(1)}(0),$ (S.114)

where $d_{G}$ is the maximum degree of the graph, and the matrix $\mathcal{M}$
has nonzero elements only for $d_{i,j}=1$ with $\mathcal{M}_{i,j}=1$. As shown
in Ref. [98, 67], the matrix $e^{\bar{J}\mathcal{M}t}$ satisfies

$\displaystyle[e^{\bar{J}\mathcal{M}t}]_{i,j}\leq
Ce^{\tilde{v}_{0}t-d_{i,j}},$
$\displaystyle\tilde{v}_{0}=\chi\bar{J}\Delta,\quad\chi\approx
3.59,\quad\Delta=\frac{\|M\|}{2}\leq d_{G}/2,\quad
C=\frac{2\chi^{2}}{\chi-1}\approx 10.$ (S.115)

Then, because $M^{(1)}_{i}(0)\leq n_{0}$ for $i\notin X$, the upper bound of
$M^{(1)}_{i}(t)$ is given by

$\displaystyle M^{(1)}_{i}(t)$ $\displaystyle\leq
e^{\bar{J}d_{G}t}\sum_{j\in\Lambda}Ce^{\tilde{v}_{0}t-d_{i,j}}M^{(1)}_{j}(0)$
$\displaystyle\leq
10e^{\tilde{v}_{0}t+\bar{J}d_{G}t}\left(N_{X}e^{-d_{i,X}}+n_{0}\sum_{j\in\Lambda}e^{-d_{i,j}}\right)=:10\tilde{N}(d_{i,X})e^{(\tilde{v}_{0}+\bar{J}d_{G})t}\leq
10\tilde{N}(d_{i,X})e^{3\bar{J}d_{G}t},$ (S.116)

where we define $\tilde{N}(d_{i,X})=N_{X}e^{-d_{i,X}}+n_{0}\lambda_{0}$ [see
Eq. (S.5) for the definition of $\lambda_{0}$].

[ End of Proof of Lemma 7]

For general $s$, we use a similar approach to obtain the upper bound. By using
the relation $[b_{i},\hat{n}_{i}]=b_{i}$ (or
$[b^{\dagger}_{i},\hat{n}_{i}]=-b_{i}^{\dagger}$), we first prove the
following lemma.

######  Lemma 8.

For an arbitrary function $f(x)$, the commutator $[b_{i},f(\hat{n}_{i})]$ is
given by

$\displaystyle[b_{i},f(\hat{n}_{i})]=[f(\hat{n}_{i}+1)-f(\hat{n}_{i})]b_{i}$
(S.117)

or

$\displaystyle[b_{i},f(\hat{n}_{i})]=b_{i}[f(\hat{n}_{i})-f(\hat{n}_{i}-1)].$
(S.118)

####  S.VI.2.2 Proof of Lemma 8

For the proof, let us define $|q,r\rangle$ as an eigenstate of $n_{i}$ as
$n_{i}|q,r\rangle=q|q,r\rangle$, where the index $r$ indicates the degenerate
eigenstate. Then, we have

$\displaystyle[b_{i},f(\hat{n}_{i})]|q,r\rangle$
$\displaystyle=b_{i}f(q)|q,r\rangle-f(\hat{n}_{i})b_{i}|q,r\rangle$
$\displaystyle=f(q)\sqrt{q}|q-1,r\rangle-f(q-1)\sqrt{q}|q-1,r\rangle$
$\displaystyle=[f(q)-f(q-1)]\sqrt{q}|q-1,r\rangle$
$\displaystyle=[f(\hat{n}_{i}+1)-f(\hat{n}_{i})]\sqrt{q}|q-1,r\rangle=[f(\hat{n}_{i}+1)-f(\hat{n}_{i})]b_{i}|q,r\rangle.$
(S.119)

This equation holds for arbitrary eigenstates $|q,r\rangle$, and we obtain Eq.
(S.117). For the proof of Eq. (S.118), we take the same approach:

$\displaystyle[b_{i},f(\hat{n}_{i})]|q,r\rangle$
$\displaystyle=[f(q)-f(q-1)]\sqrt{q}|q-1,r\rangle$
$\displaystyle=b_{i}[f(q)-f(q-1)]|q,r\rangle=b_{i}[f(\hat{n}_{i})-f(\hat{n}_{i}-1)]|q,r\rangle.$
(S.120)

We thus prove Lemma 8. $\square$

[ End of Proof of Lemma 8]

To obtain the upper bounds for a higher-order moment $M^{(s)}_{i}(t)$, we need
to consider $[b_{i},\hat{n}_{i}^{s}]$ in Eq. (S.108). Using Lemma 8, we obtain

$\displaystyle[b_{i},\hat{n}_{i}^{s}]=\sum_{s_{1}=0}^{s_{1}-1}\binom{s}{s_{1}}\hat{n}_{i}^{s_{1}}b_{i}.$
(S.121)

From this equation, the $s$th order moment $M^{(s)}_{i}(t)$ depends on the
moments with lower degrees, $M^{(s^{\prime})}_{i}(t)$ ($s^{\prime}<s$). This
point complicates the analyses significantly.

To overcome this difficulty, let us define an $s$th-order function $f_{s}(x)$
that satisfies

$\displaystyle f_{s}(x+1)-f_{s}(x)=sx^{s-1},\quad f_{s}(0)=0.$ (S.122)

We can always find such a function by iteratively determining the coefficient
for $x^{s_{1}}$ ($s_{1}\leq s$) in $f_{s}(x)$. From Lemma 8, the function
$f_{s}(x)$ satisfies

$\displaystyle[b_{i},f_{s}(\hat{n}_{i})]=s\hat{n}_{i}^{s-1}b_{i}.$ (S.123)

Although the explicit form of $f_{s}(x)$ is not simple, as

$\displaystyle f_{1}(x)=x,\quad f_{2}(x)=x^{2}-x,\quad
f_{3}(x)=x^{3}-\frac{3}{2}x^{2}+\frac{1}{2}x,\quad
f_{4}(x)=x^{4}-2x^{3}+x^{2},\quad
f_{5}(x)=x^{5}-\frac{5}{2}x^{4}+\frac{5}{3}x^{3}-\frac{1}{6}x^{2},$
$\displaystyle f_{6}(x)=x^{6}-3x^{5}+\frac{5}{2}x^{4}-\frac{1}{2}x^{2},\quad
f_{7}(x)=x^{7}-\frac{7}{2}x^{6}+\frac{7}{2}x^{5}-\frac{7}{6}x^{3}+\frac{1}{6}x,\quad\cdots,$
(S.124)

we can prove the following lemma on the properties of the function $f_{s}(x)$.

######  Lemma 9.

For an arbitrary positive integer $m\in\mathbb{N}$, we prove

$\displaystyle(m-1)^{s}\leq f_{s}(m)\leq m^{s}.$ (S.125)

In addition, the following inequality holds:

$\displaystyle f_{s}(m)+\frac{s^{s}}{4}\geq\frac{m^{s}}{4}.$ (S.126)

####  S.VI.2.3 Proof of Lemma 9

For the proof, we use Eq. (S.122). First, for $m=1$, we have $f_{s}(1)=0$
because

$\displaystyle f_{s}(1)=f_{s}(0)+s\cdot 0^{s-1}=0,$ (S.127)

where we use $f_{s}(0)=0$. For $m=2$, from $f_{s}(x+1)-f_{s}(x)=sx^{s-1}$, we
have

$\displaystyle f_{s}(2)=f_{s}(1)+s=s.$ (S.128)

Similarly, for $m=3$ and $m=4$, we have

$\displaystyle f_{s}(3)=f_{s}(2)+s\cdot 2^{s-1}=s(1+2^{s-1})$ (S.129)

and

$\displaystyle f_{s}(4)=f_{s}(3)+s\cdot 3^{s-1}=s(1+2^{s-1}+3^{s-1}).$ (S.130)

By repeating this procedure, we obtain

$\displaystyle f_{s}(m)=s\sum_{j=1}^{m-1}j^{s-1}.$ (S.131)

For $s\geq 1$, when we use

$\displaystyle\sum_{j=1}^{m-1}j^{s-1}\leq\int_{1}^{m}x^{s-1}dx=\frac{m^{s}-1}{s},$
$\displaystyle\sum_{j=1}^{m-1}j^{s-1}\geq\int_{0}^{m-1}x^{s-1}dx=\frac{(m-1)^{s}}{s},$
(S.132)

the function $f_{s}(m)$ is bounded from above/below by

$\displaystyle f_{s}(m)\leq m^{s}-1\leq m^{s}$ (S.133)

and

$\displaystyle f_{s}(m)\geq(m-1)^{s}.$ (S.134)

We thus prove the first part (S.125) of the lemma.

Next, we prove Ineq. (S.126). By using Ineq. (S.134), we obtain the following
for arbitrary $\kappa$:

$\displaystyle f_{s}(m)+\kappa\geq(m-1)^{s}+\kappa.$ (S.135)

We then prove that for $\kappa=s^{s}/4$, the inequality

$\displaystyle(m-1)^{s}+\frac{s^{s}}{4}\geq\frac{m^{s}}{4}$ (S.136)

holds. After proving the inequality (S.136), we can obtain the main inequality
(S.126) by using Ineq. (S.135) as $f_{s}(m)+s^{s}/4\geq(m-1)^{s}+s^{s}/4\geq
m^{s}/4$.

In order to prove the inequality (S.136), we first note that the inequality
(S.136) holds trivially for $s\geq m$. We thus consider the case of $m\geq s$.
The cases of $m\geq 1$ and $s=1$ are trivial; hence, we need to consider the
case of $m\geq s\geq 2$. Ineq. (S.136) reduces to

$\displaystyle\frac{s^{s}}{4m^{s}}\geq\frac{1}{4}-(1-1/m)^{s}.$ (S.137)

Because $(1-1/m)^{s}\geq 1/4$ for $m\geq s\geq 2$, the RHS of the above
inequality becomes negative, and hence it always holds for $m\geq s\geq 2$.

This completes the proof of the lemma. $\square$

[ End of Proof of Lemma 9]

In the following, we consider the time evolution of $f_{s}(\hat{n}_{i})$
instead of $\hat{n}_{i}^{s}$, which we define as

$\displaystyle F^{(s)}_{i}(t)={\rm
tr}\left[f_{s}(\hat{n}_{i})\tilde{\rho}(t)\right].$ (S.138)

Then, using Eq. (S.123), we have

$\displaystyle[H_{0},f_{s}(\hat{n}_{i})]$
$\displaystyle=\sum_{j:d_{i,j}=1}sJ_{i,j}[b_{i}b_{j}^{\dagger}+{\rm
h.c.},f_{s}(\hat{n}_{i})]=\sum_{j:d_{i,j}=1}sJ_{i,j}\left(\hat{n}_{i}^{s-1}b_{i}b_{j}^{\dagger}-{\rm
h.c.}\right),$ (S.139)

where we use

$\displaystyle[b_{i}^{\dagger},f_{s}(\hat{n}_{i})]=-\left([b_{i},f_{s}(\hat{n}_{i})]\right)^{\dagger}=-sb_{i}^{\dagger}\hat{n}_{i}^{s-1}.$
(S.140)

To estimate the upper bound of

$\displaystyle\left|\frac{d}{dt}F^{(s)}_{i}(t)\right|$
$\displaystyle=\left|{\rm
tr}\bigl{(}[H,f_{s}(\hat{n}_{i})]\tilde{\rho}(t)\bigr{)}\right|\leq\sum_{j:d_{i,j}=1}|J_{i,j}|\cdot
2s\left|{\rm
tr}\bigl{(}\hat{n}_{i}^{s-1}b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\bigr{)}\right|,$
(S.141)

we need to obtain the upper bound of

$\displaystyle\left|{\rm
tr}\left(\hat{n}_{i}^{s-1}b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\right)\right|.$
(S.142)

Here we derive the following lemma.

######  Lemma 10.

For arbitrary integers $s$ and $s_{1}$ such that $s_{1}\leq s$, we obtain the
upper bound as

$\displaystyle\left|{\rm
tr}\left(\hat{n}_{i}^{s-s_{1}}b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\right)\right|\leq\left(1-\frac{1}{2(s-s_{1}+1)}\right)M_{i}^{(s-s_{1}+1)}(t)+\frac{1}{2(s-s_{1}+1)}M_{j}^{(s-s_{1}+1)}(t).$
(S.143)

####  S.VI.2.4 Proof of Lemma 10

From the Cauchy–Schwarz inequality, we have

$\displaystyle\left|{\rm
tr}\left(\hat{n}_{i}^{s-s_{1}}b_{i}b_{j}^{\dagger}\tilde{\rho}(t)\right)\right|$
$\displaystyle=\left|{\rm
tr}\left(\hat{n}_{i}^{(s-s_{1})/2}b_{i}\tilde{\rho}(t)b_{j}^{\dagger}\hat{n}_{i}^{(s-s_{1})/2}\right)\right|$
$\displaystyle\leq\sqrt{{\rm
tr}\left(\hat{n}_{i}^{(s-s_{1})/2}b_{j}\tilde{\rho}(t)b_{j}^{\dagger}\hat{n}_{i}^{(s-s_{1})/2}\right){\rm
tr}\left(\hat{n}_{i}^{(s-s_{1})/2}b_{i}\tilde{\rho}(t)b_{i}^{\dagger}\hat{n}_{i}^{(s-s_{1})/2}\right)}$
$\displaystyle\leq\frac{1}{2}{\rm
tr}\left[\hat{n}_{i}^{s-s_{1}}\hat{n}_{j}\tilde{\rho}(t)\right]+\frac{1}{2}{\rm
tr}\left(b_{i}^{\dagger}\hat{n}_{i}^{s-s_{1}}b_{i}\tilde{\rho}(t)\right).$
(S.144)

For the first term, we use the Hölder inequality to derive

$\displaystyle{\rm
tr}\left[\hat{n}_{i}^{s-s_{1}}\hat{n}_{j}\tilde{\rho}(t)\right]$
$\displaystyle\leq\left({\rm
tr}\left[\hat{n}_{i}^{p(s-s_{1})}\tilde{\rho}(t)\right]\right)^{1/p}\left({\rm
tr}\left[\hat{n}_{j}^{q}\tilde{\rho}(t)\right]\right)^{1/q}$
$\displaystyle=\left({\rm
tr}\left[\hat{n}_{i}^{s-s_{1}+1}\tilde{\rho}(t)\right]\right)^{(s-s_{1})/(s-s_{1}+1)}\left({\rm
tr}\left[\hat{n}_{j}^{s-s_{1}+1}\tilde{\rho}(t)\right]\right)^{1/(s-s_{1}+1)}$
$\displaystyle\leq\frac{s-s_{1}}{s-s_{1}+1}M_{i}^{(s-s_{1}+1)}(t)+\frac{1}{s-s_{1}+1}M_{j}^{(s-s_{1}+1)}(t),$
(S.145)

where we choose $p$ and $q$ such that $1/p=(s-s_{1})/(s-s_{1}+1)$ and
$1/q=1/(s-s_{1}+1)$, and we use the inequality of arithmetic and geometric
means as $x^{t}y^{1-t}\leq tx+(1-t)y$ for $0\leq t\leq 1$ and $x,y\geq 0$.

To estimate the second term, we consider the spectral decomposition of
$\tilde{\rho}(t)$ as

$\displaystyle\tilde{\rho}(t)=\sum_{q,r}\tilde{p}_{q,r}(t)|q,r\rangle\langle
q,r|,$ (S.146)

where $|q,r\rangle$ is an eigenstate of $\hat{n}_{i}$ satisfying
$\hat{n}_{i}|q,r\rangle=q|q,r\rangle$. We then obtain

$\displaystyle{\rm
tr}\left(b_{i}^{\dagger}\hat{n}_{i}^{s-s_{1}}b_{i}\tilde{\rho}(t)\right)$
$\displaystyle=\sum_{q,r|q\geq
1}\tilde{p}_{q,r}(t)q(q-1)^{s-s_{1}}|q,r\rangle\langle q,r|$
$\displaystyle\leq\sum_{q,r}\tilde{p}_{q,r}(t)q^{s-s_{1}+1}|q,r\rangle\langle
q,r|=M_{i}^{(s-s_{1}+1)}(t).$ (S.147)

By applying Ineqs. (S.145) and (S.147) to Eq. (S.144), we prove Ineq. (S.143).
$\square$

[ End of Proof of Lemma 10]

By applying Lemma 10 with $s_{1}=1$ to Eq. (S.142), we reduce Ineq. (S.141) to

$\displaystyle\left|\frac{d}{dt}F^{(s)}_{i}(t)\right|$ $\displaystyle\leq
2\bar{J}\sum_{j:d_{i,j}=1}\left(\frac{2s-1}{2}M_{i}^{(s)}(t)+\frac{1}{2}M_{j}^{(s)}(t)\right).$
(S.148)

So that the above inequality includes only
$\\{F^{(s)}_{i}(t)\\}_{i\in\Lambda}$, we use Ineq. (S.126) in Lemma 9:

$\displaystyle F^{(s)}_{i}(t)={\rm
tr}\left[f_{s}(\hat{n}_{i})\tilde{\rho}(t)\right]\geq{\rm
tr}\left[\frac{\hat{n}_{i}^{s}-s^{s}}{4}\tilde{\rho}(t)\right]=\frac{M_{i}^{(s)}(t)-s^{s}}{4},$
(S.149)

which yields

$\displaystyle M_{i}^{(s)}(t)\leq 4F^{(s)}_{i}(t)+s^{s}.$ (S.150)

By using the above bound, we reduce Ineq. (S.148) to

$\displaystyle\left|\frac{d}{dt}F^{(s)}_{i}(t)\right|$ $\displaystyle\leq
2\bar{J}\sum_{j:d_{i,j}=1}\left(\frac{2s-1}{2}[4F^{(s)}_{i}(t)+s^{s}]+\frac{1}{2}[4F^{(s)}_{k}(t)+s^{s}]\right).$
(S.151)

By defining $\tilde{F}^{(s)}_{i}(t)=F_{i}^{(s)}(t)+s^{s}/4$, we obtain

$\displaystyle\left|\frac{d}{dt}\tilde{F}^{(s)}_{i}(t)\right|$
$\displaystyle\leq
4\bar{J}\sum_{j:d_{i,j}=1}\left[(2s-1)\tilde{F}^{(s)}_{i}(t)+\tilde{F}^{(s)}_{k}(t)\right],$
(S.152)

where we use $\frac{d}{dt}\tilde{F}^{(s)}_{i}(t)=\frac{d}{dt}F^{(s)}_{i}(t)$.

Then, we can derive an inequality similar to Eq. (S.114):

$\displaystyle|\vec{\tilde{F}}^{(s)}(t)|\leq
e^{4(2s-1)\bar{J}d_{G}t}e^{4\bar{J}\mathcal{M}t}\vec{\tilde{F}}^{(s)}(0),$
(S.153)

where the matrix $\mathcal{M}$ has been defined in (S.114). We also obtain a
bound similar to Eq. (S.VI.2.1):

$\displaystyle[e^{4\bar{J}\mathcal{M}t}]_{i,j}\leq Ce^{v_{0}t-d_{i,j}},$
$\displaystyle v_{0}=4\chi\bar{J}\Delta,\quad\chi\approx
3.59,\quad\Delta=\frac{\|M\|}{2}\leq d_{G}/2,\quad
C=\frac{2\chi^{2}}{\chi-1}<10.$ (S.154)

We introduce the following upper bounds:

$\displaystyle\tilde{F}^{(s)}_{i}(0)\leq\tilde{F}^{(s)}_{0}\quad\textrm{for}\quad
i\in X^{\rm c},$ $\displaystyle\sum_{i\in
X}\tilde{F}^{(s)}_{i}(0)\leq\tilde{F}^{(s)}_{X}.$ (S.155)

We will calculate $\tilde{F}^{(s)}_{0}$ and $\tilde{F}^{(s)}_{X}$ explicitly
later. Using the above inequality, we obtain

$\displaystyle\tilde{F}^{(s)}_{i}(t)$ $\displaystyle\leq
10e^{4(2s-1)\bar{J}d_{G}t}\sum_{j\in\Lambda}e^{v_{0}t-d_{i,j}}\tilde{F}^{(s)}_{j}(0)$
$\displaystyle\leq
10e^{4(2s+1)\bar{J}d_{G}t}\left(\tilde{F}^{(s)}_{X}e^{-d_{i,X}}+\tilde{F}^{(s)}_{0}\sum_{j\in\Lambda}e^{-d_{i,j}}\right)\leq
10\tilde{\mathcal{F}}^{(s)}(d_{i,X})e^{4(2s+1)\bar{J}d_{G}t},$ (S.156)

where we use $v_{0}<8\bar{J}d_{G}$, and we define
$\tilde{\mathcal{F}}^{(s)}(d_{i,X})$ as

$\displaystyle\tilde{\mathcal{F}}^{(s)}(d_{i,X}):=\tilde{F}^{(s)}_{X}e^{-d_{i,X}}+\lambda_{0}\tilde{F}^{(s)}_{0}.$
(S.157)

Note that $\lambda_{0}$ has been defined in (S.5).

Because $\tilde{F}^{(s)}_{i}(t):=F_{i}^{(s)}(t)+s^{s}/4$ and
$F_{i}^{(s)}(t)\geq M_{i}^{(s)}(t)/4-s^{s}/4$ from Lemma 9, Ineq. (S.VI.2.4)
reduces to

$\displaystyle M_{i}^{(s)}(t)\leq
40\tilde{\mathcal{F}}^{(s)}(d_{i,X})e^{4(2s+1)\bar{J}d_{G}t}.$ (S.158)

Finally, we obtain the explicit forms of $\tilde{F}^{(s)}_{0}$ and
$\tilde{F}^{(s)}_{X}$ in Eq. (S.VI.2.4). First, from $F_{i}^{(s)}(0)\leq
M_{i}^{(s)}(0)$ in Lemma 9, we have

$\displaystyle\tilde{F}^{(s)}_{i}(0)=F_{i}^{(s)}(0)+\frac{s^{s}}{4}\leq
M^{(s)}_{i}(0)+\frac{s^{s}}{4}=\tilde{F}^{(s)}_{0}.$ (S.159)

Second, we have

$\displaystyle\sum_{i\in X}M^{(s)}_{i}(0)=\sum_{i\in X}{\rm
tr}(\hat{n}_{i}^{s}\tilde{\rho})\leq{\rm tr}\left[\left(\sum_{i\in
X}\hat{n}_{i}\right)^{s}\tilde{\rho}\right]=M^{(s)}_{X}(0),$ (S.160)

where we define $M^{(s)}_{X}(0)$ as

$\displaystyle M^{(s)}_{X}(0):={\rm
tr}\left(\hat{n}_{X}^{s}\tilde{\rho}\right).$ (S.161)

Using $F_{i}^{(s)}(0)\leq M_{i}^{(s)}(0)$ in Lemma 9, we obtain

$\displaystyle\sum_{i\in X}\tilde{F}^{(s)}_{i}(0)\leq\sum_{i\in
X}\left(M^{(s)}_{i}(0)+\frac{s^{s}}{4}\right)\leq
M^{(s)}_{X}(0)+\frac{|X|s^{s}}{4}=\tilde{F}^{(s)}_{X}(0).$ (S.162)

In summary, we obtained

$\displaystyle\tilde{F}^{(s)}_{0}=M^{(s)}_{i}(0)+\frac{s^{s}}{4},\quad\tilde{F}^{(s)}_{X}(0)=M^{(s)}_{X}(0)+\frac{|X|s^{s}}{4}.$
(S.163)

Therefore, to upper-bound the quantity $\tilde{\mathcal{F}}^{(s)}(d_{i,X})$ in
Eq. (S.157), we need to upper-bound $M^{(s)}_{i}(0)$ and $M^{(s)}_{X}(0)$. We
can prove the following lemma.

######  Lemma 11.

Let $\tilde{\rho}$ be defined as in Eq. (S.37), that is,
$\tilde{\rho}=O_{X}\rho O_{X}^{\dagger}$. Then, under the condition given in
Eq. (S.13), $M_{i}^{(s)}(0)$ ($i\in X^{\rm c}$) is upper-bounded by

$\displaystyle M_{i}^{(s)}(0)={\rm
tr}(\hat{n}_{i}^{s}\tilde{\rho})\leq\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s!c_{0}^{-s-1}.$
(S.164)

Next, for $M_{X}^{(s)}(0)$, we obtain the upper bound of

$\displaystyle M_{X}^{(s)}(0)\leq
4\zeta_{0}^{2}(|X|/c_{0})^{s+3}s!e^{c_{0}(\bar{q}+1+q_{0}/|X|)}.$ (S.165)

####  S.VI.2.5 Proof of Lemma 11

Let $\Pi_{i,q}$ be a projection onto the eigenspace of $\hat{n}_{i}$ with
eigenvalue $q$:

$\displaystyle\hat{n}_{i}\Pi_{i,q}=q\Pi_{i,q}.$ (S.166)

From Ineq. (S.13), we obtain

$\displaystyle{\rm
tr}(\Pi_{i,q}e^{c_{0}(\hat{n}_{i}-\bar{q})}\tilde{\rho}\Pi_{i,q})={\rm
tr}(O_{X}\Pi_{i,q}e^{c_{0}(\hat{n}_{i}-\bar{q})}\rho_{0}\Pi_{i,q}O_{X}^{\dagger})\leq\|O_{X}\|^{2}{\rm
tr}(e^{c_{0}(\hat{n}_{i}-\bar{q})}\rho_{0})\leq\zeta_{0}^{2},$
$\displaystyle{\rm
tr}(\Pi_{i,q}e^{c_{0}(\hat{n}_{i}-\bar{q})}\tilde{\rho}\Pi_{i,q})={\rm
tr}(e^{c_{0}(q-\bar{q})}\Pi_{i,q}\tilde{\rho}\Pi_{i,q}),$ (S.167)

where we use $[\Pi_{i,q},O_{X}]=0$ for $i\in X^{\rm c}$ and
$\|O_{X}\|=\|U_{X}^{\dagger}O_{X_{0}}U_{X}\|=\|O_{X_{0}}\|$ from Eq. (S.38).
By combining these two inequalities, we obtain

$\displaystyle{\rm
tr}(\Pi_{i,q}\tilde{\rho}\Pi_{i,q})\leq\zeta_{0}^{2}e^{-c_{0}(q-\bar{q})}.$
(S.168)

From this inequality, the $s$th order moment is bounded from above by

$\displaystyle{\rm tr}(\hat{n}_{i}^{s}\tilde{\rho})$
$\displaystyle=\sum_{q=1}^{\infty}{\rm
tr}(\hat{n}_{i}^{s}\Pi_{i,q}\tilde{\rho}\Pi_{i,q})\leq\zeta_{0}^{2}e^{c_{0}\bar{q}}\sum_{q=1}^{\infty}q^{s}e^{-c_{0}q}\leq\zeta_{0}^{2}e^{c_{0}\bar{q}}\int_{0}^{\infty}x^{s}e^{-c_{0}(x-1)}dx=\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s!c_{0}^{-s-1}.$
(S.169)

We thus prove Ineq. (S.164).

Next, we prove Ineq. (S.165). By using the condition in Eq. (S.17) and
$[U_{X},\hat{n}_{X}]=0$, we obtain

$\displaystyle{\rm tr}(\Pi_{X,q}\tilde{\rho}\Pi_{X,q})$
$\displaystyle=\|\Pi_{X,q}U_{X}^{\dagger}O_{X_{0}}U_{X}\rho_{0}U_{X}^{\dagger}u_{X_{0}}^{\dagger}U_{X}\Pi_{X,q}\|_{1}$
$\displaystyle=\|\Pi_{X,q}O_{X_{0}}\Pi_{X,\geq
q-q_{0}}U_{X}\rho_{0}U_{X}^{\dagger}\Pi_{X,\geq
q-q_{0}}O_{X_{0}}^{\dagger}\Pi_{X,q}\|_{1},$ (S.170)

where we define $\Pi_{X,\geq q-q_{0}}:=\sum_{s\geq q-q_{0}}\Pi_{X,s}$, and
$\|\cdot\|_{1}$ is the trace norm. Here, we use
$\Pi_{X,q}O_{X_{0}}\Pi_{X,<q-q_{0}}=0$. The above equation yields

$\displaystyle{\rm
tr}(\Pi_{X,q}\tilde{\rho}\Pi_{X,q})\leq\|O_{X_{0}}\|^{2}{\rm tr}(\Pi_{X,\geq
q-q_{0}}U_{X}^{\dagger}\rho_{0}U_{X}\Pi_{X,\geq q-q_{0}})=\zeta_{0}^{2}{\rm
tr}(\Pi_{X,\geq q-q_{0}}\rho_{0}\Pi_{X,\geq q-q_{0}})$ (S.171)

because $\|\Pi_{X,q}O_{X_{0}}\|\leq\|O_{X_{0}}\|$,
$\|O_{1}O_{2}O_{3}\|_{1}\leq\|O_{1}\|\cdot\|O_{3}\|\cdot\|O_{2}\|_{1}$, and
$[U_{X},\Pi_{X,\geq q-q_{0}}]=0$.

Next, we estimate the following probability for the projection $\Pi_{X,q}$
with respect to $\rho_{0}$:

$\displaystyle{\rm tr}(\Pi_{X,q^{\prime}}\rho_{0}\Pi_{X,q^{\prime}}).$ (S.172)

For an arbitrary eigenstate $|q^{\prime}\rangle$ such that
$\hat{n}_{X}|q^{\prime}\rangle=q^{\prime}|q^{\prime}\rangle$, at least one
site has more than $(q^{\prime}/|X|)$ bosons, and hence

$\displaystyle{\rm tr}(\Pi_{X,q^{\prime}}\rho_{0}\Pi_{X,q^{\prime}})$
$\displaystyle\leq\sum_{i\in X}\sum_{s\geq q^{\prime}/|X|}{\rm
tr}(\Pi_{i,s}\rho_{0}\Pi_{i,s})$ $\displaystyle\leq|X|\sum_{s\geq
q^{\prime}/|X|}e^{-c_{0}(s-\bar{q})}\leq|X|\frac{e^{c_{0}\bar{q}}}{1-e^{-c_{0}}}e^{-c_{0}q^{\prime}/|X|},$
(S.173)

where we use Ineq. (S.I.3) in the second inequality. Thus, from Eq. (S.171),
we obtain

$\displaystyle{\rm tr}(\Pi_{X,q}\tilde{\rho}\Pi_{X,q})$
$\displaystyle\leq\zeta_{0}^{2}|X|\frac{e^{c_{0}\bar{q}}}{1-e^{-c_{0}}}\sum_{q^{\prime}\geq
q-q_{0}}e^{-c_{0}q^{\prime}/|X|}\leq\frac{\zeta_{0}^{2}|X|e^{c_{0}\bar{q}}}{(1-e^{-c_{0}})(1-e^{-c_{0}/|X|})}e^{-\frac{c_{0}(q-q_{0})}{|X|}}.$
(S.174)

Using the above upper bound, we obtain an inequality similar to (S.169):

$\displaystyle{\rm tr}(\hat{n}_{X}^{s}\tilde{\rho})$
$\displaystyle\leq\frac{\zeta_{0}^{2}|X|e^{c_{0}\bar{q}}}{(1-e^{-c_{0}})(1-e^{-c_{0}/|X|})}\sum_{q=1}^{\infty}q^{s}e^{-\frac{c_{0}(q-q_{0})}{|X|}}$
$\displaystyle\leq\frac{\zeta_{0}^{2}|X|e^{c_{0}(\bar{q}+1+q_{0}/|X|)}}{(1-e^{-c_{0}})(1-e^{-c_{0}/|X|})}s!(|X|/c_{0})^{s+1}\leq
4\zeta_{0}^{2}(|X|/c_{0})^{s+3}s!e^{c_{0}(\bar{q}+1+q_{0}/|X|)},$ (S.175)

where we use $c_{0}\leq 1$ and the inequality $1/(1-e^{-x})\leq 2/x$ for
$x\leq 1$. This completes the proof. $\square$

[ End of Proof of Lemma 11]

By applying the lemma to Eq. (S.163), we obtain

$\displaystyle\tilde{F}^{(s)}_{0}(0)=\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s!c_{0}^{-s-1}+s^{s}/4,$
$\displaystyle\tilde{F}^{(s)}_{X}(0)=4\zeta_{0}^{2}(|X|/c_{0})^{s+3}s!e^{c_{0}(\bar{q}+1+q_{0}/|X|)}+|X|s^{s}/4.$
(S.176)

By combining the above inequalities with Eq. (S.157), we obtain the inequality

$\displaystyle\tilde{\mathcal{F}}^{(s)}(d_{i,X})$
$\displaystyle\leq\left[4\zeta_{0}^{2}(|X|/c_{0})^{s+3}s!e^{c_{0}(\bar{q}+1+q_{0}/|X|)}+|X|s^{s}/4\right]e^{-d_{i,X}}+\lambda_{0}(\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s!c_{0}^{-s-1}+s^{s}/4)$
$\displaystyle\leq
8\zeta_{0}^{2}(|X|/c_{0})^{s+3}s^{s}e^{c_{0}(\bar{q}+1+q_{0}/|X|)-d_{i,X}}+2\lambda_{0}\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s^{s}c_{0}^{-s-1},$
(S.177)

where we use $c_{0}\leq 1$, $\zeta_{0}\geq 1$, and $s!\leq s^{s}$. We thus
reduce Ineq. (S.158) to

$\displaystyle M_{i}^{(s)}(t)$ $\displaystyle\leq
320\zeta_{0}^{2}(|X|/c_{0})^{s+3}s^{s}e^{c_{0}(\bar{q}+1+q_{0}/|X|)-d_{i,X}}e^{4(2s+1)\bar{J}d_{G}t}+80\lambda_{0}\zeta_{0}^{2}e^{c_{0}(\bar{q}+1)}s^{s}c_{0}^{-s-1}e^{4(2s+1)\bar{J}d_{G}t}$
$\displaystyle\leq
320\zeta_{0}^{2}c_{0}^{-3}e^{4\bar{J}d_{G}t+c_{0}(\bar{q}+1+q_{0}/|X|)}|X|^{3}\left(\frac{s|X|e^{8\bar{J}d_{G}t}}{c_{0}}\right)^{s}e^{-d_{i,X}}+80\lambda_{0}\zeta_{0}^{2}c_{0}^{-1}e^{4\bar{J}d_{G}t+c_{0}(\bar{q}+1)}\left(\frac{se^{8\bar{J}d_{G}t}}{c_{0}}\right)^{s}.$

The RHS of this inequality increases monotonically with $t$. Hence, when
$c_{1}$, $c^{\prime}_{1}$, and $c^{\prime\prime}_{1}$ are defined as in Eq.
(S.VI.1), the above inequality reduces to the main inequality (S.105) for
$t\leq t_{0}$. This completes the proof of Proposition 3. $\square$

##  S.VII Proof of Proposition 5: effective Hamiltonian

###  S.VII.1 Restatement

Proposition 5 Let us choose $\ell_{0}$ so that it satisfies

$\displaystyle\ell_{0}\geq c_{2}\log^{2}(r).$ (S.178)

Recall that $r$ has been defined in Eq. (S.37). Then, there exists a constant
$\eta$ that gives

$\displaystyle\left\|\left[O_{X}(t)-O_{X}(\tilde{H}[\tilde{L},\eta\ell_{0}],t)\right]\rho_{0}(-t)\right\|_{1}\leq\frac{1}{2}e^{c_{0}\bar{q}}\zeta_{0}e^{-2\ell_{0}/\log(r)}$
(S.179)

for $t\leq t_{0}$. Here, $c_{2}$ and $\eta$ are $\mathcal{O}(1)$ constants
which do not depend on $\bar{q}$.

###  S.VII.2 Proof

First, we choose $\ell_{0}$ such that the condition (S.69) in Corollary 4
holds for $d_{i,X}=\ell_{0}$:

$\displaystyle\ell_{0}\geq
2\log\left(\frac{\gamma^{3}c^{\prime}_{1}}{c^{\prime\prime}_{1}}\right)+6D\log(r),$
(S.180)

where $r$ has been defined as $X\subseteq i_{0}[r]$ for an appropriate
$i_{0}\in\Lambda$. Then, for arbitrary $i\in L_{1}^{\rm c}$
($L_{1}=X[\ell_{0}]$ as in Eq. (S.77)), we obtain

$\displaystyle P_{i,\geq z_{0}}^{(t)}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}\left(\frac{\tilde{c}_{1}\ell_{0}}{z_{0}}\right)^{\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)}.$
(S.181)

For technical reasons, we adopt the following conditions on $\ell_{0}$ in
addition to Eq. (S.180):

$\displaystyle\ell_{0}\geq\frac{6\log
r}{\tilde{c}^{\prime}_{1}},\quad\ell_{0}\geq\log^{2}(r).$ (S.182)

We notice that the first condition is equivalent to
$\frac{1}{2}\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)\geq 3$. The conditions
(S.180) and (S.182) are satisfied by Ineq. (S.178) if we choose the parameter
$c_{2}$ appropriately.

We first notice that if we truncate the boson number using $\bar{\Pi}_{L,q}$
in Eq. (S.79), the quantum state $\tilde{\rho}(-t)$ is almost the same, i.e.,
$\bar{\Pi}_{L,q}\tilde{\rho}(-t)\approx\tilde{\rho}(-t)$, if $q$ is
sufficiently large. This point is rigorously justified by the following lemma.

######  Lemma 12.

For the time-evolved operator $O_{X}(t)$ with $t\leq t_{0}$, we obtain the
following inequality if $d_{X,L}\geq\ell_{0}$ with Eq. (S.180):

$\displaystyle\|(\bar{\Pi}_{L,q}-1)O_{X}(t)\rho_{0}(-t)\|_{1}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}|L|\left(\frac{\tilde{c}_{1}\ell_{0}}{q+1}\right)^{\frac{1}{2}\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}|L|\left(\frac{\tilde{c}_{1}\ell_{0}}{q}\right)^{\frac{1}{2}\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)}.$
(S.183)

####  S.VII.2.1 Proof of Lemma 12

We start with the Cauchy–Schwartz inequality as follows:

$\displaystyle\|(\bar{\Pi}_{L,q}-1)O_{X}(t)\rho_{0}(-t)\|_{1}$
$\displaystyle=\|(\bar{\Pi}_{L,q}-1)O_{X}(t)\sqrt{\rho_{0}(-t)}\sqrt{\rho_{0}(-t)}\|_{1}$
$\displaystyle\leq\sqrt{{\rm
tr}\left[(\bar{\Pi}_{L,q}-1)O_{X}(t)\rho_{0}(-t)O_{X}(t)^{\dagger}(\bar{\Pi}_{L,q}-1)\right]}\sqrt{{\rm
tr}[\rho_{0}(-t)]}$ $\displaystyle=\sqrt{{\rm
tr}\left[(1-\bar{\Pi}_{L,q})\tilde{\rho}(-t)\right]},$ (S.184)

where we use $\tilde{\rho}(-t)=O_{X}(t)\rho_{0}(-t)O_{X}(t)^{\dagger}$ from
Eq. (S.64), and $(1-\bar{\Pi}_{L,q})^{2}=1-\bar{\Pi}_{L,q}$. For an arbitrary
quantum state $\sigma$, we have

$\displaystyle{\rm tr}(\sigma\bar{\Pi}_{L,q})\geq 1-\sum_{i\in L}{\rm
tr}(\sigma\Pi_{i,>q}),$ (S.185)

and hence we obtain

$\displaystyle 1-{\rm tr}(\bar{\Pi}_{L,q}\tilde{\rho}(-t))$
$\displaystyle\leq\sum_{i\in L}{\rm tr}[\tilde{\rho}(-t)\Pi_{i,>q}]=\sum_{i\in
L}P_{i,>q}^{(t)}\leq
2c^{\prime\prime}_{1}e^{c_{0}\bar{q}}\zeta_{0}^{2}|L|\left(\frac{\tilde{c}_{1}\ell_{0}}{q+1}\right)^{\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)},$
(S.186)

where we use Eq. (S.181) in the last inequality, with $z_{0}=q$. By combining
Ineqs. (S.184) and (S.186), we prove the main inequality (S.183), where we use
$(2c^{\prime\prime}_{1}|L|)^{1/2}\leq 2c^{\prime\prime}_{1}|L|$ from
$c^{\prime\prime}_{1}\geq 1$ [see Eq. (3)]. This completes the proof.
$\square$

[ End of Proof of Lemma 12]

From Lemma 12, if the boson truncation $q$ is sufficiently large, we expect to
be able to approximate the time evolution $e^{-iHt}$ using the effective
Hamiltonian $\tilde{H}[L,q]$ defined in Eq. (S.V.2):

$\displaystyle\tilde{H}[L,q]:=\bar{\Pi}_{L,q}H\bar{\Pi}_{L,q}.$ (S.187)

A key ingredient is the following lemma, which specifies the approximation
error between $\tilde{\rho}(t)$ and $\tilde{\rho}(\tilde{H}[L,q],t)$.

######  Lemma 13.

For $t\leq t_{0}$, the difference between the time evolutions by $H$ and
$\tilde{H}[L,q]$ is bounded from above by

$\displaystyle\|O_{X}(t)\rho_{0}(-t)-O_{X}(\tilde{H}[L,q],t)\rho_{0}(-t)\|_{1}\leq
8\sqrt{2}e^{c_{0}\bar{q}}\zeta_{0}t_{0}c^{\prime\prime}_{1}d_{G}\bar{J}q|L|(2|L|+q)\left(\frac{\tilde{c}_{1}\ell_{0}}{q}\right)^{\frac{1}{2}\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)}.$
(S.188)

####  S.VII.2.2 Proof of Lemma 13

We start with the inequality

$\displaystyle\|O_{X}(t)\rho_{0}(-t)-O_{X}(\tilde{H}[L,q],t)\rho_{0}(-t)\|_{1}$
$\displaystyle\leq\left\|e^{iHt}O_{X}e^{-iHt}\rho_{0}(-t)-e^{i\tilde{H}[L,q]t}O_{X}e^{-iHt}\rho_{0}(-t)\right\|_{1}+\left\|e^{i\tilde{H}[L,q]t}O_{X}(e^{-iHt}-e^{-i\tilde{H}[L,q]t})\rho_{0}(-t)\right\|_{1}$
$\displaystyle\leq\left\|e^{iHt}O_{X}\rho_{0}-e^{i\tilde{H}[L,q]t}O_{X}\rho_{0}\right\|_{1}+\zeta_{0}\left\|\rho_{0}-e^{-i\tilde{H}[L,q]t}e^{iHt}\rho_{0}\right\|_{1}$
$\displaystyle\leq\left\|e^{iHt}O_{X}\rho_{0}-e^{i\tilde{H}[L,q]t}O_{X}\rho_{0}\right\|_{1}+\zeta_{0}\left\|e^{i\tilde{H}[L,q]t}\rho_{0}-e^{iHt}\rho_{0}\right\|_{1},$
(S.189)

where we use $\|O_{X}\|\leq\zeta_{0}$ in the second inequality.

To estimate the RHS of the above inequality, we need to upper-bound

$\displaystyle\|e^{iHt}O_{X}\rho_{0}-e^{i\tilde{H}[L,q]t}O_{X}\rho_{0}\|_{1}$
(S.190)

for general $O_{X}$. The second term on the RHS of Eq. (S.VII.2.2) also
reduces to the above form when $O_{X}=\hat{1}$ is chosen. Let us decompose the
time to $m_{0}$ pieces (i.e., $t=m_{0}dt$) and take the limit $dt\to 0$. To
estimate the norm
$\|e^{iHt}O_{X}\rho_{0}-e^{i\tilde{H}[L,q]t}O_{X}\rho_{0}\|$, we start with
the following identical equation:

$\displaystyle e^{iHt}O_{X}\rho_{0}-e^{i\tilde{H}[L,q]t}O_{X}\rho_{0}$
$\displaystyle=\sum_{m=0}^{m_{0}-1}\Bigl{(}e^{i\tilde{H}[L,q]mdt}e^{iH(m_{0}-m)dt}O_{X}-e^{i\tilde{H}[L,q](m+1)dt}e^{iH(m_{0}-1-m)dt}O_{X}\Bigr{)}\rho_{0}.$
(S.191)

To upper-bound the norm of the above operator, we would like to calculate the
norm of

$\displaystyle\left(e^{i\tilde{H}[L,q]t_{1}}e^{iHt_{2}}O_{X}-e^{i\tilde{H}[L,q](t_{1}+dt)}e^{iH(t_{2}-dt)}O_{X}\right)\rho_{0}$
(S.192)

for arbitrary $t_{1}$ and $t_{2}$ such that $t_{1}+t_{2}=t$.

From $e^{i\tilde{H}[L,q]t_{1}}=e^{i\tilde{H}[L,q]t_{1}}\bar{\Pi}_{L,q}$, we
first obtain

$\displaystyle e^{i\tilde{H}[L,q]t_{1}}e^{iHt_{2}}=$ $\displaystyle
e^{i\tilde{H}[L,q]t_{1}}\bar{\Pi}_{L,q}e^{iHdt}e^{iH(t_{2}-dt)}$
$\displaystyle=$ $\displaystyle
e^{i\tilde{H}[L,q]t_{1}}\left[\bar{\Pi}_{L,q}(1+iH_{0,L^{\rm
c}}+iV)\bar{\Pi}_{L,q}dt+\bar{\Pi}_{L,q}i(H_{0,L}+H_{0,\partial
L})dt\right]e^{iH(t_{2}-dt)}+\mathcal{O}(dt^{2})$ $\displaystyle=$
$\displaystyle
e^{i\tilde{H}[L,q]t_{1}}\left[\bar{\Pi}_{L,q}(1+iH)\bar{\Pi}_{L,q}dt+\bar{\Pi}_{L,q}i(H_{0,L}+H_{0,\partial
L})(1-\bar{\Pi}_{L,q})dt\right]e^{iH(t_{2}-dt)}+\mathcal{O}(dt^{2})$
$\displaystyle=$ $\displaystyle
e^{i\tilde{H}[L,q](t_{1}+dt)}e^{iH(t_{2}-dt)}+e^{i\tilde{H}[L,q]t_{1}}\left[\bar{\Pi}_{L,q}i(H_{0,L}+H_{0,\partial
L})(1-\bar{\Pi}_{L,q})dt\right]e^{iH(t_{2}-dt)}+\mathcal{O}(dt^{2}),$ (S.193)

where we use $[H_{0,L^{\rm c}},\bar{\Pi}_{L,q}]=0$, $[V,\bar{\Pi}_{L,q}]=0$,
and $\bar{\Pi}_{L,q}^{2}=\bar{\Pi}_{L,q}$. Then, the upper bound of the norm
of Eq. (S.192) is given by

$\displaystyle\left\|e^{i\tilde{H}[L,q]t_{1}}\left[\bar{\Pi}_{L,q}i(H_{0,L}+H_{0,\partial
L})(1-\bar{\Pi}_{L,q})\right]e^{iH(t_{2}-dt)}O_{X}\rho_{0}\right\|_{1}dt+\mathcal{O}(dt^{2})$
$\displaystyle\leq$
$\displaystyle\left\|\left[\bar{\Pi}_{L,q}(H_{0,L}+H_{0,\partial
L})(1-\bar{\Pi}_{L,q})\right]O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\right\|_{1}dt+\mathcal{O}(dt^{2})$
$\displaystyle\leq$
$\displaystyle\|\bar{\Pi}_{L,q}H_{0,L}\|\cdot\|(1-\bar{\Pi}_{L,q})O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\|_{1}dt$
$\displaystyle+\left\|\bar{\Pi}_{L,q}H_{0,\partial
L}(1-\bar{\Pi}_{L,q})O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\right\|_{1}dt+\mathcal{O}(dt^{2}),$
(S.194)

where we use
$e^{iH(t_{2}-dt)}O_{X}\rho_{0}=O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)e^{-iH(t_{2}-dt)}$
in the first inequality. By using the definition (S.79), we have

$\displaystyle\|\bar{\Pi}_{L,q}b_{i}b_{j}^{\dagger}\|\leq\sqrt{q(q+1)}\leq\sqrt{2}q$
(S.195)

for arbitrary $i,j\in L$, where we use $\|b_{i}\Pi_{i,z}\|\leq\sqrt{z}$ and
$\|b_{i}^{\dagger}\Pi_{i,z}\|\leq\sqrt{z+1}$. This inequality yields

$\displaystyle\|\bar{\Pi}_{L,q}H_{0,L}\|\leq\sqrt{2}q\sum_{\langle
i,j\rangle:i,j\in L}2|J_{i,j}|\leq 2\sqrt{2}d_{G}\bar{J}q|L|.$ (S.196)

Recall that $\bar{J}$ is the upper bound of $J_{i,j}$, and $d_{G}$ is the
maximum degree of the lattice. Hence, the first term in Eq. (S.VII.2.2) is
upper-bounded by

$\displaystyle\|\bar{\Pi}_{L,q}H_{0,L}\|\cdot\|(1-\bar{\Pi}_{L,q})O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\|_{1}dt\leq
4\sqrt{2}e^{c_{0}\bar{q}}c^{\prime\prime}_{1}\zeta_{0}d_{G}\bar{J}q|L|^{2}\left(\frac{\tilde{c}_{1}\ell_{0}}{q}\right)^{\frac{1}{2}\tilde{c}^{\prime}_{1}\ell_{0}/\log(r)},$
(S.197)

where we use Lemma 12 to upper-bound
$\|(1-\bar{\Pi}_{L,q})O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\|_{1}$ ($t_{2}\leq
t_{0}$).

We next consider the upper bound of the second term in Eq. (S.VII.2.2). Here,
the Hamiltonian $H_{0,\partial L}$ includes boundary interactions such as
$J_{i,j}b_{i}b_{j}^{\dagger}$ with $i\in L$ and $j\notin L$. Because the norm
of $\bar{\Pi}_{L,q}b_{i}b_{j}^{\dagger}$ ($i\in L$, $j\notin L$) is unbounded,
we need to consider

$\displaystyle\left\|\bar{\Pi}_{L,q}b_{i}b_{j}^{\dagger}(1-\bar{\Pi}_{L,q})O_{X}(t_{2}-dt)\rho_{0}(-t_{2}+dt)\right\|_{1}$