# Multi-particle quantum walks in one-dimensional lattice

Daer Feng  School of Physics, Sun Yat-Sen University, Guangzhou, Guangdong
510275, China  Shengshi Pang<EMAIL_ADDRESS>School of Physics, Sun
Yat-Sen University, Guangzhou, Guangdong 510275, China

(August 28, 2024)

###### Abstract

Quantum walk is a counterpart of classical random walk in the quantum regime
that exhibits non-classical behaviors and outperforms classical random walk in
various aspects. It has been known that the spatial probability distribution
of a single-particle quantum walk can expand quadratically in time while a
single-particle classical random walk can do only linearly. In this paper, we
analytically study the discrete-time quantum walk of non-interacting multiple
particles in a one-dimensional infinite lattice, and investigate the role of
entanglement and exchange symmetry in the position distribution of the
particles during the quantum walk. To analyze the position distribution of
multi-particle quantum walk, we consider the relative distance between
particles, and study how it changes with the number of walk steps. We compute
the relative distance asymptotically for a large number of walk steps and find
that the distance increases quadratically with the number of walk steps. We
also study the extremal relative distances between the particles, and show the
role of the exchange symmetry of the initial state in the distribution of the
particles. Our study further shows the dependence of two-particle
correlations, two-particle position distributions on the exchange symmetry,
and find exponential decrement of the entanglement of the extremal state with
the number of particles.

##  I Introduction

Quantum walk is an interesting quantum protocol with distinct properties such
as quadratic speed up in the position distribution, exponentially fast hitting
[1, 2], etc., compared to its classical counterpart. These non-classical
properties have motivated intense research interest in exploring the intrinsic
properties of quantum walks and exploiting its potential applications
accordingly [3, 4]. And it has been shown to be a powerful tool in various
quantum computing and quantum information tasks such as spatial search [5, 6,
7], graph isomorphism testing [8, 9, 10, 11], quantum simulation [12],
universal quantum computation [13, 14, 15] and teleportation [16, 17].

Two types of formulations have been developed for quantum walks, the discrete-
time quantum walk (DTQW) and the continuous-time quantum walk (CTQW). The CTQW
was proposed by Farhi and co-workers in 1998 as a powerful computational model
to penetrate a decision tree in polynomial time to outperform its classical
counterpart which requires exponential time to penetrate some specific types
of trees [18]. In parallel with the continuous-time quantum walk, the
discrete-time model of quantum walk was introduced by Aharonov et al. and by
Ambainis et al. as another potential quantum-enhanced computational tool [19].
The difference lies in that the walkers evolve with discrete-time steps and
extra degrees of freedom—coin, which is usually a two-dimensional quantum
system, govern the shift direction of a walker at each step. The coin flip
operator can be position-dependent [20], which leads to the quantum-walk-based
search algorithms [5]. A walker can possess more than one coins and the coins
can be entangled [21, 22, 23]. In two or higher dimensional situation, the
evolution of quantum walk can be either constituted of alternating coin flips
and position displacements in each dimension [24] or that the conditional
displacement only occurs in one dimension at each step depending on the coin
state [25]. It was conjectured that the existence of the coin space enlarges
the total Hilbert space of the walkers compared with the continuous-time walk,
and thus it enhances the distinguishing power of the quantum walks on non-
isomorphic graphs [11].

The coin flip has substantial influence on the evolution of quantum walks. The
Hadamard coin is regarded as a balanced coin which flips both the upside and
the down side of coins to the superpositions of the two sides with equal
probabilities. It has been analytically shown that the position deviation of a
single-particle quantum walker on a one-dimensional lattice driven by a
Hadamard coin is $O(t)$ after $t$ steps[26], implying the quadratic speed up
of the spreading of the particles on the lattice. In addition to the Hadamard
coin, the Grover coin has also been studied in quantum walk, which enables the
quantum-walk-based search algorithm to reach the time complexity $O(\sqrt{N})$
for a database of size $N$ [7], in analogy to the well-known Grover search
algorithm [27, 28]. It has been widely believed that there is close connection
between the increased position spreading of quantum walks and the quadratic
speed-up of quantum walk-based algorithms [3]. Moreover, Tregenna et al. have
analytically investigated the most general two-dimensional unitary coin flips
of single-particle quantum walks and showed the effect of the coin flip on the
asymmetry of the walker position distribution [29]. And the effect of
decoherence in the coins of quantum walk has been studied, showing that
decoherence could slow down the spatial spreading of the particles in quantum
walk and lead to approximately linear growth of the position distribution
which is typically a classical behavior [30].

The exchange symmetry of the initial state of the particles can also influence
the quantum walk significantly. The work of Omar et al. [31] elucidates how
the entanglement and the relative phases in the initial coin states influence
the evolution of two-particle quantum walk, which shows that the probability
of finding at least one particle in a particular position in a one-dimensional
lattice is dependent on the entanglement of the particles and can be
drastically different from that of two uncorrelated particles. It is shown in
[32] that the presence of the correlations between the coins can induce the
development of spatial correlations between the walkers, and the conditional
shift generates entanglement between the coin and position spaces [33]. The
entanglement production is also reported in multi-particle quantum walks [34].
If we treat the position space of the walkers as an open system with the coins
as the environment, entanglement and entropy are produced during the quantum
walk [35].

The interactions between particles bring a variety of interesting properties
to quantum walks. In the continuous-time quantum walk, as one increases the
nearest-neighbor interaction of three bosons or fermions, the ballistic
expansion of the density distribution gradually fades, and the density-density
correlation indicates the occurrence of co-walking or independent walking
particles with different exchange statistics [36, 37]. Bloch oscillations are
observed in continuous-time quantum walk with nearest-neighbor interaction
[36]. In the discrete-time situation, Berry et al. [11] introduced two types
of interactions—$\mathbb{\mathds{1}}$ interaction and $\pi$-phase interaction,
leading to similar graph-distinguishing powers. A kind of quantum Hash scheme
is developed based on such interactions to improve the security of the Hash
scheme [38]. Similar gaslike interaction is introduced in [39], and the
position distribution and entanglement of the walkers are analyzed. A long-
range Coulomb-like interaction has also been studied for discrete-time quantum
walk, which is effectively to append an auxillary interaction right after the
original non-interacting evolution, and the spectrum of the system is
investigated [40]. Moreover, the interaction in quantum walk have been found
useful in some graph algorithms. In the graph isomorphism testing problem, the
Hubbard interaction prevents the evolution operator from depending only on the
strongly regular graph family parameters and thus enables the quantum walk to
distinguish non-isomorphic strongly regular graphs within the same family
[10].

Quantum walk can simulate physical phenomena in condensed matter physics such
as Anderson localization and Bloch oscillations. Anderson localization [41],
the absence of diffusion of a quantum mechanical state, can be produced in
quantum walk when the walk is subjected to disorder [42, 43] or when the walk
is inhomogeneous in the transition of a walker from one site to the others
dependent on both the coin state and the position state [44, 45]. The relation
between the localization and the degeneracy of the eigenstates of quantum walk
is studied in [46, 47]. Quan Lin et al. simulated the non-Hermitian
topological Anderson insulator experimentally by disordered photonic quantum
walks [48].

Quantum walk can be realized experimentally in a vast variety of systems [49]
such as superconducting qubits [50, 51], optical lattices [52], photons [53,
54, 55, 12, 56, 57, 58, 59, 60, 61, 62], Bose-Einstein condensates [63, 64]
and trapped ions [65, 66, 67, 68], etc. Experimental realizations of quantum
walk via photons [69, 70] and particularly Hong-Ou-Mandel effect [71, 72]
exhibit different dynamics depending on the photon indistinguishability and
statistics. Vertex search in graphs and graph isomorphism testing have been
implemented on silicon photonic quantum walk processors [73].

The dynamics of the entangled particles in multi-particle quantum walks is
interesting, and much effort has been devoted to multi-particle quantum walk
with interactions in a vast literature as reviewed above. However, the effect
of exchange statistics and entanglement between non-interacting particles on
the quantum walk is still left unkown so far. In this paper, we focus on the
evolution of non-interacting multi-particle discrete-time quantum walk in one-
dimensional lattice to fill this gap. We will study the general evolution of
quantum walk with non-ineracting particles, and obtain the asymptotic relative
distance between the particles for arbitrary initial states when the number of
walk steps is sufficiently large. The extremal relative distances will be
derived when all the particles evolve from the same position initially, and an
interesting result is that the initial states that reach the the extremal
relative distances have specific exchange symmetry between the particles, as
will be illustrated by graphs. We will also investigate how the evolution of
the position distribution and entanglement of the particles are affected by
the symmetry of the initial state.

The paper is structured as follows. In Sec. II, we provide preliminaries for
multi-particle discrete-time quantum walk. Sec. III studies the general
evolution of multi-particle quantum walk and obtains the asymptotic average
distance to characterize the spreading property of multi-particle quantum
walk. In Sec. IV, we apply the results to the case where the coins of the
initial state are entangled, and derive the upper and lower bound of the
average relative distance between the particles. In Sec. V, we show that the
exchange symmetries of the eigenstates of the relative distance operator are
critical to the average distance, and explore the entanglement and position
distribution properties of the walkers.

##  II Preliminaries

In this section, we introduce preliminaries of multi-particle discrete-time
quantum walks and notations that will be used in this paper.

For discrete-time quantum walks on an infinite one-dimensional lattice,
particles are located on discrete sites and possess extra degrees of freedom
called coins to determine the shift direction of the particles at each walk
step. The Hilbert space of single-particle quantum walk on a one-dimensional
lattice can be decomposed as
$\mathscr{H}=\mathscr{H}_{position}\otimes\mathscr{H}_{coin}$, in which a
single-particle state at position $x$ with the coin upwards or downwards can
be written as $|x\rangle\otimes|\uparrow\rangle$ or
$|x\rangle\otimes|\downarrow\rangle$, respectively. The state of the particle
evolves according to the unitary operator
$\hat{U}=\hat{S}\cdot(\hat{I}\otimes\hat{C})$, where $\hat{C}$ is a coin
operator that flips the coin and $\hat{S}$ is the shift operator that changes
the position of the particle dependent on its coin state. $\hat{S}$ usually
takes the following form,

$\hat{S}=\hat{Q}\otimes|\uparrow\rangle\langle\uparrow|+\hat{Q}^{\dagger}\otimes|\downarrow\rangle\langle\downarrow|,$
(1)

where $\hat{Q}$ and $\hat{Q}^{\dagger}$ are the shift operators in the
positive and negative $x$ directions respectively,

$\hat{Q}=\sum_{x}|x+1\rangle\langle
x|,\;\hat{Q}^{\dagger}=\sum_{x}|x-1\rangle\langle x|.$ (2)

We generalize the quantum walks on a one-dimensional lattice to non-
interacting multi-particle cases. The evolution operator is the tensor product
of single-particle evolution operators, i.e., $\hat{U}_{n}=\hat{U}^{\otimes
n}$ for $n$ particles. If we define the initial state as $|\psi_{0}\rangle$,
the final state after $t$ steps of evolution is
$|\psi_{t}\rangle=\text{$\hat{U}_{n}^{t}$}|\psi_{0}\rangle.$ A general
$n$-particle initial state for discrete-time quantum walk can be written as

$|\psi_{0}\rangle=\sum
a_{x_{1}s_{1},\cdots,x_{n}s_{n}}|x_{1}s_{1}\rangle\otimes\cdots\otimes|x_{n}s_{n}\rangle.$
(3)

After $t$ steps of evolution, the final state of the particles can be obtained
by performing single-particle quantum walk on each particle individually,

$|\psi_{t}\rangle=\sum
a_{x_{1}s_{1},\cdots,x_{n}s_{n}}\hat{U}^{t}|x_{1}s_{1}\rangle\otimes\cdots\otimes\hat{U}^{t}|x_{n}s_{n}\rangle.$
(4)

The quantum walk with two particles is studied in [31], and the bosonic or
fermionic nature of the two particles is shown to have significant influence
on the evolution. While for more particles, Chandrashekar et al. [74]
considered such multi-particle quantum walk for distinguishable particles.
When the particles are uncorrelated initially, the final state is
$\hat{U}^{t}|\phi_{1}\rangle\otimes\hat{U}^{t}|\phi_{2}\rangle\otimes\cdots\otimes\hat{U}^{t}|\phi_{n}\rangle$,
where $|\phi_{i}\rangle$ is the initial single-particle state of the $i$th
particle.

The coin operator can be generalized by a weight parameter. A general SU(2)
coin up to a global phase is

$\hat{C}=\left[\begin{array}[]{cc}\sqrt{\rho}&\sqrt{1-\rho}e^{i\theta}\\\
\sqrt{1-\rho}e^{i\varphi}&-\sqrt{\rho}e^{i(\varphi+\theta)}\end{array}\right],$
(5)

where $0\leq\rho\leq 1$ and $0\leq\theta,\varphi\leq\pi$. We choose the
Hadamard operator as the coin operator throughout this paper.

In a classical multi-particle random walk, the average distance between the
walkers grows linearly with the number of walk steps (see Appendix A1). In
this paper, we will show that for a general multi-particle quantum walk, the
average distance grows quadratically with the number of walk steps. This has
been well known for single-particle quantum walk [19], and we will prove
similar result for multiple particles, particularly considering the the
evolution of different particles may interfere each other when they are
initially entangled.

##  III Distance of multi-particle quantum walks

In this section, we study the relative distance between particles and
calculate the average distance of a general $n$-particle state for a
sufficiently large number of walk steps.

The relative distance between particles can be defined as

$\hat{D}=\sum_{i}\left[\hat{x}_{i}-\frac{1}{n}\sum_{j}\hat{x}_{j}\right]^{2}.$
(6)

where $\hat{x}_{i}$ represents the position of the $i$th particle, and
$\hat{D}$ can be decomposed into two parts,

$\hat{D}=\frac{n-1}{n}\sum_{i}\hat{x}_{i}^{2}+\frac{1}{n}\sum_{j,j^{\prime}}\hat{x}_{j}\hat{x}_{j^{\prime}},$
(7)

consisting of single-particle operators and two-particle operators. In this
paper, we mainly concern with how the average relative distance with respect
to the evolved multi-particle state changes with the evolution of the walkers
on the lattice.

###  III.1 Diagonalization of the evolution operator

Diagonalizing the evolution operator in the momentum space of the particles
will help simplify the calculation of the evolution. By following the idea of
discrete-time Fourier transform [30, 75], one can transform the state of the
particles from the position space into the momentum space, the basis of which
is given by

$|K\rangle=\frac{1}{\sqrt{2\pi}}\sum_{x}e^{iKx}|x\rangle,$ (8)

where $|K\rangle$ preserves the orthonormality $\langle
K^{\prime}|K\rangle=\delta(K^{\prime}-K)$. It can be verified that $|K\rangle$
is the eigenstate of the shift operator $\hat{Q}$ (2) with eigenvalue
$e^{-iK}$ and the eigenstate of $\hat{Q}^{\dagger}$ with eigenvalue $e^{iK}$.
Thus, applying the conditional shift $\hat{S}$ (1) to $|K\rangle$ merely
induces extra phases in the coin space, and the evolution of the particles in
the joint Hilbert space of the position and coin can be reduced to one in the
coin space with an arbitrary fixed momentum $K$,

$\hat{U}|K\rangle\otimes|s\rangle=|K\rangle\otimes\left(e^{-iK}|\uparrow\rangle\langle\uparrow|+e^{iK}|\downarrow\rangle\langle\downarrow|\right)\hat{C}|m\rangle.$
(9)

Denote
$\hat{U}_{K}\equiv(e^{-iK}|\uparrow\rangle\langle\uparrow|+e^{iK}|\downarrow\rangle\langle\downarrow|)\hat{C}$,
and choose the Hardamard operator $\hat{H}$ to be the coin flip operator,

$\hat{H}=\frac{1}{\sqrt{2}}\left(|\uparrow\rangle\langle\uparrow|+|\uparrow\rangle\langle\downarrow|+|\downarrow\rangle\langle\uparrow|\right)-\frac{1}{\sqrt{2}}|\downarrow\rangle\langle\downarrow|,$
(10)

then $\hat{U}_{K}$ can be written as

