Sharp asymptotics on the compression of

two-layer neural networks

The work of M.H.A. and R.P. was performed while at ISTA. M.H.A., S.B., M.M., and R.P. were partially supported by the 2019 Lopez-Loreta Prize. The work of S.R. was supported by the MOST grant 110-2221-E-A49 -052.

Mohammad Hossein Amani,
Simone Bombari,
Marco Mondelli,
Rattana Pukdee, and
Stefano Rini

EPFL, Switzerland, email:

ISTA, Austria, email:

CMU, USA, email:

NYCU, Taiwan, email:

To appear at the IEEE Information Theory Workshop (ITW) 2022, Mumbai, India.

In this paper, we study the compression of a target two-layer neural network with $N$ nodes into a compressed network with $M<N$ nodes.

More precisely, we consider the setting in which the weights of the target network are i.i.d. sub-Gaussian, and we minimize the population $L_2$ loss between the outputs of the target and of the compressed network, under the assumption of Gaussian inputs.

By using tools from high-dimensional probability, we show that this non-convex problem can be simplified when the target network is sufficiently over-parameterized, and provide the error rate of this approximation as a function of the input dimension and $N$.

In this mean-field limit,
the simplified objective, as well as the optimal weights of the compressed network, does not depend on the realization of the target network, but only on expected scaling factors.

Furthermore, for networks with ReLU activation, we conjecture that the optimum of the simplified optimization problem is achieved by taking weights on the Equiangular Tight Frame (ETF), while the scaling of the weights and the orientation of the ETF depend on the parameters of the target network.

Numerical evidence is provided to support this conjecture.

§ INTRODUCTION

The compression of Deep Neural Networks (DNNs) plays a crucial role in many modern-day applications in which inference has to be performed at minimal computational cost.

This occurs, for instance, in the IoT paradigm: here, classic DNN tasks (e.g., image classification, speech recognition, and anomaly detection) have to be implemented over diverse platforms with limited computational capabilities.

In this scenario, large DNNs trained over vast datasets have to be severely compressed to fit the current hardware limitations.

Several empirical studies have shown that over-parameterized DNNs are often highly redundant and, hence,
smaller sub-networks can be identified which perform close to the original DNN [1].

This observation has motivated the emergence of a host of different approaches to identify such sub-networks

Despite the very promising empirical results, the ultimate performance of DNN compression is yet to be characterized.

This problem can be summarized as the following question:

What is the minimum performance loss that can be guaranteed when compressing any network of $N$ nodes into a smaller network of $M<N$ nodes?

Main contributions. In this paper, we address the above question
by investigating the

limiting performance when approximating
a two-layer target network with $N$ neurons via a two-layer compressed network with $M$ neurons.

More precisely, we assume that (i) the input features are Gaussian, and that (ii) the network weights are i.i.d.

Then, we consider the problem of determining the parameters of the compressed network which minimize the expected $L_2$ distance (population loss) between the outputs of the target and compressed network.

This optimization can be divided into the two sub-problems: that of determining the amplitudes and the directions of the weights.

While finding the amplitudes boils down to linear regression, finding the directions

is  challenging.

Thus, by using tools from high-dimensional probability, we show that this latter problem can be simplified in the asymptotic regime in which $N, d\gg 1$.

We highlight that the simplified objective and the optimal parameters of the pruned network do not depend on the realization of the target network.

Furthermore, if the activation function is ReLU, we formulate the Equiangular Tight Frame (ETF) conjecture, that is, we conjecture that the optimal choice of weights on the compressed network forms an ETF (see Fig. <ref>), while  the orientation of the ETF and the scaling of the neuronal outputs depends on the
parameters of the target network. We also provide numerical evidence supporting our conjecture.

Related work.

Early work used the second order Taylor approximation of the cost function to determine the relevance of the network weights [3, 4].

Using this change as a measure of salience, weights corresponding to low salience are pruned from the network.

The magnitude of the weights
has also been considered as a measure of saliency [5, 6, 7].

The main advantage of these methods relies on their
computational efficiency.

The so called “lottery ticket hypothesis” suggests that pruning should be performed in a structured manner [8],

and sparsification techniques which promote a structured form of sparsity have been considered e.g. in [9].
In [10], the authors consider a “greedy forward selection” strategy in which
the sub-network is obtained by successively adding (with replacement) neurons to an empty network. In [11], pruning is performed via

“neuron merging”, i.e., similar neurons are identified by the cosine similarity of their weights and they are successively merged.

We  adopt the short-hands  $[m:n] \triangleq \{m, \ldots, n\}$

and  $[n] \triangleq [1:n]$.
We let $\one$ be the all-1 vector, and $\Iv$ the identity matrix.

Given a vector $v$, $\|v\|_2$ denotes its $L_2$ norm.
Given a matrix $M$, $\|M\|_{\rm op}$ denotes its operator norm, $\|M\|_{\Fsf}$ the Frobenius norm,  and $M^+$ its
Moore-Penrose pseudo-inverse.

The unit sphere of dimension $d$ is indicated as $\Scal^{d-1}$. The Gaussian distribution with mean $\muv$ and covariance matrix $\Cv$ is denoted as
$\Ncal(\muv,\Cv)$. Given a sub-Gaussian random variable $w$, its sub-Gaussian norm is defined as $\|w\|_{\psi_2} = \inf \{ t>0 \,\,: \,\,\mathbb E[\exp(w^2/t^2)] \le 2 \}$. Given a random vector $\wv\in\mathbb R^d$, its sub-Gaussian norm is defined as $\|\wv\|_{\psi_2} := \sup_{\uv \,:\, \|\uv\|_2=1} \|\uv^\intercal \wv\|_{\psi_2}$.

[tdplot_main_coords, scale = 2.25]

[ball color = lightgray,
opacity = 0.5
] (0,0,0) circle (1cm);

(P1) at  (-0.81649658, 0.40824829, 0.40824829);
(P2) at  (0.40824829, - 0.81649658, 0.40824829);
(P3) at  (0.40824829, 0.40824829,  -0.81649658);

[fill = vamsiblue, vamsiblue ] (P1) circle (0.75pt);
[fill = vamsiblue,vamsiblue] (P2) circle (0.75pt);
[fill = vamsiblue,vamsiblue] (P3) circle (0.75pt);

[vamsiblue, line cap=round, line width = 1mm, line join=round, ] (0,0,0) – (P1);
[vamsiblue, line cap=round, line width = 1mm ] (0,0,0) – (P2);
[vamsiblue, line cap=round, line width = 1mm  ] (0,0,0) – (P3);

[thin, dashed,vamsiblue, opacity=0.5 ] (-0.81649658, 0.40824829, 0) –  (-0.816 49658, 0, 0);
[thin, dashed,vamsiblue, opacity=0.5] (-0.81649658, 0.40824829, 0) –  (0, 0.40824829, 0);
[thin, dashed,vamsiblue, opacity=0.5] (-0.81649658, 0.40824829, 0) – (P1);

[thin, dashed,vamsiblue, opacity=0.5 ] (0.40824829, - 0.81649658, 0) –  (0.40824829, 0, 0);
[thin, dashed,vamsiblue, opacity=0.5] (0.40824829, - 0.81649658, 0) –  (0, - 0.81649658, 0);
[thin, dashed,vamsiblue, opacity=0.5] (0.40824829, - 0.81649658, 0) – (P2);

[thin, dashed,vamsiblue, opacity=0.5 ] (0.40824829, 0.40824829, 0) –  (0.40824829, 0, 0);
[thin, dashed,vamsiblue, opacity=0.5] (0.40824829, 0.40824829, 0) –
(0, 0.40824829, 0);
[thin, dashed,vamsiblue, opacity=0.5] (0.40824829, 0.40824829, 0) – (P3);

[-stealth] (0,0,0) – (1.80,0,0);
[-stealth] (0,0,0) – (0,1.50,0);
[-stealth] (0,0,0) – (0,0,1.30);
[dashed, gray] (0,0,0) – (-1,0,0);
[dashed, gray] (0,0,0) – (0,-1,0);
[dashed, gray] (0,0,0) – (0,0,-1);

A representation of the equiangular tight frame  (ETF).

§ PROBLEM FORMULATION

Let the target network $\Wsf$ be specified by (i) a set of weight vectors $\{ \wv_n \}_{n\in [N]}$  with $\wv_n \in \Scal^{d-1}$, and (ii) a set of scaling coefficients $\{a_n\}_{n\in [N]}$ with $a_n \in \Rbb$.

Then, the output $f_{\Wsf}(\xv)$ is obtained as

f_() = 1 N ∑_n ∈[N] a_n σ(

⟨_n, ⟩),

where $\sigma$ is the activation function.

The compressed network $\Vsf$ is specified by the weights $\{\vv_m \}_{m\in [M]}$ with $\vv_m \in \Scal^{d-1}$ and scaling coefficients  $\{b_m \}_{m\in [M]}$, so that the output $f_{\Vsf}(\xv)$ is obtained as

f_() = 1 M ∑_m ∈[M] b_m σ(

⟨_m, ⟩).

As we allow for the coefficients $\{a_n\}_{n\in [N]}$ and $\{b_m \}_{m\in [M]}$, the assumption that $\wv_n, \vv_m \in \Scal^{d-1}$ is not particularly restrictive. Furthermore, if we assume that the activation function is ReLU ($\sigma(x)=\max(x, 0)$) then, by using the homogeneity of $\sigma$, one can incorporate the norm of the weights $\wv_n$ and $\vv_m$ into the coefficients $a_n$ and $b_m$, respectively.

Our goal is to compress (<ref>) into (<ref>) so that the expected $L_2$ loss between the two outputs when the input is Gaussian, i.e.,

L(,) =  _∼(, )    f_()-f_() _2^2  ,

is minimized over the choice of $\{b_m, \vv_m \}_{m\in [M]}$. We remark that, in the theoretical literature, it is common to assume that the input data is Gaussian, see e.g. <cit.>.
Furthermore, the average $L_2$ distance of two functions on a Gaussian input, as in (<ref>), can be regarded as a proxy for understanding how close they behave on normalized (i.e. whitened) data.

A compact formulation of the expression in (<ref>) using vector notation can be achieved by defining

g(,)    = _∼(, ) [σ(^⊺ ) σ(^⊺ )],

G_WW _ij     = g(_i, _j),   G_VV _ij  = g(_i, _j) ,

G_VW _ij      = g (_i, _j).

Similarly, we define the column vectors $\av$, $\bv$ as  $[\av]_n=a_n$ $[\bv]_m=b_m$.

Finally, defining
$\sv := \frac{1}{N} G_{VW} \av$,
(<ref>) becomes

L(,)    = 1 N^2 ^⊺ G_WW - 2 M ^⊺

+ 1 M^2 ^⊺ G_VV .

The entries of $G_{VV}, G_{WW}$, and $G_{VW}$ can be interpreted as correlations between neurons, and $[\sv]_i$ can be interpreted as the average correlation between the $i^{\rm th}$ neuron in the compressed network and all neurons in the target network.

In this problem definition, we use the term “compression” with the following acceptation.

The network $\Wsf$ has been trained to perform an inference task over a given dataset.

The problem at hand is to produce a smaller network $\Vsf$  approximating $\Wsf$ using only the knowledge of the network parameters.

In the remainder of the paper, we consider the following set of assumptions:

* Assumption (A-$\xv$): the network input $\Xv$  is Gaussian with zero mean and identity covariance.

* Assumption (A-$\wv$): the target network weights $\{ \wv_n \}_{n\in [N]}$  are i.i.d. zero-mean, unit-norm (i.e., $\wv_n\in \Scal^{d-1}$) sub-Gaussian random vectors, with a sub-Gaussian norm $\|\wv_n\|_{\psi_2}\le\sigma_W/\sqrt{d}$.[The scaling $1/\sqrt{d}$ ensures that the weights belong to the unit sphere $\Scal^{d-1}$.]

* Assumption (A-$\av$): the scaling coefficients of the target network $\{a_n\}_{n\in [N]}$ are i.i.d. bounded random variables, i.e. $|a_n| \leq A$, and they are independent of the weights $\{ \wv_n \}_{n\in [N]}$. We also assume $\Ebb[a_n] = \mu_a \neq 0$.

* Assumption (A-$\bv$): the scaling coefficients $\{b_m\}_{m\in [M]}$ are optimized in a compact set, i.e. $|b_m| \leq B$.

A few remarks on the assumptions are in order. First, the quantities $\sigma_W, \mu_a, A$ and $B$ are all numerical constants (independent of $d, N, M$). Second, the assumption that $\{ \wv_n \}_{n\in [N]}$  are i.i.d. is compatible with gradient descent training in the mean-field regime [13, 14, 15, 16]. In fact, in this regime, the weights obtained from stochastic gradient descent are close to $N$ i.i.d. particles distributed according to the solution of a certain non-linear partial differential equation [13, 17]. This perspective has been used to prove a number of properties of neural networks, including convergence [13, 15, 18], mode connectivity [19] and implicit bias [20].

§ OPTIMIZATION OF SCALING COEFFICIENTS

We begin
by arguing that the minimization over $\{b_m\}_{m\in [M]}$ in (<ref>)