$\hat{U}_{K}=\frac{e^{-iK}}{\sqrt{2}}\left(|\uparrow\rangle\langle\uparrow|+|\uparrow\rangle\langle\downarrow|\right)+\frac{e^{iK}}{\sqrt{2}}\left(|\downarrow\rangle\langle\uparrow|-|\downarrow\rangle\langle\downarrow|\right).$
(11)

Further diagonalization of the evolution $\hat{U}_{K}$ in the coin space gives
the eigenstates of $\hat{U}_{K}$, $|d_{j}(K)\rangle(j=1,2)$,

$\displaystyle|d_{1}(K)\rangle=-\frac{e^{-iK}}{\sqrt{2N(K)}}|\uparrow\rangle+\frac{1}{\sqrt{2N(\pi-K)}}|\downarrow\rangle,$
(12)
$\displaystyle|d_{2}(K)\rangle=\frac{e^{-iK}}{\sqrt{2N(\pi-K)}}|\uparrow\rangle+\frac{1}{\sqrt{2N(K)}}|\downarrow\rangle,$

with $\lambda_{j}\,(j=1,2)$ as the corresponding eigenvalues,

$\lambda_{j}(K)=(-1)^{j}\frac{1}{2}\sqrt{3+2\cos(2K)}-\frac{i\sin(K)}{\sqrt{2}},$
(13)

where

$N(K)=(1+\cos^{2}K)+\cos K\sqrt{1+\cos^{2}K}.$ (14)

Finally, we arrive at the diagonalized representation of the evolution
operator,

$\hat{U}=\int{\rm d}K|K\rangle\langle
K|\otimes\left(\lambda_{1}|d_{1}(K)\rangle\langle
d_{1}(K)|+\lambda_{2}|d_{2}(K)\rangle\langle d_{2}(K)|\right),$ (15)

where we apply momentum-dependent evolution in the coin space.

###  III.2 Average distance of the particles

Now, we move on to the average distance between particles. As the exact
distance between particles is complex generally, we will focus on the
asymptotic results for a suffciently large number of walk steps. For a general
$n$-particle quantum walk with initial state given in Eq. (3), the average
distance of the particles after $t$ steps of evolution is
$\langle\psi_{0}|(\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}|\psi_{0}\rangle$,
and it can be broken down into the sum of following terms

$\displaystyle\langle x_{1}^{\prime}s_{1}^{\prime}|\cdots\langle
x_{n}^{\prime}s_{n}^{\prime}|(\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}|x_{1}s_{1}\rangle\cdots|x_{n}s_{n}\rangle$
(16) $\displaystyle=$ $\displaystyle\frac{n-1}{n}\sum_{i}\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle\prod_{\gamma=1,\gamma\neq
i}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}$
$\displaystyle-\frac{1}{n}\sum_{j,k,j\neq k}\langle
x_{j}^{\prime}s_{j}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{j}\hat{U}^{t}|x_{j}s_{j}\rangle\langle
x_{k}^{\prime}s_{k}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{k}\hat{U}^{t}|x_{k}s_{k}\rangle$
$\displaystyle\cdot\prod_{\gamma=1,\gamma\neq
j,k}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}.$

$\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}$ is $1$
when $x_{\gamma}^{\prime}=x_{\gamma}$, $s_{\gamma}^{\prime}=s_{\gamma}$ and
$0$ otherwise, implying the summands in the second line can contribute to the
average distance only when the states of the particles other than the $i$th
particle in the bra $\langle x_{1}^{\prime}s_{1}^{\prime}|\cdots\langle
x_{n}^{\prime}s_{n}^{\prime}|$ and the ket
$|x_{1}s_{1}\rangle\cdots|x_{n}s_{n}\rangle$ match. And similar constraint
applies to the summands in the third and fourth lines.

To work out the matrix element $\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{j}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle$,
note that Eq. (15) is the spectral decomposition of the evolution operator and
it brings each basis state $|K\rangle\otimes|d\rangle$ just the corresponding
eigenvalue, so it preserves the eigenstate and applying $\hat{U}$
consecutively yields the power of the same eigenvalue. By the completeness of
the basis, we have

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle$
(17) $\displaystyle=$ $\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm d}K\sum_{d,d^{\prime}}\langle
x_{i}^{\prime}s_{i}^{\prime}|K^{\prime}d^{\prime}\rangle\langle
K^{\prime}d^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|Kd\rangle\langle
Kd|x_{i}s_{i}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}\langle x_{i}^{\prime}|K^{\prime}\rangle\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle K|x_{i}\rangle\langle
d|s_{i}\rangle\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle$
$\displaystyle\cdot\langle K^{\prime}|\hat{x}_{i}^{2}|K\rangle\langle
d^{\prime}|d\rangle\langle d|\hat{U}_{K}^{t}|d\rangle.$

As transforming $x^{m}$ into the momentum space gives the periodic extension
of the derivative of Dirac delta function with period $2\pi$, it gives

$\displaystyle\langle K^{\prime}|\hat{x}_{i}^{2}|K\rangle=$
$\displaystyle\sum_{x=-\infty}^{\infty}\langle
K^{\prime}|\hat{x}_{i}^{2}|x\rangle\langle
x|K\rangle=-\delta^{(2)}(K^{\prime}-K).$ (18)

So, the matrix element (17) becomes

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle$
(19) $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}e^{iK^{\prime}x_{i}^{\prime}}e^{-ikx_{i}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle d|s_{i}\rangle$
$\displaystyle\cdot\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle\delta^{(2)}(K^{\prime}-K)\langle
d^{\prime}|d\rangle\langle d|\hat{U}_{K}^{t}|d\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}[t^{2}I_{A2}(x_{i}^{\prime}-x_{i},s_{i}^{\prime},s_{i})+tI_{A1}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})$
$\displaystyle+I_{AC}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})+I_{Ao}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i};t)].$

Similarly,

$\displaystyle\langle
x_{j}^{\prime}s_{j}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{j}\hat{U}^{t}|x_{j}s_{j}\rangle$
(20) $\displaystyle=$
$\displaystyle\frac{i}{2\pi}[tI_{B}(x_{j}^{\prime}-x_{j},s_{j}^{\prime},s_{j})+I_{B_{1}}(x_{j}^{\prime},x_{j},s_{j}^{\prime},s_{j})$
$\displaystyle+I_{Bo}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i};t)].$

In the above equations, $I_{A2}$, $I_{A1}$, $I_{AC}$, $I_{B}$, $I_{B_{1}}$ are
coefficients for different powers of $t$ and can be calculated using the
residue theorem. The integrals in the term $t^{2}$ are

$\displaystyle I_{A2}(x,s^{\prime},s)=$
$\displaystyle\begin{cases}(1-\sqrt{2})f(0)&x=0,s^{\prime}=s;\\\ f(x)&x\neq
0,s^{\prime}=s;\\\ 0&s^{\prime}\neq s,\end{cases}$ (21)

$\displaystyle I_{B}(x,s^{\prime},s)=$
$\displaystyle\begin{cases}(-1)^{s}iI_{A2}(x,s^{\prime},s)&s^{\prime}=s;\\\
\frac{i}{2}[f(x)+f(x-(-1)^{s}\cdot 2)]&s^{\prime}\neq s,\end{cases}$ (22)

where $f(x)=\sqrt{2}\pi\left(\sqrt{2}-1\right)^{|x|}\cos\left(\pi x/2\right)$.
When $s$ is $\uparrow$, it act as number 1, and $\downarrow$ as number $0$,
i.e. when $s=\downarrow$, $(-1)^{s}=1$. The detailed expressions of integrals
$I_{A1}$, $I_{B_{1}}$, $I_{AC}$, $I_{Ao}$ and $I_{Bo}$ are given in Appendix
A3.

$I_{Ao}$ and $I_{Bo}$ are integrals with rapidly varying phases when $t$ is
large, which can be calculated using the stationary phase approximation [76].
Eq. (19) and Eq. (20) are the functions of the position and coin states
$|x_{i}s_{i}\rangle$, $|x_{i}^{\prime}s_{i}^{\prime}\rangle$ as well as the
number of time steps $t$, with the highest order of $t$ being $t^{2}$. After a
long time evolution, the $t^{2}$ term becomes dominant, so the other lower-
order terms can be dropped. By plugging Eq. (19) and (20) into Eq. (16) and
summing up all the terms in Eq. (16), it yields the average distance for a
given initial state in Eq. (3).

##  IV Maximization of particle distance

In this section, we consider the case that the walkers are initially
uncorrelated on the spatial lattice but entangled in the coin space. We
concern how the relative distance is dependent on the initial coin state, and
obtain the upper and lower bounds of the distance when the number of walk
steps is large.

###  IV.1 The upper and lower bounds of the coefficient

To simplify the representation of the coin state, we regard $\uparrow$ as
number $1$ and $\downarrow$ as number $0$ as mentioned before and denote the
coin state of the $i$th particle as $|s_{i}\rangle$. Then the vector
$|s_{i}\rangle$ takes the value of $|1\rangle$ or $|0\rangle$, and the
combination of the coin states of all particles can be treated as a binary
number and transformed into a decimal number, i.e., $|s_{1}s_{2}\cdots
s_{n}\rangle=|2^{n-1}\times s_{1}+2^{n-2}\times
s_{2}+\cdots+2s_{n-1}+s_{n}\rangle=|\xi\rangle$, where $\xi$ is the decimal
representation of $s_{1}s_{2}\cdots s_{n}$. We assume the spatial state and
the coin state of the particles are not entangled initially, and the initial
state can be written as

$|\Psi\rangle=|x_{1}\rangle|x_{2}\rangle\cdots|x_{n}\rangle\otimes\sum_{\xi=0}^{2^{n}-1}a_{\xi}|\xi\rangle.$
(23)

The average distance after $t$ steps of evolution
$\langle\Psi|(\hat{U}^{\dagger})^{t}\hat{D}\hat{U}^{t}|\Psi\rangle$ can broken
down into the power series of $t$. For a large $t$, we only consider the
coefficient of $t^{2}$ in the average distance, denoted as $c_{2}$, and
rewrite it in a matrix form,

$c_{2}=[a_{0}^{*},a_{1}^{*},\cdots,a_{2^{n-1}}^{*}]M[a_{0},a_{1},\cdots,a_{2^{n-1}}]^{T}.$
(24)

The matrix $M\in\mathbb{R}^{2^{n}\times 2^{n}}$ is obtained in Appendix A4.
The diagonal elements are

$\displaystyle M_{ii}=$ $\displaystyle(n-1)\left(1-\frac{1}{\sqrt{2}}\right)$
(25)
$\displaystyle+\frac{1}{n}[n-(2W(i)-n)^{2}]\left(1-\frac{1}{\sqrt{2}}\right)^{2},$

and the off-diagonal elements are

$\displaystyle M_{jk}=$
$\displaystyle\frac{2}{n}\left(1-\frac{1}{\sqrt{2}}\right)^{2}[\delta_{d(j,k),1}(n-1$
(26) $\displaystyle-2\min\\{W(j),W(k)\\})-\delta_{d(j,k),2}],$

where $i,j,k\in\\{0,1,2,\cdots,2^{n}-1\\}$ and $j\neq k$, and $W(i)$ is the
Hamming weight of the binary representation of integer $i$, i.e. the number of
“$1$”s in the $n$ bit binary representation of integer $i$. $d(j,k)$ is the
Hamming distance between integer $j$ and $k$, i.e., the number of bit
positions where $i$ differs from $j$. $\delta$ is the Kronecker delta
function. Eq. (24) can be considered as a quadratic form of the paramters
$a_{0},a_{1},\cdots,a_{2^{n-1}}$ and their conjugates, and the eigenvectors of
the matrix $M$ are the extremal points for the eigenvalues of $M$.

It can be proven that the eigenvalues of the matrix $M$ are

$\eta_{k}=\left(1-\frac{1}{\sqrt{2}}\right)^{2}\left[\frac{8W(k)(n-W(k))}{n}+(n-1)\sqrt{2}\right],$
(27)

where $k\in\\{0,2,\cdots,2^{n}-2\\}$. The mathematical detail of the
derivation is lengthy and left in the Appendix A5. Therefore, $0\leq W(k)\leq
n-1$, $W(k)\in\mathbb{Z}$ and there are $C_{n-1}^{W(k)}$ different values of
$k$ that correspond to the same value of $W(k)$. Notice that for different
values of $k$, e.g., $k$ and $k^{\prime}$, $\eta_{k}$ and $\eta_{k^{\prime}}$
take the same value when $W(k)=n-W(k^{\prime})$, which should to be taken into
account in calculating the degeneracies of the eigenvalues. The corresponding
eigenvectors for $\eta_{k}$ turn out to be $P_{k}$ and $P_{k+1}$. They are
unnormalized eigenstates, which are initial coin states of the walk. And
$P_{k}$ is the $k$th column of matrix $P$ and $P$ consist of $2\times 2$ sub-
matrices

$P_{ik}=(-1)^{c(i,k)}D^{d(i,k)}=(-1)^{c(i,k)}\left[\begin{array}[]{cc}0&1\\\
1&2\end{array}\right]^{d(i,k)},$ (28)

where $i,k\in\\{0,2,\cdots,2^{n}-2\\}$, and they represent the position of the
upper left corner element of the sub-matrices in the matrix $P$. $c(i,j)$
represents the number of bits in the binary form on which $i$ and $j$ have a
common $1$. For instance, when $n=2$,

$P=\begin{array}[]{c}\begin{array}[]{cc}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\overbrace{\rule{31.29802pt}{0.0pt}}^{P_{00}}}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\overbrace{\rule{31.29802pt}{0.0pt}}^{P_{02}}}\end{array}\\\\[-2.0pt]
\left[\begin{array}[]{cc|cc}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1\begin{array}[]{c}\end{array}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\begin{array}[]{c}\hphantom{}\end{array}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0\begin{array}[]{c}\end{array}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\begin{array}[]{c}\hphantom{}\end{array}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}\\\
\hline\cr{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}0}\begin{array}[]{c}\end{array}&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}1}\begin{array}[]{c}\hphantom{}\end{array}&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}-1}&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}0}\\\
{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}1}\begin{array}[]{c}\end{array}&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}2}\begin{array}[]{c}\hphantom{}\end{array}&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}0}&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}-1}\end{array}\right]\\\\[-7.0pt]
\begin{array}[]{cc}{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underbrace{\rule{31.29802pt}{0.0pt}}_{P_{20}}}&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}\underbrace{\rule{31.29802pt}{0.0pt}}_{P_{22}}}\end{array}\\\
\end{array}$ (29)

The eigenvalues of the matrix $M$ give the range of $c_{2}$, $\eta_{min}\leq
c_{2}\leq\eta_{max}$. According to Eq. (27), the largest eigenvalue for a
given number of particles $n$ is

$\eta_{max}(n)=\begin{cases}\left(1-\frac{1}{\sqrt{2}}\right)^{2}[(n-1)\sqrt{2}+2n],&n\
even,\\\
\left(1-\frac{1}{\sqrt{2}}\right)^{2}[(n-1)\sqrt{2}+2n-\frac{2}{n}],&n\
odd,\end{cases}$ (30)

while the smallest eigenvalue is

$\eta_{min}(n)=\left(1-\frac{1}{\sqrt{2}}\right)^{2}\sqrt{2}(n-1).$ (31)

The relation of the number of particles to the eigenvalues of $M$ and their
degeneracies is explicitly shown in Fig. 1. As the number of particles $n$
approaches infinity,

$\eta_{max}(n)/\eta_{min}(n)\rightarrow 1+\sqrt{2},$ (32)

and the lower bound and upper bound of the eigenvalues grow (approximately)
linearly.

Figure 1: Eigenvalues and their degeneracies for different number $n$ of
particles. The y-coordinates of these points correspond to the magnitudes of
the eigenvalues, and x-coordinates the numbers of particles. The numbers next
to the points are the corresponding degeneracies of the eigenvalues. Figure 2:
Vertex partitions for $n=5$ particle eigenstates. The numbers near the
vertices mark different particles that number $i$ denotes the $i$th bit of
$k$. From left to right and then top to bottom, each cell of the table
represents complete graph(s) corresponding to a different even number $k$
increasing from $0$ to $2^{n}-2$ respectively. When $k$ is odd, the vertex
partition is the same as the situation of $k-1$. The same background colors is
assigned to graphs of the same structure, which also indicates the same
eigenvalue.

##  V Symmetry and entanglement in coin space

In this section, we show that the coin states for the maximum and minimum
distance between the particles is related to the exchange symmetry of the coin
states.

###  V.1 Partial exchange symmetry of the eigenstates

To study the relation between the bounds of the particle distance and the
symmetry of the initial state, we define “partial exchange symmetry” here. We
say that the eigenstates of matrix $M$ possess partial exchange symmetry if
the state is unchanged after arbitrary exchange among a given subset of the
particles. The number of possible ways $p_{k}$ to exchange two particles that
preserve the state is related to the eigenvalue $\eta_{k}$ by
$4(C_{n}^{2}-p_{k})=\eta_{k}$.

If we visualize the partial exchange symmetry by assigning the particles to
the vertices of an undirected graph and connecting a pair of particles with an
edge if the state is unchanged after swapping those two particles, it turns
out that the graph may be composed of two disconnected subgraphs, each of
which is a complete graph respectively. It can be proven that (see Appendix
A6) if we denote $n_{\uparrow}$ and $n_{\downarrow}$ as the number of “$1$”s
and “$0$”s in the binary form of $k$ ($k$ even) respectively, then
$p_{k}=\frac{n_{\uparrow}(n_{\uparrow}-1)}{2}+\frac{n_{\downarrow}(n_{\downarrow}-1)}{2}$,
$n_{\uparrow}+n_{\downarrow}=n$, and that the two disconnected subgraphs have
$n_{\uparrow}$ and $n_{\downarrow}$ vertices respectively, as shown in Fig. 2.
We can then deduce that the larger the $p_{k}$ is, the smaller the eigenvalue
$\eta_{k}$ will be, and thus a greater degree of partial exchange symmetry of
the eigenstate (the initial state) results in a smaller spreading distance of
the quantum walk, which is essentially induced by a greater extent of
destructive interference between the walking routes of the particles and more
particles oriented in the same direction and thus closer to each other.

When $k=0$, all of the $n$ particles are assigned to the same subgraph, which
corresponds to the greatest value of $p_{k}$ and thus the smallest $\eta_{k}$.
The particles are closest to each other in this case and it possess the most
exchange symmetry. When the particles are distributed as uniformly as possible
between the two subgraphs, we obtain the smallest value of $p_{k}$ and the
particles turn out to be the most distant. Note that while the eigenstates of
the matrix $M$ possess partial exchange symmetry, there is no anti-symmetry in
the eigenstates when $n>2$, in contrast to the case of two particles [31].

It should be clarified that the bipartite structures of the graphs in Fig. 2
only represent the partial exchange symmetry between particles, but does mean
that the two parts of the graph are not correlated. Actually, the bipartite
entanglement between the two disconnected subgraphs is non-zero. For the
eigenstate associated with the $k$th eigenvector, if $k\neq 0$, the
entanglement between the two parts of the graphs can be quantified by the von
Neumann entropy,

$S(\rho_{\uparrow})=-\sum_{i=1,2}\nu_{i}{\rm\log}\nu_{i},$ (33)

in which $\nu_{1}=[1+(3-2\sqrt{2})^{n+2}]^{-1}$,
$\nu_{2}=[1+(3+2\sqrt{2})^{n+2}]^{-1}$ for an even $k$ and
$\nu_{1}=[1+(3-2\sqrt{2})^{n}]^{-1}$, $\nu_{2}=[1+(3+2\sqrt{2})^{n}]^{-1}$ for
an odd $k$. $\rho_{\uparrow}$ denotes the reduced density matrix of the
particles associated with $n_{\uparrow}$. Whether $k$ is odd or even, the
expressions of $\nu_{1}$ and $\nu_{2}$ are independent of $k$. For detailed
derivation, we refer the readers to AppendixA7. When $n$ is sufficiently
large, $S(\rho_{\uparrow})$ decreases exponentially to $0$ with respect to
$n$.

Figure 3: The average of $\hat{x}_{i}^{2}$ and two-particle correlations of
evolved eigenstate corresponding to $k=(1001010)_{2}$, $n=7$, $t=100$. (a)
Average value of $\hat{x}_{i}^{2}$ of each particle. Particle label $i$
corresponds to the $i$th bit in $k$. Bits $1,4,6$ are $1$, and $2,3,4,7$ are
$0$. (b) Two-particle correlations. There are three types of two-particle
combinations: both two particles from the same subgraphs (there are two
disconnected subgraphs and thus two different situations) and one particle
from each of the two subgraphs. Therefore, the 3D bar graph is divided into
four regions by two vertical planes. In the current model, the evolution does
not change the exchange symmetry of particles in the subgraphs. Thus the
correlations between arbitrary two particles from the same subgraph are
identical. Figure 4: Bipartite entanglement of coins over $t=50$ steps of
evolution. Numbers label the particles and “$|$” represents the partition,
i.e. in the case of $2|134567$, the von Neumann entropy is calculated by
reducing the state to the density matrix of particle $2$ to obtain the
bipartite entanglement between particle 2 and particle 1,3,4,5,6,7.

###  V.2 Dynamics of the eigenstates

In the final part of this section, we try to dig a bit deeper into the
dynamics of the particles in the two subgraphs. The initial positions of all
particles are assumed to be 0.

For the case of $n=7$, the average of $\hat{x}_{i}^{2}$ and two-particle
position correlations of the eigenstate after an evoluton of $t=100$ steps are
shown in Fig. 3. The binary number for the label of the eigenstate is
$k=(1001010)_{2}$. They are essentially the components of the average distance
operator in Eq. (7). It can be seen that the particles in the same subgraph
posssess the same average value of $\hat{x}_{i}^{2}$, since these particles
are identical due to the partial exchange symmetry and the evolution does not
break this symmetry. The position correlations of two particles are positive
if the particles are from the same subgraph and negative if from different
subgraphs, because the particles from the two subgraphs are initially placed
in the position $0$ but oriented in different directions by their coins. The
two-particle position correlations always keep the same if either particle is
replaced by another particle in the same subgraph, since picking any particle
from the same subgraph does not make difference due to the partial exchange
symmetry.

Figure 5: (a) Two-particle position distribution of particle $1$ and 2, 1 and
4, 3 and 5, in the eigenstate with $k=(1001010)_{2}$, evolving for
$t=10,20,30$ steps respectively. For particle 1 and 2, as they belong to
different subgraphs, they are more likely to be found in different positions,
particularly particle 1 in negative positions and particle 2 in positive
positions. For two particles in the same subgraph, i.e., 1 and 4, 3 and 5, on
account of the exchange symmetry among particles within the same subgraph,
they are close to each other. (b) Two-particle position distribution in the
eigenstate associated with $k=(0110100)_{2}$ after $t=30$ steps of evolution.

A more insightful observation of the position distribution over time for
$k=(1001010)_{2}$ is presented in Fig. 5(a). Particles in the same subgraph
are more likely to be close to each other. But particles from different
subgraphs, such as particle 1 and particle 2, are more distant. As is shown in
Fig. 5(a), particle 1 is more likely to be found in negative positions while
particle 2 is more likely to be found in positive positions, because we can
see from the reduced density matrix of the initial coin states that the coins
are in average oriented in the opposite directions. As a comparison, Fig. 5(b)
is the position distribution of the evolved eigenstate with $k=(0110100)_{2}$
after $t=30$ steps of evolution. The Hamming weight of $(1001010)_{2}$ and
$(0110100)_{2}$ are identical, and $(0110100)_{2}$ can be perceived as
relabeling particle $1,4,6$ as $2,3,5$ and $2,3,5$ as $1,4,6$ in
$(1001010)_{2}$. Referring to the expression of $P_{ik}$ (28), the two
eigenstates are essentially the same after the relabeling and therefore have
the same average distance, particularly the same second order term of $t$.

The bipartite entanglement of the particles of the initial state
$k=(1001010)_{2}$ is small compared with that of the evolved state, as shown
in Fig. 4. We use the von Neumann entropy as the entanglement measure. The
entanglement increases rapidly at the beginning and then starts to oscillate.

##  VI Conclusions

In this paper, we mainly studied the discrete-time multi-particle quantum
random walks driven by the Hadamard coin. We defined the distance operator to
quantify how distant the particles are apart from each other, and obtain an
asymptotic result of the distance for a large number of evolution steps.

To explore the influence of the initial state on the distance after the
evolution, we considered a situation where the initial state may be entangled
in the coin space but separable in the position space. We focused on the
leading quadratic term of $t$ in the average distance between particles, since
the $t^{2}$ term is dominant when $t$ is large enough. It turns out that when
there are two or more particles, the lower bounds and upper bounds for the
quadratic term are non-zero and increase with the number of particles, and
that the eigenstates of the evolved distance operator, which are essentially
the local extremas of the quadratic terms of $t$, possess partial exchange
symmetry. The number of different swaps between particles that preserve the
state is related to the eigenvalue.

It is found that the particles can be grouped into two subsets for each
eigenstate of the relative distance such that the eigenstate has partial
exchange symmetry with respect to either subset, i.e., permutations within
either subset of the particles do not change the eigenstate. But it does not
mean that the two subsets of particles are not correlated. We studied the
entanglement between the two subsets quantified by the von Neumann entropy,
which is identical among different eigenstates for the same number of
particles but is variable for different numbers of particles. As the
eigenvalues are all positive, there is no localization phenomena in this
scenario. Owing to the symmetry of the evolution, the particles in the same
subset remain identical during the evolution. We numerically illustrated the
two-particle correlations, two-particle position distribution and entanglement
between the two subgraphs of the eigenstates to reveal their dynamical
properties due to the symmetry of the states.

###### Acknowledgements.

The authors acknowledge the helpful discussions with Litong Deng from Tsinghua
University and Prof. Xiaoyong Fu from Sun Yat-sen University. This work is
supported by the National Natural Science Foundation of China (Grant No.
12075323) and (national) college students innovation and entrepreneurship
training program, Sun Yat-sen University.