for a given $\{\vv_m\}_{m\in [M]}$ has a rather simple closed-form expression (Lemma <ref>).

Next, we show that the optimization of the compressed weights $\{\vv_m\}_{m \in [M]}$
can be approached using the Hermite expansion of $\sigma$, when the input is Gaussian.

The optimal value of $\bv$ in  (<ref>) for fixed $\{ \vv_m \}_{m\in [M]}$ is

^*(G_VV)= M N G^+_VV G_VW =M G^+_VV .

Furthermore, when replacing the solution of (<ref>) into  (<ref>), the loss is given by

L(,) _=(<ref>)

= 1 N^2 ^⊺ G_WW -

G^+_VV .

Given any set of unit-norm vectors $\{ \vv_m \}_{m\in [M]}$ in (<ref>), the loss is a quadratic (hence convex) program in $\{ b_m \}_{m\in [M]}$.
Thus, the optimal value can be found by setting the gradient $\nabla_{\bv} L$ to $0$, which readily gives (<ref>).

By plugging (<ref>) into (<ref>), we also obtain (<ref>).

We note that, since by definition the loss is $\ge 0$, $\sv$ always lies in the image of $G_{VV}$, justifying the use of the Moore-Penrose pseudo-inverse.
In fact, if there exists $\kv \in \text{kernel}(G_{VV})$ such that $\langle \kv, \sv \rangle \neq 0$, then the loss evaluated at $\bv = \eta \kv$ is