## References

  * Childs  _et al._ [2003] A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman, Proceedings of the thirty-fifth ACM symposium on Theory of computing - STOC ’03 , 59 (2003).
  * Tang  _et al._ [2018a] H. Tang, C. Di Franco, Z.-Y. Shi, T.-S. He, Z. Feng, J. Gao, K. Sun, Z.-M. Li, Z.-Q. Jiao, T.-Y. Wang, M. S. Kim, and X.-M. Jin, Nature Photonics 12, 754 (2018a).
  * Venegas-Andraca [2012] S. E. Venegas-Andraca, Quantum Information Processing 11, 1015 (2012).
  * Kadian  _et al._ [2021] K. Kadian, S. Garhwal, and A. Kumar, Computer Science Review 41, 100419 (2021).
  * Portugal [2013] R. Portugal,  _Quantum Walks and Search Algorithms_ (Springer New York, New York, NY, 2013).
  * Wong  _et al._ [2016] T. G. Wong, L. Tarrataca, and N. Nahimov, Quantum Information Processing 15, 4029 (2016).
  * Shenvi  _et al._ [2003] N. Shenvi, J. Kempe, and K. B. Whaley, Physical Review A 67, 052307 (2003).
  * Rudinger  _et al._ [2012] K. Rudinger, J. K. Gamble, M. Wellons, E. Bach, M. Friesen, R. Joynt, and S. N. Coppersmith, Physical Review A 86, 022334 (2012).
  * Douglas and Wang [2008] B. L. Douglas and J. B. Wang, Journal of Physics A: Mathematical and Theoretical 41, 075303 (2008).
  * Gamble  _et al._ [2010] J. K. Gamble, M. Friesen, D. Zhou, R. Joynt, and S. N. Coppersmith, Physical Review A 81, 052313 (2010).
  * Berry and Wang [2011] S. D. Berry and J. B. Wang, Physical Review A 83, 042317 (2011).
  * Aspuru-Guzik and Walther [2012] A. Aspuru-Guzik and P. Walther, Nature Physics 8, 285 (2012).
  * Childs  _et al._ [2013] A. M. Childs, D. Gosset, and Z. Webb, Science 339, 791 (2013).
  * Singh  _et al._ [2021] S. Singh, P. Chawla, A. Sarkar, and C. M. Chandrashekar, Scientific Reports 11, 11551 (2021).
  * Childs [2009] A. M. Childs, Physical Review Letters 102, 180501 (2009).
  * Yan  _et al._ [2019a] Y.-G. Yan, S.-N. Cao, W.-F. Cao, D. Li, Y.-H. Zhou, and W.-M. Shi, Modern Physics Letters B 33, 1950070 (2019a).
  * Li  _et al._ [2019] H.-J. Li, X.-B. Chen, Y.-L. Wang, Y.-Y. Hou, and J. Li, Quantum Information Processing 18, 266 (2019).
  * Farhi and Gutmann [1998] E. Farhi and S. Gutmann, Physical Review A 58, 915 (1998).
  * Kempe [2003] J. Kempe, Contemporary Physics 44, 307 (2003).
  * Ahmad  _et al._ [2020] R. Ahmad, U. Sajjad, and M. Sajid, Communications in Theoretical Physics 72, 065101 (2020).
  * Liu and Petulante [2009] C. Liu and N. Petulante, Physical Review A 79, 032312 (2009).
  * Venegas-Andraca [2018] S. E. Venegas-Andraca, in  _2018 IEEE Globecom Workshops (GC Wkshps)_ (2018) pp. 1–6.
  * Venegas-Andraca  _et al._ [2005] S. E. Venegas-Andraca, J. L. Ball, K. Burnett, and S. Bose, New Journal of Physics 7, 221 (2005).
  * Roldán  _et al._ [2013] E. Roldán, C. Di Franco, F. Silva, and G. J. de Valcárcel, Physical Review A 87, 022336 (2013).
  * Annabestani  _et al._ [2010] M. Annabestani, M. R. Abolhasani, and G. Abal, Journal of Physics A: Mathematical and Theoretical 43, 075301 (2010).
  * Konno [2002] N. Konno, Quantum Information Processing 1, 345 (2002).
  * Grover [1997] L. K. Grover, Physical Review Letters 79, 325 (1997).
  * Grover [1996] L. K. Grover, in  _Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing - STOC ’96_ (ACM Press, Philadelphia, Pennsylvania, United States, 1996) pp. 212–219.
  * Tregenna  _et al._ [2003] B. Tregenna, W. Flanagan, R. Maile, and V. Kendon, New Journal of Physics 5, 83 (2003).
  * Brun  _et al._ [2003] T. A. Brun, H. A. Carteret, and A. Ambainis, Physical Review A 67, 032304 (2003).
  * Omar  _et al._ [2006] Y. Omar, N. Paunkovic, L. Sheridan, and S. Bose, Physical Review A 74, 042304 (2006).
  * Orthey and Angelo [2019] A. C. Orthey and R. M. Angelo, Physical Review A 100, 042110 (2019).
  * Abal  _et al._ [2006] G. Abal, R. Siri, A. Romanelli, and R. Donangelo, Physical Review A 73, 042302 (2006).
  * Goyal and Chandrashekar [2010] S. K. Goyal and C. M. Chandrashekar, Journal of Physics A: Mathematical and Theoretical 43, 235303 (2010).
  * Vallejo  _et al._ [2019] A. Vallejo, A. Romanelli, R. Donangelo, and R. Portugal, Physical Review A 99, 032319 (2019).
  * Cai  _et al._ [2021] X. Cai, H. Yang, H.-L. Shi, C. Lee, N. Andrei, and X.-W. Guan, Physical Review Letters 127, 100406 (2021).
  * Qin  _et al._ [2014] X. Qin, Y. Ke, X. Guan, Z. Li, N. Andrei, and C. Lee, Physical Review A 90, 062301 (2014).
  * Li  _et al._ [2013] D. Li, J. Zhang, F.-Z. Guo, W. Huang, Q.-Y. Wen, and H. Chen, Quantum Information Processing 12, 1501 (2013).
  * Costa  _et al._ [2019] P. C. S. Costa, F. de Melo, and R. Portugal, Physical Review A 100, 042320 (2019).
  * Alonso-Lobo  _et al._ [2018] C. Alonso-Lobo, M. Martinez-Quesada, M. Hinarejos, G. J. de Valcarcel, and E. Roldan, Journal of Physics a-Mathematical and Theoretical 51, 455301 (2018).
  * Anderson [1958] P. W. Anderson, Physical Review 109, 1492 (1958).
  * Crespi  _et al._ [2013] A. Crespi, R. Osellame, R. Ramponi, V. Giovannetti, R. Fazio, L. Sansoni, F. De Nicola, F. Sciarrino, and P. Mataloni, Nature Photonics 7, 322 (2013).
  * Mendes  _et al._ [2019] C. V. C. Mendes, G. M. A. Almeida, M. L. Lyra, and F. A. B. F. de Moura, Physical Review E 99, 022117 (2019).
  * Shikano and Katsura [2010] Y. Shikano and H. Katsura, Physical Review E 82, 031122 (2010).
  * Konno [2010] N. Konno, Quantum Information Processing 9, 405 (2010).
  * Inui and Konno [2005] N. Inui and N. Konno, Physica A: Statistical Mechanics and its Applications 353, 133 (2005).
  * Inui  _et al._ [2004] N. Inui, Y. Konishi, and N. Konno, Physical Review A 69, 052323 (2004).
  * Lin  _et al._ [2022] Q. Lin, T. Li, L. Xiao, K. Wang, W. Yi, and P. Xue, Nature Communications 13, 3229 (2022).
  * Manouchehri and Wang [2014] K. Manouchehri and J. Wang,  _Physical Implementation of Quantum Walks_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2014).
  * Gong  _et al._ [2021] M. Gong, S. Wang, C. Zha, M.-C. Chen, H.-L. Huang, Y. Wu, Q. Zhu, Y. Zhao, S. Li, S. Guo, H. Qian, Y. Ye, F. Chen, C. Ying, J. Yu, D. Fan, D. Wu, H. Su, H. Deng, H. Rong, K. Zhang, S. Cao, J. Lin, Y. Xu, L. Sun, C. Guo, N. Li, F. Liang, V. M. Bastidas, K. Nemoto, W. J. Munro, Y.-H. Huo, C.-Y. Lu, C.-Z. Peng, X. Zhu, and J.-W. Pan, Science 372, 948 (2021).
  * Yan  _et al._ [2019b] 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, Science 364, 753 (2019b).
  * Preiss  _et al._ [2015] P. M. Preiss, R. Ma, M. E. Tai, A. Lukin, M. Rispoli, P. Zupancic, Y. Lahini, R. Islam, and M. Greiner, Science 347, 1229 (2015).
  * Schreiber  _et al._ [2010] A. Schreiber, K. N. Cassemiro, V. Potoček, A. Gábris, P. J. Mosley, E. Andersson, I. Jex, and C. Silberhorn, Physical Review Letters 104, 050502 (2010).
  * Broome  _et al._ [2010] M. A. Broome, A. Fedrizzi, B. P. Lanyon, I. Kassal, A. Aspuru-Guzik, and A. G. White, Physical Review Letters 104, 153602 (2010).
  * Schreiber  _et al._ [2011] A. Schreiber, K. N. Cassemiro, V. Potoček, A. Gábris, I. Jex, and C. Silberhorn, Physical Review Letters 106, 180403 (2011).
  * Peruzzo  _et al._ [2010] A. Peruzzo, M. Lobino, J. C. F. Matthews, N. Matsuda, A. Politi, K. Poulios, X.-Q. Zhou, Y. Lahini, N. Ismail, K. Worhoff, Y. Bromberg, Y. Silberberg, M. G. Thompson, and J. L. OBrien, Science 329, 1500 (2010).
  * Schreiber  _et al._ [2012] A. Schreiber, A. Gábris, P. P. Rohde, K. Laiho, M. Štefaňák, V. Potoček, C. Hamilton, I. Jex, and C. Silberhorn, Science 336, 55 (2012).
  * Tang  _et al._ [2018b] H. Tang, X.-F. Lin, Z. Feng, J.-Y. Chen, J. Gao, K. Sun, C.-Y. Wang, P.-C. Lai, X.-Y. Xu, Y. Wang, L.-F. Qiao, A.-L. Yang, and X.-M. Jin, Science Advances 4, eaat3174 (2018b).
  * Cardano  _et al._ [2015] F. Cardano, F. Massa, H. Qassim, E. Karimi, S. Slussarenko, D. Paparo, C. de Lisio, F. Sciarrino, E. Santamato, R. W. Boyd, and L. Marrucci, Science Advances 1, e1500087 (2015).
  * Defienne  _et al._ [2016] H. Defienne, M. Barbieri, I. A. Walmsley, B. J. Smith, and S. Gigan, Science Advances 2, e1501054 (2016).
  * Qiang  _et al._ [2016] X. Qiang, T. Loke, A. Montanaro, K. Aungskunsiri, X. Zhou, J. L. O’Brien, J. B. Wang, and J. C. F. Matthews, Nature Communications 7, 11511 (2016).
  * Jeong  _et al._ [2013] Y.-C. Jeong, C. Di Franco, H.-T. Lim, M. Kim, and Y.-H. Kim, Nature Communications 4, 2471 (2013).
  * Xie  _et al._ [2020] D. Xie, T.-S. Deng, T. Xiao, W. Gou, T. Chen, W. Yi, and B. Yan, Physical Review Letters 124, 050502 (2020).
  * Dadras  _et al._ [2018] S. Dadras, A. Gresch, C. Groiseau, S. Wimberger, and G. S. Summy, Physical Review Letters 121, 070402 (2018).
  * Huerta Alderete  _et al._ [2020] C. Huerta Alderete, S. Singh, N. H. Nguyen, D. Zhu, R. Balu, C. Monroe, C. M. Chandrashekar, and N. M. Linke, Nature Communications 11, 3720 (2020).
  * Schmitz  _et al._ [2009] H. Schmitz, R. Matjeschk, C. Schneider, J. Glueckert, M. Enderlein, T. Huber, and T. Schaetz, Physical Review Letters 103, 090504 (2009).
  * Zähringer  _et al._ [2010] F. Zähringer, G. Kirchmair, R. Gerritsma, E. Solano, R. Blatt, and C. F. Roos, Physical Review Letters 104, 100503 (2010).
  * Xue  _et al._ [2009] P. Xue, B. C. Sanders, and D. Leibfried, Physical Review Letters 103, 183602 (2009).
  * Pathak and Agarwal [2007] P. K. Pathak and G. S. Agarwal, Physical Review A 75, 032351 (2007).
  * Sansoni  _et al._ [2012] L. Sansoni, F. Sciarrino, G. Vallone, P. Mataloni, A. Crespi, R. Ramponi, and R. Osellame, Physical Review Letters 108, 010502 (2012).
  * Karczewski  _et al._ [2019] M. Karczewski, R. Pisarczyk, and P. Kurzyński, Physical Review A 99, 042102 (2019).
  * Simon  _et al._ [2020] D. S. Simon, S. Osawa, and A. V. Sergienko, Physical Review A 101, 032118 (2020).
  * Qiang  _et al._ [2021] X. Qiang, Y. Wang, S. Xue, R. Ge, L. Chen, Y. Liu, A. Huang, X. Fu, P. Xu, T. Yi, F. Xu, M. Deng, J. B. Wang, J. D. A. Meinecke, J. C. F. Matthews, X. Cai, X. Yang, and J. Wu, Science Advances 7, eabb8375 (2021).
  * Chandrashekar and Busch [2012] C. M. Chandrashekar and T. Busch, Quantum Information Processing 11, 1287 (2012).
  * Ambainis  _et al._ [2001] A. Ambainis, E. Bach, A. Nayak, A. Vishwanath, and J. Watrous, in  _Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing - STOC ’01_ (ACM Press, Hersonissos, Greece, 2001) pp. 37–49.
  * Bleistein and Handelsman [1975] N. Bleistein and R. A. Handelsman, Asymptotic expansions of integrals (1975).

## Appendix: Details of our derivation

In this appendix, we give more details about the derivation. Sec. A1
calculates the distance of multi-walker classical random walk. Sec. A2
illustrates the preliminaries of discrete time multi-particle quantum walks
and process of diagonalizing the evolution operator. Sec. A3 derives the
average distance of the walkers as a function of the evolution time $t$ and
initial state. We consider a situation where the particles are entangled in
the coin space, and express the coefficient $c_{2}$ of $t^{2}$ into matrix
multiplication in Sec. A4. Sec. A5 provides an approach to derive the
eigenvalues and eigenvectors of matrix $M$. We analyse and proof the partial
exchange symmetry of the particles of the eigenstates corresponding to the
eigenvectors in Sec. A6, and then in Sec. A7 we calculate the bipartite
entanglement between the two subsets of particles that the particles within
the same subset possess exchange symmetry.

###  A1 Distance of multi-walker classical random walk

Considering a classical random walk with $n$ independent walkers. We define
the average distance of the particles as

$\bar{D}=\overline{\sum_{i=1}^{n}\left[x_{i}-\frac{1}{n}\sum_{j=1}^{n}x_{j}\right]^{2}}=\frac{n-1}{n}\sum_{i}\overline{x_{i}^{2}}-\frac{1}{n}\sum_{j,k,j\neq
k}\bar{x}_{j}\bar{x}_{k}$ (S1)

$x_{i}=x_{i0}+x_{i1}+\cdots+x_{it}$ is the final position of the $i$th
particle. And $x_{it}\in\\{-1,1\\}$ is the displacement of the $i$th particle
at step $t$. And the probability of $x_{it}$ to be $-1$ or $1$ is equal.
$x_{i0}$ is the initial position of the $i$th particle. Therefore
$\bar{x}_{it}=0$ and $\overline{x_{it}^{2}}=1$.

$\displaystyle\bar{D}=$
$\displaystyle\frac{n-1}{n}\sum_{i}\overline{(x_{i0}+x_{i1}+x_{i2}+\cdots+x_{it})^{2}}-\frac{1}{n}\sum_{j,k,j\neq
k}\bar{x}_{j}\bar{x}_{k}$ (S2) $\displaystyle=$
$\displaystyle\frac{n-1}{n}\sum_{i}\overline{\left(\sum_{\gamma=0}^{t}\sum_{\gamma^{\prime}=0}^{t}x_{i\gamma}x_{i\gamma^{\prime}}\right)}-\frac{1}{n}\sum_{j,k,j\neq
k}\bar{x}_{j}\bar{x}_{k}$ $\displaystyle=$
$\displaystyle\frac{n-1}{n}\sum_{i}\left[x_{i0}^{2}+2x_{i0}\sum_{\gamma=1}^{t}\bar{x}_{i\gamma}+\sum_{\gamma=1}^{t}\sum_{\gamma^{\prime}=1}^{t}\bar{x}_{i\gamma}\bar{x}_{i\gamma^{\prime}}+\sum_{\gamma=1}^{t}\overline{x_{i\gamma}^{2}}\right]-\frac{1}{n}\sum_{j,k,j\neq
k}x_{j0}x_{k0}$ $\displaystyle=$
$\displaystyle\frac{n-1}{n}\sum_{i}x_{i0}^{2}-\frac{1}{n}\sum_{j,k,j\neq
k}x_{j0}x_{k0}+(n-1)t$

$\bar{D}$ grows linearly as the function of $t$.

###  A2 Diagonalization of the evolution operator

The $n$-particle evolution operator is

$\hat{U}_{n}=\hat{U}^{\otimes n},$ (S3)

and $\hat{U}=\hat{S}\cdot(\hat{I}\otimes\hat{C})$, where $\hat{C}$ is the coin
operator that flips the coin and $\hat{S}$ is the shift operator that changes
the position of the particle according to its spin. $\hat{S}$ takes on the
following form $\hat{S}=\sum_{x}|x+1\rangle\langle
x|\otimes|\uparrow\rangle\langle\uparrow|+\sum_{x}|x-1\rangle\langle
x|\otimes|\downarrow\rangle\langle\downarrow|$.

We perform discrete time Fourier transform

$|K\rangle=\frac{1}{\sqrt{2\pi}}\sum_{x}e^{iKx}|x\rangle,$ (S4)

to diagonalize the evolution operator on the position space that

$\hat{U}|K\rangle\otimes|s\rangle=|K\rangle\otimes(e^{-iK}|\uparrow\rangle\langle\uparrow|+e^{iK}|\downarrow\rangle\langle\downarrow|)\hat{C}|m\rangle.$
(S5)

The calculation of the evolution can be converted into applying the operator
$\hat{U}_{K}\equiv(e^{-iK}|\uparrow\rangle\langle\uparrow|+e^{iK}|\downarrow\rangle\langle\downarrow|)\hat{C}$
on the coin space. As a specific case, we choose the Hardamard operator
$\hat{H}$ as the coin operator

$\hat{H}=\frac{1}{\sqrt{2}}\left(|\uparrow\rangle\langle\uparrow|+|\uparrow\rangle\langle\downarrow|+|\downarrow\rangle\langle\uparrow|\right)-\frac{1}{\sqrt{2}}|\downarrow\rangle\langle\downarrow|.$
(S6)

Then $\hat{U}_{K}$ can be written as,

$\hat{U}_{K}=\frac{e^{-iK}}{\sqrt{2}}\left(|\uparrow\rangle\langle\uparrow|+|\uparrow\rangle\langle\downarrow|\right)+\frac{e^{iK}}{\sqrt{2}}\left(|\downarrow\rangle\langle\uparrow|-|\downarrow\rangle\langle\downarrow|\right).$
(S7)

By diagonalizing it in the coin space further on, denoting $|d_{i}(K)\rangle$
as the eigenvectors of the new vector space, the diagonal representation of
operator $\hat{U}_{K}$ is shown in Eq. (S8), where $|d_{1}\rangle$ and
$|d_{2}\rangle$ are the two eigenvectors.

$\displaystyle\hat{U}_{K}=$
$\displaystyle\frac{e^{-iK}\left(-\sqrt{1+6e^{2iK}+e^{4iK}}-e^{2iK}+1\right)}{2\sqrt{2}}|d_{1}\rangle\langle
d_{1}|+\frac{e^{-iK}\left(\sqrt{1+6e^{2iK}+e^{4iK}}-e^{2iK}+1\right)}{2\sqrt{2}}|d_{2}\rangle\langle
d_{2}|$ (S8) $\displaystyle=$ $\displaystyle\lambda_{1}|d_{1}\rangle\langle
d_{1}|+\lambda_{2}|d_{2}\rangle\langle d_{2}|$

The actual values of $\lambda_{1}$ and $\lambda_{2}$ varies if we choose
different branches of the square-root function. Here we define

$\sqrt{1+6e^{2iK}+e^{4iK}}=\sqrt{2}\cos(K)\sqrt{\cos(2K)+3}+i\sqrt{2}\sin(K)\sqrt{\cos(2K)+3},$
(S9)

then

$\displaystyle\lambda_{1}(K)=$
$\displaystyle-\frac{1}{2}\sqrt{3+2\cos(2K)}-\frac{i\sin(K)}{\sqrt{2}},$ (S10)
$\displaystyle\lambda_{2}(K)=$
$\displaystyle+\frac{1}{2}\sqrt{3+2\cos(2K)}-\frac{i\sin(K)}{\sqrt{2}},$

$|d_{1}\rangle=-\frac{e^{-iK}}{\sqrt{2N(K)}}|\uparrow\rangle+\frac{1}{\sqrt{2N(\pi-K)}}|\downarrow\rangle,$
(S11)

$|d_{2}\rangle=\frac{e^{-iK}}{\sqrt{2N(\pi-K)}}|\uparrow\rangle+\frac{1}{\sqrt{2N(K)}}|\downarrow\rangle,$
(S12)

where $N(K)=(1+\cos^{2}K)+\cos K\sqrt{1+\cos^{2}K}=[\sqrt{1+\cos^{2}K}+\cos
K]\sqrt{1+\cos^{2}K}$.

###  A3 Calculation of distance

As is mentioned in the article, the initial state of an $n$-particle quantum
walk can be written as

$|\psi_{0}\rangle=\sum_{x_{1}s_{1}\cdots x_{n}s_{n}}a_{x_{1}s_{1}\cdots
x_{n}s_{n}}|x_{1}s_{1}\rangle|x_{2}s_{2}\rangle\cdots|x_{n}s_{n}\rangle,$
(S13)

where $x_{i}$ and $s_{i}$ are the position and spin of the $i$th particle.

And the distance operator is defined as

$\hat{D}=\frac{n-1}{n}\sum_{i}\hat{x}_{i}^{2}+\frac{1}{n}\sum_{j,j^{\prime}}\hat{x}_{j}\hat{x}_{j^{\prime}}.$
(S14)

Similar to evolving operators in Heisenberg picture, the evolved distance
operator can be decomposed into the the summation of single-particle operators
and two particle-operators.

$\begin{aligned}
(\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}=&\frac{n-1}{n}\sum_{i}(\hat{U}_{i}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}_{i}^{t}-\frac{1}{n}\sum_{j,k,j\neq
k}\end{aligned}(\hat{U}_{j}^{\dagger})^{t}\hat{x_{j}}\hat{U}_{j}^{t}\otimes(\hat{U}_{k}^{\dagger})^{t}\hat{x_{k}}\hat{U}_{k}^{t}.$
(S15)

The subscripts $i,j,k$ represent which particle the operators apply on. The
average distance after $t$ steps of evolution is given by

$\displaystyle\langle\psi_{0}|(\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}|\psi_{0}\rangle=$
$\displaystyle\sum_{x_{1}^{\prime}s_{1}^{\prime}\cdots
x_{n}^{\prime}s_{n}^{\prime}}\sum_{x_{1}s_{1}\cdots
x_{n}s_{n}}a_{x_{1}^{\prime}s_{1}^{\prime}\cdots
x_{n}^{\prime}s_{n}^{\prime}}^{*}a_{x_{1}s_{1}\cdots
x_{n}s_{n}}[\frac{n-1}{n}\sum_{i}\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle\prod_{\gamma=1,\gamma\neq
i}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}$ (S16)
$\displaystyle-\frac{1}{n}\sum_{j,k,j\neq k}\langle
x_{j}^{\prime}s_{j}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{j}\hat{U}^{t}|x_{j}s_{j}\rangle\langle
x_{k}^{\prime}s_{k}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{k}\hat{U}^{t}|x_{k}s_{k}\rangle\prod_{\gamma=1,\gamma\neq
j,k}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}].$

$\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}$ is $1$
if $x_{\gamma}^{\prime}=x_{\gamma}$ and $s_{\gamma}^{\prime}=s_{\gamma}$, and
is $0$ otherwise.

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle$
(S17) $\displaystyle=$ $\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm d}K\sum_{d,d^{\prime}}\langle
x_{i}^{\prime}s_{i}^{\prime}|K^{\prime}d^{\prime}\rangle\langle
K^{\prime}d^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|Kd\rangle\langle
Kd|x_{i}s_{i}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}\langle x_{i}^{\prime}|K^{\prime}\rangle\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle K|x_{i}\rangle\langle
d|s_{i}\rangle\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle\langle
K^{\prime}|\hat{x}_{i}^{2}|K\rangle\langle d^{\prime}|d\rangle\langle
d|\hat{U}_{K}^{t}|d\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}e^{iK^{\prime}x_{i}^{\prime}}e^{-iKx_{i}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle d|s_{i}\rangle\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle\delta^{(2)}(K^{\prime}-K)\langle
d^{\prime}|d\rangle\langle d|\hat{U}_{K}^{t}|d\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d,d^{\prime}}\partial_{K^{\prime}}^{2}[e^{i(K^{\prime}x_{i}^{\prime}-Kx_{i})}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle
d|s_{i}\rangle\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}\lambda_{d}(K)^{t}\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}.$

We used

$\displaystyle\langle K^{\prime}|\hat{x}_{i}^{2}|K\rangle=$
$\displaystyle\sum_{x=-\infty}^{\infty}\langle
K^{\prime}|\hat{x}_{i}^{2}|x\rangle\langle
x|K\rangle=\frac{1}{2\pi}\sum_{x=-\infty}^{\infty}x^{2}e^{-i(K^{\prime}-K)x}=-\delta^{(2)}(K^{\prime}-K),$
(S18)

where we utilized the relation

$\frac{1}{2\pi}\sum_{x=-\infty}^{\infty}x^{m}e^{-i(K^{\prime}-K)x}=(i)^{m}\sum_{l=-\infty}^{\infty}\delta^{(m)}(K^{\prime}-K-2\pi
l),$ (S19)

which is the periodic extension of the derivative of Dirac delta function with
period $2\pi$. We can omit the summation on the right side because the
integration region is $2\pi$, so there is one term remaining.

The order of $t$ is relative to the order of derivatives. If we expand the
derivatives,

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}^{2}\hat{U}^{t}|x_{i}s_{i}\rangle$
(S20) $\displaystyle=$ $\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d,d^{\prime}}\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{k^{\prime}=K}e^{-iKx_{i}}\langle
d|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d,d^{\prime}}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]\partial_{K^{\prime}}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]$
$\displaystyle+e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle\partial_{K^{\prime}}^{2}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]\\}|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]\lambda_{d}^{*}(K^{\prime})^{t}+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]\partial_{K^{\prime}}[\lambda_{d}^{*}(K^{\prime})^{t}]$
$\displaystyle+e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle\partial_{K^{\prime}}^{2}[\lambda_{d}^{*}(K^{\prime})^{t}]\\}|_{|K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle-$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d^{\prime}}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]\partial_{K^{\prime}}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]$
$\displaystyle+e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle\partial_{K^{\prime}}^{2}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]\\}|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}\\}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle$
$\displaystyle+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}t\lambda_{d}^{*}(K)^{-1}\partial_{K}[\lambda_{d}^{*}(K)]e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\\}$ $\displaystyle+e^{iK(x_{i}^{\prime}-x_{i})}\langle
s_{i}^{\prime}|d(k)\rangle t(t-1)\lambda_{d}^{*}(K)^{-2}\langle
d(K)|s_{i}\rangle(\partial_{K}[\lambda_{d}^{*}(K)])^{2}+e^{iK(x_{i}^{\prime}-x_{i})}\langle
s_{i}^{\prime}|d(K)\rangle\langle d(K)|s_{i}\rangle
t\lambda_{d}^{*}(K)^{-1}\partial_{K}^{2}[\lambda_{d}^{*}(K)]\\}\\}$
$\displaystyle-$ $\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d\neq
d^{\prime}}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{k^{\prime}=K}t\lambda_{d^{\prime}}^{*}(K)^{-1}\partial_{K}[\lambda_{d^{\prime}}^{*}(K)]$
$\displaystyle+e^{iKx_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K)\rangle\langle
d^{\prime}(K)|d(K)\rangle\\{t(t-1)\lambda_{d^{\prime}}^{*}(K)^{t-2}(\partial_{K}[\lambda_{d^{\prime}}^{*}(K)])^{2}+t\lambda_{d^{\prime}}^{*}(K)^{t-1}\partial_{K}^{2}[\lambda_{d^{\prime}}^{*}(K)]\\}\\}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d^{\prime}}^{*}(K)^{t}\lambda_{d}(K)^{t}\\}$
$\displaystyle=$
$\displaystyle-\frac{1}{2\pi}[t^{2}I_{A2}(x_{i}^{\prime}-x_{i},s_{i}^{\prime},s_{i})+tI_{A1}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})+I_{AC}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})+I_{Ao}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i};t)].$

Denote

$I_{A2}(x_{i}^{\prime}-x_{i},s_{i}^{\prime},s_{i})=\int_{-\pi}^{\pi}{\rm
d}Ke^{iK(x_{i}^{\prime}-x_{i})}\sum_{d}\langle
s_{i}^{\prime}|d(k)\rangle\langle
d(K)|s_{i}\rangle\lambda_{d}^{*}(K)^{-2}(\partial_{K}[\lambda_{d}^{*}(K)])^{2},$
(S21)

$\displaystyle I_{A1}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})=$
$\displaystyle\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}\lambda_{d}^{*}(K)^{-1}\partial_{K}[\lambda_{d}^{*}(K)]e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle$ (S22) $\displaystyle+e^{iK(x_{i}^{\prime}-x_{i})}\langle
s_{i}^{\prime}|d(K)\rangle\langle
d(K)|s_{i}\rangle\\{-\lambda_{d}^{*}(K)^{-2}(\partial_{K}[\lambda_{d}^{*}(K)])^{2}+\lambda_{d}^{*}(K)^{-1}\partial_{K}^{2}[\lambda_{d}^{*}(K)]\\}\\},$

$\displaystyle I_{AC}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})=$
$\displaystyle\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle,$ (S23)

$\displaystyle I_{Ao}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})=$
$\displaystyle\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d^{\prime}}\\{\partial_{K^{\prime}}^{2}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}$ (S24)
$\displaystyle+2\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}t\lambda_{d^{\prime}}^{*}(K)^{-1}\partial_{K}[\lambda_{d^{\prime}}^{*}(K)]\\}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle(\lambda_{d^{\prime}}^{*}(K)\lambda_{d}(K))^{t}.$

And for the another term,

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}\hat{U}^{t}|x_{i}s_{i}\rangle$
(S25) $\displaystyle=$ $\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm d}K\sum_{d,d^{\prime}}\langle
x_{i}^{\prime}s_{i}^{\prime}|K^{\prime}d^{\prime}\rangle\langle
K^{\prime}d^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}\hat{U}^{t}|Kd\rangle\langle
Kd|x_{i}s_{i}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}\langle x_{i}^{\prime}|K^{\prime}\rangle\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle K|x_{i}\rangle\langle
d|s_{i}\rangle\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle\langle
K^{\prime}|\hat{x}_{i}|K\rangle\langle d^{\prime}|d\rangle\langle
d|\hat{U}_{K}^{t}|d\rangle$ $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}{\rm
d}K^{\prime}{\rm
d}K\sum_{d,d^{\prime}}e^{iK^{\prime}x_{i}^{\prime}}e^{-iKx_{i}}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle d|s_{i}\rangle\langle
d^{\prime}|(\hat{U}_{K}^{\dagger})^{t}|d^{\prime}\rangle\delta^{(1)}(K^{\prime}-K)\langle
d^{\prime}|d\rangle\langle d|\hat{U}_{K}^{t}|d\rangle$ $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d,d^{\prime}}\partial_{K^{\prime}}[e^{i(K^{\prime}x_{i}^{\prime}-Kx_{i})}\langle
s_{i}^{\prime}|d^{\prime}\rangle\langle
d|s_{i}\rangle\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}\lambda_{d}(K)^{t}\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}.$

Similarly,

$\displaystyle\langle
x_{i}^{\prime}s_{i}^{\prime}|(\hat{U}^{\dagger})^{t}\hat{x}_{i}\hat{U}^{t}|x_{i}s_{i}\rangle$
(S26) $\displaystyle=$ $\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d,d^{\prime}}\\{\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}$
$\displaystyle+e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle\partial_{K^{\prime}}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]\\}|_{k^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]\lambda_{d}^{*}(K^{\prime})^{t}$
$\displaystyle+e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle\partial_{K^{\prime}}[\lambda_{d}^{*}(K^{\prime})^{t}]\\}|_{k^{\prime}=K}e^{-ikx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle+$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d}\\{\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}\lambda_{d^{\prime}}^{*}(K)^{t}$
$\displaystyle+e^{iKx_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K)\rangle\langle
d^{\prime}(K)|d(K)\rangle\partial_{K^{\prime}}[\lambda_{d^{\prime}}^{*}(K^{\prime})^{t}]|_{k^{\prime}=K}\\}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}\lambda_{d}^{*}(K)^{t}$
$\displaystyle+e^{iKx_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K)\rangle\partial_{K^{\prime}}[\lambda_{d}^{*}(K^{\prime})^{t}]|_{k^{\prime}=K}\\}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle\lambda_{d}(K)^{t}$ $\displaystyle+$
$\displaystyle\frac{i}{2\pi}\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d}\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{k^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle(\lambda_{d^{\prime}}^{*}(K)\lambda_{d}(K))^{t}$
$\displaystyle=$
$\displaystyle\frac{i}{2\pi}[tI_{B}(x_{j}^{\prime}-x_{j},s_{j}^{\prime},s_{j})+I_{B_{1}}(x_{j}^{\prime},x_{j},s_{j}^{\prime},s_{j})+I_{Bo}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i};t)].$

$I_{B}(x_{i}^{\prime}-x_{i},s_{i}^{\prime},s_{i})=\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}e^{iK(x_{i}^{\prime}-x_{i})}\langle
s_{i}^{\prime}|d(K)\rangle\langle
d(K)|s_{i}\rangle\lambda_{d}^{*}(K)^{-1}\partial_{k}[\lambda_{d}^{*}(K)],$
(S27)

$\displaystyle I_{B_{1}}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})=$
$\displaystyle\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle$ (S28) $\displaystyle=$ $\displaystyle\int_{-\pi}^{\pi}{\rm
d}K\sum_{d}\\{\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K^{\prime})\rangle]|_{K^{\prime}=K}+e^{iKx_{i}^{\prime}}\langle
s_{i}^{\prime}|d(K)\rangle\partial_{K^{\prime}}[\langle
d(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}\\}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle,$

$\displaystyle I_{Bo}(x_{i}^{\prime},x_{i},s_{i}^{\prime},s_{i})=$
$\displaystyle\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d}\partial_{K^{\prime}}[e^{iK^{\prime}x_{i}^{\prime}}\langle
s_{i}^{\prime}|d^{\prime}(K^{\prime})\rangle\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}e^{-iKx_{i}}\langle
d(K)|s_{i}\rangle(\lambda_{d^{\prime}}^{*}(K)\lambda_{d}(K))^{t}$ (S29)
$\displaystyle=$ $\displaystyle\int_{-\pi}^{\pi}{\rm d}K\sum_{d\neq
d}e^{iK(x_{i}^{\prime}-x_{i})}\langle
s_{i}^{\prime}|d^{\prime}(K)\rangle\langle
d(K)|s_{i}\rangle\partial_{K^{\prime}}[\langle
d^{\prime}(K^{\prime})|d(K)\rangle]|_{K^{\prime}=K}(\lambda_{d^{\prime}}^{*}(K)\lambda_{d}(K))^{t}.$

The above expressions apply to any coin operator of SU(2). If we use Hadamard
coin,

$\displaystyle I_{A2}(x^{\prime}-x,s^{\prime},s)=$
$\displaystyle\begin{cases}(\sqrt{2}-2)\pi&x^{\prime}-x=0,s^{\prime}=s;\\\
\sqrt{2}\pi\left(\sqrt{2}-1\right)^{|x^{\prime}-x|}\cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)&x^{\prime}-x\neq
0,s^{\prime}=s;\\\ 0&s^{\prime}\neq s.\end{cases}$ (S30) $\displaystyle=$
$\displaystyle\begin{cases}(1-\sqrt{2})f(0)&x^{\prime}-x=0,s^{\prime}=s;\\\
f(x^{\prime}-x)&x^{\prime}-x\neq 0,s^{\prime}=s;\\\ 0&s^{\prime}\neq
s,\end{cases}$

where $f(x)=\sqrt{2}\pi\left(\sqrt{2}-1\right)^{|x|}\cos\left(\pi x/2\right)$.

$I_{A1}(x^{\prime},x,s^{\prime},s)=\begin{cases}\begin{aligned}
16x^{\prime}[I_{a}(\Delta x-2)+I_{a}(\Delta x+2)]+28x^{\prime}I_{a}(\Delta
x)\\\ +2x^{\prime}[I_{a}(\Delta x+4)+I_{a}(\Delta x-4)]&+8[I_{a}(\Delta
x+2)-I_{a}(\Delta x-2)]\end{aligned}&s^{\prime}=s=\uparrow;\\\ \begin{aligned}
-\\{16x^{\prime}[I_{a}(\Delta x-2)+I_{a}(\Delta x+2)]+28x^{\prime}I_{a}(\Delta
x)\\\ +2x^{\prime}[I_{a}(\Delta x+4)+I_{a}(\Delta x-4)]&+8[I_{a}(\Delta
x+2)-I_{a}(\Delta x-2)]\\}\end{aligned}&s^{\prime}=s=\downarrow;\\\
\begin{aligned} (4x^{\prime}-4)I_{a}(\Delta x+2)+(28x^{\prime}-8)I_{a}(\Delta
x)+(28x^{\prime}-20)I_{a}(\Delta x-2)+4x^{\prime}I_{a}(\Delta
x-4)\end{aligned}&s^{\prime}=\uparrow,s=\downarrow;\\\ \begin{aligned}
(4x^{\prime}+4)I_{a}(\Delta x-2)+(28x^{\prime}+8)I_{a}(\Delta
x)&+(28x^{\prime}+20)I_{a}(\Delta x+2)+4x^{\prime}I_{a}(\Delta
x+4)\end{aligned}&s^{\prime}=\downarrow,s=\uparrow.\end{cases}$ (S31)