\[
\frac{1}{N^2} \av^{\intercal} G_{WW} \av - \eta \frac{2}{M} \kv^{\intercal} \sv + \frac{\eta^2}{M^2} \kv^{\intercal} G_{VV} \kv = C - C^{'} \eta,
\]

which can be taken to $\pm \infty$ by taking $\eta$ to $\pm \infty$ depending on the sign of $C^{'}$.

The function $g(\uv,\vv)$

in (<ref>) is a function of $\alpha=\langle \uv, \vv \rangle$, which can be
expressed as

g(,)  = g(α)= 1 2π √(1-α^2)+ (π- arccos(α).

Additionally, $g(\alpha)$ admits the well-defined Taylor expansion

g(α)     =  1/2 π +
1/4 α+

∑_k ≥1 1/2 π ((2k-3)!!)^2/(2k)! α^2k

The expression in (<ref>) follows from
[We note that the expression in (<ref>) has a missing factor of $1/2$, which we corrected.]

The expression in (<ref>) is derived through some rather straightforward manipulations. Here, the Taylor coefficients of $g$ are the squared Hermite coefficients of $\sigma$, see <cit.>.

We note that, for any $\sigma\in L^2(e^{-w^2/2})$[$\sigma\in L^2(e^{-w^2/2})$ means that the integral $\int_{\mathbb R} \sigma^2(w)e^{-w^2/2}{\rm d}t$ is finite, which ensures that the Taylor expansion of $g$ is well-defined.], $g(\uv,\vv)$ is a function of $\alpha=\langle \uv, \vv \rangle$, as long as the input distribution is invariant under rotations.

LIMIT LOSS AND LIMIT OBJECTIVE

By solving the minimization over $\{b_m\}_{m\in [M]}$, we have reduced the problem to optimizing
(<ref>) over $\{ \vv_m \}_{m\in [M]}$. Note that the vectors $\{ \vv_m \}_{m\in [M]}$ appear both in $G^{+}_{VV}$ and in $\sv$. Next, we show that, as $N, d\gg 1$, the dependence of $\sv$ on $\{ \vv_m \}_{m\in [M]}$ can be essentially dropped. More specifically, Theorem <ref> proves that the loss (<ref>) is close to the limit loss

\begin{equation}
\tilde{L}(\Vsf,\Wsf) =
\f 1 {N^2} \av^{\intercal} G_{WW} \av- \f 2 {M}\f {\mu_a}{2\pi} \bv^{\intercal} \one
%\nonumber \\
%& \quad
+ \f 1 {M^2} \bv^{\intercal} G_{VV} \bv,
\label{eq:limit loss}
\end{equation}

obtained by plugging $\sv = \frac{\mu_a}{2 \pi}\one$ into (<ref>).
To prove the claim, the idea is to exploit the Taylor expansion (<ref>) of $g$, and show that only the zero-th order term has an impact on $\sv$ (as $N, d\gg 1$). To do so, we begin by deriving a bound on the magnitude of the linear term coming from the expansion in (<ref>).

Let the assumptions (A-$\wv$)-(A-$\av$) hold. Then, we have that

\[
\sup_{\vv \in \Scal^{d-1}} \labs \frac{1}{N} \sum_{n \in [N]}  a_n \langle \vv, \wv_n \rangle\rabs \leq \frac{t A\sigma_W}{\sqrt{N}},
\]

with probability at least $1 - 2 \exp(-t^2/16)$.

We first note that

\[
\sup_{\vv \in \Scal^{d-1}} \labs \frac{1}{N} \sum_{n\in [N]} a_n \langle \vv, \wv_n \rangle \rabs = \frac{1}{N} \left\| \dsum_{n\in [N]} a_n \wv_n \right\|_2.
\]

Since $\{\wv_n\}_{n\in [N]}$ are zero-mean sub-Gaussian vectors with norm $\sigma_W/\sqrt{d}$, $\wv_{\rm avg} := \f 1 N \sum_{n\in [N]} a_n \wv_n$ is also a sub-Gaussian random vector with sub-Gaussian norm $|| \wv_{\rm avg}||_{\psi_2} \leq A \sigma_W/(\sqrt{N d})$.

This can be seen by the following bound on the moment generating function: for any $\uv \in \Scal^{d-1}$, it holds that, for any $\lambda\in\mathbb R$,

e^λ^⊺ _avg     = ∏_n∈[N] _a_n, _n e^λ1 N a_n ^⊺ _n

≤∏_n∈[N] _a_n e^λ^2a^2_n N^2 σ^2_W d   ≤e^λ^2A^2 N σ^2_W d,

where in the first inequality we use $|| \wv_n||_{\psi_2} \leq \sigma_W/\sqrt{d}$ and Proposition 2.5.2 in [23], and in the second inequality we use that $|a_n|\le A$.
By <cit.>, we have that

|| _avg||_2 ≥tAσ_W √(N) ≤2 e^-t^2/16,

which gives the desired result.

Next, we bound the quadratic term of the expansion in (<ref>).

Let the assumptions (A-$\wv$)-(A-$\av$) hold. Then, we have that

\begin{equation}
\begin{split}
\sup_{\vv \in \Scal^{d-1}}& \f 1 N \sum_{n\in [N]} |a_n| \langle \vv, \wv_n \rangle^{2} \le  C A \lb\f {\sigma_W^2} {d} + \f 1 N +\f {1} {\sqrt{Nd}}\rb,
\label{eq:sum squared}
\end{split}
\end{equation}

with probability at least $1 - 2 \exp(-d)$, where $C$ is a numerical constant.

Manipulating the LHS of (<ref>), we have

\begin{equation}
\begin{split}
&\sup_{\vv \in \Scal^{d-1}}\f 1 N \sum_{n\in [N]} |a_n| \langle \vv, \wv_n \rangle^{2} \le \sup_{\vv \in \Scal^{d-1}}\f 1 N \sum_{n\in [N]} A \langle \vv, \wv_n \rangle^{2} \\
& \hspace{1em}\leq\sup_{\vv \in \Scal^{d-1}} A \labs \f 1 N \vv^{\intercal} \sum_{n\in [N]}\wv_n \wv_n^{\intercal} \vv \rabs   \leq \f A d \left\|\f 1 N\Wv\Wv^\intercal\right\|_{\rm op},
\end{split}
\end{equation}

where in the first inequality we use that $|a_n|\le A$, and in the last line we have defined $\Wv := [\sqrt{d}\wv_1, \ldots, \sqrt{d}\wv_N]$.

To bound the operator norm of $\f 1 N \Wv\Wv^\intercal$, let $\boldsymbol \Sigma$ be the second moment matrix of $\sqrt{d}\wv_n$, i.e., $\boldsymbol \Sigma:= d\mathbb E[\wv_n\wv_n^\intercal]$. Then, by  <cit.>, we have that

\begin{equation}\label{eq:bds0}
\left\|\f 1 N\Wv\Wv^\intercal\right\|_{\rm op}\le \left\|\boldsymbol \Sigma\right\|_{\rm op}+\max(\delta, \delta^2),
\end{equation}

with probability at least $1-2\exp(-c_1 t^2)$, where $\delta= C_1\sqrt{d/N}+t/\sqrt{N}$ and $c_1, C_1$ are numerical constants. To bound $\left\|\boldsymbol \Sigma\right\|_{\rm op}$, let $\bbu$ be the unitary eigenvector associated to the maximum eigenvalue of $\boldsymbol \Sigma$. Then,

\begin{equation}\label{eq:int1}
\left\|\boldsymbol \Sigma\right\|_{\rm op} = \bbu^\intercal \boldsymbol \Sigma \bbu = d\,\mathbb E \left[ \bbu^\intercal \wv_n\wv_n^\intercal \bbu \right]=d\,\mathbb E \left[ (\bbu^\intercal \wv_n)^2 \right].
\end{equation}

Furthermore, we have that

\begin{equation}\label{eq:int2}
\begin{split}
\|\wv_n\|_{\psi_2} &:= \sup_{\bbu' \,:\, \|\bbu'\|_2 = 1} \|(\bbu')^\intercal \wv_n\|_{\psi_2}  \\
&\geq \| \bbu^\intercal \wv_n\|_{\psi_2}\geq \frac{1}{C_2} \sqrt{ \mathbb E \left[ (\bbu^\top \wv_n)^2 \right]},
\end{split}
\end{equation}

where $C_2$ is a numerical constant, and the last inequality comes from Eq. (2.15) of [23]. By combining (<ref>) and (<ref>) with the fact that $\|\wv_n\|_{\psi_2}\le \sigma_W/\sqrt{d}$, we conclude that

\begin{equation}\label{eq:bds}
\left\|\boldsymbol \Sigma\right\|_{\rm op}\le C_2^2\sigma_W^2.
\end{equation}

Hence, the desired result follows by taking $t=d/\sqrt{c_1}$ in (<ref>) and using (<ref>).

At this point, we are ready to prove our main result.

Let the assumptions (A-$\wv$)-(A-$\av$) hold. Then, we have that

\begin{equation}
\begin{split}
%\sup_{i\in [M]}\sup_{\vv_i \in \Scal^{d-1}}
\sup_{\substack{\{\vv_m\}_{m \in [M]} \\ \vv_m \in \Scal^{d-1}}}
\labs [\sv]_i - \f {\mu_a} {2 \pi} \rabs
&\le {\sf err}(N, d, A, \sigma_W, t),
\label{eq:approx sv}
\end{split}
\end{equation}

with probability at least $1 -  4 e^{-t^2/16}-2e^{-d}$, where

\begin{equation}
\begin{split}
&{\sf err}(N, d, A, \sigma_W, t)
:=  C A \bigg(\f {\sigma_W^2} {d} + \f {t(1+\sigma_W)}{\sqrt{N}} \bigg),
\label{eq:approx l}
\end{split}
\end{equation}

and $C$ is a numerical constant. Furthermore, by assuming also (A-$\bv$), we have that, with the same probability,

sup_{_m}_m ∈[M]
_m ∈^d-1
|L(,)-L̃(,)| ≤B M err(N, d, A, σ_W, t),

where $L(\Vsf,\Wsf)$ and $\tilde{L}(\Vsf,\Wsf)$ are defined in (<ref>) and (<ref>), respectively.

By plugging the Taylor expansion of $g$ (cf.  (<ref>)) into the definition $\sv := \frac{1}{N} G_{VW} \av$, we have that, for $i\in [M]$,

\begin{equation}
\begin{split}
[\sv]_i =  \f 1 N \sum_{n \in [N]} \f {a_n} {2 \pi} &+ \f 1 N \sum_{n \in [N]} \f 1 4 a_n \langle \vv_i, \wv_n \rangle \\
\f 1 N \sum_{n \in [N]} \sum_{k\ge 1} a_n c_k\langle \vv_i, \wv_n \rangle^{2k},
\end{split}
\end{equation}

where $c_k = \frac{1}{2 \pi} \frac{((2k-3)!!)^2}{k!}$.

\begin{equation}
\begin{split}
\labs [\sv]_i - \f {\mu_a} {2 \pi} \rabs &\leq   \labs \f 1 N \sum_{n \in [N]} \f {a_n} {2 \pi}  -  \f {\mu_a} {2 \pi} \rabs
+ \labs\f 1 N \sum_{n \in [N]} \f 1 4 a_n \langle \vv_i, \wv_n \rangle \rabs\\
&\hspace{-1em} +
\labs \f 1 N \sum_{n \in [N]} \sum_{k \ge 1} a_n c_k\langle \vv_i, \wv_n \rangle^{2k} \rabs:= T_1 +T_2+T_3.
\end{split}
\end{equation}

Using Hoeffding inequality for bounded variables, we get

\begin{equation}\label{eq:bd1}
\Prob \lsb T_1 \geq \frac{tA}{2\pi\sqrt{N}} \rsb \leq 2 e^{-\f {t^2} {2}}.
\end{equation}

Furthermore, using Lemma <ref>, we have

\begin{equation}\label{eq:bd2}
\Prob \lsb \sup_{\vv_i \in \Scal^{d-1}} T_2 \geq \frac{tA\sigma_W}{4\sqrt{N}} \rsb \leq 2 e^{-\f {t^2} {16}}.
\end{equation}

For the last term, we first note that

\begin{equation}\label{eq:bd3}
T_3 \le    \f 1 N \sum_{n \in [N]} \sum_{k \ge 1} |a_n| c_k\langle \vv_i, \wv_n \rangle^{2} \le \f C N \sum_{n \in [N]} |a_n| \langle \vv_i, \wv_n \rangle^{2},
\end{equation}

where $C$ denotes a numerical constant. Here, in the first inequality we use that $c_k\ge 0$ for all $k\ge 1$ and that $|\langle \vv_i, \wv_n \rangle|\le 1$ since $\vv_i, \wv_n \in \Scal^{d-1}$; and in the second equality we use that $\sum_{k \ge 1} c_k$ is finite, since $\sigma\in L^2(e^{-w^2/2})$. Thus, we can use Lemma <ref> to give a uniform bound[It is uniform in the sense that it upper bounds $\sup_{\vv_i \in \Scal^{d-1}} T_3$.] on the RHS of (<ref>). By combining this with (<ref>) and (<ref>), after some simplifications we obtain (<ref>).

Recall from assumption (A-$\bv$) that $|b_m|\le B$. Then, an immediate consequence of (<ref>) is that

\begin{equation}
%_{\vv_1, \ldots, \vv_M \in \Scal^{d-1}}
\sup_{\substack{\{\vv_m\}_{m \in [M]} \\ \vv_m \in \Scal^{d-1}}}
\labs \bv^{\intercal} \sv - \bv^{\intercal}  \f {\mu_a} {2 \pi}\one\rabs\le M B{\sf err}(N, d, A, \sigma_W, t),
\end{equation}

which readily gives (<ref>).

Theorem <ref> shows that, if we optimize the limit loss (<ref>) instead of the original loss (<ref>), we suffer at most an additive error which scales as $1/\sqrt{N}+1/d$.

The same result can be obtained for any $\sigma\in L^2(e^{-w^2/2})$, as long as the zero-th, first and second Hermite coefficients are non-zero. If at least one of these coefficients is zero, we still expect convergence to the limit loss, but with a different additive error term.

By the same argument of Lemma <ref>, the optimal value of $\bv$ in (<ref>) for a fixed $G_{VV}$ is

^*(G_VV)= Mμ_a2π G^+_VV .

Furthermore, when replacing the solution of (<ref>) into  (<ref>), the loss is given by

(,) _=(<ref>)

= 1 N^2 ^⊺ G_WW -

G^+_VV .

This is equivalent to maximizing

\begin{equation}
\label{eqn:simpleobj}
\tilde \Rsf(\Vsf) =
%\arg \max
%\sup_{\substack{\{\vv_m\}_{m \in [M]} \\ \vv_m \in \Scal^{d-1}}}
%%_{\vv_1, \ldots, \vv_M\in \Scal^{d-1}}
%_{\{\vv_m\}_{m \in [M]}, \  \vv_m \in \Scal^{d-1}}
\one^{\intercal} G_{VV}^{+} \one.
\end{equation}

We refer
to  (<ref>)
as the limit loss and
to $\tilde \Rsf(\Vsf)$ as the limit objective.

Note that, in (<ref>), we have removed the dependency on the target network $\Wsf$, as the limit objective no longer depend on the target network.

We also remark that our bound for the approximation gap in (<ref>) worsens as $M$ increases. This is to be expected since, when $M$ scales e.g. linearly with $N$, the weights of the compressed network should depend on the realization of the target network.

§ EQUIANGULAR TIGHT FRAME (ETF) CONJECTURE

We conjecture that the limit objective (<ref>) is maximized when the weights of $\Vsf$ form an ETF, i.e.,
$\vv_i^\intercal\vv_j =-1/(M-1)$ for all
$i \neq j$.

More precisely, we conjecture that any ETF achieves the maximum in (<ref>):

optimizing  the orientation of the ETF results in a reduction of the compression loss bounded as in (<ref>) – with high probability.

We indicate the value of (<ref>) when $\{\vv_m\}$ is an ETF as $\tilde \Rsf_{\rm ETF}$ and the corresponding compressed network as $\Vsf_{\rm ETF}$. We note that the assumption (A-$\bv$) is satisfied when $\bv$ is given by (<ref>) and the weights of $\Vsf$ form an ETF. We compute the value of $\tilde \Rsf_{\rm ETF}$ in Lemma <ref>, and prove the conjecture in the $M=2$ case in Lemma <ref>.

When  the weights of $\Vsf$ form an ETF, we have that

() = _ETF

= M/g(1) + (M - 1) g-1/M-1 ,

where $g(\alpha)$ is given by (<ref>).

Define $c_1=g(-1/(M-1))$ and $c_2= g(1)-g(-1/(M-1))$. Then, $g(\vv_i^\intercal\vv_j)=c_1$ if $i\neq j$, and $g(\vv_i^\intercal\vv_j)=c_1+c_2$ otherwise. Hence, $G_{VV} = c_1 \one\one^\intercal + c_2 \Iv$ and, by using the Sherman-Morrison formula,

For $M=2$, an ETF is the optimal solution of the optimization in (<ref>).

To see this, let us rewrite (<ref>) as

\begin{multline}
\one^{\intercal}
\lsb \p{
1 & g_{12}\\
g_{12} & 1 }
\rsb^{-1}
\one =
% \one^{\intercal}
% \frac{1}{1 - a^2}
% \begin{pmatrix}
% 1 & -a\\
% -a & 1
% \end{pmatrix}^{+}
%\one =
\frac{4 (1 - g_{12})}{1 - g_{12}^2} = \frac{4}{1 + g_{12}},
\end{multline}

where $g_{12} = 2  g(\langle \vv_1, \vv_2 \rangle)$.

Since $g$ in (<ref>)  is an increasing function in $\al$, the maximum of the limit objective
happens when $\vv_1 = -\vv_2$, which is an ETF.

In the rest of the section, we provide numerical results giving a rather strong indication that the ETF maximizes

First, in Fig. <ref>, we compare (i) the value of the limit objective attained by the ETF (cf. Lemma <ref>) and (ii) the limit objective obtained after running $10^3$ iteration of gradient descent (GD), as a function of $M$.

These two values coincide up to the $12^{\rm th}$ decimal place for $M$ between $5$ and $30$.

height=3.9 cm,
width=9 cm,
x label style=at=(axis description cs:0.5,-0.1),anchor=north,
legend style=at=(0.65,0.25),anchor=west,
ylabel=limit objective,
(top) at (axis cs:1,/pgfplots/ymax);

+[mark=+,draw=mycolor1, line width=1.0pt]
table[col sep=comma]./data/SGDmaximas.txt;

+[mark=otimes,dashed, draw=mycolor2, line width=1.0pt]
table[col sep=comma]./data/simplexresults.txt;

GD solution, ETF solution

Limit objective as estimated through GD and as computed at ETF.

height=4.5 cm,
width=9 cm,

x label style=at=(axis description cs:0.5,-0.1),anchor=north,

xlabel=iteration number,
ylabel= $d \lb \Vsf_{\sf ETF},\Vsf_{\sf GD} \rb $

(top) at (axis cs:1,/pgfplots/ymax);

+[ no marks,draw=mycolor1,line width=1.0pt]
table[col sep=comma]./data/traindyn_frobdists_M=5_iter0.txt;

+[ no marks,draw=mycolor2, line width=1.0pt]
table[col sep=comma]./data/traindyn_frobdists_M=10_iter0.txt;

+[ no marks,draw=mycolor3, line width=1.0pt]
table[col sep=comma]./data/traindyn_frobdists_M=15_iter0.txt;

+[ no marks,draw=mycolor4, line width=1.0pt]
table[col sep=comma]./data/traindyn_frobdists_M=30_iter0.txt;

Distance between the GD solution and the ETF solution, as in (<ref>).

To further investigate the GD solution, let us define the distance between two networks $\Vsf$ and $\Vsf'$ as

d (, ')=

G_VV -G_V'V'

In Fig. <ref>, we plot $d \lb \Vsf_{\sf ETF},\Vsf_{\sf GD} \rb $, i.e., the distance between the ETF compressed network, $\Vsf_{\sf ETF}$, and  the network compressed using gradient descent, $\Vsf_{\sf GD}$, as a function of the iteration number and for various values of $M$.

We note that
the convergence of GD to the ETF solution is rather stable and constant with the number of iterations.

We  also observe that the convergence rate of GD  decreases as $M$ grows large.

height=3.9 cm,
width=9 cm,

legend style=at=(0.75,0.25),anchor=west,
xlabel=iteration number,
ylabel=limit objective
(top) at (axis cs:1,/pgfplots/ymax);

+[name path=A,draw=mycolor1!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/min_m=2.txt;

+[name path=B, draw=mycolor1!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/max_m=2.txt;

[mycolor1!10, forget plot] fill between[of=A and B];

+[draw=mycolor1, no marks,line legend, line width=1.0pt]
table[x expr=, y=q]./data/minmax/avg_m=2.txt;

+[name path=A,draw=mycolor2!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/min_m=4.txt;

+[name path=B, draw=mycolor2!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/max_m=4.txt;

[mycolor2!10, forget plot] fill between[of=A and B];

+[draw=mycolor2, no marks,line legend, line width=1.0pt]
table[x expr=, y=q]./data/minmax/avg_m=4.txt;

+[name path=A,draw=mycolor3!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/min_m=8.txt;

+[name path=B, draw=mycolor3!50,dotted, no marks,forget plot]
table[x expr=, y=q]./data/minmax/max_m=8.txt;

[mycolor3!10, forget plot] fill between[of=A and B];

+[draw=mycolor3, no marks,line legend, line width=1.0pt]
table[x expr=, y=q]./data/minmax/avg_m=8.txt;

Convergence of
GD as a function of the iteration number.
The boundaries of the shaded area indicate the minimum and maximum of the objective across 10 independent initializations.

Finally, in Fig. <ref>, we plot the variability of the convergence behaviour across multiple GD initializations.

§ CONCLUSION

In this paper, the compression of a two-layer neural network with $N$ neurons to a network with $M$ neurons is investigated.

An asymptotic expression for the compression loss is obtained for the case in which the input features are Gaussian and the target weights are i.i.d. sub-Gaussian.

This mean-field view allow to significantly simplify the optimization problem, as the limit objective depends only on $M$ (and not on $N\gg M$).

Numerical simulations suggest that, for a ReLU activation, the optimal structure of the compressed network is the one in which weights lie on an equiangular tight frame,
while the orientation and the scaling of the neurons is obtained from the parameters of the original network.
Proving this claim is left as an open problem. Interesting avenues for future research also consist in understanding the class of activations for which the ETF conjecture holds and extending the current results to deep networks.

[1]
T. Gale, E. Elsen, and S. Hooker, “The state of sparsity in deep neural
networks,” arXiv preprint arXiv:1902.09574, 2019.

[2]
T. Hoefler, D. Alistarh, T. Ben-Nun, N. Dryden, and A. Peste, “Sparsity in
deep learning: Pruning and growth for efficient inference and training in
neural networks,” Journal of Machine Learning Research, vol. 22, no.
241, pp. 1–124, 2021.

[3]
Y. LeCun, J. Denker, and S. Solla, “Optimal brain damage,” in Advances
in Neural Information Processing Systems, vol. 2, 1989.

[4]
B. Hassibi and D. Stork, “Second order derivatives for network pruning:
Optimal brain surgeon,” in Advances in Neural Information Processing
Systems, vol. 5, 1992.

[5]
S. Han, J. Pool, J. Tran, and W. Dally, “Learning both weights and connections
for efficient neural network,” Advances in Neural Information
Processing Systems, vol. 28, 2015.

[6]
A. See, M.-T. Luong, and C. D. Manning, “Compression of neural machine
translation models via pruning,” in Proceedings of The 20th SIGNLL
Conference on Computational Natural Language Learning, 2016, pp. 291–301.

[7]
S. Narang, E. Elsen, G. Diamos, and S. Sengupta, “Exploring sparsity in
recurrent neural networks,” arXiv preprint arXiv:1704.05119, 2017.

[8]
J. Frankle and M. Carbin, “The lottery ticket hypothesis: Finding sparse,
trainable neural networks,” in International Conference on Learning
Representations, 2018.

[9]
S. Han, X. Liu, H. Mao, J. Pu, A. Pedram, M. A. Horowitz, and W. J. Dally,
“Eie: Efficient inference engine on compressed deep neural network,”
ACM SIGARCH Computer Architecture News, vol. 44, no. 3, pp. 243–254,

[10]
M. Ye, C. Gong, L. Nie, D. Zhou, A. Klivans, and Q. Liu, “Good subnetworks
provably exist: Pruning via greedy forward selection,” in
International Conference on Machine Learning, 2020, pp.
10 820–10 830.

[11]
W. Kim, S. Kim, M. Park, and G. Jeon, “Neuron merging: Compensating for pruned
neurons,” Advances in Neural Information Processing Systems, vol. 33,
pp. 585–595, 2020.

[12]
M. Mondelli and A. Montanari, “On the connection between learning two-layer
neural networks and tensor decomposition,” in The 22nd International
Conference on Artificial Intelligence and Statistics, 2019, pp. 1051–1060.

[13]
S. Mei, A. Montanari, and P.-M. Nguyen, “A mean field view of the landscape of
two-layer neural networks,” Proceedings of the National Academy of
Sciences, vol. 115, no. 33, pp. E7665–E7671, 2018.

[14]
G. M. Rotskoff and E. Vanden-Eijnden, “Neural networks as interacting particle
systems: Asymptotic convexity of the loss landscape and universal scaling of
the approximation error,” in Advances in Neural Information Processing
Systems, vol. 32, 2018.

[15]
L. Chizat and F. Bach, “On the global convergence of gradient descent for
over-parameterized models using optimal transport,” in Advances in
Neural Information Processing Systems, vol. 32, 2018.

[16]
J. Sirignano and K. Spiliopoulos, “Mean field analysis of neural networks: A
law of large numbers,” SIAM Journal on Applied Mathematics, vol. 80,
no. 2, pp. 725–752, 2020.

[17]
S. Mei, T. Misiakiewicz, and A. Montanari, “Mean-field theory of two-layers
neural networks: Dimension-free bounds and kernel limit,” in
Conference on Learning Theory, 2019, pp. 2388–2464.

[18]
A. Javanmard, M. Mondelli, and A. Montanari, “Analysis of a two-layer neural
network via displacement convexity,” The Annals of Statistics,
vol. 48, no. 6, pp. 3619–3642, 2020.

[19]
A. Shevchenko and M. Mondelli, “Landscape connectivity and dropout stability
of sgd solutions for over-parameterized neural networks,” in
International Conference on Machine Learning, 2020, pp. 8773–8784.

[20]
A. Shevchenko, V. Kungurtsev, and M. Mondelli, “Mean-field analysis of
piecewise linear solutions for wide ReLU networks,” Journal of
Machine Learning Research, 2022.

[21]
V. Dutordoir, J. Hensman, M. van der Wilk, C. H. Ek, Z. Ghahramani, and
N. Durrande, “Deep neural networks as point estimates for deep Gaussian
processes,” Advances in Neural Information Processing Systems,
vol. 34, 2021.

[22]
R. O'Donnell, Analysis of boolean functions.1em plus 0.5em minus
0.4emCambridge University Press, 2014.

[23]
R. Vershynin, High-dimensional probability: An introduction with
applications in data science.1em plus 0.5em minus 0.4emCambridge university press, 2018.

[24]
C. Jin, P. Netrapalli, R. Ge, S. M. Kakade, and M. I. Jordan, “A short note on
concentration inequalities for random vectors with subgaussian norm,”
arXiv preprint arXiv:1902.03736, 2019.

[25]
R. Vershynin, “Introduction to the non-asymptotic analysis of random
matrices,” arXiv preprint arXiv:1011.3027, 2010.


# Interaction-induced multiparticle bound states in the continuum

Boning Huang 1,2,3 Yongguan Ke 1,2,3,4 Corresponding author. Email:
<EMAIL_ADDRESS>Honghua Zhong 5 Yuri S. Kivshar6 Chaohong Lee 1,2,3,7
Corresponding author. Email<EMAIL_ADDRESS>1Institute of Quantum
Precision Measurement, State Key Laboratory of Radio Frequency Heterogeneous
Integration, Shenzhen University, Shenzhen 518060, China  2College of Physics
and Optoelectronic Engineering, Shenzhen University, Shenzhen 518060, China
3Laboratory of Quantum Engineering and Quantum Metrology, School of Physics
and Astronomy, Sun Yat-Sen University (Zhuhai Campus), Zhuhai 519082, China
4International Quantum Academy, Shenzhen 518048, China  5Institute of
Mathematics and Physics, Central South University of Forestry and Technology,
Changsha, 410004, China  6Nonlinear Physics Center, Australian National
University, Canberra ACT 2601, Australia  7Quantum Science Center of
Guangdong-Hongkong-Macao Greater Bay Area (Guangdong), Shenzhen 518045, China

###### Abstract

Bound states in the continuum (BICs) are localized modes residing in the
radiation continuum. They were first predicted for single-particle states, and
became a general feature of many wave systems. In many-body quantum physics,
it is still unclear what would be a close analog of BICs, and whether
interparticle interaction may induce BICs. Here, we predict a novel type of
multiparticle states in the interaction-modulated Bose-Hubbard model that can
be associated with the BIC concept. Under periodic boundary conditions, a so-
called quasi-BIC appears as a bounded pair residing in a standing wave formed
by the third particle. Under open boundary conditions, such a hybrid state
becomes an eigenstate of the system. We demonstrate that the Thouless pumping
of the quasi-BICs can be realized by modulating the on-site interactions in
space and time. Surprisingly, while the center-of-mass of the quasi-BIC is
shifted by a unit cell in one cycle, the bounded pair moves into the opposite
direction with the standing wave.

Bound states in continuum (BICs) are known as spatially localized states
residing in the continuum spectrum of radiative states [1, 2, 3, 4, 5]. In
practice, an ideal BIC being a dark mode with infinite lifetime always turns
into a quasi-BIC with finite lifetime [6, 7, 8]. The study of BIC and quasi-
BIC in photonic systems has attracted intense interest because of the
fundamental mechanism they provide for many problems of light-matter
interaction [9, 10, 11] and important applications such as lasers [12, 13, 14,
15, 16], high-harmonic generation [17, 18, 19, 20], sensing and imaging [21,
22, 23]. The mechanics of generating BIC vary from symmetry protection [24,
25], separable subspace [26], destructive interference [27, 28], Fabry-Pérot
cavity [29, 27], time-periodic modulation [30, 31, 32, 33], to inverse
construction [34, 35]. Single-particle BICs can be protected by a topological
charge being robust to a variation of the system parameters [36, 37].
Tunability and dynamic control of BICs can be achieved via optical pumping
[38, 39], the use of phase-change materials [40, 41], which are of interest
due to potential applications.

Beyond conventional BICs for classical waves and single-particle quantum
systems, in the presence of impurities, multiparticle quantum BICs may exist
in one-dimensional Hubbard systems [42, 43, 44, 45] and may prevent the system
from thermalization [45]. Even without impurity, two-particle BICs appear in a
Hubbard lattice with anyonic statistics [46] or an intense oscillating
electric field [47]. Up to now, it is still unclear whether particle-particle
interactions may induce BICs and there are no studies of dynamical control of
multiparticle BICs. On the other hand, although it has been theoretically
predicted and experimentally demonstrated the topological Thouless pumping of
conventional bound states [3, 49, 50], it is still challenging and appealing
to achieve the topological Thouless pumping of multiparticle BICs.

Figure 1: Schematics for (a) bound state in continuum under open boundary
condition and (b) quasi-bound state in continuum under periodic boundary
condition.

In this Letter, we predict the existence of few-particle BICs and quasi-BICs
in a Bose-Hubbard model without impurity and propose how to achieve
topological pumping of quasi-BICs by modulating particle-particle interactions
in space and time. Under open boundary condition, multiparticle BIC of three
particles appears as a bound pair localized by the standing wave of the other
particle in a finite lattice; see Fig. 1(a). Under periodic boundary
condition, we develop a new method to construct multiparticle maximally
localized Wannier states (MLWS) via the projected position operator. We find
that quasi-BIC, constructed from three-particle MLWS, appears as a localized
bound pair in the standing wave of the other particle; see Fig. 1(b). In the
pumping process, the transport direction of the standing wave is opposite to
the one of the bounded pair. While the topological pumping of BIC breaks down,
the mean position of quasi-BIC is shifted anticlockwise by one unit cell. The
displacement is directly related to the Chern number of the multiparticle band
which is uniformly populated by the quasi-BIC. Our work paves the way for
systematic generation and topological control of the multiparticle BICs.

Depicted in Fig. 1(a) under open boundary condition and Fig. 1(b) under
periodic boundary condition, we consider interacting bosons in a lattice with
spatial modulation of onsite interaction,

$\hat{H}=J\sum_{j}(\hat{a}_{j+1}^{{\dagger}}\hat{a}_{j}+h.c.)+\frac{1}{2}\sum_{j}U_{j}(\phi)\hat{n}_{j}(\hat{n}_{j}-1).$
(1)

Here, $\hat{a}_{j}^{{\dagger}}$ ($\hat{a}_{j}$) and
$\hat{n}_{j}=\hat{a}_{j}^{{\dagger}}\hat{a}_{j}$ are the operators of bosonic
creation (annihilation) and particle number at the $j$th site, respectively.
The system size is $M$. $J$ is the nearest neighboring hopping strength.
$U_{j}=U_{0}+\delta\cos(2\pi\beta j+\phi)$ is the on-site interaction
strength, where $\delta$, $\beta=p/q$ ($p$, $q$ being co-prime numbers), and
$\phi$ are the modulation strength, spatial frequency and modulation phase,
respectively. $U_{j}$ can be tuned by varying the modulation phase as
$\phi=\omega t$, where $\omega$ is the temporally modulation frequency with
modulation period $T=2\pi/\omega$. In this work, we concentrate on the strong
interaction regime ($|U|>>|J|$), and set $\beta=1/3$ without loss of
generality. Such a model could be readily realized in an array of
superconducting circuits [50]. The onsite interaction can be independently and
dynamically tuned by controlling the anharmonicity of individual qubits. If
the two particles are separately located at different sites, they can tunnel
independently without onsite energy. If two more particles are at the same
site, due to the strong onsite interaction, they form bound states which
tunnel as a whole [51] and feel a periodic onsite energy with period $q=3$,
described by an effective Aubry-André-Harper model [52]. Thus, the
spatiotemperal modulation only breaks time reversal symmetry of bound states
but plays no role for a single particle.

Figure 2: BIC under open boundary condition and quasi-BIC under periodic
boundary condition. (a) Energy spectrum at $\phi=0$ under open boundary
condition, and the inset shows the bound state in continuum given by the
type-(ii) eigenstates with energy $\epsilon=321.9826$(left) and
$\epsilon=320.8269$(right), which are indicated by the orange circle and cross
in the spectrum respectively. (b) Multiparticle Bloch bands at $\phi=0$, and
the inset shows the quasi-bound state in continuum given by the maximally
localized Wannier state of the highest (left) and middle (right) type-(ii)
bands at $t=0$, which are indicated by the orange asterisk and plus sign in
the bands respectively. (c) Decay ratio
$D=\langle\hat{n}_{j}\rangle(t)/\langle\hat{n}_{j}\rangle(0)$ as a function of
time, where the bounded two bosons locating at the $j$-th site. The green,
blue solid lines, green dashed line and blue dotted line correspond to the
four (quasi-)BICs in the insets of Fig.2(a) and (b), respectively. The
parameters are $N=3,M=30,U_{0}=300,\delta=20,J=1$.

First, we first consider three-boson eigenvalues and eigenstates at $\phi=0$
under the open boundary condition and periodic boundary condition; see the
left and middle panels of Fig. 2, respectively. Under the periodic boundary
condition the system has co-translational symmetry, that is, the energy is
conserved when all particles are shifted as a whole by unit cells. Therefore,
the center of mass momentum $\kappa$ is a good quantum number, and the energy
spectrum as a function of $\kappa$ forms multiparticle energy bands. Due to
the strong repulsive interaction ($U/J=300$), the energy spectrum [Fig. 2 (a)]
and the multiparticle energy band [Fig. 2(b)] are separated into three
clusters. The other parameters are chosen as $J=1$, $\delta=20$, and $M=30$.
From the bottom to the top, the three clusters (marked by red, blue, and
green) correspond to (i) three independent bosons, (ii) two bound bosons and
an independent bosons, and (iii) three bound bosons. Since the type-(i,iii)
states are well known [52], we focus on the type-(ii) states. We notice that
in the conventional Bose-Hubbard model both the bound pair and independent
bosons in the type-(ii) states are always extended over the lattice.
Surprisingly, in our systems under open boundary we find all type-(ii) states
are a localized bound pair by a standing wave; see insets of Fig. 2(a) for two
typical eigenstates with energies $\epsilon=321.9826$ (left) and
$\epsilon=320.8269$ (right) and their higher-order correlation functions in
the Supplemental Materials [52]. Considering the cluster of type-(ii) bound
states forms the continuum spectrum as the system size increases, we term the
localized type-(ii) states under the open boundary condition as multiparticle
bound state in continuum. Since three-paticle BICs can be decomposed to an
extended single-particle state and a localized two-particle state, we can
obtain an effective Hamiltonian for the localized part of some BICs [52]

$\displaystyle\hat{H}_{\rm
eff}=\hat{H}+\sum_{j}2U_{j}(\phi){|\varphi_{j}|}^{2}\hat{n}_{j},$ (2)

where $\varphi_{j}$ is the amplitude of the extended single-particle state,
which provides background potential. The idea of background potential provided
by the particle itself can explain a variety of localization phenomena [53,
54, 55, 56, 57, 58, 59, 60, 61, 62].

There is no multiparticle BIC under periodic boundary condition because the
systems are translation invariant. However, we find that the multiparticle
energy bands for the type-(ii) states are flat because of the strong
interaction. The flat energy band means that the particles are almost
localized as the group velocities of the multiparticle Bloch states
$|\psi_{m}(\kappa)\rangle$ ($m\in$ the type-(ii) states) almost vanish. We can
construct multiparticle maximally localized Wannier states (MLWS) of the flat
band which are approximately the eigenstates. To this end, we first make use
of the projected position operator $\hat{P}\hat{x}\hat{P}$,

$\hat{P}=\sum_{m\in\mathcal{M},\kappa}|\psi_{m}(\kappa)\rangle\langle\psi_{m}(\kappa)|,$
(3)

where position operator is given by $\hat{x}=N^{-1}\sum_{j}j\hat{n}_{j}$, and
$m$ takes targeted $\mathcal{M}$ multiparticle Bloch bands. Similarly to
single-particle MLWS [63], we choose the multiparticle Wannier state
$|W_{m}(0)\rangle$ as the eigenstate of the projected position operator with
eigenvalue $X_{0m}$, that is,

$\displaystyle\langle W_{n}(R)|\hat{x}|W_{m}(0)\rangle$ $\displaystyle=\langle
W_{n}(R)|\hat{P}\hat{x}\hat{P}|W_{m}(0)\rangle$ (4)
$\displaystyle=X_{0m}\delta_{R,0}\delta_{m,n},$

where $m,n\in\mathcal{M}$. We can prove that $|W_{m}(0)\rangle$ is the
multiparticle MLWS [52]. As the system size increases, the energy gaps between
these states vanish and the spectrum becomes continuum. We term the
multiparticle MLWS of the type-(ii) band as multiparticle quasi-bound state in
continuum (quasi-BIC). Insets of Fig. 2(b) show the density distribution
$\langle n_{j}\rangle$ of two quasi-BICs in the highest (left) and middle
(right) of the (ii) type cluster bands; see higher order correlation functions
in the Supplemental Materials [52]. These multiparticle quasi-BICs are a bound
pair in the standing wave of the other particle.

Theoretically, BICs are eigenstates which have infinite lifetimes, and quasi-
BICs are approximately eigenstates which have long lifetimes. By setting the
initial states as the BIC and quasi-BIC in insets of Figs. 2(a) and (b), the
difference between BIC and quasi-BIC can be found by tracing the ratio
$D=n_{j}(t)/n_{j}(0)$ in a long time evolution under the static Hamiltonian
with $\phi=0$, where the bound pair is located at the $j$-th site; see Fig.
2(c). We find $D=1$ preserves in all time for BIC and $D\approx 1$ even at the
time scale of $Jt\approx 10^{4}$ for the right quasi-BIC in Fig. 2(b). The
left quasi-BIC has a much smaller decay rate than the right quasi-BIC, because
it comes from the energy band which is more flat, and the BICs have no decay.

Figure 3: Pumping dynamics of BIC (left) and quasi-BIC (right). (a,d)
Schematics for the dynamical localization of BIC and topological pumping of
quasi-BIC. Colors denote the density distribution. (b,e) Density distribution
as a function of time for (a,d). The white dashed line denotes the c.m.
position of the bound pair. (c,f) Density distribution of three instantaneous
states around transition time $t=1/3T$ for BIC and $t=1/6T$ and quasi-BIC. (d)
Schematics for the adiabatic evolution of BIC under open boundary condition in
one cycle in a ring lattice with $12$ sites. Colors denote the density
distribution of the initial quasi-BIC (left) and the final quasi-BIC (right).
(e) Density distribution as a function of time for (d). The white dashed line
denotes the c.m. position of the bound pair. (f) Density distribution of two
instantaneous states around the transition moment in (e). The parameters are
set as $U_{0}=300,\delta=20,J=1,\omega=0.0005$.

Based on multiparticle Thouless pumping theory [3], if an initial state
$|\psi(0)\rangle$ uniformly fills and sweeps a multiparticle energy band with
nontrivial Chern number,

$C_{m}=\frac{i}{2\pi}\int d\kappa\int
dt(\langle\partial_{t}\psi_{m}|\partial_{\kappa}\psi_{m}\rangle-\langle\partial_{\kappa}\psi_{m}|\partial_{t}\psi_{m}\rangle),$
(5)

the position shift of the multiparticle Wannier state in one cycle
($T=2\pi/\omega$) is equal to the Chern number multiples of unit cell, $\Delta
x(T)=\langle\psi(T)|\hat{x}|\psi(T)\rangle-\langle\psi(0)|\hat{x}|\psi(0)\rangle=C_{m}q$.
The Chern numbers of the highest and lowest type-(ii) band are $-1$, and there
are finite energy gaps that separate other energy bands, which make it
possible to perform topological pumping.

We consider the particle transport for BIC and quasi-BIC under a
spatiotemperal modulation of interaction. Under open boundary condition, the
initial state is set as the multiparticle BIC given by the type-(ii)
eigenstate with energy $\epsilon=321.8212$, and it evolves with time as
$|\psi(t)\rangle=\mathcal{T}\exp[-i\int_{0}^{t}\hat{H}(\tau)d\tau]|\psi(0)\rangle$,
where $\mathcal{T}$ is time ordering operator. The numerical studies of the
dynamics require a lot of computational resources for large $M$ and $N$ and
discrete time sequences. However, strong on-site interaction supports large
energy gaps between the types (ii) and (i,iii) states. Because of negligible
coupling between states of different types under slow modulation, we can only
consider the Hamiltonian in the subspace
$\mathcal{V}=\\{|0,....,n_{i}=2,...,n_{j}=1,...\rangle\\}$, that reduces the
dimension of the Hilbert space from $(M+2)(M+1)M/6$ to $M(M-1)$. The method of
subspace projection is beneficial to numerical simulations of many-body
systems [52]. By choosing parameters as
$U_{0}=300,\delta=20,J=1,\omega=0.0005$ and $M=12$, we calculate the density
distribution
$\langle\hat{n}_{j}(t)\rangle=\langle\psi(t)|\hat{n}_{j}|\psi(t)\rangle$; see
Fig. 3(a) and 3(b) for the schematic and time evolution, respectively. At the
$6$th site, the bound pair tunnels to the $7$th site and returns to the $6$th
site in one cycle; see a dashed white line for the mean position of the bound
pair. Transitions occur mainly at $t=1/3T$ and $t=5/6T$. Clearly, BIC is
dynamically localized at multiples of the pumping cycle, and the topological
pumping of the multiparticle BIC breaks down because there is no uniform band
occupation. To show the motion of the other particle, we calculate the density
distribution before, at, and after the transition time $t=1/3T$; see left,
middle, and right panels in Fig. 3(c), respectively. We find that the standing
wave gradually gathers around the right-hand side of the bound pair before the
transition, it mixes with the bound pair at the transition, and then it
tunnels to the left-hand side, after the transition. The transport direction
of the standing wave is opposite to that of the bound pair.

However, employing periodic boundary condition, we can realize topological
pumping of a quasi-BIC filling uniformly the highest band of the type (ii),
see Fig. 3(d) for schematic and 3(e) for the time evolution of the density
distribution. The initial state is composed of a two-particle bound state
localized at the $6$th site and a standing wave of the other particle centered
at the $12$th site. Under slow modulation, the initial state sweeps the
highest band of the type (ii), and the global Berry curvatures play an
important role. Unidirectional tunnelings of the bound pair occur mainly at
$t=T/6,~{}T/2,~{}5T/6$. We find that in one cycle the mean position of the
bound pair is shifted from the $6$th to the $3$th sites; see a white dashed
line in Fig. 3(e). The left and right panels of Fig. 3(d) show almost the same
profiles of the initial and final states, except that the relative distance
between the bound pair and the standing wave is shifted by $0.0103$ unit cell.
The mean position is anticlockwise being shifted by $-0.9896$ unit cell, that
is consistent with the Chern number $C=-1$ of the filled band. To understand
the corresponding microscopic picture, we show the density distribution around
the transition time at $t=1/6T$; see Fig. 3(f). We find that the standing wave
first gathers around the left side of the bound pair and then mixes with the
two particles at $t=T/6$, and crosses to the right side of the bound pair. The
mixed state at $t=T/6$ is presented approximately by
$(|1_{5},2_{6}\rangle+|2_{5},1_{6}\rangle)/\sqrt{2}$, which can be prepared by
local excitation. This means that the standing wave moves clockwise in the
pumping process, which is opposite to the anticlockwise direction of the
bound-pair transport. This mechanics differs substantially the topological
pumping of bound states and topological resonant tunnelings, where all
particles move in the same direction as a whole or one by one [3, 49, 50]. We
notice that although the energy gaps between different quasi-BICs are small,
the topological pumping of quasi-BICs is relatively robust to disorder in the
onsite energies [52].

In summary, we have suggested a systematic method for generating BIC and
quasi-BICs, that can be extended to a variety of many-body systems, involving
fermionic particles, spin chain, and hybrid quantum systems. We notice that,
as the next step, it is interesting to construct maximally localized
multiparticle Wannier states in two dimensions, which can support richer many-
body quasi-BICs. While the Thouless pumping provides a simple way to control
the spatial freedom of quasi-BICs in a topologically quantized way, some
further studies are required to explore the possibility of fractional
topological pumping of many-body quasi-BICs. Many-body BICs have been proven
to break ergodicity in many-body systems [45]. In our case, strong interaction
is required to stabilize quasi-BICs. However, when the interaction strength
decreases, a small group of quasi-BICs of the same band at different sites
will couple with each other. Our system may provide a new possibility to
realize the Hilbert space fragmentation, many-body scar, and ergodicity
breaking.

###### Acknowledgements.

The authors acknowledge useful discussions with Alexander Poddubny, Wenjie
Liu, Ling Lin, Li Zhang, and Na Zhang. This work was supported by the National
Key Research and Development Program of China (Grant No. 2022YFA1404104), the
National Natural Science Foundation of China (Grants Nos. 12025509, 12275365,
and 12175315), the Key-Area Research and Development Program of Guangdong
Province (Grant No. 2019B030330001), the Natural Science Foundation of
Guangdong (Grant No. 2023A1515012099), and the Australian Research Council
(Grant No. DP200101168).

## References

  * Neumann and Wigner [1929] J. v. Neumann and E. Wigner, Über merkwürdige diskrete eigenwerte, Phys. Z 30, 291 (1929).
  * Hsu  _et al._ [2016] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Bound states in the continuum, Nat. Rev. Mater. 1, 1 (2016).
  * Sadrieva  _et al._ [2019] Z. Sadrieva, K. Frizyuk, M. Petrov, Y. Kivshar, and A. Bogdanov, Multipolar origin of bound states in the continuum, Phys. Rev. B 100, 115303 (2019).
  * Koshelev  _et al._ [2020a] K. Koshelev, A. Bogdanov, and Y. Kivshar, Engineering with bound states in the continuum, Opt. Photon. News 31, 38 (2020a).
  * Azzam and Kildishev [2021] S. I. Azzam and A. V. Kildishev, Photonic bound states in the continuum: from basics to applications, Adv. Opt. Mater. 9, 2001469 (2021).
  * Sadrieva  _et al._ [2017] Z. F. Sadrieva, I. S. Sinev, K. L. Koshelev, A. Samusev, I. V. Iorsh, O. Takayama, R. Malureanu, A. A. Bogdanov, and A. V. Lavrinenko, Transition from optical bound states in the continuum to leaky resonances: role of substrate and roughness, ACS Photonics 4, 723 (2017).
  * Doeleman  _et al._ [2018] H. M. Doeleman, F. Monticone, W. den Hollander, A. Alù, and A. F. Koenderink, Experimental observation of a polarization vortex at an optical bound state in the continuum, Nat. Photonics 12, 397 (2018).
  * Koshelev  _et al._ [2018a] K. Koshelev, S. Lepeshov, M. Liu, A. Bogdanov, and Y. Kivshar, Asymmetric metasurfaces with high-$q$ resonances governed by bound states in the continuum, Phys. Rev. Lett. 121, 193903 (2018a).
  * Koshelev  _et al._ [2018b] K. L. Koshelev, S. K. Sychev, Z. F. Sadrieva, A. A. Bogdanov, and I. V. Iorsh, Strong coupling between excitons in transition metal dichalcogenides and optical bound states in the continuum, Phys. Rev. B 98, 161113 (2018b).
  * Cao  _et al._ [2020] S. Cao, H. Dong, J. He, E. Forsberg, Y. Jin, and S. He, Normal-incidence-excited strong coupling between excitons and symmetry-protected quasi-bound states in the continuum in silicon nitride–ws2 heterostructures at room temperature, J. Phys. Chem. Lett. 11, 4631 (2020).
  * Kravtsov  _et al._ [2020] V. Kravtsov, E. Khestanova, F. A. Benimetskiy, T. Ivanova, A. K. Samusev, I. S. Sinev, D. Pidgayko, A. M. Mozharov, I. S. Mukhin, M. S. Lozhkin, _et al._ , Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum, Light Sci. Appl. 9, 56 (2020).
  * Meier  _et al._ [1999] M. Meier, A. Mekis, A. Dodabalapur, A. Timko, R. Slusher, J. Joannopoulos, and O. Nalamasu, Laser action from two-dimensional distributed feedback in photonic crystals, Appl. Phys. Lett. 74, 7 (1999).
  * Imada  _et al._ [1999] M. Imada, S. Noda, A. Chutinan, T. Tokuda, M. Murata, and G. Sasaki, Coherent two-dimensional lasing action in surface-emitting laser with triangular-lattice photonic crystal structure, Appl. Phys. Lett. 75, 316 (1999).
  * Kodigala  _et al._ [2017] A. Kodigala, T. Lepetit, Q. Gu, B. Bahari, Y. Fainman, and B. Kanté, Lasing action from photonic bound states in continuum, Nature 541, 196 (2017).
  * Rybin and Kivshar [2017] M. Rybin and Y. Kivshar, Supercavity lasing, Nature 541, 164 (2017).
  * Huang  _et al._ [2020] C. Huang, C. Zhang, S. Xiao, Y. Wang, Y. Fan, Y. Liu, N. Zhang, G. Qu, H. Ji, J. Han, _et al._ , Ultrafast control of vortex microlasers, Science 367, 1018 (2020).
  * Carletti  _et al._ [2018] L. Carletti, K. Koshelev, C. De Angelis, and Y. Kivshar, Giant nonlinear response at the nanoscale driven by bound states in the continuum, Phys. Rev. Lett. 121, 033903 (2018).
  * Koshelev  _et al._ [2019] K. Koshelev, Y. Tang, K. Li, D.-Y. Choi, G. Li, and Y. Kivshar, Nonlinear metasurfaces governed by bound states in the continuum, ACS Photonics 6, 1639 (2019).
  * Koshelev  _et al._ [2020b] K. Koshelev, S. Kruk, E. Melik-Gaykazyan, J.-H. Choi, A. Bogdanov, H.-G. Park, and Y. Kivshar, Subwavelength dielectric resonators for nonlinear nanophotonics, Science 367, 288 (2020b).
  * Zograf  _et al._ [2020] G. Zograf, A. Zalogina, K. Koshelev, D.-Y. Choi, V. Korolev, R. Hollinger, D. Kartashov, M. Zürch, C. Spielmann, S. Makarov, B. Luther-Davies, S. Kruk, and Y. Kivshar, High-harmonic generation in dielectric metasurfaces empowered by bound states in the continuum, in  _CLEO_ (Optica Publishing Group, 2020) p. FTh1C.5.
  * Liu  _et al._ [2017] Y. Liu, W. Zhou, and Y. Sun, Optical refractive index sensing based on high-q bound states in the continuum in free-space coupled photonic crystal slabs, Sensors 17, 1861 (2017).
  * Leitis  _et al._ [2019] A. Leitis, A. Tittl, M. Liu, B. H. Lee, M. B. Gu, Y. S. Kivshar, and H. Altug, Angle-multiplexed all-dielectric metasurfaces for broadband molecular fingerprint retrieval, Sci. Adv. 5, eaaw2871 (2019).
  * Yesilkoy  _et al._ [2019] F. Yesilkoy, E. R. Arvelo, Y. Jahani, M. Liu, A. Tittl, V. Cevher, Y. Kivshar, and H. Altug, Ultrasensitive hyperspectral imaging and biodetection enabled by dielectric metasurfaces, Nat. Photonics 13, 390 (2019).
  * Moiseyev [2009] N. Moiseyev, Suppression of feshbach resonance widths in two-dimensional waveguides and quantum dots: A lower bound for the number of bound states in the continuum, Phys. Rev. Lett. 102, 167404 (2009).
  * Plotnik  _et al._ [2011] Y. Plotnik, O. Peleg, F. Dreisow, M. Heinrich, S. Nolte, A. Szameit, and M. Segev, Experimental observation of optical bound states in the continuum, Phys. Rev. Lett. 107, 183901 (2011).
  * Rivera  _et al._ [2016] N. Rivera, C. W. Hsu, B. Zhen, H. Buljan, J. D. Joannopoulos, and M. Soljačić, Controlling directionality and dimensionality of radiation by perturbing separable bound states in the continuum, Sci. Rep. 6, 33394 (2016).
  * Bulgakov and Sadreev [2008] E. N. Bulgakov and A. F. Sadreev, Bound states in the continuum in photonic waveguides inspired by defects, Phys. Rev. B 78, 075105 (2008).
  * Gao  _et al._ [2016] X. Gao, C. W. Hsu, B. Zhen, X. Lin, J. D. Joannopoulos, M. Soljačić, and H. Chen, Formation mechanism of guided resonances and bound states in the continuum in photonic crystal slabs, Sci. Rep. 6, 31908 (2016).
  * Marinica  _et al._ [2008] D. C. Marinica, A. G. Borisov, and S. V. Shabanov, Bound states in the continuum in photonics, Phys. Rev. Lett. 100, 183902 (2008).
  * González-Santander  _et al._ [2013] C. González-Santander, P. Orellana, and F. Dominguez-Adame, Bound states in the continuum driven by ac fields, EPL 102, 17012 (2013).
  * Longhi and Valle [2013] S. Longhi and G. D. Valle, Floquet bound states in the continuum, Sci. Rep. 3, 2219 (2013).
  * Agarwala and Sen [2017] A. Agarwala and D. Sen, Effects of local periodic driving on transport and generation of bound states, Phys. Rev. B 96, 104309 (2017).
  * Zhu  _et al._ [2020] B. Zhu, Y. Ke, W. Liu, Z. Zhou, and H. Zhong, Floquet-surface bound states in the continuum in a resonantly driven one-dimensional tilted defect-free lattice, Phys. Rev. A 102, 023303 (2020).
  * Molina  _et al._ [2012] M. I. Molina, A. E. Miroshnichenko, and Y. S. Kivshar, Surface bound states in the continuum, Phys. Rev. Lett. 108, 070401 (2012).
  * Gallo and Molina [2014] N. Gallo and M. Molina, Bulk and surface bound states in the continuum, J. Phys. A Math. Theor. 48, 045302 (2014).
  * Yang  _et al._ [2013] B.-J. Yang, M. Saeed Bahramy, and N. Nagaosa, Topological protection of bound states against the hybridization, Nat. Commun. 4, 1524 (2013).
  * Zhen  _et al._ [2014] B. Zhen, C. W. Hsu, L. Lu, A. D. Stone, and M. Soljačić, Topological nature of optical bound states in the continuum, Phys. Rev. Lett. 113, 257401 (2014).
  * Fan  _et al._ [2019] K. Fan, I. V. Shadrivov, and W. J. Padilla, Dynamic bound states in the continuum, Optica 6, 169 (2019).
  * Han  _et al._ [2019] S. Han, L. Cong, Y. K. Srivastava, B. Qiang, M. V. Rybin, A. Kumar, R. Jain, W. X. Lim, V. G. Achanta, S. S. Prabhu, _et al._ , All-dielectric active terahertz photonics driven by bound states in the continuum, Adv. Mater. 31, 1901921 (2019).
  * Mikheeva  _et al._ [2019] E. Mikheeva, K. Koshelev, D.-Y. Choi, S. Kruk, J. Lumeau, R. Abdeddaim, I. Voznyuk, S. Enoch, and Y. Kivshar, Photosensitive chalcogenide metasurfaces supporting bound states in the continuum, Opt. Express 27, 33847 (2019).
  * Chen and Fan [2020] X. Chen and W. Fan, Tunable bound states in the continuum in all-dielectric terahertz metasurfaces, Nanomaterials 10, 623 (2020).
  * Zhang  _et al._ [2012] J. M. Zhang, D. Braak, and M. Kollar, Bound states in the continuum realized in the one-dimensional two-particle hubbard model with an impurity, Phys. Rev. Lett. 109, 116405 (2012).
  * Zhang  _et al._ [2013] J. M. Zhang, D. Braak, and M. Kollar, Bound states in the one-dimensional two-particle hubbard model with an impurity, Phys. Rev. A 87, 023613 (2013).
  * Zhong  _et al._ [2017] H.-H. Zhong, Z. Zhou, B. Zhu, Y.-G. Ke, and C.-H. Lee, Floquet bound states in a driven two-particle Bose-Hubbard model with an impurity, CPL 34, 070304 (2017).
  * Sugimoto  _et al._ [2023] S. Sugimoto, Y. Ashida, and M. Ueda, Many-body bound states in the continuum (2023), arXiv:2307.05456 [quant-ph] .
  * Zhang  _et al._ [2023a] W. Zhang, L. Qian, H. Sun, and X. Zhang, Anyonic bound states in the continuum, Commun. Phys. 6, 139 (2023a).
  * Della Valle and Longhi [2014] G. Della Valle and S. Longhi, Floquet-hubbard bound states in the continuum, Phys. Rev. B 89, 115118 (2014).
  * Ke  _et al._ [2017] Y. Ke, X. Qin, Y. S. Kivshar, and C. Lee, Multiparticle wannier states and thouless pumping of interacting bosons, Phys. Rev. A 95, 063630 (2017).
  * Walter  _et al._ [2023] A.-S. Walter, Z. Zhu, M. Gächter, J. Minguzzi, S. Roschinski, K. Sandholzer, K. Viebahn, and T. Esslinger, Quantization and its breakdown in a $h$ubbard–thouless pump, Nat. Physics , 1 (2023).
  * Tao  _et al._ [2023] Z. Tao, W. Huang, J. Niu, L. Zhang, Y. Ke, X. Gu, L. Lin, J. Qiu, X. Sun, X. Yang, J. Zhang, J. Zhang, S. Zhao, Y. Zhou, X. Deng, C. Hu, L. Hu, J. Li, Y. Liu, D. Tan, Y. Xu, T. Yan, Y. Chen, C. Lee, Y. Zhong, S. Liu, and D. Yu, Interaction-induced topological pumping in a solid-state quantum system (2023), arXiv:2303.04582 [quant-ph] .
  * Winkler  _et al._ [2006] K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. Hecker Denschlag, A. Daley, A. Kantian, H. Büchler, and P. Zoller, Repulsively bound atom pairs in an optical lattice, Nature 441, 853 (2006).
  * [52] See Supplemental Material for details of (S1) Derivation of effective Hamiltonian for bound states; (S2) Dynamics of the type-(i,iii) states for three particles; (S3) higher order correlation of the BIC and quasi-BIC; (S4) Demonstration of multiparticle MLWS as eigenstate of projected position operator; (S5) Derivation of effective Hamiltonian for the BIC; (S6) Topological pumping of quasi-BIC in subspace; (S7) Robustness against disorder for topological pumping of the quasi-BIC, which includes Refs.[48,62,64,65].
  * Kagan and Maksimov [1984] Y. Kagan and L. Maksimov, Localization in a system of interacting particles diffusing in a regular crystal, Zh. Eksp. Teor. Fiz. 87, 348 (1984).
  * Grover and Fisher [2014] T. Grover and M. P. A. Fisher, Quantum disentangled liquids, J. Stat. Mech. Theory Exp. 2014, P10010 (2014).
  * De Roeck and Huveneers [2014] W. De Roeck and F. m. c. Huveneers, Scenario for delocalization in translation-invariant systems, Phys. Rev. B 90, 165137 (2014).
  * Yao  _et al._ [2016] N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin, and J. E. Moore, Quasi-many-body localization in translation-invariant systems, Phys. Rev. Lett. 117, 240601 (2016).
  * Zhong  _et al._ [2020] J. Zhong, N. A. Olekhno, Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Photon-mediated localization in two-level qubit arrays, Phys. Rev. Lett. 124, 093604 (2020).
  * Zhong and Poddubny [2021] J. Zhong and A. N. Poddubny, Classification of three-photon states in waveguide quantum electrodynamics, Phys. Rev. A 103, 023720 (2021).
  * Lee [2021] C. H. Lee, Many-body topological and skin states without open boundaries, Phys. Rev. B 104, 195102 (2021).
  * Kirsch  _et al._ [2021] M. S. Kirsch, Y. Zhang, M. Kremer, L. J. Maczewsky, S. K. Ivanov, Y. V. Kartashov, L. Torner, D. Bauer, A. Szameit, and M. Heinrich, Nonlinear second-order photonic topological insulators, Nat. Physics 17, 995 (2021).
  * Poddubny [2023] A. N. Poddubny, Interaction-induced analog of a non-hermitian skin effect in a lattice two-body problem, Phys. Rev. B 107, 045131 (2023).
  * Zhang  _et al._ [2023b] N. Zhang, Y. Ke, L. Lin, L. Zhang, and C. Lee, Stable interaction-induced anderson-like localization embedded in standing waves, New J. Phys. 25, 043021 (2023b).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
  * Takahashi [1977] M. Takahashi, Half-filled hubbard model at low temperature, J. Phys. C: Solid State Phys. 10, 1289 (1977).
  * Qin  _et al._ [2014] X. Qin, Y. Ke, X. Guan, Z. Li, N. Andrei, and C. Lee, Statistics-dependent quantum co-walking of two particles in one-dimensional lattices with nearest-neighbor interactions, Phys. Rev. A 90, 062301 (2014).

Supplemental Materials:

Topological inverse band theory in waveguide quantum electrodynamics

Boning Huang1,2,3, Yongguan Ke1,2,3,4,∗, Honghua Zhong5, Yuri S. Kivshar6,
Chaohong Lee1,2,3,7†

1 _Institute of Quantum Precision Measurement, State Key Laboratory of Radio
Frequency Heterogeneous Integration, Shenzhen University, Shenzhen 518060,
China_
2 _College of Physics and Optoelectronic Engineering, Shenzhen University,
Shenzhen 518060, China_
3 _Laboratory of Quantum Engineering and Quantum Metrology, School of Physics
and Astronomy, Sun Yat-Sen University (Zhuhai Campus), Zhuhai 519082, China_
4 _International Quantum Academy, Shenzhen 518048, China_
5 _Institute of Mathematics and Physics, Central South University of Forestry
and Technology, Changsha, 410004, China_
6 _Nonlinear Physics Center, Australian National University, Canberra ACT
2601, Australia_
7 _Quantum Science Center of Guangdong-Hongkong-Macao Greater Bay Area
(Guangdong), Shenzhen 518045, China_

##  S1 Derivation of effective Hamiltonian for bound states

In this section, we show how to derive effective Hamiltonian for three-
particle bound states through perturbation theory [1, 2]. When
$U_{0}\gg(\delta,J)$, the term noted as

$\displaystyle\hat{H}^{\prime}=$ $\displaystyle
J\sum_{j}(\hat{a}^{{\dagger}}_{j+1}\hat{a}_{j}+h.c.)+\frac{1}{2}\sum_{j}\delta\cos(2\pi\beta
j+{\varphi})n_{j}(n_{j}-1)$ (S1)

can be treated as a perturbation to the term noted as

$\hat{H}_{0}=\frac{1}{2}\sum_{j}U_{0}n_{j}(n_{j}-1).$ (S2)

When there are three bosons, the Hilbert space is divided into the subspace
$u$ expanded with degenerate eigenstates $|3\rangle_{j}$ of $\hat{H}_{0}$ with
eigenvalues $E_{j}=3U_{0}$ and its complement subspace $\nu$ expanded by
degenerate eigenstates $|2\rangle_{j}|1\rangle_{k}(j\neq k)$ with eigenvalues
$E_{j,k}=U_{0}$ and degenerate eigenstates
$|1\rangle_{j}|1\rangle_{k}|1\rangle_{l}(j\neq k\neq l)$ with eigenvalues
$E_{j,k,l}=0$. Here, $|3\rangle_{j}$, $|2\rangle_{j}|1\rangle_{k}$, and
$|1\rangle_{j}|1\rangle_{k}|1\rangle_{l}$ are short for Fock states
$|0,...,n_{j}=3,...,0\rangle$, $|0,...,n_{j}=2,...,n_{k}=1,...,0\rangle$, and
$|0,...,n_{j}=1,...,n_{k}=1,...,n_{l}=1,...,0\rangle$, respectively. The
projection operators upon $u,\nu$ are respectively defined as

$\displaystyle\hat{P}=$ $\displaystyle\sum_{j}|3\rangle_{j}\langle 3|_{j},$
(S3) $\displaystyle\hat{S}=\sum_{j\neq
k}\frac{1}{E_{j}-E_{j,k}}|2\rangle_{j}|1\rangle_{k}\langle 1|_{k}\langle
2|_{j}$ $\displaystyle+\sum_{j\neq k,j\neq l,k\neq
l}\frac{1}{E_{j}-E_{j,k,l}}|1\rangle_{j}|1\rangle_{k}|1\rangle_{l}\langle
1|_{l}\langle 1|_{k}\langle 1|_{j}$

Applying the degenerate perturbation theory up to the third order, the
effective Hamiltonian of subspace $u$ is given by

$\hat{H}_{\rm
eff}=\hat{P}\hat{H}\hat{P}+\hat{P}\hat{H}^{\prime}\hat{S}\hat{H}^{\prime}\hat{P}+\hat{P}\hat{H}^{\prime}\hat{S}\hat{H}^{\prime}\hat{S}\hat{H}^{\prime}\hat{P}-\hat{P}\hat{H}^{\prime}\hat{S}^{2}\hat{H}^{\prime}\hat{P}\hat{H}^{\prime}\hat{P}.$
(S4)

Substituting Eq. S3 into Eq. S4, the effective Hamiltonian is given by

$\hat{H}_{\rm
eff}=\sum_{j}3\left[U_{0}+(1-\frac{J^{2}}{U_{0}^{2}})\delta\cos(2\pi\beta
j+{\varphi})\right]\hat{c}_{j}^{{\dagger}}\hat{c}_{j}+\frac{3J^{3}}{2U_{0}^{2}}\sum_{j}(\hat{c}_{j+1}^{{\dagger}}\hat{c}_{j}+h.c.).$
(S5)

Here,
$\hat{c}_{j}^{{\dagger}}=\frac{1}{\sqrt{6}}\hat{a}_{j}^{{\dagger}}\hat{a}_{j}^{{\dagger}}\hat{a}_{j}^{{\dagger}}$
is the creation operator of three bosons as a whole at the $j$th site, and a
uniform on-site energy shift $3J^{2}/U_{0}$ is neglected. Similarly, the
effective model can be extended to the case of $N$-particle bound state with
the general form

$\hat{H}_{\rm eff}=\sum_{j}U_{\rm
eff}\hat{d}_{j}^{{\dagger}}\hat{d}_{j}+J_{\rm
eff}\sum_{j}(\hat{d}_{j+1}^{{\dagger}}\hat{d}_{j}+h.c.),$ (S6)

with
$\hat{d}_{j}^{{\dagger}}=\frac{1}{\sqrt{N!}}(\hat{a}_{j}^{{\dagger}})^{N}$.
The effective onsite energy $U_{\rm eff}$ and hopping strength $J_{\rm eff}$
can be obtained in a similar way. In the general case, the effective
Hamiltonian is the well-known Aubry-André-Harper model of a quasiparticle as a
$N$-particle bound state.

##  S2 Dynamics of the type-(i,iii) states for three particles

Figure S1: Density distribution as a function of time for a adiabatic
evolution of (a) bound state and (b) scattering state and (c) evolution of
same initial state under static Hamiltonian. The red dashed line in (a)
denotes mean position of the bounded three bosons. The parameters are
$M=12,U_{0}=30,J=3,\delta=2$, and $\omega=0.0001$.

The type-(i,iii) states are the three-particle bound state and scattering
state. With the initial state given by the multiparticle MLWS of the highest
multiparticle Bloch band with Chern number $-1$, Fig. S1(a) shows the density
distribution $\langle\hat{n}_{j}(t)\rangle$ as a function of time in one cycle
of time-dependent adiabatic evolution, where the red dashed line shows its
mean position as a function of time. The three-particle bound state is
unidirectionally shifted by $2.9759$ sites during a pumping cycle, consistent
with the corresponding Chern number. For the three-particle scattering states,
The behaviors are quite different whether modulation is turned on or off. Fig.
S1(b) corresponds to the initial state given by the eigenstate of projected
position operator of the highest scattering-state band under periodic boundary
condition, and it is noticed that the state is extended. Fig. S1(c) shows a
time-independent evolution for the same initial state in (b). The dynamics is
governed by
$|\psi(t)\rangle=\mathcal{T}\exp[-i\int\hat{H}(t)dt]|\psi(0)\rangle$ for (a)
and (b), and $|\psi(t)\rangle=\exp[-i\int\hat{H}(0)t]|\psi(0)\rangle$ for (c).
The parameters are $M=12,U_{0}=30,\delta=2,J=3,\omega=0.0001$. We find that
particles in scattering state also have opportunity to interact, so the
modulation of interaction can take effect on the three-particle scattering
states.

##  S3 Higher order correlation of the BIC and quasi-BIC

We have shown the density distribution $\langle\hat{n}_{j}\rangle$ of two BICs
under open boundary condition at ${\varphi}=0$ with eigenvalues
$\epsilon=321.9826$ and $\epsilon=320.8269$ and quasi-BICs assisted with the
type-(ii) cluster bands at $t=0$ with the highest energy and middle energy
under periodic boundary condition in the main text.

Figure S2: (a,c) the second-correlation and (b,d) the third correlation of two
BICs at ${\varphi}=0$ under open boundery condition with eigenvalues
$\epsilon=321.9826$ and $\epsilon=320.8269$ respectively. (e,g) the second-
correlation and (f,h) the third correlation of two quasi-BICs given by the the
maximally localized Wannier state of the highest and middle type-(ii) bands at
$t=0$ under periodic boundery condition respectively. The parameters are
$M=30,U_{0}=300,\delta=20,J=1$.

The density can be regarded as the first-order correlation of the states.
Here, to get a full picture of the two states, we respectively calculate the
second- and third-order correlation functions,

$\displaystyle
C_{ij}^{(2)}=\langle\psi|\hat{a}_{i}^{{\dagger}}\hat{a}_{j}^{{\dagger}}\hat{a}_{j}\hat{a}_{i}|\psi\rangle,$
(S7) $\displaystyle
C_{ijk}^{(3)}=\langle\psi|\hat{a}_{i}^{{\dagger}}\hat{a}_{j}^{{\dagger}}\hat{a}_{k}^{{\dagger}}\hat{a}_{k}\hat{a}_{j}\hat{a}_{i}|\psi\rangle.$

Figs. S2(a,c) and (b,d) show the second- and third-correlation for the two
BICs in the main text, and (e,g) and (f,h) for quasi-BICs in the main text,
respectively. Parameters set as $M=30,U_{0}=300,\delta=20,J=1$. We can clearly
see that the BICs are bound pair by the single-particle standing wave, while
quasi-BICs are bound pair in the standing wave of the other particle.

##  S4 Demonstration of multiparticle MLWS as eigenstate of projected position
operator

A multiparticle MLWS can be obtained by minimizing its spread functional
[PhysRevA.95.063630S]. A spread functional of multiparticle Wannier states in
the $\mathcal{M}$ cluster band is written as

$\Omega=\sum_{m\in\mathcal{M}}\langle\hat{x}^{2}\rangle_{m}-\langle\hat{x}\rangle^{2}_{m},$
(S8)

with $\langle\hat{x}^{2}\rangle_{m}=\langle
W_{m}(0)|\hat{x}^{2}|W_{m}(0)\rangle$ and $\langle\hat{x}\rangle_{m}=\langle
W_{m}(0)|\hat{x}|W_{m}(0)\rangle$. The spread functional can be decomposed to
two terms conveniently,

$\Omega=\Omega_{I}+\Omega_{V},$ (S9)

where

$\displaystyle\Omega_{I}=\sum_{m\in\mathcal{M},n\notin\mathcal{M},R}|\langle
W_{n}(R)|\hat{x}|W_{m}(0)\rangle|^{2},$ (S10)
$\displaystyle\Omega_{V}=\sum_{m,n\in\mathcal{M}}\sum_{Rn\neq 0m}|\langle
W_{n}(R)|\hat{x}|W_{m}(0)\rangle|^{2}.$

Here, It can be proved that $\Omega_{I}$ is gauge invariant under the unitary
transformation
$|u_{m}(\kappa)\rangle\rightarrow\sum_{n\in\mathcal{M}}U_{mn}(\kappa)|u_{n}(\kappa)\rangle$.
So, the purpose of minimizing the spread functional is equal to minimizing
$\Omega_{V}$. In one dimensional systems, a convenient method is to use the
projected position operator $\hat{P}\hat{x}\hat{P}$, where

$\hat{P}=\sum_{m\in\mathcal{M},R}|W_{m}(R)\rangle\langle
W_{m}(R)|=\sum_{m\in\mathcal{M},\kappa}|u_{m}(\kappa)\rangle\langle
u_{m}(\kappa)|$ (S11)

is the projected operator for the cluster of targeted multiparticle Bloch
bands. Here, $|W_{m}(R)\rangle$ and $|u_{m}(\kappa)\rangle$ represent a
multiparticle Wannier state and multiparticle Bloch state, respectively.
Briefly, It can be seen that $\Omega_{V}$ vanishes when choosing the
multiparticle Wannier states $|W_{m}(0)\rangle$ to be the eigenstate of the
projected position operator $\hat{P}\hat{x}\hat{P}$ with an eigenvalue
$X_{0m}$, that is

$\displaystyle\langle W_{n}(R)|\hat{x}|W_{m}(0)\rangle$ $\displaystyle=\langle
W_{n}(R)|\hat{P}\hat{x}\hat{P}|W_{m}(0)\rangle$ (S12)
$\displaystyle=X_{0m}\delta_{R,0}\delta_{m,n},$

where $m,n\in\mathcal{M}$. So, by this way one can obtain the multiparticle
MLWS.

##  S5 Derivation of effective Hamiltonian for the BIC

In this section, we show how to derive effective Hamiltonian for the BICs
under open boundary condition, where singular value decomposition (SVD) and
some semi-analytical derivation are used [4]. A three-particle BIC
$|\psi\rangle$ with a basis of Fock states can be reshaped to
$|\psi\rangle=\psi_{i,j,k}|\psi_{i,j,k}\rangle$, with the index
$i,j,k\in[1,M]$ for the three particles. Due to the symmetry of the bosonic
particles, any swap of $i,j,k$ does not change the amplitude $\psi_{i,j,k}$.
So, we have a tensor whose elements are $\psi_{i,j,k}$, which can be reshaped
from a $M\times M\times M$ tensor to a $M\times M^{2}$ matrix. After that, we
denote the elements of the matrix as $\tilde{\psi}_{i,r}$, where
$r=(j-1)\times M+k$. Performing SVD, there are two largest singular values
that play the dominant role. Then, $\tilde{\psi}_{i,r}$ can be represented by
two terms

$\tilde{\psi}_{i,r}\approx D_{11}S_{i1}W_{1r}+D_{22}S_{i2}W_{2r}.$ (S13)

$D_{11}$ and $D_{22}$ are the two largest singular values, $S_{i1},S_{i2}$
represent a single-particle state, and $W_{1r},W_{2r}$ represent two-particle
state. We denote localized single-particle state denoted as $S^{(l)}$, and the
other extended state as $S^{(f)}$. For the two-particle state, we denote the
two-particle localized state as $W^{(l,l)}$, and the state of one localized
particle and one extended particle as $W^{(l,f)}$. Then, $\tilde{\psi}_{i,r}$
can be represented as

$\tilde{\psi}_{i,r}\approx
D_{11}S^{(l)}_{i}W^{(l,f)}_{r}+D_{22}S^{(f)}_{i}W^{(l,l)}_{r}.$ (S14)

We next reshape and perform SVD on $W^{(l,f)}$, and obtain a localized single-
particle state $\mu^{(l)}$ and an extended state $\mu^{(f)}$,

$W^{(l,f)}_{j,k}\approx\mu^{(l)}_{j}\mu^{(f)}_{k}+\mu^{(l)}_{k}\mu^{(f)}_{j}.$
(S15)

We denote $S^{(f)}\approx\mu^{(f)}$ as $\varphi^{(f)}$, the rest localized
part as $\chi^{(l)}$, and consider the exchange symmetry of bosons, the
amplitudes of the there-particle BIC take the form of

$\psi_{i,j,k}=\varphi^{(f)}_{i}\chi^{(l)}_{j,k}+\varphi^{(f)}_{j}\chi^{(l)}_{i,k}+\varphi^{(f)}_{k}\chi^{(l)}_{i,j}.$
(S16)

Motivated by this form, the three-particle Hamiltonian can be represented by

$\hat{H}=\hat{H}^{1}\otimes I^{2}\otimes I^{3}+I^{1}\otimes\hat{H}^{2}\otimes
I^{3}+I^{1}\otimes I^{2}\otimes H^{3}+U,$ (S17)

where $\hat{H}$ with the superscript $1,2,3$ is linear hopping term, which
works on the subspace of three individual particles, respectively. $I$ is the
identity operator, and $U$ is the interacting term of $\hat{H}$. Submitting
this form to stationary Schrödinger equation, in the form of matrix operation,
one can obtain

$\sum_{i_{1},i_{2},i_{3}}(\hat{H}^{1}\otimes I^{2}\otimes
I^{3}+I^{1}\otimes\hat{H}^{2}\otimes I^{3}+I^{1}\otimes I^{2}\otimes
H^{3}+U)_{i_{1},i_{2},i_{3}}^{j_{1},j_{2},j_{3}}\psi_{i_{1},i_{2},i_{3}}=E\psi_{j_{1},j_{2},j_{3}}.$
(S18)

Here, $i_{1},i_{2},i_{3}$ and $j_{1},j_{2},j_{3}$ are the matrix indices of
the Hamiltonian for the first, second and third particles. The elements of the
matrix of $U$ are given by

$U_{i_{1},i_{2},i_{3}}^{j_{1},j_{2},j_{3}}=\delta_{i_{1},j_{1}}\delta_{i_{2},j_{2}}\delta_{i_{3},j_{3}}(U_{i_{1}}\delta_{i_{1},i_{2}}+U_{i_{3}}\delta_{i_{1},i_{3}}+U_{i_{2}}\delta_{i_{2},i_{3}}),$
(S19)

where $U_{i}=U_{0}+\delta\cos(2\pi\beta i+{\varphi})$. Numerically, we find
the single-particle extended state $\varphi^{(f)}$ is close to the eigenstate
of $\hat{H}^{1},\hat{H}^{2},\hat{H}^{3}$, that is

$\hat{H}^{1(2,3)}\varphi\approx\epsilon_{n}\varphi,$ (S20)

where $\epsilon_{n}$ is corresponding eigenvalue. Hereafter, we omit the
superscript of $\varphi$ and $\chi$ for brevity. Then, we can obtain

$\displaystyle\epsilon_{n}\varphi_{j_{1}}\chi_{j_{2},j_{3}}+\hat{H}^{1}_{i_{1},j_{1}}\varphi_{j_{2}}\chi_{i_{1},j_{3}}+\hat{H}^{1}_{i_{1},j_{1}}\varphi_{j_{3}}\chi_{i_{1},j_{2}}+$
(S21)
$\displaystyle\hat{H}_{i_{2},j_{2}}^{2}\varphi_{j_{1}}\chi_{i_{2},j_{3}}+\epsilon_{n}\varphi_{j_{2}}\chi_{j_{1},j_{3}}+\hat{H}_{i_{2},j_{2}}^{2}\varphi_{j_{3}}\chi_{j_{1},i_{2}}+$
$\displaystyle\hat{H}_{i_{3},j_{3}}^{3}\varphi_{j_{1}}\chi_{j_{2},i_{3}}+\hat{H}_{i_{3},j_{3}}^{3}\varphi_{j_{2}}\chi_{j_{1},i_{3}}+\epsilon_{n}\varphi_{j_{3}}\chi_{j_{1},j_{2}}+$
$\displaystyle(\varphi_{j_{1}}\chi_{j_{2},j_{3}}+\varphi_{j_{2}}\chi_{j_{1},j_{3}}+\varphi_{j_{3}}\chi_{j_{1},j_{2}})(U_{j_{1}}\delta_{j_{1},j_{2}}+U_{j_{3}}\delta_{j_{1},j_{3}}+U_{j_{2}}\delta_{j_{2},j_{3}})$
$\displaystyle=E(\varphi_{j_{1}}\chi_{j_{2},j_{3}}+\varphi_{j_{2}}\chi_{j_{1},j_{3}}+\varphi_{j_{3}}\chi_{j_{1},j_{2}}),$

The same subscript of $i_{1},i_{2},i_{3}$ with $\hat{H}^{1,2,3}$ needs to be
summed over this index. Next, multiplying the equation by
$\varphi_{j_{1}}^{*}$ and summing over $j_{1}$, we have

$\epsilon_{n}\chi_{j_{2},j_{3}}+\hat{H}^{2}_{i_{2},j_{2}}\chi_{i_{2},j_{3}}+\hat{H}^{3}_{i_{3},j_{3}}\chi_{j_{2},i_{3}}+(2U_{j_{2}}|\varphi_{j_{2}}|^{2}\chi_{j_{2},j_{3}}+2U_{j_{3}}|\varphi_{j_{3}}|^{2}\chi_{j_{2},j_{3}}+U_{j_{2}}\chi_{j_{2},j_{3}}\delta_{j_{2},j_{3}})=E\chi_{j_{2},j_{3}}.$
(S22)

In the above derivation, we assume that $\varphi_{j_{1}}\chi_{j_{2},j_{3}}=0$
if $j_{1}=j_{2}$, or $j_{1}=j_{3}$, or $j_{1}=j_{2}=j_{3}$. This is because
$\chi$ should have little overlap to $\varphi$ for the BICs, and the energies
of type-(ii) states have a large gap to those of type-(iii) states. From the
above equation, we can extract the effective Hamiltonian for $\chi$,

$\hat{H}_{\rm
eff}=J\sum_{j}(\hat{a}_{j+1}^{{\dagger}}\hat{a}_{j}+h.c.)+\sum_{j}2U_{j}{|\varphi_{j}^{(f)}|}^{2}\hat{n}_{j}+\frac{1}{2}\sum_{j}U_{j}\hat{a}_{j}^{{\dagger}}\hat{a}_{j}^{{\dagger}}\hat{a}_{j}\hat{a}_{j}.$
(S23)

Here, $U_{j}=U_{0}+\delta$cos$(2\pi\beta j+{\varphi})$, and
$2U_{j}{|\varphi_{j}^{(f)}|}^{2}$ is the effective modulated on-site
potential.

##  S6 Topological pumping of quasi-BIC in subspace

Figure S3: Comparison of energy spectrum (a,b) and dynamics (c,d) between the
total Hilbert space and the reduced Hilbert space expanded by type-(ii) Fock
states. The dashed line denotes the c.m. position of the bound pair. The
parameters are $M=12,U_{0}=300,J=1,\delta=20$ and $\omega=0.0005$.

As the system size and particle number increase, the numerical simulation of
dynamics costs enormous computational resources. For strong interaction
leading to small correlation between the different types of states, we can
simulate the dynamics in a subspace of a certain type which the initial state
belongs. The subspace of type-(ii) states is expanded by the basis
$\mathcal{V}=\\{|0,....,n_{i}=2,...,n_{j}=1,...\rangle\\}$, where two bosons
in the same site are bounded and the other boson at a different site. The
dimension of this subspace is $M(M-1)$ for a $M$-sites lattice system, which
is much less than $(M+2)(M+1)M/6$ of the total Hilbert space. Figs. S3(a) and
(b) show the energy spectrum of type-(ii) states as a function of time
obtained in the full Hilbert space and type-(ii) subspace, respectively. Figs.
S3(c) and (d) respectively show the density distribution
$\langle\hat{n}_{j}\rangle$ as a function of time during a pumping cycle and
in the total Hilbert space and type-(ii) subspace, where the initial state is
the quasi-BIC that fills the highest type-(ii) state band. The parameters are
set $M=12,U_{0}=300,J=1,\delta=20$, and $\omega=0.0005$. Both the energy
spectrum and the dynamics are almost the same, indicating validation of the
projected subspace method. This method can also apply to the type-(i,iii)
subspaces, and it is obviously beneficial to the numerical calculation of
many-body systems.

##  S7 Robustness against disorder for topological pumping of the quasi-BIC

In this section, we show that topological pumping of the quasi-BIC can be
robust to disorder to some extent. We add an extra disordered onsite energies
to the Hamiltonian,

$\hat{H}=J\sum_{j}(\hat{a}_{j+1}^{{\dagger}}\hat{a}_{j}+h.c.)+\frac{1}{2}\sum_{j}U_{j}({\varphi})\hat{n}_{j}(\hat{n}_{j}-1)+\sum_{j}FV_{j}\hat{n}_{j},$
(S24)

where $V_{j}$ are random numbers in the range of $(0,1)$, and $F$ is the
strength of disorder. Without loss of generality, we take a group of random
values of $V_{j}$ in numerical calculation,

$\\{V_{j}\\}=\\{0.8308,0.5853,0.5497,0.9172,0.2858,0.7572,0.7537,0.3804,0.5678,0.0759,0.0540,0.5308\\}.$
(S25)

Setting the quasi-BIC filling the highest type-(ii) band without disorder as
initial state, the time-dependent adiabatic evolution governed by
$|\psi(t)\rangle=\mathcal{T}\exp[-i\int\hat{H}(t)dt]|\psi(0)\rangle$ is shown
in Fig. S4, where $F=1$ in (a), $F=3$ in (b) and $F=5$ in (c). The other
parameters are set as $M=12,U_{0}=300,\delta=20,J=1,\omega=0.0005$, and
periodic boundary condition is adopted. As can be seen, although it is a
common trend that increasing disordered strength will break the topological
pumping, with disorder strength to some extent such as $F=1$, the topological
pumping is robust and persists well.

Figure S4: Density distribution as a function of time in a time-dependent
adiabatic evolution with different disordered strengths $F=1,3,5$ for (a,b,c),
respectively. The parameters are $M=12,U_{0}=300,\delta=20,J=1,\omega=0.0005$,
and $F=1$ for (a), $F=3$ for (b), $F=5$ for (c).

## References

  * Qin  _et al._ [2014] X. Qin, Y. Ke, X. Guan, Z. Li, N. Andrei, and C. Lee, Statistics-dependent quantum co-walking of two particles in one-dimensional lattices with nearest-neighbor interactions, Phys. Rev. A 90, 062301 (2014).
  * Takahashi [1977] M. Takahashi, Half-filled hubbard model at low temperature, J. Phys. C: Solid State Phys. 10, 1289 (1977).
  * Ke  _et al._ [2017] Y. Ke, X. Qin, Y. S. Kivshar, and C. Lee, Multiparticle wannier states and thouless pumping of interacting bosons, Phys. Rev. A 95, 063630 (2017).
  * Zhang  _et al._ [2023] N. Zhang, Y. Ke, L. Lin, L. Zhang, and C. Lee, Stable interaction-induced anderson-like localization embedded in standing waves, New J. Phys. 25, 043021 (2023).