$\displaystyle I_{AC}(x^{\prime},x,s^{\prime},s)=$ $\displaystyle
x_{i}^{\prime 2}I_{a}(x^{\prime}-x-4)-4(1+x^{\prime}-3x^{\prime
2})I_{a}(x^{\prime}-x-2)$ (S32) $\displaystyle+(16-24x^{\prime}+38x^{\prime
2})I_{a}(x^{\prime}-x)+(4-4x^{\prime}+12x^{\prime 2})I_{a}(x^{\prime}-x+2)$
$\displaystyle+x^{\prime 2}I_{a}(x^{\prime}-x+4),$

where $\Delta x=x^{\prime}-x$,
$I_{a}(x)=-\frac{\pi}{64}(\sqrt{2}-1)^{|x|}(3\sqrt{2}+2|x|)\cos(\frac{\pi
x}{2})$.

$\displaystyle I_{Ao}(x^{\prime},x,s^{\prime},s)$ (S33) $\displaystyle=$
$\displaystyle\begin{cases}(-1)^{\frac{3}{4}}\sqrt{\frac{\pi}{4t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left((-1+(1+i)x^{\prime})(e^{i\pi(t+x)}+e^{i\pi(3t+x^{\prime})})+((1+i)x^{\prime}-i)(e^{i\pi(2t+x)}+e^{i\pi
x^{\prime}})\right)&s^{\prime}=s=\uparrow\\\
-(-1)^{\frac{3}{4}}\sqrt{\frac{\pi}{4t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left((1+(1+i)x^{\prime})(e^{i\pi(t+x)}+e^{i\pi(3t+x^{\prime})})+((1+i)x^{\prime}+i)(e^{i\pi(2t+x)}+e^{i\pi
x^{\prime}})\right)&s^{\prime}=s=\downarrow\\\
(-1)^{\frac{1}{4}}\sqrt{\frac{\pi}{4t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left(((1+i)x^{\prime}-i)(e^{i\pi(2t+x)}+e^{i\pi
x^{\prime}})+(1-(1+i)x^{\prime})(e^{i\pi(t+x)}+e^{i\pi(3t+x^{\prime})})\right)&s^{\prime}=\uparrow,s=\downarrow\\\
(-1)^{\frac{3}{4}}\sqrt{\frac{\pi}{4t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left((1+(1-i)x^{\prime})(e^{i\pi(2t+x)}+e^{i\pi
x^{\prime}})+i(1+(1+i)x^{\prime})(e^{i\pi(t+x)}+e^{i\pi(3t+x^{\prime})})\right)&s^{\prime}=\downarrow,s=\uparrow\end{cases}+o(\sqrt{\frac{1}{t}}).$

$\displaystyle I_{B}(x^{\prime}-x,s^{\prime},s)=$
$\displaystyle\begin{cases}-iI_{A2}(x^{\prime}-x,s^{\prime},s)&s^{\prime}=s=\uparrow;\\\
iI_{A2}(x^{\prime}-x,s^{\prime},s)&s^{\prime}=s=\downarrow;\\\
\frac{i\pi}{\sqrt{2}}\cos(\frac{\pi(x^{\prime}-x)}{2})[(\sqrt{2}-1)^{|x^{\prime}-x|}-(\sqrt{2}-1)^{|x^{\prime}-x-2|}]&s^{\prime}=\uparrow,s=\downarrow;\\\
\frac{i\pi}{\sqrt{2}}\cos(\frac{\pi(x^{\prime}-x)}{2})[(\sqrt{2}-1)^{|x^{\prime}-x|}-(\sqrt{2}-1)^{|x^{\prime}-x+2|}]&s^{\prime}=\downarrow,s=\uparrow;\end{cases}$
(S34) $\displaystyle=$
$\displaystyle\begin{cases}(-1)^{s}iI_{A2}(x^{\prime}-x,s^{\prime},s)&s^{\prime}=s;\\\
\frac{i}{2}[f(x^{\prime}-x)+f(x^{\prime}-x-(-1)^{s}\cdot 2)]&s^{\prime}\neq
s.\end{cases}$

In the above equation, when $s$ is $\uparrow$, it act as number $1$, and
$\downarrow$ as number $0$, i.e. when $s=\uparrow$, $(-1)^{s}=-1$, and when
$s=\downarrow$, $(-1)^{s}=1$.

$\displaystyle I_{B_{1}}(x^{\prime},x,s^{\prime},s)=$
$\displaystyle\begin{cases}2i\pi
x\delta_{x^{\prime}-x,0}-\frac{i\pi\left(\sqrt{2}-1\right)^{|x^{\prime}-x|}\cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)}{\sqrt{2}}&s^{\prime}=s=\uparrow;\\\
2i\pi
x\delta_{x^{\prime}-x,0}+\frac{i\pi\left(\sqrt{2}-1\right)^{|x^{\prime}-x|}\cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)}{\sqrt{2}}&s^{\prime}=s=\downarrow;\\\
\frac{i\pi\left(\sqrt{2}-1\right)^{|x^{\prime}-x|}\cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)}{\sqrt{2}}&s^{\prime}\neq
s.\end{cases}$ (S35)

$\begin{aligned}
I_{Bo}(x^{\prime},x,s^{\prime},s)=&\begin{cases}\frac{1}{2}i\sqrt{\frac{\pi}{2}}\left(1+e^{i\pi
t}\right)\sqrt{\frac{1}{t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left(-e^{i\pi(t+x^{\prime})}+e^{i\pi(2t+x^{\prime})}+e^{i\pi(t+x)}+e^{i\pi
x^{\prime}}\right)&s^{\prime}=s=\uparrow\\\
-\frac{1}{2}i\sqrt{\frac{\pi}{2}}\left(1+e^{i\pi
t}\right)\sqrt{\frac{1}{t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left(-e^{i\pi(t+x^{\prime})}+e^{i\pi(2t+x^{\prime})}+e^{i\pi(t+x)}+e^{i\pi
x^{\prime}}\right)&s^{\prime}=s=\downarrow\\\
-\frac{1}{2}\sqrt{\frac{\pi}{2}}\left(-1+e^{i\pi
t}\right)\sqrt{\frac{1}{t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left(e^{i\pi(t+x^{\prime})}+e^{i\pi(2t+x^{\prime})}-e^{i\pi(t+x)}+e^{i\pi
x^{\prime}}\right)&s^{\prime}=\uparrow,s=\downarrow\\\
-\frac{1}{2}\sqrt{\frac{\pi}{2}}\left(-1+e^{i\pi
t}\right)\sqrt{\frac{1}{t}}e^{-\frac{1}{2}i\pi(3t+x^{\prime}+x)}\left(e^{i\pi(t+x^{\prime})}+e^{i\pi(2t+x^{\prime})}-e^{i\pi(t+x)}+e^{i\pi
x^{\prime}}\right)&s^{\prime}=\downarrow,s=\uparrow\end{cases}\end{aligned}+o(\sqrt{\frac{1}{t}}).$
(S36)

$I_{Ao}$ and $I_{Bo}$ are integrals with rapidly varying phase when $t$ is
large. The asymptotic formula can be calculate by stationary phase
approximation. Letting $\Sigma$ denote the set of critical points of the
function $f(x)$,

$\int_{R}g(x)e^{itf(x)}{\rm
d}x=\sum_{x_{0}\in\Sigma}g(x_{0})g^{itf(x_{0})+sign(f^{\prime\prime}(x_{0}))\frac{i\pi}{4}}(\frac{2\pi}{t|f^{\prime\prime}(x_{0})|})^{\frac{1}{2}}+o(t^{-\frac{1}{2}}).$
(S37)

To illustrate how we calculate the integrals like $I_{A2}$, $I_{A1}$,
$I_{AC}$, $I_{B}$, $I_{B_{1}}$ using residue theorem, we give an example of
calculating the integral

$I_{A2}(x^{\prime}-x,s^{\prime},s)=\int_{-\pi}^{\pi}{\rm
d}Ke^{iK(x^{\prime}-x)}\sum_{d}\langle s^{\prime}|d(k)\rangle\langle
d(K)|s\rangle\lambda_{d}^{*}(K)^{-2}(\partial_{K}[\lambda_{d}^{*}(K)])^{2}.$
(S38)

When $s^{\prime}=s$,

$\displaystyle I_{A2}(x^{\prime}-x,s^{\prime},s)=$
$\displaystyle\int_{-\pi}^{\pi}[-\frac{2e^{iK(x^{\prime}-x)}\cos^{2}(K)}{cos(2K)+3}]{\rm
d}K$ (S39) $\displaystyle\overset{z=e^{iK}}{=}$
$\displaystyle\oint_{|z|=1}[i\frac{z^{(x^{\prime}-x-1)}(z+z^{-1})^{2}}{((z+z^{-1})^{2}+4)}]dz$
$\displaystyle=$
$\displaystyle\oint_{|z|=1}i[z^{(x^{\prime}-x-1)}-\frac{4z^{(x^{\prime}-x+1)}}{z^{4}+6z^{2}+1}]dz.$

Roots of equation $z^{4}+6z^{2}+1=0$ are
$-i(\sqrt{2}-1),i(\sqrt{2}-1),-i(\sqrt{2}+1),i(\sqrt{2}+1)$, among which
$-i(\sqrt{2}-1),i(\sqrt{2}-1)$ are inside the unit circle on the complex plane
and $-i(\sqrt{2}+1),i(\sqrt{2}+1)$ outside. They are the first-order poles of
the integrand. In addition, $z^{(x^{\prime}-x-1)}$ may produce residue for
$x^{\prime}-x-1=-1$. Then

$\displaystyle I_{A2}(x^{\prime}-x,s^{\prime},s)=$
$\displaystyle\begin{cases}-2\pi
i\\{Res[f(-i\sqrt{3+2\sqrt{2}})]+Res[f(i\sqrt{3+2\sqrt{2}})]\\}&x^{\prime}-x<0;\\\
+2\pi i\\{Res[f(-i\sqrt{3-2\sqrt{2}})]+Res[f(i\sqrt{3-2\sqrt{2}})]\\}-2\pi
iRes[f(\infty)]&x^{\prime}-x=0;\\\ +2\pi
i\\{Res[f(-i\sqrt{3-2\sqrt{2}})]+Res[f(i\sqrt{3-2\sqrt{2}})]\\}&x^{\prime}-x>0;\end{cases}$
(S40) $\displaystyle=$
$\displaystyle\begin{cases}\sqrt{2}\pi\left(3+2\sqrt{2}\right)^{(x^{\prime}-x)/2}cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)&x^{\prime}-x<0;\\\
(\sqrt{2}-2)\pi&x^{\prime}-x=0;\\\
\sqrt{2}\pi\left(3-2\sqrt{2}\right)^{(x^{\prime}-x)/2}cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)&x^{\prime}-x>0;\end{cases}$
$\displaystyle=$ $\displaystyle\begin{cases}(\sqrt{2}-2)\pi&x^{\prime}-x=0;\\\
\sqrt{2}\pi\left(3-2\sqrt{2}\right)^{|x^{\prime}-x|/2}cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)&x^{\prime}-x\text{$\neq$}0.\end{cases}$

The notation $Res[f(a)]$ represents the residue of function $f(z)$ at point
$a$. When $s^{\prime}\neq s$,

$I_{A2}(x^{\prime}-x,s^{\prime},s)=0.$ (S41)

To sum up,

$\displaystyle I_{A2}(x^{\prime}-x,s^{\prime},s)=$
$\displaystyle\begin{cases}(\sqrt{2}-2)\pi&x^{\prime}-x=0,s^{\prime}=s;\\\
\sqrt{2}\pi\left(\sqrt{2}-1\right)^{|x^{\prime}-x|}\cos\left(\frac{\pi(x^{\prime}-x)}{2}\right)&x^{\prime}-x\neq
0,s^{\prime}=s;\\\ 0&s^{\prime}\neq s.\end{cases}$ (S42)

###  A4 Optimization problem of coefficient

The coefficient $c_{2}$ of $t^{2}$ in the expression of the average distance
after $t$ steps of evolution
$\langle\Psi|(\hat{U}^{\dagger})^{t}\hat{D}\hat{U}^{t}|\Psi\rangle$ can be
calculated from the combination of inner products of $\langle
x_{1}s_{1}^{\prime}|\cdots\langle
x_{n}s_{n}^{\prime}|(\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}|x_{1}s_{1}\rangle\cdots|x_{n}s_{n}\rangle$,
see Eq. (S43). For simplicity, we regard $\uparrow$ as number $1$ and
$\downarrow$ as number $0$ and denote the spin of the $i$th particle as
$|s_{i}\rangle$. Then the vector $|s_{i}\rangle$ takes on the value of
$|1\rangle$ or $|0\rangle$, and the combination of the spin notation of these
particles can be treated as a binary number, that is, $|s_{1}s_{2}\cdots
s_{n}\rangle=|2^{n-1}\times s_{1}+2^{n-2}\times
s_{2}+\cdots+2s_{n-1}+s_{n}\rangle=|\xi\rangle$, the real value of which is
$\xi$. We can rewrite the state vectors as
$|x_{1}\rangle|x_{2}\rangle\cdots|x_{n}\rangle|\xi\rangle$.

The average distance is given by

$\displaystyle\langle\Psi|(\hat{U}^{\dagger})\hat{D}\hat{U}^{t}|\Psi\rangle=$
$\displaystyle[\langle x_{1}|\langle x_{2}|\cdots\langle
x_{n}|\sum_{\xi^{\prime}=0}^{2^{n}-1}a_{\xi^{\prime}}^{*}\langle\xi^{\prime}|](\hat{U}_{n}^{\dagger})^{t}\hat{D}\hat{U}_{n}^{t}[|x_{1}\rangle|x_{2}\rangle\cdots|x_{n}\rangle\sum_{\xi=0}^{2^{n}-1}a_{\xi}|\xi\rangle],$
(S43)

in which the coefficient of $t^{2}$ is

$\displaystyle c_{2}=$
$\displaystyle\sum_{\xi^{\prime}\xi}a_{\xi^{\prime}}^{*}a_{\xi}m_{\xi^{\prime}\xi}$
(S44) $\displaystyle=$
$\displaystyle[a_{0}^{*},a_{1}^{*},\cdots,a_{2^{n-1}}^{*}]M\left[\begin{array}[]{c}a_{0}\\\
a_{1}\\\ \vdots\\\ a_{2^{n-1}}\end{array}\right].$

The matrix elements of $M$ are given by

$\displaystyle M_{\xi^{\prime}\xi}=$
$\displaystyle-\frac{n-1}{n}\frac{1}{2\pi}\sum_{i}I_{A2}(0,s_{i}^{\prime},s_{i})\prod_{\gamma=1,\gamma\neq
i}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}$ (S45)
$\displaystyle+\frac{1}{n}(\frac{1}{2\pi})^{2}\sum_{j,k,j\neq
k}I_{B}(0,s_{j}^{\prime},s_{j})I_{B}(0,s_{k}^{\prime},s_{k})\prod_{\gamma=1,\gamma\neq
j,k}\delta_{x_{\gamma}^{\prime}s_{\gamma}^{\prime},x_{\gamma}s_{\gamma}}.$

In $I_{A2}(0,s_{i}^{\prime},s_{i})$, if $s_{i}^{\prime}\neq s_{i}$,
$I_{A2}(0,s_{i}^{\prime},s_{i})=0$. Thus the first part in
$m_{\xi^{\prime}\xi}$ is zero if $\xi^{\prime}\neq\xi$. If $\xi^{\prime}=\xi$,
it takes on the value of $(n-1)(1-\frac{1}{\sqrt{2}})$.

For the second part in $m_{\xi^{\prime}\xi}$, integrals
$I_{B}(0,\uparrow,\uparrow)$, $I_{B}(0,\uparrow,\downarrow)$ and
$I_{B}(0,\downarrow,\uparrow)$ are all $i(2-\sqrt{2})\pi$. While
$I_{B}(0,\downarrow,\downarrow)=-i(2-\sqrt{2})\pi$. So the values of
$I_{B}(0,s_{j}^{\prime},s_{j})I_{B}(0,s_{k}^{\prime},s_{k})$ corresponding to
different combinations of $s_{j}^{\prime}$,$s_{j}$,$s_{k}^{\prime}$ and
$s_{k}$ are listed below:

$(s_{j}^{\prime},s_{j},s_{k}^{\prime},s_{k})$ | value
---|---
$\begin{array}[]{c}(\uparrow,\uparrow,\uparrow,\uparrow),(\uparrow,\uparrow,\uparrow,\downarrow),(\uparrow,\uparrow,\downarrow,\uparrow),(\uparrow,\uparrow,\downarrow,\downarrow),(\uparrow,\downarrow,\uparrow,\uparrow),\\\ (\uparrow,\downarrow,\downarrow,\uparrow),(\downarrow,\uparrow,\uparrow,\uparrow),(\downarrow,\uparrow,\uparrow,\downarrow),(\downarrow,\downarrow,\uparrow,\uparrow),(\downarrow,\downarrow,\downarrow,\downarrow)\end{array}$ | $-(2-\sqrt{2})^{2}\pi^{2}$
$(\uparrow,\downarrow,\uparrow,\downarrow)$,$(\uparrow,\downarrow,\downarrow,\downarrow)$,$(\downarrow,\uparrow,\downarrow,\uparrow)$,$(\downarrow,\uparrow,\downarrow,\downarrow)$,$(\downarrow,\downarrow,\uparrow,\downarrow)$,$(\downarrow,\downarrow,\downarrow,\uparrow)$ | $(2-\sqrt{2})^{2}\pi^{2}$

When $\xi^{\prime}=\xi$, i.e. the left statevector is the same as the right
statevector, denote the number of particles that are spinning up as
$n_{\uparrow}$, spinning down as $n_{\downarrow}$, then
$n_{\uparrow}+n_{\downarrow}=n$.

$\displaystyle m_{\xi^{\prime}\xi}=$
$\displaystyle(n-1)\left(1-\frac{1}{\sqrt{2}}\right)+\frac{1}{n}\sum_{j,k,j\neq
k}(-1)^{\delta_{s_{j},s_{k}}}\left(1-\frac{1}{\sqrt{2}}\right)^{2}$ (S46)
$\displaystyle=$
$\displaystyle(n-1)\left(1-\frac{1}{\sqrt{2}}\right)+\frac{2}{n}[-(C_{n_{\uparrow}}^{2}+C_{n_{\downarrow}}^{2})+n_{\uparrow}n_{\downarrow}]\left(1-\frac{1}{\sqrt{2}}\right)^{2}$
$\displaystyle=$
$\displaystyle(n-1)\left(1-\frac{1}{\sqrt{2}}\right)-\frac{1}{n}[((n_{\uparrow}-n_{\downarrow})^{2}-n)]\left(1-\frac{1}{\sqrt{2}}\right)^{2}.$

When $d(\xi^{\prime},\xi)>2$, $m_{\xi^{\prime}\xi}=0$, and when
$d(\xi^{\prime},\xi)=1$, there is an mismatch in one of the spins pairs.
Denote the number of spinning up pairs in the matched $n-1$ pairs as
$n_{\uparrow}^{\prime}$ and spinning down as $n_{\downarrow}^{\prime}$, the
matrix elements are

$m_{\xi^{\prime}\xi}=\frac{2}{n}[-n_{\uparrow}^{\prime}+n_{\downarrow}^{\prime}]\left(1-\frac{1}{\sqrt{2}}\right)^{2}.$
(S47)

When $d(\xi^{\prime},\xi)=2$,

$m_{\xi^{\prime}\xi}=-\frac{2}{n}\left(1-\frac{1}{\sqrt{2}}\right)^{2}.$ (S48)

If we rewrite the subscript of matrix $M$ as $i,j,k$,

$\displaystyle m_{ii}=$
$\displaystyle(n-1)\left(1-\frac{1}{\sqrt{2}}\right)+\frac{1}{n}[n-(2W(i)-n)^{2}]\left(1-\frac{1}{\sqrt{2}}\right)^{2},$
(S49)

$\displaystyle M_{jk}=$
$\displaystyle\frac{2}{n}\left(1-\frac{1}{\sqrt{2}}\right)^{2}[\delta_{d(j,k),1}(n-1-2\min\\{W(j),W(k)\\})-\delta_{d(j,k),2}],(j\neq
k),$ (S50)

where $i,j,k=0,1,2,\cdots,2^{n}-1$, $W(i)$ represents the Hamming weight of
integer $i$, i.e. the number of “$1$”s in the $n$ bit binary form of integer
$i$. $d(j,k)$ represents the Hamming distance between integer $j$ and $k$,
which is the number of bit positions where the corresponding two bits in the
binary form of $i$ and $j$ are different. $\delta$ is the Kronecker delta.

The maximum and minimum values of $c_{2}$ as a function of the initial coin
state are the maximum and minimum eigenvalues of $M$ respectively.

###  A5 Solving eigenvalues and eigenvectors

Given that the number of particles is $n$,

$\begin{cases}m_{ii}=&(n-1)\left(1-\frac{1}{\sqrt{2}}\right)+\frac{1}{n}[n-(2W(i)-n)^{2}]\left(1-\frac{1}{\sqrt{2}}\right)^{2}\\\
m_{jk}=&\frac{2}{n}\left(1-\frac{1}{\sqrt{2}}\right)^{2}[\delta_{d(j,k),1}(n-1-2\min\\{W(j),W(k)\\})-\delta_{d(j,k),2}],(j\neq
k),\end{cases}$ (S51)

where $i,j,k=0,1,\cdots,2^{n}-1$. Define

$B_{n}=\left(1-\frac{1}{\sqrt{2}}\right)^{-2}\frac{n}{2}M_{n}-\frac{n(n-1)}{\sqrt{2}}I,$
(S52)

the matrix elements of which are

$\begin{cases}b_{ii}=n(n-1)+\frac{1}{2}[n-(2W(i)-n)^{2}]\\\
b_{jk}=\delta_{d(j,k),1}(n-1-2\min\\{W(j),W(k)\\})-\delta_{d(j,k),2},(j\neq
k).\end{cases}$ (S53)

When $n\geq 2$, define the 2 by 2 matrix $P_{ij}=(-1)^{c(i,j)}D^{d(i,j)}$,
$2|i$, $2|j$ and $i,j\in\\{0,1,\cdots,2^{n}-1\\}$. $c(i,j)$ represents the
number of bits in the binary form on which $i$ and $j$ have a common $1$.

######  Theorem 1.

The submatrices $P_{ij}$ make up of the matrix $P_{n}$, such that
$P_{n}^{-1}B_{n}P_{n}$ is diagonal matrix $(\mu_{k}I_{2}),$ where $k$ in even
number and $I_{2}$ is $2\times 2$ identity matrix.

###### Proof.

When $n=2$,

$P_{n}^{-1}B_{n}P_{n}=\left[\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&0&4&0\\\ 0&0&0&4\end{array}\right].$ (S54)

Assume that the theorem 1 holds true for $n$. Then for $n+1$, let
$B_{n}=(B_{ij})$ is the $2\times 2$ submatrix partition and similarly
$B_{n+1}=(B_{ij}^{\prime}),$ $P_{n}=(P_{ij})$, $P_{n+1}=(P_{ij}^{\prime})$,
where $i,j$ are even numbers. It’s to be proved that $\forall i,k<2^{n+1}$,
and that $i,k$ are odd numbers, $\exists\mu$ s.t.

$\sum_{\begin{subarray}{c}0\leq j<2^{n+1}\\\
2|j\end{subarray}}B_{ij}^{\prime}P_{jk}^{\prime}=\mu^{\prime}P_{ik}^{\prime}.$
(S55)

Let’s proof it in the case of $i<2^{n},k<2^{n}$.

Let $j<2^{n}$,

$B_{ij}^{\prime}=B_{ij}+\left[\begin{array}[]{cc}\delta_{d(i,j),1}&\delta_{d(i,j+1),1}\\\
\delta_{d(i,j+1),1}&\delta_{d(i+1,j+1),1}\end{array}\right]=B_{ij}+\delta_{d(i,j),1}I_{2}(i\neq
j),$ (S56)

$B_{ii}^{\prime}=B_{ii}+\left[\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right],$ (S57)

$B_{i,j+2^{n}}^{\prime}=\left[\begin{array}[]{cc}-\delta_{d(i,j+2^{n}),2}&-\delta_{d(i,j+2^{n}+1),2}\\\
-\delta_{d(i+1,j+2^{n}),2}&-\delta_{d(i+1,j+2^{n}+1),2}\end{array}\right]=-\delta_{d(i,j),1}I_{2}(i\neq
j),$ (S58)

$B_{i,i+2^{n}}^{\prime}=\left[\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right].$ (S59)

Hence,

$\displaystyle\sum_{\begin{subarray}{c}0\leq j<2^{n+1}\\\
2|j\end{subarray}}B_{ij}^{\prime}P_{jk}^{\prime}=$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}(B_{ij}+\delta_{d(i,j),1}I)P_{jk}+[B_{ii}+\left(\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right)]P_{ik}$ (S60) $\displaystyle+$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}[-\delta_{d(i,j),1}DP_{jk}+\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)DP_{ik}].$

Considering $\exists\mu$, s.t. $\sum_{\begin{subarray}{c}j<2^{n}\\\
2|j\end{subarray}}B_{ij}P_{jk}=\mu P_{ik}$, then

$\displaystyle\sum_{\begin{subarray}{c}0\leq j<2^{n+1}\\\
2|j\end{subarray}}B_{ij}^{\prime}P_{jk}^{\prime}=$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}\delta_{d(i,j),1}P_{jk}(I-D)+[\left(\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right)+\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)D]P_{ik}+\mu P_{ik}.$ (S61)

Considering $\sum_{\begin{subarray}{c}j<2^{n}\\\
2|j,d(i,j)=1\end{subarray}}P_{jk}P_{ik}^{-1}=rI+yD$, where
$r=2\cdot\\#\\{(i,k)=(0,1)\\}-2\cdot\\#\\{(i,k)=(1,0)\\}$, $y=\\#\\{(i,k)=(0\
or\ 1,0)\\}-\\#\\{(i,k)=(0\ or\ 1,1)\\}$. $\\#\\{(i,k)=(0,1)\\}$ represents
the number of $j$ for $0\leq j<2^{n},2|j$ s.t. that the binary form of $i$ is
0 and $j$ is 1 on the bit where the binary form of $i$ differs from $j$. One
can deduce that $r=2(W(k)-W(i))$, $y=n-2W(k)-1$. Taking
$\mu^{\prime}=\mu+4W(k)$, it can be verified that

$\displaystyle(rI+yD)(I-D)P_{ik}+[\left(\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right)+\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)D]P_{ik}+\mu P_{ik}$ (S62) $\displaystyle=$
$\displaystyle[(rI+yD)(I-D)+(n+2W(i)-1)I+(n-2W(i)-1)D+\mu
I]P_{ik}=\mu^{\prime}P_{ik}=\mu^{\prime}P_{ik}^{\prime}.$

Notice that $P_{ik}$ commutes with $D$.

When $i<2^{n}$, $k\geq 2^{n}$,

$\displaystyle\sum_{\begin{subarray}{c}\text{0$\leq$j<$2^{n+1}$}\\\
2|J\end{subarray}}B_{ij}^{\prime}P_{jk}^{\prime}=$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}(B_{ij}+\delta_{d(i,j),1}I)DP_{jk^{\prime}}+[B_{ii}+\left(\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right)]DP_{ik^{\prime}}$ (S63) $\displaystyle+$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}\delta_{d(i,j),1}P_{jk^{\prime}}-\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)P_{ik^{\prime}}.$

Hence,

$\displaystyle\sum_{\begin{subarray}{c}\text{0$\leq$j<$2^{n+1}$}\\\
2|J\end{subarray}}B_{ij}^{\prime}P_{jk}^{\prime}=$
$\displaystyle(rI+yD)(D+I)P_{ik^{\prime}}+[\left(\begin{array}[]{cc}n+2W(i)&1\\\
1&n+2+2W(i)\end{array}\right)D-\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)]P_{ik^{\prime}}+\mu DP_{ik^{\prime}}$ (S64)
$\displaystyle=$
$\displaystyle[4(n-W(k^{\prime}))+\mu]P_{ik}=\mu^{\prime}P_{ik}.$

When $i\geq 2^{n}$, $k<2^{n}$, let $j<2^{n}$,

$B_{i+2^{n},j}^{\prime}=B_{j,i+2^{n}}^{\prime}=\left(\begin{array}[]{cc}-\delta_{d(j,i+2^{n}),2}&-\delta_{d(j,i+2^{n}+1),2}\\\
-\delta_{d(j+1,i+2^{n}),2}&-\delta_{d(j+1,i+2^{n}+1),2}\end{array}\right)=-\delta_{d(i,j),1}(i\neq
j),$ (S65)

$B_{i+2^{n},1}^{\prime}=\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right),$ (S66)

$B_{i+2^{n},j+2^{n}}^{\prime}=\left(\begin{array}[]{cc}-2\delta_{d(i,j),1}&-2\delta_{d(i,j+1),1}\\\
-2\delta_{d(i+1,j),1}&-2\delta_{d(i+1,j+1),1}\end{array}\right)+B_{ij}^{\prime}=B_{ij}-\delta_{d(i,j),1}I_{2}(i\neq
j),$ (S67)

$B_{i+2^{n},i+2^{n}}^{\prime}=\left(\begin{array}[]{cc}-4W(i)+2n&-2\\\
-2&-4-4W(i)+2n\end{array}\right)+B_{ii}^{\prime}=B_{ii}+\left(\begin{array}[]{cc}3n-2W(i)&-1\\\
-1&3n-2W(i)-2\end{array}\right),$ (S68)

$\displaystyle\sum_{\begin{subarray}{c}0\leq j<2^{n+!}\\\
2|j\end{subarray}}B_{i+2^{n},j}^{\prime}P_{jk}^{\prime}=$
$\displaystyle\sum_{\begin{subarray}{c}j<2\\\ 2|j\neq
i\end{subarray}}(-\delta_{d(i,j),1})P_{jk}+\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)P_{ik}$ (S69) $\displaystyle+$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}(B_{ij}-\delta_{d(i,j),1})DP_{jk}+[B_{ii}+\left(\begin{array}[]{cc}3n-2W(i)&-1\\\
-1&3n-2W(i)-2\end{array}\right)]DP_{ik}$ $\displaystyle=$
$\displaystyle[4W(k)+\mu]P_{i+2^{n},k}^{\prime}.$

When $i^{\prime}\geq 2^{n}$, $k^{\prime}\geq 2^{n}$,

$\displaystyle\sum_{\begin{subarray}{c}0\leq j<2^{n+1}\\\
2|j\end{subarray}}B_{i+2^{n},j}^{\prime}P_{j,k+2^{n}}=$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\
2|j\end{subarray}}(-\delta_{d(i,j),1})DP_{jk}+\left(\begin{array}[]{cc}n-2W(i)&-1\\\
-1&n-2-2W(i)\end{array}\right)DP_{ik}$ (S70) $\displaystyle-$
$\displaystyle\sum_{\begin{subarray}{c}j<2^{n}\\\ 2|j\neq
i\end{subarray}}(B_{ij}-\delta_{d(i,j),1})P_{jk}-[B_{ii}+\left(\begin{array}[]{cc}3n-2W(i)&-1\\\
-1&3n-2W(i)-2\end{array}\right)]P_{ik}$ $\displaystyle=$
$\displaystyle[4n-4W(k)+\mu]P_{i+2^{n},k+2^{n}}^{\prime}.$

∎

According to theorem 1, the eigenvalues corresponding to the column subscript
$k$, is

$\mu_{k}^{n}=\begin{cases}\mu_{k}^{n-1}+4W(k)&k<2^{n-1};\\\
4(n-1)-4W(k-2^{n-1})+\mu_{k-2^{n-1}}^{n-1}&k\geq 2^{n-1}.\end{cases}$ (S71)

$\mu_{0}^{2}=0$, $\mu_{2}^{2}=4$,

$\Rightarrow\mu_{k}=4W(k)W(2^{n}-1-k).$ (S72)

And each $\mu_{k}$($0\leq k<2^{n}$, $2|k$) correspond to 2 eigenvectors. So
eigenvalues of $M_{n}$ are
$\left(1-\frac{1}{\sqrt{2}}\right)^{2}[\frac{2\mu_{k}}{n}+(n-1)\sqrt{2}]$.

###  A6 Proof of partial exchange symmetry

Using the relation $D^{2}=2D+I$, it can be proved that

$D^{\alpha}=\left[\begin{array}[]{cc}F(\alpha-1)&F(\alpha)\\\
F(\alpha)&F(\alpha+1)\end{array}\right].$ (S73)

$\alpha\in\mathbb{Z},\alpha\geq-1$ and

$F(\alpha)=\frac{(1+\sqrt{2})^{\alpha}-(1-\sqrt{2})^{\alpha}}{2\sqrt{2}}.$
(S74)

So

$P_{ij}=(-1)^{c(i,j)}D^{d(i,j)}=(-1)^{c(i,j)}\left[\begin{array}[]{cc}F(d(i,j)-1)&F(d(i,j))\\\
F(d(i,j))&F(d(i,j)+1)\end{array}\right].$ (S75)

To analyze the exchange symmetry, we look at he binary form of subscript $i$
and $j$, that $i,j\in\\{0,1,\cdots,2^{n}-1\\}$. Each column of the matrix
$(P_{ij})$ is an unnormalized eigenvector of matrix $M$. In the matrix
representation of a state vector, $i$ is the corresponding row subscript of
the matrix elements and the binary form of $i$ correspond to a product state.
Each bit in the binary form of $i$ corresponds to a particle. We can say that
the eigenstate possess exchange symmetry between two particles if the state
vector is unchanged when the basis of the state vector are changed as we swap
corresponding two bits. We call this exchange symmetry in a subset of
particles as partial exchange symmetry.

######  Lemma 1.

If swap in the binary form of $i$ happens between two of the first $n-1$ bits
where the corresponding bits in $j$ are both $0$ or $1$, $c(i,j)$ and $d(i,j)$
are unchanged.

######  Lemma 2.

If swap in the binary form of $i$ happens between two of the first $n-1$ bits
where the corresponding bits in $j$ are $1$ and $0$ respectively, there is
always an $i$ where the $c(i,j)$ or $d(i,j)$ is changed.

######  Lemma 3.

When we swap particles corresponding to the first $n-1$ bits in $i$, iff. we
swap the particles in the way described in lemma 1, the state is unchanged.

###### Proof.

The adequacy is evident. The following is the proof of necessity. If one does
not perform the swap as described in lemma 1, i.e. performing a swap described
in lemma 2, there exist an $i$, s.t. $c(i,j)$ and $d(i,j$) are changed, so the
sign of the elements are changed. And since every column in $P_{n}$ can’t be
zero, the eigenvector is changed.∎

######  Lemma 4.

When swapping the last bit in $i$ with one of the first $n-1$ bits, iff. the
last bit $s_{n}$ swaps with the bit where the corresponding bit in $j$ is $0$,
the state is unchanged.

###### Proof.

The proof of the sufficiency: If $s_{n}=0$, $c(i,j)$ is unchanged. When the
bit in $i$ which swaps with the last bit is $0$, then the swap does not change
$i$, i.e. the bases are unchanged, so the state vector is unchanged. But when
the bit is $1$, $d(i,j)$ decreases by $1$ (if we don’t consider the last bit),
and the corresponding submatrix is

$\mathfrak{P}_{ij}=(-1)^{c(i,j)}D^{d(i,j)-1}=(-1)^{c(i,j)}\left[\begin{array}[]{cc}F(d(i,j)-2)&F(d(i,j)-1)\\\
F(d(i,j)-1)&F(d(i,j))\end{array}\right].$ (S76)

While in the original matrix,

$P_{ij}=(-1)^{c(i,j)}\left[\begin{array}[]{cc}F(d(i,j)-1)&F(d(i,j))\\\
F(d(i,j))&F(d(i,j)+1)\end{array}\right].$ (S77)

We can see that the first row of $P_{ij}$ is identical to the second row of
$\mathfrak{P}_{ij}$, so if we modify a basis where the last bits are
$s_{n}=0$, it turns out that it maps into another basis where the
corresponding elements in $P_{ij}$ are the same.

If $s_{n}=1$, the bases are unchanged if the bit which swaps with $s_{n}$ is
$1$. If the bit is 0, then $d(i,j)$ increases with $1$ (if we don’t consider
the last bit). That

$\mathfrak{P}_{ij}=(-1)^{c(i,j)}\left[\begin{array}[]{cc}F(d(i,j))&F(d(i,j)+1)\\\
F(d(i,j)+1)&F(d(i,j)+2)\end{array}\right].$ (S78)

The second row of $P_{ij}$ is the same as the first row of
$\mathfrak{P}_{ij}$.

The proof of the necessity: if the last bit $s_{n}$ of $i$ swaps with one of
the first $n-1$ bits in $i$ where the corresponding bit in $j$ is $1$,
$c(i,j)$ decreases by $1$ and $d(i,j)$ increases by $1$ when $s_{n}=0$, and
reversely $c(i,j)$ increases by $1$ and $d(i,j)$ decreases by $1$ when
$s_{n}=1$. If we ignore the change in sign factor $(-1)^{c(i,j)}$, the
corresponding rows are the same. However, changing $c(i,j)$ by $\pm 1$ will
induce a factor $(-1)$, and because there always exist $i$ that the elements
on the basis is not $0$, the state would change by a sign $(-1)$ on some of
the basis.∎

To sum up, according to lemma 1$\sim$4, iff. the swap happens among the first
$n-1$ particles where the corresponding bits in $j$ are all $1$ or $0$, or
between the last particle and one of the first $n-1$ particles where the
corresponding bit in $j$ is 0, the state vector is unchanged. So we can group
the particles where the first $n-1$ bits in $j$ are $1$ into a set, and the
other particles into another. Then we connect each pair of particles with an
undirected edge if the state is unchanged if they are swapped. Each of the two
sets (when $j\neq 0$) of particles consists of a complete graph. If we denote
$n_{\uparrow}$ and $n_{\downarrow}$ as the number of “$1$”s and “$0$”s in the
binary form of $j$ respectively, then
$p_{k}=\frac{n_{\uparrow}(n_{\uparrow}-1)}{2}+\frac{n_{\downarrow}(n_{\downarrow}-1)}{2}$,
$n_{\uparrow}+n_{\downarrow}=n$, then the two disconnected complete graphs
have $n_{\uparrow}$ and $n_{\downarrow}$ vertices respectively. When $j=0$,
all of the bits in $j$ are $0$, so there is only one complete graph and all of
the particles are vertices of this graph.

###  A7 Entanglement between complete graphs

Since the particles in each of the two complete graphs possess exchange
symmetry. For simplicity, the eigenstates can be represent by the direct
product of two Fock states. When $k$ is even number,

$\displaystyle|\psi_{k}\rangle=$ $\displaystyle
C\sum_{\begin{subarray}{c}0\leq w_{1}\leq W(k)\\\ 0\leq w_{0}\leq
n-W(k)\end{subarray}}(-1)^{w_{1}}F(W(k)-w_{1}+w_{0}-1)\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}$
(S79)
$\displaystyle\cdot|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle\otimes|1_{w_{1}}0_{W(k)-w_{1}}\rangle$
$\displaystyle=$ $\displaystyle
C[\frac{(1+\sqrt{2})^{w(k)-1}}{2\sqrt{2}}\sum_{0\leq w_{0}\leq
n-W(k)}(1+\sqrt{2})^{w_{0}}\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle$
$\displaystyle\otimes\sum_{0\leq w_{1}\leq
W(k)}(-1)^{w_{1}}(1+\sqrt{2})^{-w_{1}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}|1_{w_{1}}0_{W(k)-w_{1}}\rangle$
$\displaystyle-\frac{(1-\sqrt{2})^{w(k)-1}}{2\sqrt{2}}\sum_{0\leq w_{0}\leq
n-W(k)}(1-\sqrt{2})^{w_{0}}\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle$
$\displaystyle\otimes\sum_{0\leq w_{1}\leq
W(k)}(-1)^{w_{1}}(1-\sqrt{2})^{-w_{1}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}|1_{w_{1}}0_{W(k)-w_{1}}\rangle].$

$|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle$($0\leq w_{0}\leq n-W(k)$) and
$|1_{w_{1}}0_{W(k)-w_{1}}\rangle$($0\leq w_{1}\leq W(k)$) are two sets of
orthonormal vectors. $C$ is the normalization factor.

Denote

$|\tilde{\varphi}_{A}^{1}\rangle\equiv\sum_{0\leq w_{0}\leq
n-W(k)}(1+\sqrt{2})^{w_{0}}\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle,$
(S80) $|\tilde{\varphi}_{B}^{1}\rangle\equiv\sum_{0\leq w_{1}\leq
W(k)}(-1)^{w_{1}}(1+\sqrt{2})^{-w_{1}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}|1_{w_{1}}0_{W(k)-w_{1}}\rangle,$
(S81) $|\tilde{\varphi}_{A}^{2}\rangle\equiv\sum_{0\leq w_{0}\leq
n-W(k)}(1-\sqrt{2})^{w_{0}}\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle,$
(S82) $|\tilde{\varphi}_{B}^{2}\rangle\equiv\sum_{0\leq w_{1}\leq
W(k)}(-1)^{w_{1}}(1-\sqrt{2})^{-w_{1}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}|1_{w_{1}}0_{W(k)-w_{1}}\rangle.$
(S83)

and normalize them using the following equations,

$\sum_{0\leq w_{0}\leq
n-W(k)}(1+\sqrt{2})^{2w_{0}}C_{n-W(k)}^{w_{0}}=(4+2\sqrt{2})^{n-W(k)},$ (S84)

$\sum_{0\leq w_{0}\leq
n-W(k)}(1+\sqrt{2})^{-2w_{0}}C_{n-W(k)}^{w_{0}}=(4-2\sqrt{2})^{W(k)},$ (S85)

and $(1-\sqrt{2})^{-1}=-(1+\sqrt{2})$. We remove the tilde notation for the
normalized vectors, the eigenstate turns into

$\displaystyle|\psi_{k}\rangle=$ $\displaystyle
C[\frac{(1+\sqrt{2})^{W(k)-1}}{2\sqrt{2}}\sqrt{(4+2\sqrt{2})^{n-W(k)}}|\varphi_{A}^{1}\rangle\otimes\sqrt{(4-2\sqrt{2})^{W(k)}}|\varphi_{B}^{1}\rangle$
(S86)
$\displaystyle-\frac{(1-\sqrt{2})^{W(k)-1}}{2\sqrt{2}}\sqrt{(4-2\sqrt{2})^{n-W(k)}}|\varphi_{A}^{2}\rangle\otimes\sqrt{(4+2\sqrt{2})^{W(k)}}|\varphi_{B}^{2}\rangle].$

The density matrix of part $A$ is

$\displaystyle\rho_{A}=\rho_{\uparrow}=$ $\displaystyle
tr_{B}(|\psi_{k}\rangle\langle\psi_{k}|)$ (S87) $\displaystyle=$
$\displaystyle|\varphi_{A}^{1}\rangle\langle\varphi_{A}^{1}|C^{2}\frac{(1+\sqrt{2})^{2(W(k)-1)}}{8}(4+2\sqrt{2})^{n-W(k)}(4-2\sqrt{2})^{W(k)}$
$\displaystyle+$
$\displaystyle|\varphi_{A}^{2}\rangle\langle\varphi_{A}^{2}|C^{2}\frac{(1-\sqrt{2})^{2(W(k)-1)}}{8}(4-2\sqrt{2})^{n-W(k)}(4+2\sqrt{2})^{W(k)},$

where $\langle\varphi_{A}^{1}|\varphi_{A}^{2}\rangle=\sum_{0\leq w_{0}\leq
n-W(k)}(-1)^{w_{0}}C_{n-W(k)}^{w_{0}}=0$.
$C^{2}=\\{2^{-3+n}(\sqrt{2}-1)^{2}[(2+\sqrt{2})^{n}+(2-\sqrt{2})^{n}(17+12\sqrt{2})]\\}^{-1}$.
Finally, we obtain

$\begin{aligned}
\rho_{\uparrow}=\rho_{A}=&|\varphi_{A}^{1}\rangle\langle\varphi_{A}^{1}|[1+(3-2\sqrt{2})^{n}(17+12\sqrt{2})]^{-1}+|\varphi_{A}^{2}\rangle\langle\varphi_{A}^{2}|[1+(3+2\sqrt{2})^{n}/(17+12\sqrt{2})]^{-1}\\\
=&|\varphi_{A}^{1}\rangle\langle\varphi_{A}^{1}|[1+(3-2\sqrt{2})^{n-2}]^{-1}+|\varphi_{A}^{2}\rangle\langle\varphi_{A}^{2}|[1+(3+2\sqrt{2})^{n-2}]^{-1}.\end{aligned}.$
(S88)

So the von Neumann entropy is given by

$S(\rho_{\uparrow})=-\sum_{i=1,2}\nu_{i}{\rm\log}\nu_{i},$ (S89)

in which $\nu_{1}=[1+(3-2\sqrt{2})^{n+2}]^{-1}$,
$\nu_{2}=[1+(3+2\sqrt{2})^{n+2}]^{-1}$.

Similarly, when $k$ is odd number,

$\displaystyle|\psi_{k}\rangle=$ $\displaystyle
C\sum_{\begin{subarray}{c}0\leq w_{1}\leq W(k)\\\ 0\leq w_{0}\leq
n-W(k)\end{subarray}}(-1)^{w_{1}}F(W(k)-w_{1}+w_{0})\sqrt{\frac{(n-W(k))!}{w_{0}!(n-W(k)-w_{0})!}}\sqrt{\frac{W(k)!}{w_{1}!(W(k)-w_{1})!}}$
(S90)
$\displaystyle\cdot|1_{w_{0}}0_{n-W(k)-w_{0}}\rangle\otimes|1_{w_{1}}0_{W(k)-w_{1}}\rangle,$

where
$C^{2}=\\{2^{-3+n}(\sqrt{2}-1)^{2}[(2+\sqrt{2})^{n}+(2-\sqrt{2})^{n}(17+12\sqrt{2})]\\}^{-1}$.

$\displaystyle\rho_{A}=\rho_{\uparrow}=$ $\displaystyle
tr_{B}(|\psi_{k}\rangle\langle\psi_{k}|)$ (S91) $\displaystyle=$
$\displaystyle|\varphi_{A}^{1}\rangle\langle\varphi_{A}^{1}|C^{2}\frac{(1+\sqrt{2})^{2W(k)}}{8}(4+2\sqrt{2})^{n-W(k)}(4-2\sqrt{2})^{W(k)}$
$\displaystyle+$
$\displaystyle|\varphi_{A}^{2}\rangle\langle\varphi_{A}^{2}|C^{2}\frac{(1-\sqrt{2})^{2W(k)}}{8}(4-2\sqrt{2})^{n-W(k)}(4+2\sqrt{2})^{W(k)}$
$\displaystyle=$
$\displaystyle|\varphi_{A}^{1}\rangle\langle\varphi_{A}^{1}|[1+(3-2\sqrt{2})^{n}]^{-1}+|\varphi_{A}^{2}\rangle\langle\varphi_{A}^{2}|[1+(3+2\sqrt{2})^{n}]^{-1}.$

So in Eq. (S89), $\nu_{1}=[1+(3-2\sqrt{2})^{n}]^{-1}$,
$\nu_{2}=[1+(3+2\sqrt{2})^{n}]^{-1}$.