These authors contributed equally to this work.

These authors contributed equally to this work.

These authors contributed equally to this work.

[2] Aleksandra Maluckov

1]Institute of Physics, University of Rostock, Rostock, 18051, Germany

2]Vinča Institute of Nuclear Sciences, National Institute of the Republic of
Serbia, University of Belgrade, P.O.B. 522, Belgrade, 11001, Serbia

3]Departamento de Física and Millenium Institute for Research in Optics–MIRO,
Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Santiago,
8370448, Chile

# Supplementary information for “Observation of topologically protected
compact edge states in flux-dressed photonic graphene ribbon”

Gabriel Cáceres-Aravena<EMAIL_ADDRESS>Milica Nedić
<EMAIL_ADDRESS>Paloma Vildoso<EMAIL_ADDRESS>Goran
Gligorić<EMAIL_ADDRESS>Jovana Petrovic<EMAIL_ADDRESS>Rodrigo
A. Vicencio<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [

1\. Localized modes in ribbon

Compact localized modes

The compact localized states (CLS) are steady states of the $\pi$-flux
Hamiltonian. Analytically to derive the CLSs we start from model equation Eq.
(2) (main text) transformed to the configuration (lattice) space

$\displaystyle i\frac{da_{n}}{dz}$ $\displaystyle=$ $\displaystyle
b_{n}Ve^{-i\phi}+f_{n}Ve^{i\phi}+D_{n-1}V$ $\displaystyle i\frac{db_{n}}{dz}$
$\displaystyle=$ $\displaystyle a_{n}Ve^{i\phi}+c_{n}Ve^{-i\phi}$
$\displaystyle i\frac{dc_{n}}{dz}$ $\displaystyle=$ $\displaystyle
b_{n}Ve^{i\phi}+d_{n}Ve^{-i\phi}$ $\displaystyle i\frac{dd_{n}}{dz}$
$\displaystyle=$ $\displaystyle c_{n}Ve^{i\phi}+e_{n}Ve^{-i\phi}+A_{n+1}$ (S1)
$\displaystyle i\frac{de_{n}}{dz}$ $\displaystyle=$ $\displaystyle
d_{n}Ve^{i\phi}+f_{n}Ve^{-i\phi}$ $\displaystyle i\frac{df_{n}}{dz}$
$\displaystyle=$ $\displaystyle a_{n}Ve^{-i\phi}+e_{n}Ve^{i\phi}$

where $n$ is the index of unit cell, $\psi_{n}=(a\,b\,c\,d\,e\,f)_{n}$ is the
six-component vector, $V$ is the coupling strength (Fig. 1 in the main text)
and $\phi=\Phi/6$. The presence of flux is modeled by complex coupling
coefficients inside the unit-cells.

Being a steady state of $\Phi=0$ and $\Phi=\pi$-flux ribbon, the CLS can be
presented in a form $(a_{n},b_{n},c_{n},d_{n},e_{n},f_{n})(z)=exp(-i\lambda
z)(A_{n},B_{n},C_{n},D_{n},E_{n},F_{n})$. Therefore the set of eigenvalue
equations ($V=1$) is:

$\displaystyle\lambda A_{n}$
$\displaystyle=B_{n}e^{-i\phi}+F_{n}e^{i\phi}+D_{n-1}$ (S2)
$\displaystyle\lambda B_{n}$ $\displaystyle=A_{n}e^{i\phi}+C_{n}e^{-i\phi}$
$\displaystyle\lambda C_{n}$ $\displaystyle=B_{n}e^{i\phi}+D_{n}e^{-i\phi}$
$\displaystyle\lambda D_{n}$
$\displaystyle=C_{n}e^{i\phi}+E_{n}e^{-i\phi}+A_{n+1}$ $\displaystyle\lambda
E_{n}$ $\displaystyle=D_{n}e^{i\phi}+F_{n}e^{-i\phi}$ $\displaystyle\lambda
F_{n}$ $\displaystyle=A_{n}e^{-i\phi}+E_{n}e^{i\phi}$

where $\lambda$ is eigenenergy.

After simple algebraic procedure the CLS amplitudes can be estimated from
expressions:

$\displaystyle\\{A,\,B,\,C,\,D,\,E,\,F\\}_{n}$
$\displaystyle=\kappa\\{0,\,1,\,-\lambda
e^{i\phi},\,(\lambda^{2}-1)e^{2i\phi},\,-\lambda e^{3i\phi},\,e^{4i\phi}\\}$
(S3) $\displaystyle\\{A,\,B,\,C,\,D,\,E,\,F\\}_{n+1}$
$\displaystyle=\frac{\kappa(3-\lambda^{2})}{1-\lambda^{2}}\\{\lambda(1-\lambda^{2})e^{2i\phi},\,\lambda^{2}e^{3i\phi},\,-\lambda
e^{4i\phi},\,0,\lambda,\,\lambda^{2}e^{i\phi}\\}$

where $\kappa$ is an arbitrary constant ($\kappa=1$ due to convenience). Two
FBs at fluxless ribbon are singular due to the crossing of FBs with dispersion
bands, while at $\Phi=\pi$ the spectrum consists of six isolated FBs. The
fundamental CLSs hosted by two unit-cells, in the last case form the compact
nonorthogonal basis (Fig. S1).

Figure S1: The fundamental compact localized modes in $\pi$-flux graphene
ribbon. The CLSs corresponding to the same absolute value of $\lambda$ differ
by phase distribution among the nonzero amplitude sites, Eqs. S3.

Edge modes

The edge modes are caged to the lattice left/right edge, Figure S2 (a) and
(c). Their topological properties are determined in this study. The edge mode
amplitudes can be analytically estimated from recursive relations obtained
from Eq. (1) for $\lambda=0,\pm\sqrt{3}$ assuming the $D$ sites (left ribbon
edge) or $A$ sites (right ribbon edge) are empty. Here are presented
expressions for the left edge mode:

$\displaystyle\\{A,\,B,\,C,\,D,\,E,\,F\\}_{n}$ $\displaystyle=$ $\displaystyle
A_{n}\left\\{1,\,\frac{\lambda\,e^{i\phi}}{(1-\lambda^{2})},\,\frac{e^{i2\phi}}{\lambda^{2}-1},\,0,\,\frac{e^{-i2\phi}}{\lambda^{2}-1},\,\frac{\lambda\,e^{-i\phi}}{1-\lambda^{2}}\right\\}$
$\displaystyle A_{n+1}$ $\displaystyle=$
$\displaystyle-\frac{2\cos{(3\phi)}}{\lambda^{2}-1}A_{n}.$ (S4)

Taking arbitrary value of $A_{1}=c$ we obtain $A_{n\neq
1}=(-2\cos(3\phi)/(\lambda^{2}-1))^{n-1}c$, which in the finite lattice case
with $N$-sites provides condition:

$\left(-\frac{2\cos(3\phi)}{\lambda^{2}-1}\right)^{N}c=0\ .$ (S5)

Only in the fully FB ribbon, $\Phi=\pi$ ($\phi=\pi/6$) Eq. (S5) is satisfied
for all $N$, while in the other flux configurations the edge modes posses
$N$-dependent tails, as confirmed by the calculation of the participation
ratio:

$PR=\frac{\mathcal{P}^{2}}{2N}\sum_{i}^{N}\left(|\Psi_{i}|^{4}\right)^{-1},\quad\Psi_{i}=(A\,B\,C\,D\,E\,F)_{i}^{T}$

where $\mathcal{P}=\sum_{i}\Psi_{i}^{\dagger}\Psi_{i}$ is the total power,
Fig. S2 (b) and (d). The participation ratio gives information regarding the
fraction of strongly excited lattice sites. Dynamical simulations confirm
extreme robustness of edge modes.

Figure S2:  (a) The schematic view of the amplitude distribution of zero-edge
double degenerated mode ($\lambda=0$) and (c) double degenerated edge mode
($\lambda=\sqrt{3}$) in graphene ribbon vs. flux in the area
$\Phi=[2\pi/3,4\pi/3]$. (b) and (d) The participation ratio in the logarithmic
scale for modes shown in (a) and (c), respectively.

2\. Symmetries of armchair graphene-like ribbon

Following the Altland-Zirnbauer symmetry classification of 1D systems, the
graphene armchair ribbon belongs to AIII symmetry class [1, 2]. It is
characterized by bipartite symmetry regardless of the flux value. As a
consequence, the energy band spectrum is symmetric with respect to zero
energy.

The bipartite symmetric Hamiltonian $\hat{H}$ is related to a chirality
operator $\hat{G}$, $\hat{H}=\hat{G}\hat{H}\hat{G}=-\hat{H}$, which can be
expressed in a symbolic form:

$\quad\hat{G}=\hat{B}+\hat{D}+\hat{F}-(\hat{A}+\hat{C}+\hat{E})\ ,$ (S6)

with $\hat{A}=\sum_{n}\ket{A_{n}}\bra{A_{n}}$ being the sublattice projector
for the $A$ sites, $\hat{B}$ for $B$ sites, etc. Chirality operator provides
the possibility to create the topologically non-trivial phases in the system.

In addition, the armchair fluxless ribbon [3, 4] is symmetric with respect to
the horizontal axis for any number of hexagons and with respect to the
vertical axis for even number of hexagons. The introduction of nonzero flux
asymmetrizes the ribbon, which is associated with opening of gaps between the
periphery bands at the boundaries of the first Brillouin zone (BZ).

The gap between the central bands is closed and reopened at the flux values
$\Phi=2\pi/3$ and $\Phi=4\pi/3$ at the center of first BZ. At these particular
flux values, the C3 symmetry inside of unit cells (hexagons) is reestablished.
Hence, the opening/closing of the central gap, which is followed by the
creation/destruction of the chiral zero-edge states, can be related to
establishing/breaking of this symmetry. For $\Phi=\pi$ zero edge mode is
hosted by only three sites inside the boundary unit cell.

3\. Determination of topological phase: bulk-edge correspondence

The nontrivial topological phase of periodic systems (lattices, chains) is
described by quantized topological invariant derived from Berry (geometric)
phase after the adiabatic evolution of the system [5, 6]. The role of the
topological invariant in 1D lattice is played by the winding number or related
to it, Zak-phase, which describes how a particular function (connected with
Bloch states) winds around the 1st Brillouin zone. These two quantities in 1D
are proportional: $Z_{p}=W\pi$, where $W$ is winding number and $Z_{p}$ the
Zak phase. In 1D lattices integer values of Zak-phase indicate the presence of
topologically protected edge states which is in pursuance of the bulk-edge
correspondence.

The SSH lattice is a paradigm of 1D topological insulators formed of two-
sites’ unit cells. Peculiarity of the SSH lattice is topological
nontriviallity in the parametric area $v/w>1$ described by $Z_{p}=\pm 1$,
where parameters $v,w$ are the intra and inter site coupling strengths,
respectively. Strictly, the threshold value is dependent on the number of
cells. The band spectrum is formed of two bands which are connected by the
topologically protected edge (zero) mode in topologically nontrivial phase.
There are a few mathematical methods to derive the value of $Z_{p}$, the
geometrical, mean chiral displacement and projector method.

The first attempts regarding the topology of multi-band periodic systems are
developed on the effective SSH-like equivalent two-band systems [7]. Analyzing
our multi FB system we theoretically and experimentally prove the nontrivial
topological properties of central bands and extremely robust compact zero-edge
modes. Moreover, in the flux range $[2\pi/3,4\pi/3]$, the topological non-
triviality is numerically confirmed too. The topological nontriviality in the
rest of parametric space of flux-dressed ribbon (see the main text) is just
indicated by utilizing the projector method. Final confirmation of these
finding requires additional study because the other two methods presently
could not be applied to periphery bands.

The geometric method: reduction to SSH-like Hamiltonian

The governing idea of the geometric method is based on the reduction of the
Hamiltonian to the SSH-type (two-band) one, which can be expressed as:

$\hat{H}=\left(\begin{array}[]{cc}\hat{0}&\hat{h}\\\
\hat{h}^{\dagger}&\hat{0}\\\ \end{array}\right).$ (S7)

In the $N$-dimensional SSH Hamiltonian the block matrices $\hat{h}$ are
$(N/2)\times(N/2)$ matrices.

We found that the following unitary matrix $S$ ($S^{\dagger}S=I$, therefore
$S^{-1}=S^{\dagger}$):

$\displaystyle S=\left(\begin{array}[]{cccccc}1&0&0&0&0&0\\\ 0&0&0&0&0&1\\\
0&1&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&0&1&0\\\
\end{array}\right)\ ,$ (S14)

transforms $H(k)$ (Eq. (2) in main text) into :

$\displaystyle
S^{-1}HS=\left(\begin{array}[]{cccccc}0&0&0&e^{-ik_{x}}&e^{i\phi}&e^{-i\phi}\\\
0&0&0&e^{-i\phi}&0&e^{i\phi}\\\ 0&0&0&e^{i\phi}&e^{-i\phi}&0\\\
e^{ik_{x}}&e^{i\phi}&e^{-i\phi}&0&0&0\\\ e^{-i\phi}&0&e^{i\phi}&0&0&0\\\
e^{i\phi}&e^{-i\phi}&0&0&0&0\\\ \end{array}\right)\ ,$ (S21)

which can be rewritten in the form (S7) in terms of off-diagonal matrix $h$

$\displaystyle h=\left(\begin{array}[]{ccc}e^{-ik_{x}}&e^{i\phi}&e^{-i\phi}\\\
e^{-i\phi}&0&e^{i\phi}\\\ e^{i\phi}&e^{-i\phi}&0\\\ \end{array}\right).$ (S25)

By calculating the determinant of the off-diagonal sub-matrix we obtained
$\det(h)=-\cos k_{x}+2\cos(\Phi/2)+i\sin k_{x}$, where $\Phi=6\phi$. In the
complex plane this determinant can be presented by unit circle which position
depends on the parameter $\Phi$. The winding number can be defined by how many
types loops encircle the origin. Here the origin is circled ones for the
fluxes $\Phi=[2\pi/3,4\pi/3]$. From this follows that only in this parameter
region the winding number is $1$, Fig. S3.

Figure S3: The determinant of $h$ from Eq.(S25) for
$\Phi=\\{0,\pi/3,2\pi/3,\pi,4\pi/3,2\pi\\}$.

In addition, the value of the winding number $w$ can be determined from
expression:

$\displaystyle w=\frac{1}{2\pi
i}\int_{0}^{2\pi}dk_{x}\;\frac{\partial}{\partial
k_{x}}\log\left(\det(h(k_{x}))\right),$ (S26)

showing that it is unity in the range of fluxes $\Phi=[2\pi/3,4\pi/3]$, Fig.
S4. Thus, the topological nontriviallity of the central bands and via the
bulk-edge correspondence of the compact edge modes is confirmed.

Figure S4: (a) The calculation of the winding number using equation (S26) with
$h$ for the red line and $h^{\dagger}$ the black line. (b) Numerically
obtained values of MCD by initially exciting single A (blue circle), C (red
squares), and E (cyan triangles) site belonging to the first sublattice. (c)
The same as (b) just by initially exciting single D (blue circle), B (red
squares), and F (cyan triangles) site belonging to the second sublattice. The
black dashed lines are the average values of the MCD single-site excitations
in the corresponding sublattice.

The mean field/chirality displacement

In the root of the Mean Field Displacement (MFD) method and the projector
methods is the concept of geometric, i.e. the Berry phase which describes how
a global phase accumulates as a complex Bloch eigenvector is carried around a
closed loop in the $\vec{k}$ space inside the 1st Brillouin zone. In the 1D
periodic system in the basis of Wannier functions, the Fourier-transform
partners to the Bloch functions, the Berry phase is proportional to the
displacement of the center of mass of basis vectors [5]. In the literature,
the MFD method has been widely used in the 1D two-band systems to derive
numerically and experimentally the value of the topological invariant of
isolated band, which is determined by the mean displacement of the center of
propagating wavepacket initiated by single site lattice excitation [8, 9, 10]:

$M_{FD}(L_{z})=\frac{1}{L_{z}}\int_{z=0}^{z=L_{z}}\overline{P}(z)\;dz=\frac{1}{L_{z}}\int_{z=0}^{z=L_{z}}\frac{\sum_{n=0}^{N-1}(x_{n}(z)-x_{n}(0))|\psi_{n}(z)|^{2}}{\sum_{n=0}^{N-1}|\psi_{n}(z)|^{2}}\;dz,$
(S27)

where $L_{z}$ is the propagation distance in the $z$-axis direction and $x$ is
position in ribbon, and $\Psi=[\psi_{1}\,\psi_{2}...\psi_{N}]^{T}$ is the
wavepacket. The Zak-phase is given by the expression: $Z_{p}[\pi]=2\pi
M_{FD}$. The integer value of $Z_{p}$ in units of $\pi$ is related with the
total number of edge states formed in the gaps between the isolated bands,
just after the topological phase transition. This is a confirmation of the
bulk-edge correspondence [11, 12].

The MFD method is experimentally utilized to clarify the theoretical
predictions in diverse 1D and quasi 1D lattices [8, 9, 10]. It has been proved
that the value $2M_{FD}$ approaches asymptotically the winding number [9]
(Fig. S4), or related Zak phase (with a $\pi$ factor) [8] for large values of
the propagation distance.

Strictly speaking the above described procedure is fully developed for the 1D
lattice systems with two bands [13, 8, 9, 10]. However, the MFD approach can
be applied to the lattices’ with multilevel spectra which geometry and related
symmetries allow the reduction to the effective two-level spectrum [14, 15,
7]. In other circumstances problem is model-dependent.

Further implementations of this method for more complex lattices had been made
taking advantage of the chiral symmetry of the bipartite lattices [13, 14,
15]. The main idea is to utilize the sublattice related peculiarities of such
systems. Thus in the chiral lattice with $4$ sites ($2$ from each sublattice)
[14], the MFD is calculated using two different initial excitations of sites
within the same unit cell, each meaning the excitation of just one of
sublattices, are followed. The corresponding MFD are summed and their average
value is classified as a Mean Chiral Displacement (MCD)[13], Fig. S4 (b). The
winding number (Zak-phase) is then $Z_{p}=2MCD/\pi$. Here we adopted the same
idea to our chiral ribbon which unit cell consists of six sites, three from
each of sublattice. Chiral symmetry is associted with the chirality operator
$G$. We excite each of sites inside the middle unit cell hosting the eigen-
states of the chiral operator, $A,C,E$ or $B,D,F$ (the eigenvalue of $G$ are
respectivelly, $1$ and $-1$). The averaged value of three MFD is taken to
determine the Zak-phase of central ribbon bands. The shape of the MCD vs.
$\Phi$ fully corresponds to the $W$ vs. $\Phi$ derived by geometric method,
Fig. S4.

Calculation of Zak-phase: projector method

The projector method is established in the basis of Bloch vectors [5, 16]. The
central feature is the projection operator of eigenmodes
$\hat{P}(k)=|\tilde{\psi}(k)\rangle\langle\tilde{\psi}(k)|$ [17, 5, 16] from
which the topological invariant is derived:

$Zp=\imaginary\left(\ln(\Tr\left[\prod_{n=1}^{N}\hat{P}(k_{n})\right])\right)\
,$ (S28)

where $k_{n}=2\pi n/L$ is the discrete momentum.

The numerical procedure starts by solving the eigenvalue problem of the ribbon
Hamiltonian in $k$-space, which is now discretized in $N$-grids separated by
$2\pi/N$ distance. For each $k_{i},\,i=[1,..,N]$ within the 1st BZ,
$k=[-\pi,\pi]$ the set of $6$ eigenvalues and corresponding eigenvectors is
determined. Each of $6$ eigenvalues belongs to different bands. The $6\times
6$ projector matrix is formed as the correlation matrix between the
eigenvectors from the whole eigenvalue set. After straightforward numerical
calculation the values of Zak-phase determined by Eq. S28, for each band are
summarized in Table I.

Flux | band 1 | band 2 | band 3 | band 4 | band 5 | band6
---|---|---|---|---|---|---
$\Phi=0$ | 0 | 0 | 0 | 0 | 0 | 0
$0<\Phi<2\pi/3$ | -1 | 1 | 0 | 0 | 1 | 1
$2\pi/3<\Phi<4\pi/3$ | 1 | -1 | 1 | -1 | -1 | 1
$4\pi/3<\Phi<2\pi$ | 1 | -1 | 0 | 0 | -1 | 1
$\Phi=2\pi$ | 0 | 0 | 0 | 0 | 0 | 0
Table 1: Table I: The values of Zak -phase for each band for different values
of flux.

In contrast to other methods the projector method gives information about the
topological phase of all isolated bands, the central and periphery ones.
Assuming the validity of the bulk-surface correspondence we claim that the
edge modes between the nontrivial bands are nontrivial too, i.e. topologically
protected modes.

Theoretical vs Experimental model

The theoretical model is composed of $6$ sites per unit cell as described in
the main text and at the beginning of the supplementary material. On the other
hand, the experimental model for $\pi$ flux is composed of $12$ sites per unit
cell but can be mapped to the theoretical model. The mapping of the
experimental model without flux is simple given by assuming the sites are
waveguides with the same propagation constant. However, in the case of the
$\pi$ flux, the experimental model is composed of $12$ sites per unit cell and
the value of the phase in the plaquette is achieved by a synthetic magnetic
flux [18] owing a negative hopping from the second order photonic mode in a
middle waveguide between two hexagons.

The Fig. S5 shows the experimental model. When the couplings are all the same
and the propagation constant (on-site energy) are equal to zero, then the
model has 6 bands each with degeneracy 2. Given that each eigenvalue has
degeneracy two, the system is reducible. The reduction is possible due to the
left plaquette is the same as the right plaquette but mirrored. This mirroring
is associated to a symmetry in the system that creates a degeneracy, therefore
we can reduce the system to the half without losing information about
eigenstates. In this manner, the 12 sites per unit cell system is equivalent
to a 6 sites system after mirroring every other unit cell. Moreover, we
calculate the eigenvalues of the momentum Hamiltonian (the basis of Bloch
waves) and we find 12 values, each with degeneracy 2, therefore different 6
eigenvalues are $\\{\pm 2.17,\pm 1.48,\pm 0.31\\}$ which coincides with the 6
eigenvalues of the theoretical model when $\Phi=\pi$ and $V=1$ (main text).

Figure S5: The scheme of the experimentally realized $\pi$-flux graphene
ribbon. The gray dashed square denotes the $n$-th unit cell composed of 12
sites.

4\. Experimental characterization setups

Supercontinuum characterization Setup. The excitation of a photonic lattice is
typically achieved by focusing a polarized gaussian-like laser beam at the
input facet of a photonic lattice under study. Then, the output facet is
imaged and obtained the intensity output profile on a CCD camera. This
standard method is relatively simple to implement and it is limited to the
excitation of a single waveguide only, at a fixed wavelength. The coupling
constants of a photonic lattice depend on the system’s geometry (determined by
waveguide distances), but also on the excitation wavelength. The width of the
mode profile at each waveguide increases for a larger wavelength and,
therefore, for a fixed geometry it is possible to change the coupling
constants by using a different excitation wavelength. This is the basis of the
Supercontinuum setup, where by changing the input wavelength we modify the
coupling constants, what simultaneously alters the effective propagation
distance of the system. The dynamical variable on tight-binding-like models is
given by $``Vz^{\prime\prime}$, where $z$ is the propagation coordinate and
$V$ the coupling constant. For a fixed glass length, the coordinate $z$ is
also fixed and, therefore, a modification of the coupling constants produces a
shorter/larger effective dynamics. The experimental setup is sketched in Fig.
S6 and it relies on a Supercontinuum laser source (SC) YSL SC-5, emitting in
the wavelength range of $450-2400$ nm. The operational range is narrowed down
to $450-1450$ nm by using acoustic filtering, with a resolution of $5$ nm and
$\sim 1$ mW at each wavelength. The SC beam is first horizontally polarized
and, then, focused onto the input facet of the photonic chip (PC) using a
$10\times$ microscope objective. The output image is obtained by using a
second $10\times$ microscope objective and a beam profiler.

Figure S6: Scheme of SC setup. It consists on a Supercontinuum (SC) laser
source, a $10\times$ Microscope Objectives, and a beam profiler (mirrors, nd-
filters, iris, polarizers, XYZ stages, goniometers are not including here for
simplicity).

Image Setup. The excitation of compact edge states on a graphene ribbon
demands the generation of an input condition having a specific amplitude and
phase profile. To achieve this we use an image setup described in Figure
S7(a). We first expand and collimate a $640$ nm laser beam using a $20\times$
microscope objective (MO) and a lens with a focal length of $f=125$ mm, which
is then modulated in amplitude by a Holoeye LC2012 transmission spatial light
modulator (SLM). This enables the generation of an arbitrary amplitude pattern
projected onto the beam, facilitating the simultaneous excitation of multiple
waveguides. During this part, the intensity and size of the light spot (disk)
can be configured for each waveguide under study. Once the image amplitude
profile is formed, the phase structure is added by means of a Holoeye PLUTO
reflection SLM. An example of an input image condition, formed by amplitude
and phase modulation, is shown in Fig. S7(b1). Upon completing the modulation
process, the image is then optically transferred to the input facet of the
photonic chip by using a telescope and a $4\times$ MO. An example of an output
profile is shown in Fig. S7(b2), after the input condition has traveled along
a 70 mm graphene photonic lattice.

Figure S7: Image setup. (a) Scheme of experimental image setup based on two
SLM’s to generate arbitrary input conditions. It consists on several linear
polarizers (P), $\lambda/2$ waveplates, beam splitters (BS), several lenses,
and two CCD cameras. (b1) Amplitude and phase modulation. (b2) Imaged output
profile.

5\. Experimental analysis of trivial and topological graphene ribbon lattices.

Single Site Bulk excitation. We first explore the strong contrast in between
the dynamical properties of the trivial (SS) and topological (SP) photonic
lattices. We start by exciting the unit cell at the lattice bulk using
different input wavelengths coming from a SC laser source. Fig. S8 shows a
compilation of several output profiles after exciting sites A, B, C, D, E, and
F of a (a) trivial SS and (b) nontrivial SP lattices. For the SS case, we
observe flat band oscillation and weak transport after exciting sites B, C, E,
and F, and ballistic transport only after exciting sites A and D. On the
contrary, for the SP case, the light is mostly confined at two neighbouring
unit cells due to the excitation of a perfectly flat band spectrum only,
independent of the input site.

Figure S8: Output profiles, at different wavelengths, after single site bulk
excitations of the unit cell for (a) SS and (b) SP lattices. Yellow circles
indicate the input site. From left to right: B, C, E, F, A and D input sites,
respectively.

We analyze the obtained data by computing the Participation Ratio (PR), center
of mass (CM), and second moment ($m_{2}$), for each image shown in Fig. S8,
and collect results in Fig. S9(a). In the SS lattice case, after exciting
central A and D sites (full-blue curves) we observe that the PR and $m_{2}$
increase upon wavelength, as an indication of the ballistic transport
originated from a dispersive-only band excitation. On the other hand, after
exciting sites B,C,E or F (see full-red curves) we notice a slight oscillation
in the PR, and a slower increment of $m_{2}$. This indicates a mixture in
between the excitation of dispersive and flat bands. The CM in Fig. S9(a2)
does not exhibit a significant increment upon wavelength, as the energy
spreads out symmetrically from the lattice center. On the other hand, the
excitation of a SP nontrivial lattice [see dashed blue and dashed red data in
Fig. S8(a)] produces an oscillatory bounded behavior, where the PR, CM and
$m_{2}$ oscillate for an increasing wavelength. This is in agreement with the
caging phenomena shown in Fig. S8(b), and directly confirms the flat-only band
spectrum for a nontrivial flux.

Figure S9: Participation ratio (PR), Center of Mass (CM) and second moment
(m2) for SS (full lines) and SP (dashed lines) lattices, and for (a) bulk and
(b) edge excitation. A-site (B and F) excitation is shown in blue (red) color.

Single Site Edge excitation. Now, we study the excitation of the lattice edges
at different wavelengths. For simplicity, we compare the excitation of sites
A, B and F at the first unit cell (left-edge) for both, SS and SP, cases. Fig.
S10 compiles our results, with the output profiles for different excitation
wavelengths. For a SS lattice, we observe that an A site excitation produces
ballistic transport as a result of exciting the dispersive part of the
spectrum only. Data analysis in Fig. S9(b) (see full-blue curve) shows a clear
increment of the CM, with a weakly increment of the PR and a rather constant
$m_{2}$. However, if we compare this result with the excitation of sites B and
F, we clearly notice the competition in between the flat (tendency to trap the
energy at the excited unit cell) and the dispersive (discrete diffraction)
bands. The analysis in Fig. S9(b) (see full-red curve) shows an oscillating
PR, a slowly increment of CM, and a rather fast increment on $m_{2}$. The edge
excitation of an SP ribbon lattice, having a nontrivial flux of $\pi$, shows
remarkable differences compared to a SS lattice. We clearly observe in Fig.
S10(b), for the three excitation sites, that the light remains trapped/caged
and oscillating at the first unit cell only, as a result of the direct
excitation of the three edge states existing at this regime. We notice a quite
low background for all wavelengths, due to the fact that the light keeps
strongly localized at the edge unit cell. The analysis shown in Fig. S9(b)
(see blue and red dashed lines) indicates a small and weakly oscillating PR, a
rather constant CM, and a very small $m_{2}$, indicating absence of transport.
These results show the robustness of the edge localization for a wide range of
wavelengths.

Figure S10: Output profiles after excitation of sites A, B and F, for the
first (left-edge) unit cell, for different input wavelengths, and for (a) SS
and (b) SP lattices. Yellow circle indicates the input site.

6\. Excitation of zero energy edge states.

Robustness to the change of input phase. We demonstrate the robustness of the
edge state in $\pi$-graphene lattices by applying a three-site input condition
using the image setup at $640$ nm. In Fig. S11 we observe the flux $\Phi=\pi$
effect for the simultaneous excitation of A,C and E sites. Similar to the case
of a single-site excitation shown in Fig. S10, the light remained trapped
inside the unit cell for the SP lattice. A precise excitation of the
topological edge state is achieved around a phase difference $\Delta\phi=\pi$
in between A and C/E sites. This behavior demonstrates the robustness of the
edge state to a phase change in an SP topological lattice. In contrast, the
excited pattern in the SS lattice shows dispersive-like behavior, as shown for
$\Delta\phi=\pi$ in Fig. S11. This clearly show the absence of the edge modes
and edge localization for $0$-flux ribbon graphene lattices.

Figure S11: Output profile for the simultaneous excitation of A,C and E sites
with different phases in between the A and C/E sites. The phase is changed at
site A only, while keeping constant the phase at C/E sites. The horizontal
arrow indicates the increment of the phase difference $\Delta\phi$.

Steady edge mode propagation. Finally, we study the dynamical stability for
the edge compact topological mode along the propagation coordinate $z$ for a
topologically nontrivial SP lattice. We first excite a three (A, C, and E)
sites input pattern having a $\Delta\phi=0$ input phase and show compiled
results in Fig. S12-top. For this in-phase initial condition the light does
not remain trapped at the excitation sites and oscillates at the first unit
cell, similar to the single-site excitation shown in Fig.S10(b). An out of
phase input excitation, corresponding to a perfectly compact edge topological
state profile, shows how this mode remains perfectly compact along the
propagation coordinate, with only slight fluctuations associated with
fabrication errors (every image is obtained for a different sample [18]). This
compiled results demonstrate the robustness and stability of the edge
topological state along propagation, as different samples introduce a
different amount of disorder.

Figure S12: Output profile versus propagation coordinate $z$ for a
simultaneous excitation of sites A, C and E, with an in-phase structure
$\Delta\phi=0$ (top) or an out of phase structure $\Delta\phi=\pi$ (bottom).

## References

  * [1] Chiu, Ch.-K., Yao, H. & Ryu, Sh. Classification of topological insulators and superconductors in the presence of reflection symmetry. Phys. Rev. B 88, 075142 (2013).
  * [2] Zurita, J., Creffield, C. & Platero, G. Tunable zero modes and quantum interferences in flat-band topological insulators. Quantum 5, 591 (2021).
  * [3] Jürß, C. & Bauer, D. High-order harmonic generation in hexagonal nanoribbons. The European Physical Journal Special Topics 230, 4081–4089 (2021).
  * [4] Drüeke, H. & Bauerb, D. High-harmonic spectra of hexagonal nanoribbons from real-space time-dependent Schrödinger calculations. The European Physical Journal Special Topics 230, 4065–4070 (2021).
  * [5] Vanderbilt, D. Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators. (Cambridge: Cambridge University Press, 2018).
  * [6] Asbóth J. K., Oroszlány, L. & Pályi, A. A Short Course on Topological Insulators. (Springer, Berlin, 2016).
  * [7] Lee, C.-S., Io, I.-F. & Kao, H.-C. Winding number and Zak phase in multi-band SSH models. Chinese Journal of Physics 78, 96-110 (2022).
  * [8] Longhi, S. Probing one-dimensional topological phases in waveguide lattices with broken chiral symmetry. Opt. Lett. 43, 4639-4642 (2018).
  * [9] Jiao, Z.-Q. et al. Experimentally detecting quantized Zak phases without chiral symmetry in photonic lattices. Phys. Rev. Lett. 127, 147401 (2021).
  * [10] Cáceres-Aravena, G. et al. Edge-to-edge topological spectral transfer in diamond photonic lattices. APL Photonics 8, 080801 (2023).
  * [11] Lu, L., Joannopoulos, J. D. & Soljačić, M. Topological photonics. Nat. Photon. 8, 821 (2014).
  * [12] Ozawa, T. et al. Topological photonics. Rev. Mod. Phys. 91, 015006 (2019).
  * [13] Cardano, F. et al. Detection of Zak phases and topological invariants in a chiral quantum walk of twisted photons.Nat. Comm. 8, 15516 (2017).
  * [14] Maffei, M., Dauphin, A., Cardano, F., Lewenstein, M. & Massignan, P. Topological characterization of chiral models through their long time dynamics. New J. Phys. 20, 013023 (2018).
  * [15] D’Errico, A. et al. Bulk detection of time-dependent topological transitions in quenched chiral models. Phys. Rev. Research 2, 023119 (2020).
  * [16] Leykam, D. & Smirnova, D. A. Probing bulk topological invariants using leaky photonic lattices. Nat. Phys. 17, 632–638 (2021).
  * [17] Leykam, D. & Smirnova, D. A. Photonic quantum Hall effect. In Encyclopedia of Condensed Matter Physics (Second Edition) (eds. Chakraborty, T.) 575-586 (Academic Press, 2024).
  * [18] Cáceres-Aravena, G., Guzmán-Silva, D., Salinas, I. & Vicencio, R. A. Controlled transport based on multiorbital Aharonov-Bohm photonic caging. Phys. Rev. Lett. 128, 256602 (2022).


So suppose $Z$ is a reducing sphere for the chamber $C$ in $\mathcal{C}$, and
$C=A_{C}\cup_{T_{C}}\cup B_{C}$ is the Heegaard splitting determined by $T$.
By [FS2, Theorem 1.6], $T_{C}$ can be isotoped in $C$ so that $Z$ becomes a
reducing sphere for $T_{C}$ (that is, $Z$ intersects $T_{C}$ in a single
essential circle) and the isotopy of $T_{C}$ is unique up to eyeglass moves.
Amalgamate the Heegaard splittings of all the chambers into the Heegaard
splitting $T$ and observe that $Z$ remains a reducing sphere for the Heegaard
surface $T$. But since $M$ is irreducible, $Z$ bounds a ball in $M$. Since $Z$
is a reducing sphere for the Heegaard splitting $T$, $T$ restricts to a non-
trivial Heegaard splitting of this ball. The induced Heegaard splitting of the
ball is necessarily a bubble, via Waldhausen again. ∎

######  Proposition 11.7.

Suppose $\mathcal{C}$ and $\mbox{$\mathcal{C}$}^{\prime}$ are flagged chamber
complexes in $M$, each supporting the Heegaard splitting $M=A\cup_{T}B$.
Suppose $\mathcal{C}$ and $\mbox{$\mathcal{C}$}^{\prime}$ belong to certifying
flagged chamber complex decomposition sequences $\vec{\mbox{$\mathcal{C}$}}$
and $\vec{\mbox{$\mathcal{C}$}^{\prime}}$ respectively. Suppose $Z$ is a
certificate issued by $\vec{\mbox{$\mathcal{C}$}}$ and $Z^{\prime}$ is a
certificate issued by $\vec{\mbox{$\mathcal{C}$}^{\prime}}$. If
$\vec{\mbox{$\mathcal{C}$}}\sim\vec{\mbox{$\mathcal{C}$}^{\prime}}$, there is
a sequence of reducing spheres $K(Z)=K_{1},...,K_{m}=K(Z^{\prime})$ for $T$,
each bounding a non-trivial bubble in $T$, so that for each $1\leq i<m$,
$K_{i}$ is disjoint from $K_{i+1}$.

###### Proof.

It suffices to consider the elementary situation in which
$\vec{\mbox{$\mathcal{C}$}}$ interacts with
$\vec{\mbox{$\mathcal{C}$}^{\prime}}$. In that case, let
$\mbox{$\mathcal{C}$}_{j}$ be a flagged chamber complex in the sequence
$\vec{\mbox{$\mathcal{C}$}}$, $\mbox{$\mathcal{C}$}^{\prime}_{k}$ a flagged
chamber complex in the sequence $\vec{\mbox{$\mathcal{C}$}^{\prime}}$, each
containing chambers $C,C^{\prime}$ respectively, so that, say $C$ contains
certifying spheres for both $\mbox{$\mathcal{C}$}_{j}$ and
$\mbox{$\mathcal{C}$}^{\prime}_{k}$.

If $C$ itself is an occupied ball, so its certificate is its boundary, the
certificate for $\mbox{$\mathcal{C}$}^{\prime}_{k}$ lying within it is
disjoint. Suppose $C$ is not an occupied ball, so its certificate $Z_{C}$ is a
reducing sphere for $C$. A standard innermost (in $Z^{\prime}_{C}$) disk
argument shows that there is a reducing sphere in $C$ (hence a certificate for
$\vec{\mbox{$\mathcal{C}$}}$) that is disjoint from the certificate
$Z^{\prime}_{C}$ for $\mbox{$\mathcal{C}$}^{\prime}_{k}$. So in either case,
there are corresponding reducing spheres $K_{C},K^{\prime}_{C}$ for $T$ that
are disjoint, lie in this chamber, and certify for
$\vec{\mbox{$\mathcal{C}$}}$ and $\vec{\mbox{$\mathcal{C}$}^{\prime}}$
respectively.

It suffices now to show there is a sequence of spheres as in the Proposition,
beginning with $K(Z)$ and ending with $K_{C}$. Then, symmetrically, there is
such a sequence from $K^{\prime}_{C}$ to $K(Z^{\prime})$.

We may as well assume that the original certificate $Z$ lies in a chamber
$C^{\prime\prime}$ of some $\mbox{$\mathcal{C}$}_{i}$ with $i\leq j$, since we
can argue symmetrically, switching the roles of $Z$ and $Z_{C}$ if $i>j$.
Suppose $i=j$, so $K(Z)$ and $K_{C}$ lie in the same chamber complex. If they
lie in separate chambers then of course they are disjoint. If they both lie in
the chamber $C$, and so are reducing spheres there, apply a classic innermost
disk argument to show that there is a sequence of reducing spheres in $C$,
each disjoint from the next, beginning with $K(Z)$ and ending with $K_{C}$.

Finally, suppose $i<j$ and induct on $j-i$. It is again a classic innermost
disk argument that there is a reducing sphere $\hat{K}$ in $C$ that is
disjoint from decomposing disks $\mathcal{D}$ in
$\mbox{$\mathcal{C}$}_{i}\xrightarrow{\mbox{$\mathcal{D}$}}\mbox{$\mathcal{C}$}_{i+1}$
and so a sequence of reducing spheres for $C^{\prime\prime}$ beginning with
$K(Z)$ and ending with $\hat{K}$. Then we may as well substitute $\hat{K}$ for
$K(Z)$, and note that $\hat{K}$ survives intact as a sphere in a chamber
$\hat{C}$ of $\mbox{$\mathcal{C}$}_{i+1}$. If $\hat{K}$ is a reducing sphere
for $\hat{C}$ this completes the inductive step, and we are done.

It is possible that $\hat{K}$ is not a reducing sphere for $\hat{C}$, though.
For example, the ball that $\hat{K}$ bounds in $M$ could have contained a
handlebody chamber adjacent to $C^{\prime\prime}$ in
$\mbox{$\mathcal{C}$}_{i}$, and the handlebody chamber disappears in
$\mbox{$\mathcal{C}$}_{i+1}$ because it is decomposed by $\mathcal{D}$ into a
goneball. The bubble in $T$ that it bounds then already appears as a bubble in
the induced Heegaard splitting of $\mbox{$\mathcal{C}$}_{i+1}$. In this case,
simply make the bubble very small in $\hat{C}$, avoiding all further
decomposing disks, until the bubble arrives intact in
$\mbox{$\mathcal{C}$}_{j}$ disjoint from the sphere $K_{C}$. There is
ambiguity as to which chamber the bubble lies in as the decomposition sequence
proceeds, but recall from the comments following Corollary 4.4 that how
decomposing disks intersect the Heegaard surface are anyway ambiguous, partly
by action of the eyeglass group, but also because bubbles may slide across
decomposing disks. So just ensuring that later decomposing disks miss the
bubble adds no new ambiguity. ∎

######  Theorem 11.8.

For any Heegaard splitting $S^{3}=A\cup_{T}B$ of $S^{3}$, the Goeritz group
$G(S^{3},T)$ is generated by the eyeglass group.

###### Proof.

Let $g=genus(T)$ and suppose inductively that the theorem has been proven for
all splittings of $S^{3}$ of lower genus. Per [FS1] we may as well take $g\geq
4$, though $\geq 3$ is all that’s used. Given $\tau\in G(S^{3},T)$,
Proposition 11.5 provides certifying flagged chamber complex decomposition
sequences $\vec{\mbox{$\mathcal{C}$}},\vec{\mbox{$\mathcal{C}$}^{\prime}}$ in
$S^{3}$ so that
$\vec{\mbox{$\mathcal{C}$}}\sim\vec{\mbox{$\mathcal{C}$}^{\prime}}$, and for
any certificate $Z$ issued by $\vec{\mbox{$\mathcal{C}$}}$, $\tau(Z)$ is a
certificate issued by $\vec{\mbox{$\mathcal{C}$}^{\prime}}$.

Choose one such certificate $Z$, defining (up to action of the eyeglass group)
a reducing sphere $K$ for $T$, which naturally divides $T$ into two non-
trivial bubbles. According to Proposition 11.7 there is then a sequence of
non-trivial bubbles
$\mbox{$\mathfrak{b}$}=\mbox{$\mathfrak{b}$}_{1},...,\mbox{$\mathfrak{b}$}_{m}=\tau(\mbox{$\mathfrak{b}$})$
so that for each $1\leq i<m$, the bubbles
$\mbox{$\mathfrak{b}$}_{i},\mbox{$\mathfrak{b}$}_{i+1}$ are disjoint, and
$\mathfrak{b}$ is a bubble bounded by $K$.

Motivated by [FS1, Section 3] let $T_{g}$ be a fixed genus $g$ Heegaard
splitting of $S^{3}$, with orthogonal meridian disks in $A$ and in $B$,
labeled $a_{1},...,a_{g}$ and $b_{1},...,b_{g}$ respectively. (By orthogonal
disks we mean $|a_{i}\cap b_{j}|=\delta_{ij}$.) Two homeomorphisms
$h_{1},h_{2}:(S^{3},T)\to(S^{3},T_{g})$ are eyeglass equivalent if there is
$\alpha\in\mbox{$\mathcal{E}$}$ so that $h_{2}$ is isotopic to $h_{1}\alpha$.

Claim 1: Under our inductive assumption, a non-trivial bubble $\mathfrak{b}$
in $T$ defines an eyeglass equivalence class of homeomorphisms
$(S^{3},T)\to(S^{3},T_{g})$.

Let $a\subset A$ be a primitive meridian disk, meaning that there is a
meridian disk $b\subset B$ so that $|a\cap b|=1$. It is shown in [FS1,
Corollary 3.6] that, under our inductive assumption, a choice of such a
meridian disk $a$ determines a homeomorphism
$(S^{3},T,a)\to(S^{3},T_{g},a_{1})$ that is unique up to eyeglass equivalence.
(In [FS1] it was shown to be true up to Powell equivalence, but this stricter
outcome relied on the stricter inductive assumption that the Powell conjecture
was known for splittings of lower genus.)

Let ($\mbox{$\mathfrak{b}$},T_{\mbox{$\mathfrak{b}$}}$) denote the induced
Heegaard splitting in the bubble, of genus $\leq g-1$. Suppose $a,a^{\prime}$
are primitive meridians for $A$ lying in $\mathfrak{b}$. There is a
homeomorphism
$(\mbox{$\mathfrak{b}$},T_{\mbox{$\mathfrak{b}$}})\to(\mbox{$\mathfrak{b}$},T_{\mbox{$\mathfrak{b}$}})$
that carries $a$ to $a^{\prime}$; by inductive assumption the homeomorphism is
in the eyeglass group. So a choice of any primitive meridian for $A$ in
$\mathfrak{b}$ determines the same eyeglass equivalence class of
homeomorphisms $(S^{3},T)\to(S^{3},T_{g})$, verifying the claim.

Claim 2: Suppose $\mathfrak{b}$ and $\mbox{$\mathfrak{b}$}^{\prime}$ are
disjoint bubbles for $T$ and $h,h^{\prime}:(S^{3},T)\to(S^{3},T_{g})$ are
homeomorphisms determined by the bubbles, as described in Claim 1. Then
$h,h^{\prime}$ are eyeglass equivalent.

Pick primitive meridian disks $a,a^{\prime}$ in $\mathfrak{b}$ and
$\mbox{$\mathfrak{b}$}^{\prime}$ respectively. Consider a homeomorphism
$\hat{h}:(S^{3},T)\to(S^{3},T_{g})$ so that
$\hat{h}(a)=a_{1},\hat{h}(a^{\prime})=a_{g}$. By construction $\hat{h}$ is in
the eyeglass equivalence class defined by $\mathfrak{b}$; if we post-compose
$\hat{h}$ with the homeomorphism $\sigma:(S^{3},T_{g})\to(S^{3},T_{g})$ that
exchanges the monobubbles $(a_{1},b_{1})$ and $(a_{g},b_{g})$ we get a
homeomorphism in the eyeglass equivalence class defined by
$\mbox{$\mathfrak{b}$}^{\prime}$. But a monobubble exchange is in the eyeglass
group, so the two are eyeglass equivalent, verifying the claim.

Having now verified the claims, select a primitive $a$ from the bubble
$\mathfrak{b}$ and consider the homeoomorphisms
$h:(S^{3},T,a)\to(S^{3},T_{g},a_{1})$ and
$h\tau^{-1}:(S^{3},T,\tau(a))\to(S^{3},T_{g},a_{1})$. The former is a
homeomorphism determined by $\mathfrak{b}$ and the latter by
$\tau(\mbox{$\mathfrak{b}$})$. Following the claims, the sequence of bubbles
$\mbox{$\mathfrak{b}$}=\mbox{$\mathfrak{b}$}_{1},...,\mbox{$\mathfrak{b}$}_{m}=\tau(\mbox{$\mathfrak{b}$})$
shows that the two are eyeglass equivalent. Hence $\tau$ must be in the
eyeglass group. ∎

##  12\. Eyeglass twists are stably Powell

We briefly recall the conventions and definitions from [FS1], filling out and
reviewing the conventions introduced in the proof of Theorem 11.8: Let
$T_{g}\subset S^{3}$ be the standard genus $g$ Heegaard surface in $S^{3}$,
dividing $S^{3}$ into the genus $g$ handlebodies $A$ and $B$. Let
$\\{a_{1},...,a_{g}\\}$ and $\\{b_{1},...,b_{g}\\}$ be orthogonal sets of
disks in $A$ and $B$ respectively. We can view each pair
$\mbox{$\partial$}a_{i},\mbox{$\partial$}b_{i}$ as respectively the meridian
and longitude of one of the standard summands of $T_{g}$. Let
$\\{c_{1},...,c_{g-1}\\}$ be the disjoint separating circles on $T_{g}$ shown
in [FS1, Figure 2], with each $c_{i}$ separating the first $i$ standard
summands from the last $g-i$ standard summands. Note that each $c_{i}$ bounds
a disk in both $A$ and $B$ and so defines a reducing sphere $S_{i}$ for
$T_{g}$.

######  Definition 12.1.

Any finite composition of Powell generators, illustrated in [FS1, Fig. 1],
will be called a Powell move.

######  Proposition 12.2.

Suppose $\tau$ is an eyeglass twist whose frame $\eta$ is disjoint from
$c_{1}$. Then $\tau$ is a Powell move.

###### Proof.

Let $\\{\ell_{a},\ell_{b},v\\}$ be the frame of $\eta$. Let $u$ be an arc such
that

  * •

the ends of $u$ lie on $a_{1}$ near the point $a_{1}\cap b_{1}$ and on
$\ell_{a}$ near the point $\ell_{a}\cap v$,

  * •

$u$ crosses $c_{1}$ once and

  * •

$u$ is otherwise disjoint from $a_{1},b_{1}$, and $\eta$.

See Figure 29.

2pt $\ell_{a}$ at 110 145 $\ell_{b}$ at 240 100 $v$ at 190 145 $u$ at 265 170
$a_{1}$ at 320 190 $b_{1}$ at 435 185 $c_{1}$ at 365 150

Figure 29.

Let $\eta^{\prime}$ be the eyeglass twist whose frame is
$\\{a_{1},\ell_{b},u\cup v\\}$. We know from [FS1, Lemma 3.4] that an eyeglass
twist $\tau^{\prime}$ along $\eta^{\prime}$ is a Powell move. Let
$\ell^{\prime}_{a}$ be the band sum of $\ell_{a}$ and $a_{1}$ along $u$ and
observe, by watching the motion of $\ell_{b}$, that the composition
$\tau\tau^{\prime}$ is an eyeglass twist $\tau_{+}$ whose frame $\eta_{+}$ is
$\\{\ell^{\prime}_{a},\ell_{b},v\\}$.

2pt $\eta_{+}$ at 190 125 $\eta_{1}$ at 265 175 $b_{1}$ at 440 185

Figure 30.

Now let $\eta_{1}$ be the eyeglass whose frame is $\\{\ell_{a},b_{1},u\\}$ and
$\tau_{1}$ be the eyeglass twist along $\eta_{1}$. See Figure 30. Again, from
[FS1, Lemma 3.4], $\tau_{1}$ is a Powell move. Further observe that
$\tau^{-1}_{1}(\eta_{+})$ is a frame $\eta^{\prime}_{+}$ with lenses $a_{1}$
and $\ell_{b}$ and a bridge that intersects $c_{1}$ in a single point, so
again an eyeglass twist $\tau^{\prime}_{+}$ along $\eta^{\prime}_{+}$ is a
Powell move. See Figure 31.

Since $\tau^{-1}_{1}(\eta_{+})=\eta^{\prime}_{+}$ it follows that
$\tau_{+}=\tau_{1}\tau^{\prime}_{+}\tau^{-1}_{1}$. As a composition of Powell
moves, $\tau_{+}$ is a Powell move. But we have earlier shown that
$\tau=\tau_{+}\tau^{\prime-1}$ so, as a composition of Powell moves, $\tau$ is
a Powell move. ∎

2pt $\tau^{-1}_{1}(\eta_{+})$ at 150 470 $\eta^{\prime}_{+}$ at 265 180

Figure 31.

######  Corollary 12.3.

Suppose $\tau\in G(S,T)$ and $\tau:(S^{3},T)\to(S^{3},T)$ leaves the curve
$c_{1}$ invariant. Then $\tau$ is a Powell move.

## References

  * [Bo] F Bonahon, Cobordism of automorphisms of surfaces, Ann. Sci. École Norm. Sup.(4) 16 1983) 237–270
  * [CG] A. Casson and C. Gordon, Reducing Heegaard splittings, Topology and its applications, 27 (1987) 275-283.
  * [FS1] M. Freedman and M. Scharlemann, Powell moves and the Goeritz group, Arxiv preprint 1804.05909
  * [FS2] M. Freedman and M. Scharlemann, Uniqueness in Haken’s Theorem, Arxiv preprint 2004.07385
  * [Ha] A Hatcher, Homeomorphisms of sufficiently large $P^{2}$-irreducible $3$-manifolds, Topol- ogy 15 (1976), 343–347.
  * [JM] J. Johnson and D. McCullough, The space of Heegaard splittings J. Reine Angew. Math 679 (2013) 155-179
  * [La] M. Lackenby, An algorithm to determine the Heegaard genus of simple 3–manifolds with nonempty boundary Alg. Geom. Top. 8 (2008) 911-934.
  * [Po] J. Powell, Homeomorphisms of $S^{3}$ leaving a Heegaard surface invariant, Trans. Amer. Math. Soc. 257 (1980) 193–216.
  * [Sc1] M. Scharlemann, A Strong Haken’s Theorem , Arxiv preprint 2003.08523
  * [Sc2] M. Scharlemann, One Powell generator is redundant, to appear Proc. Amer. Math. Soc., arxiv preprint 1908.00479.
  * [Wa] F. Waldhausen, Heegaard-Zerlegungen der 3-Sphäre, Topology 7 (1968) 195–203.

## Index

  * Aligned disk set §3
  * Balanced sphere Definition 8.6
  * Blank bullseye Definition 7.1
  * Bubble §4
  * Bubble exchange item 3
  * Bubble move item 1
  * Bullseye Definition 7.1
  * Certify; certificate Definition 8.1
  * Chamber complex §2
  * Decomposition of chamber complex §2
  * Disk set §2
  * Disky §2
  * Empty chamber Definition 5.1
  * Eyeglass item 4
  * Eyeglass group Definition 4.1
  * Eyeglass twist item 4
  * Flag rubric §6
  * Flagged chamber complex Definition 5.1
  * Flip item 2
  * Goneball §2
  * Guiding sphere Definition 8.3
  * Heegaard split chamber complex §3
  * interacting flagged chamber complex decompositions Definition 8.1
  * Occupied chamber Definition 5.1
  * Powell group Definition 4.1
  * Remnant of chamber §2
  * Scars §2
  * Tiny chamber complex Definition 2.3


contributing to the entanglement entropy $S$, and $(1-\delta)$ be the
probability for the exclusive event, namely, an entangled pair staying inside
the subsystem A. Then $S$ is given by the total ‘successful’ events for a
fixed, but large, number of trials, resulting in a binomial distribution. Now,
for $L_{A}/L=1/2$ an entangled pair has equally the same probability to cross
subsystem’s boundary or stay inside A, i.e., $\delta=1/2$; it is well known
that in this case the binomial distribution reduces to a Gaussian
distribution. For small $L_{A}/L$, $\delta\approx 1$ and thus for an entangled
pair to stay in A is a rare event. This results in a heavier lower tail of the
$S$ distribution which is exponential, in the same fashion as the binomial
distribution approaches its Poissonian limit. So this picture well explains
the asymptotic behaviors of the tail bounds and their evolution with
$L_{A}/L$. This suggests that the bounds yield the form of large deviation
probability for the upper and the lower tail, respectively. That is,
inequalities (89) and (93) can be promoted to equalities, i.e., for
sufficiently large $\epsilon$,

$\displaystyle\mathbf{P}(|S-\langle
S\rangle|\geq\epsilon)=\left\\{\begin{array}[]{ll}e^{-{\epsilon^{2}\over
2\mathfrak{b}_{+}}},&{\rm for}\,S-\langle S\rangle>0\\\ e^{-{\epsilon^{2}\over
2(\mathfrak{b}_{-}+\mathfrak{c}\epsilon)}},&{\rm for}\,S-\langle
S\rangle<0\end{array}\right.,$ (152)

with the parameters $\mathfrak{b}_{\pm}\propto b_{\pm}$ and
$\mathfrak{c}\propto|c_{-}|$ and the proportionality coefficients being
absolute constants. This gives Eq. (1) in the main text, which was confirmed
by simulations (see Fig. 2(c) in the main text). In Sec. VII.2 we provide
further simulation results to show that Eq. (152) holds also for the second-
order R$\acute{\rm e}$nyi entropy $S_{2}$.

Let us further discuss semi-quantitatively how the rate of the exponential
decay $\sim e^{-{\epsilon\over
2\mathbf{c}}}\,(\epsilon\gg\mathfrak{b}_{-}/\mathfrak{c})$ of the lower tail
behaves, when $L$ approaches the infinite and the ratio $L_{A}/L$ is fixed.
First of all, Eq. (148) implies that ${\rm Var}(S)\sim L$ for sufficiently
large $L$ and fixed $L_{A}/L$. Then, expanding $dG/du$ in $u$ we obtain
$dG/du=u{\rm Var}(S)+{\cal O}(u^{2})$. So the coefficient of the leading term
is $\propto L$, and thus the leading expansion of $L^{-1}dG/du$ is well
defined in the limit $L\rightarrow\infty$. It is natural to expect that the
entire expansion of $L^{-1}dG/du$ is well defined in such limit, although this
is difficult to prove. This implies that all the coefficients of higher order
terms cannot grow faster than $L$. Next, according to Eq. (80), $\delta F(u)$
is a sum of $N\sim L$ terms. Thus we estimate $\delta F(u)\sim L$. Finally,
combining the estimations for $dG/du$ and $\delta F(u)$ with the inequality
(92), we find that $\mathfrak{c}={\cal O}(L^{0})$, i.e., depends only on the
ratio $L_{A}/L$. This is in agreement with the data in Table 2.

##  VII Numerical simulations

In this section we provide a complete description of numerical simulations and
report extended numerical results.

###  VII.1 The entanglement entropy $\boldsymbol{S}$

We first diagonalize numerically for each time $t$ the correlation matrix
$C(t)$ of Eq. (3) to obtain its eigenvalues
$\\{p_{\nu}(t)\\}_{\nu=1}^{2L_{A}}$. We then substitute them into Eq. (23) to
obtain the instantaneous $S(t)$. Upon varying $t$, the pattern $S(t)$ shown in
Figs. 1(a) and 1(b) of the main text is obtained. The time is in unit of
$\hbar/J$, where $J$ is the amplitude of hopping between nearest ${\bar{A}}$\-
and ${\bar{B}}$-site.

To obtain the statistics of the random entanglement entropy
$S(\boldsymbol{\varphi})$ we simulate an ensemble of random correlation
matrices $\tilde{C}(\boldsymbol{\varphi})$. Specifically, we generate $N$
random angles $\varphi_{m}$ ($m=0,1,\ldots,N-1$) uniformly drawn from the
interval $[0,2\pi)$, which gives a disorder realization
$\boldsymbol{\varphi}=(\varphi_{0},\varphi_{1},\ldots,\varphi_{N-1})$, and
substitute them into Eq. (15). For each $\boldsymbol{\varphi}$ we diagonalize
numerically $\tilde{C}(\boldsymbol{\varphi})$ to obtain the eigenvalue
spectrum $\\{p_{\nu}(\boldsymbol{\varphi})\\}_{\nu=1}^{2L_{A}}$. Analogous to
Eq. (23), we obtain the corresponding value of $S(\boldsymbol{\varphi})$. The
statistics of $S(\boldsymbol{\varphi})$ is obtained for an ensemble consisting
of $5\times 10^{5}$ disorder realizations. The distribution of
$S(\boldsymbol{\varphi})$ is shown by the dashed line in Fig. 1(a) of the main
text.

When the eigenfrequencies $\omega_{0},\omega_{1},\ldots,\omega_{N-1}$ are
incommensurate, generic for the energy eigenspectrum, the statistics of the
time series $S(t)$ and the random function $S(\boldsymbol{\varphi})$ are
equivalent, as shown by Eq. (45). In Fig. 1(a) of the main text, the dash line
and the histogram agree with each other, independent of the time interval used
in the sampling of $S(t)$. When $\omega_{0},\omega_{1},\ldots,\omega_{N-1}$
are commensurate, this statistical equivalence does not hold; see Sec. VIII.

###  VII.2 The $\boldsymbol{n}$th-order R${\acute{\rm{\bf e}}}$nyi entropy
$\boldsymbol{S_{n}}$

As an example for more general entanglement probes discussed in Sec. II.2, we
numerically compute the $n$th-order R$\acute{\rm e}$nyi entropy according to
Eqs. (38) and (39). We find that its behaviors are similar to those of the
entanglement entropy. As shown in Fig. 1, first, after initial growth and
damped oscillations, $S_{2}(t)$ displays quasiperiodic oscillations,
fluctuating around its average value, see (a) and (b); second, at fixed
$L_{A}$ the distribution is broadened as $L$ decreases (c), and its upper and
lower tail are sub-Gaussian and sub-Gamma, respectively (d). Furthermore, as
shown in Fig. 2, simulations confirm the relation (107) for $S_{n}$, for
$n=2,3,10$. Even more surprisingly, simulations show that the scaling law
(148) holds also for the variance ${\rm Var}(S_{n})$ (b-d). These findings
confirm that the theory developed in this work applies to general entanglement
probes, not limited to the entanglement entropy.

Figure 1: Numerical simulations of the time evolution of the second-order
R$\acute{\rm e}$nyi entropy $S_{2}(t)$ upon quench. The patterns (red) in
different time windows are shown in (a,b), where the patterns of $S(t)$ (grey,
same as Fig. 1 in the main text) are also shown for comparison. The evolution
of the distribution with $L$ at fixed $L_{A}$ is shown in (c). For different
ratios $L_{A}/L=0.1$ (yellow), $0.2$ (green) and $0.5$ (blue), numerical data
for the upper (squares) and lower (circles) tail are well fitted by the form
given by the first and the second line of Eq. (152) (dashed lines),
respectively (d). The subsystem size $L_{A}=25$ and for (a,b) the total system
size $L=124$.

###  VII.3 Numerical computation of $\boldsymbol{(b^{\pm},c^{\pm})}$

In Sec. V.1, we introduce the variance factor and the scale parameter
$(b^{\pm},c^{\pm})$ in solving the modified logarithmic Sobolev inequality
(74) for the logarithmic moment-generating function $G(u)$ of
$S(\boldsymbol{\varphi})$. The $+$ sign corresponds to the upper tail
distribution and $-$ to the lower. Here we describe the numerical method for
computing these parameters summarized in Table 2. For $b^{\pm}$, following Eq.
(76):

  * •

We realize a disorder realization $\boldsymbol{\varphi}$ and compute
$S(\boldsymbol{\varphi})$, in precisely the same way as what is described
above.

  * •

With respect to this $\boldsymbol{\varphi}$, we vary its $m$-th component
$\varphi_{m}$ over $[0,2\pi)$ while keeping other components, i.e.,
$(\varphi_{0},\ldots,\varphi_{m-1},\varphi_{m+1},\ldots,\varphi_{N-1})$ fixed,
and subsequently find the minimum (maximum) of $S(\boldsymbol{\varphi})$,
which determines $S_{m}^{\pm}$.

  * •

With the same $\boldsymbol{\varphi}$, we repeat the second step for each
component and obtain $N$ values of $S_{m}^{\pm}$, and then determine
$\sum_{m=0}^{N-1}(S-S_{m}^{\pm})^{2}$.

  * •

We repeat the previous three steps for $2\times 10^{3}$ disorder realizations
$\boldsymbol{\varphi}$ and average the outcome.

To obtain the scale parameter $c^{\pm}$, we numerically evaluate the two
quantities $\delta F(u)$ and $dG(u)/du$ for a range of $u$ around the origin.
$c_{+}$ is determined by the maximum of $\delta F/(dG/du)$ in the $u>0$
interval, whereas $c_{-}$ is determined by the corresponding minimum in the
$u<0$ interval.

Figure 2: Simulation results for the variance of $n$th-order R${\acute{\rm{\bf
e}}}$nyi entropy $S_{n}$, $n=2,3,10$. For each $S_{n}$, we vary either $L_{A}$
or $L$, while keeping the other fixed, in the same way as we generate the data
for Fig. 3 in the main text. (a) The results confirm the relation (107).
(b)-(d) They confirm the first (b) and the second (c) terms, and the overall
(d) of the scaling law (148). All theoretical predictions are presented by
dashed lines. In (d) ${\rm Var}(S_{n}),\,L,\,L_{A}$ are all rescaled in the
manner similar to what is performed in the scaling analysis of ${\rm Var}(S)$.

###  VII.4 Spectral statistics of the correlation matrix

We provide more details of the nearest-neighbor spacing statistics obtained in
Fig. 1(c) in the main text. In Fig. 3(a) we show the zoom-out view of the time
evolving spectrum of the correlation matrix $C(t)$. In Fig. 3(b) we show the
density of states of the full spectrum, with the shaded area indicating the
spectral window used to obtain the nearest-neighbor spacing distribution.

Figure 3: (a) Zoom-out view of the evolving spectrum of $C(t)$. (b) The
density of states (DOS) of $C(t)$ shown in arbitrary unit (a.u.) scale. The
spectrum lies within the interval $[0,1]$ and it is particle-hole symmetric.
The shaded area is the spectral window where we perform the level statistics
shown in Fig. 1(c) of the main text.

###  VII.5 Quench parameters

In the generic numerical simulations for the Rice-Mele (RM) model, we use the
quench parameters $(J,J^{\prime},M):(1,0.5,0.5)\rightarrow(1,1.5,1.5)$. For
TFIC we use the quench parameter $h=3\rightarrow 5$, with $h$ being the
external magnetic field. A detailed account of the dynamics of entanglement in
TFIC is given in Appendix A. In Fig. 3 of the main text, when obtaining
results of various quenches, we use:

  * •

for RM model $(J,J^{\prime},M)$ quench I:
$(1,0.7,0.3)\rightarrow(0.3,1.0,0.001)$; II:
$(0.6,1.0,0.8)\rightarrow(1.0,0.2,10)$.

  * •

for TFIC quench I: $h=3\rightarrow 2$; II: $h=3\rightarrow 5$; III:
$h=5\rightarrow 2$.

Figure 4: (a) Numerical simulations confirm that, when $\boldsymbol{\omega}$
is commensurate, the pattern of the entanglement entropy is periodic and
displays complete revival at multiple periods (yellow curves). The pattern in
the incommensurate $\boldsymbol{\omega}$ case is shown for comparison (blue
curves, same as Fig. 1 in the main text). Inset of (a) shows the energy
spectrum with commensurate energy gap $\omega_{m}=\omega_{0}+m\,r$,
$\omega_{0}=3.162$, $r=0.043$, $m=0,\ldots,62$. The dash lines are the
physical energy spectrum (159). The total system size $L=124$ and the
subsystem size $L_{A}=25$. (b) The spectrum of the corresponding evolving
correlation matrix with a zoom-in view in (c). (d) The nearest-neighbor
spacing distribution and the DOS (inset). The shaded area in the inset is the
spectral window used for statistical analysis.

##  VIII Entanglement evolution with commensurate $\boldsymbol{\omega}$

So far we have studied the dynamics of entanglement where the frequencies
governing the evolution of the correlation matrix $C(t)$ are incommensurate,
i.e., $\boldsymbol{\omega}$ is incommensurate. In this section we study the
commensurate case. In this case, by Eq. (14), $C(t)$ is periodic in time. Thus
all entanglement probes defined by Eq. (32) are also periodic in time, and
thus display complete revival to the initial values at multiple periods.

Since it is impossible to tune parameters of the physical Rice-Mele model to
obtain an energy spectrum with commensurate frequencies, we modify directly
the frequencies in the dynamical phases. Specifically, at the Bloch momentum
$k_{m}$ the frequency $\omega_{m}$ is set to $\omega_{0}+m\,r$, with $r$ being
a suitably chosen rational number and $\omega_{0}$ being the zero-momentum
eigenfrequency in the incommensurate case. In this way, the resulting spectrum
simulates as closely as the physical model (see inset of Fig. 4(a)). We assume
the eigenstate properties are not modified in any qualitative way.

In Fig. 4(a), the result for the entanglement entropy evolution $S(t)$
confirms that it is periodic and displays complete revival to the initial
value $S(0)$ at multiple periods (yellow curve). The time profile is
completely different from the incommensurate case (blue curve), although in
both cases the system has the same degrees of freedom.

We also show the evolving spectrum of $C(t)$ in Figs. 4(b) and 4(c), which
clearly show a high degree of regularity. Finally, this regularity results in
a very different nearest-neighbor level statistics shown in Fig. 4(d), in
stark contrast to the incommensurate case giving the semi-Poissonian
statistics shown in Fig. 1(c) of the main text.

##  Appendix A Entanglement evolution in transverse field Ising chain

The transverse field Ising chain Hamiltonian is given by

$\displaystyle
H=-\frac{1}{2}\sum_{i=1}^{L}(\sigma_{i}^{x}\sigma_{i+1}^{x}+h\sigma_{i}^{z}),$
(153)

with the periodic boundary condition
$\sigma_{L+1}^{\alpha}=\sigma_{1}^{\alpha}$, and the Pauli matrices
$\sigma^{x,z}_{i}$ act on the spin at site $i$ belonging to a spin chain of
length $L$. $h$ is the external magnetic field. We also carried a parallel
study of the entanglement entropy evolution of TFIC with the use of the
numerical solution of the correlation matrix. The method of finding that
solution is described in standard textbooks and we simply use analytical
expressions given in Ref. [1]. The key formulae for the numerical computations
are Eqs. (3.3) – (3.6) in that work. The expressions for finite $L$ simply
amount to converting the continuous momentum variable $\varphi$ in those
formulae back to a discrete one. That is, for a generic expression
$F(\varphi)$, we have

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}d\varphi\,F(\varphi)\,\,\,\rightarrow\,\,\,\frac{1}{L}\sum_{n=-L/2}^{L/2-1}\,F\left(\left(n+{1\over
2}\right){2\pi\over L}\right).$ (154)

Note that on the right-hand side the Bloch momenta are quantized at half-
integers. This arises from the antiperiodic boundary condition for the Jordan-
Wigner fermions in the even fermion number sector. For studies of the time-
dependent correlation matrix, we stay in the paramagnetic phase with $h>1$
when quenching the magnetic field. Moreover, the pre-quench state is the
ground state of the pre-quench Hamiltonian.

##  Appendix B Formulae for correlation matrix of the Rice-Mele model

Here we give the expressions for the elements of
$\check{\gamma}(k),\check{\alpha}(k),\check{\beta}(k)$. In the momentum space,
the Hamiltonian (1) reduces to a $2\times 2$ Bloch Hamiltonian which reads

$\displaystyle H(k)=\boldsymbol{B}(k)\cdot\boldsymbol{\sigma}$ (155)

where $\boldsymbol{B}(k)=(B_{x}(k),B_{y}(k),B_{z}(k))$ with

$\displaystyle
B_{x}(k)=-(J+J^{\prime})\cos(k/2),\,\,\,B_{y}(k)=(J-J^{\prime})\sin(k/2),\,\,\,B_{z}(k)=M,$
(156)

and $\boldsymbol{\sigma}=\\{\sigma^{x},\sigma^{y},\sigma^{z}\\}$, with
$\sigma^{x,y,z}$ being the Pauli matrices. The $L$ Bloch momenta are

$\displaystyle k=0,\pm{2\pi m\over L},\quad m=1,2,\ldots,{L-1\over 2}$ (157)

for odd $L$ and

$\displaystyle k=0,\pi,\pm{2\pi m\over L},\quad m=1,2,\ldots,{L\over 2}-1$
(158)

for even $L$.

The energy spectrum of the Bloch Hamiltonian (155) is

$\displaystyle\pm E_{k}=\pm|\boldsymbol{B}(k)|=\pm\sqrt{M^{2}+J^{2}+J^{\prime
2}+2JJ^{\prime}\cos k}.$ (159)

The $-$ ($+$) sign corresponds to the lower (upper) band. Note that this
spectrum exhibits the particle-hole symmetry. Moreover, $E_{k}$ is symmetric
with respect to $k=0$, i.e.,

$\displaystyle E_{k}=E_{-k}.$ (160)

That is, the spectrum bears the reflection symmetry. The eigenstates are
$SU(2)$ spinors. The eigenstate corresponding to $E_{k}$ is

$\displaystyle\left(\begin{array}[]{c}u^{-}_{k}\\\ v^{-}_{k}\\\
\end{array}\right)=\left(\begin{array}[]{c}\sin\frac{\theta_{k}}{2}\\\
-e^{i\phi_{k}}\cos\frac{\theta_{k}}{2}\\\ \end{array}\right)$ (165)

and to $-E_{k}$ is

$\displaystyle\left(\begin{array}[]{c}u^{+}_{k}\\\ v^{+}_{k}\\\
\end{array}\right)=\left(\begin{array}[]{c}\cos\frac{\theta_{k}}{2}\\\
e^{i\phi_{k}}\sin\frac{\theta_{k}}{2}\\\ \end{array}\right),$ (170)

where
$\cos\theta_{k}=M/E_{k},\tan\phi_{k}=B_{y}(k)/B_{x}(k)=-{J-J^{\prime}\over
J+J^{\prime}}\tan({k\over 2})$. It is obvious that

$\theta_{k}=\theta_{-k},\qquad\phi_{k}=-\phi_{-k}.$ (171)

With the help of $\Psi(t)=e^{-iH_{f}t/\hbar}\Psi(0)$ as well as Eqs. (165) and
(170), the matrix elements of the correlation matrix
$C_{i\sigma,i^{\prime}\sigma^{\prime}}(t)$ can be computed to give

$\displaystyle
C_{i{\bar{A}},i^{\prime}{\bar{A}}}(t)=\frac{1}{L}\sum_{k}\,e^{ik(i-i^{\prime})}\,|\tilde{u}_{fk}(t)|^{2}=\delta_{ii^{\prime}}-C_{i{\bar{B}},i^{\prime}{\bar{B}}}(t),$
(172) $\displaystyle
C_{i{\bar{A}},i^{\prime}{\bar{B}}}(t)=\frac{1}{L}\sum_{k}\,e^{ik(i-i^{\prime}+{1\over
2})}\,\tilde{u}^{*}_{fk}(t)\tilde{v}_{fk}(t)=C^{*}_{i^{\prime}{\bar{B}},i{\bar{A}}}(t)$
(173)

with the summation over all Bloch momenta, since we take the initial state as
the half-filling ground state of the pre-quench Hamiltonian. The time-
dependent form factors in Eqs. (172) and (173) are given by

$\displaystyle\tilde{u}_{fk}(t)=a_{k}\,u_{fk}^{-}\,e^{iE_{fk}t}+b_{k}\,u_{fk}^{+}\,e^{-iE_{fk}t}$
(174)

and

$\displaystyle\tilde{v}_{fk}(t)=a_{k}\,v_{fk}^{-}\,e^{iE_{fk}t}+b_{k}\,v_{fk}^{+}\,e^{-iE_{fk}t}.$
(175)

The coefficients $a_{k}$ ($b_{k}$) are overlaps between the half-filling
ground state of the pre-quench Hamiltonian $H_{0}$ and the lower (upper) band
eigenstates of the post-quench Hamiltonian $H_{f}$, with

$\displaystyle
a_{k}=(u^{-*}_{fk},v^{-*}_{fk})\left(\begin{array}[]{c}u^{-}_{0k}\\\
v^{-}_{0k}\end{array}\right),\,\,\,\,b_{k}=(u^{+*}_{fk},v^{+*}_{fk})\left(\begin{array}[]{c}u^{-}_{0k}\\\
v^{-}_{0k}\end{array}\right),$ (180)

where the subscript $0$ denotes the eigenstate of $H_{0}$.

For $\sigma=\sigma^{\prime}$, with the help of Eq. (174) we find that

$\displaystyle|\tilde{u}_{fk}(t)|^{2}=|a_{k}|^{2}(u_{fk}^{-})^{2}+|b_{k}|^{2}(u_{fk}^{+})^{2}+u_{fk}^{+}u_{fk}^{-}(a_{k}b_{k}^{*}+a_{k}^{*}b_{k})\cos(2E_{fk}t)+iu_{fk}^{+}u_{fk}^{-}(a_{k}b_{k}^{*}-a_{k}^{*}b_{k})\sin(2E_{fk}t).$
(181)

Inserting it into Eq. (172) we get

$\displaystyle\gamma_{{\bar{A}}{\bar{A}}}(k)=1-\gamma_{{\bar{B}}{\bar{B}}}(k)=|a_{k}|^{2}(u_{fk}^{-})^{2}+|b_{k}|^{2}(u_{fk}^{+})^{2},$
(182)
$\displaystyle\alpha_{{\bar{A}}{\bar{A}}}(k)=-\alpha_{{\bar{B}}{\bar{B}}}(k)=u_{fk}^{+}u_{fk}^{-}(a_{k}b_{k}^{*}+a_{k}^{*}b_{k}),$
(183)
$\displaystyle\beta_{{\bar{A}}{\bar{A}}}(k)=-\beta_{{\bar{B}}{\bar{B}}}(k)=iu_{fk}^{+}u_{fk}^{-}(a_{k}b_{k}^{*}-a_{k}^{*}b_{k}).$
(184)

For $\sigma\neq\sigma^{\prime}$, with the help of Eqs. (174) and (175) we find
that

$\displaystyle\tilde{u}^{*}_{fk}(t)\tilde{v}_{fk}(t)=e^{i\phi_{fk}}\biggl{(}|a_{k}|^{2}u_{fk}^{-}{v^{\prime}}_{fk}^{-}+|b_{k}|^{2}u_{fk}^{+}{v^{\prime}}_{fk}^{+}+(a_{k}^{*}b_{k}u_{fk}^{-}{v^{\prime}}_{fk}^{+}+a_{k}b_{k}^{*}u_{fk}^{+}{v^{\prime}}_{fk}^{-})\cos(2E_{fk}t)$
$\displaystyle\textrm{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
}-i(a_{k}^{*}b_{k}u_{fk}^{-}{v^{\prime}}_{fk}^{+}-a_{k}b_{k}^{*}u_{fk}^{+}{v^{\prime}}_{fk}^{-})\sin(2E_{fk}t)\,\biggr{)},$
(185)

where $v_{fk}^{\pm}=e^{i\phi_{fk}}\,{v^{\prime}}_{fk}^{\pm}$, and $\phi_{fk}$
is defined in the same way as $\phi_{k}$ but with post-quench Hamiltonian
parameters. Inserting it into Eq. (173) we get

$\displaystyle\gamma_{{\bar{A}}{\bar{B}}}(k)=\gamma_{{\bar{B}}{\bar{A}}}(k)^{*}=e^{i\delta(k)}\left(|a_{k}|^{2}u_{fk}^{-}{v^{\prime}}_{fk}^{-}+|b_{k}|^{2}u_{fk}^{+}{v^{\prime}}_{fk}^{+}\right),$
(186)
$\displaystyle\alpha_{{\bar{A}}{\bar{B}}}(k)=\alpha_{{\bar{B}}{\bar{A}}}(k)^{*}=e^{i\delta(k)}\left(a_{k}^{*}b_{k}u_{fk}^{-}{v^{\prime}}_{fk}^{+}+a_{k}b_{k}^{*}u_{fk}^{+}{v^{\prime}}_{fk}^{-}\right),$
(187)
$\displaystyle\beta_{{\bar{A}}{\bar{B}}}(k)=\beta_{{\bar{B}}{\bar{A}}}(k)^{*}=-ie^{i\delta(k)}\left(a_{k}^{*}b_{k}u_{fk}^{-}{v^{\prime}}_{fk}^{+}-a_{k}b_{k}^{*}u_{fk}^{+}{v^{\prime}}_{fk}^{-}\right),$
(188)

and $\delta(k)\equiv k/2+\phi_{fk}$. We note that $u_{fk}^{\pm}$,
${v^{\prime}}_{fk}^{\pm}$ are real, whereas $a_{k}$, $b_{k}$, and
$v_{{fk}}^{\pm}$ are generally complex except for $k=0$ when they are also
real. Furthermore, $u_{fk}^{\pm}$, ${v^{\prime}}_{fk}^{\pm}$ are even in $k$,
$\delta(k)$ is odd in $k$, and $a_{k}^{*}=a_{-k}$, $b_{k}^{*}=b_{-k}$.

##  Appendix C Some general remarks on quantum expectation value of operators
upon a quench

When a generic isolated many-body system with Hamiltonian $H_{f}$ is driven
out of equilibrium, under Schrödinger evolution it evolves from an initial
state $\Psi(0)$ to a state $\Psi(t)$ at time $t$, given by

$\displaystyle\Psi(t)=e^{-iH_{f}t}\,\Psi(0).$ (189)

In the basis of many-body eigenstates $\Psi_{\bf m}$ (labelled by m), Eq.
(189) is written as

$\displaystyle\Psi(t)=\sum_{\bf m}e^{-i\mathcal{E}_{\bf m}t}\,w_{\bf
m}\,\Psi_{\bf m},$ (190)

where $\mathcal{E}_{\bf m}$’s are the many-body eigenenergies and $w_{\bf
m}$’s are the superposition coefficients of $\Psi(0)$. Under this evolution,
the expectation value of a generic operator ${\cal A}$ at time $t$ is

$\displaystyle\langle\Psi(t)|{\cal A}|\Psi(t)\rangle=\sum_{{\bf m}{\bf
m}^{\prime}}e^{-i(\mathcal{E}_{\bf m}-\mathcal{E}_{{\bf
m}^{\prime}})t}\,w_{{\bf m}^{\prime}}^{*}\,w_{\bf m}\,\langle\Psi_{{\bf
m}^{\prime}}|{\cal A}|\Psi_{\bf m}\rangle.$ (191)

So the time parameter enters the evolution $\langle\Psi(t)|{\cal
A}|\Psi(t)\rangle$ through a large number $\sim{\cal L}^{2}$ of dynamical
phases: $(\mathcal{E}_{\bf m}-\mathcal{E}_{{\bf m}^{\prime}})t$ with
$\mathcal{E}_{\bf m}\neq\mathcal{E}_{{\bf m}^{\prime}}$. Here ${\cal L}$ is
the number of many-body eigenstates superposing the initial state.
Investigations of the emergence of statistical mechanics from such evolving
quantum expectation values were initiated in Ref. [31]. Nowadays it is
accepted [32–35] that the quantum expectation can equilibrate after long time,
when out-of-equilibrium fluctuations around the equilibrium value are small.
However, despite of some progresses it remains a difficult problem to
demonstrate under what circumstances this scenario would emerge from quantum
dynamics of an isolated system.

Compared to Eq. (191), we find that the correlation matrix Eq. (3) is
simplified substantially. The general reasons underlying this simplification
are as follows. First, the Rice-Mele model is composed of noninteracting
fermions, and many-body effects enter through the so-called exchange
interaction, namely, the particle indistinguishability [36]. Second, the
system is driven out of equilibrium by global quench. For the first reason
each many-body eigenstate ${\bf m}$ corresponds to a configuration,
$\\{n^{\zeta}_{k}\\}$, of the occupation number $n^{\zeta}_{k}=0,1$ at the
single-particle eigenstate $(u^{\zeta}_{fk},v^{\zeta}_{fk})^{T}$ of $H_{f}$,
with $\zeta=+\,(-)$ denoting the particle (hole) band and the superscript $T$
denoting the transpose (see Appendix B). Because the quench is global and the
pre-quench state $\Psi(0)$ is a half-filling ground state, the configurations
corresponding to the many-body eigenstates of $H_{f}$ superposing $\Psi(t)$ —
not all the many-body eigenstates — must satisfy

$n^{+}_{k}+n^{-}_{k}=1,\quad{\rm for\,\,every\,\,Bloch\,\,momentum\,}\,k,$
(192)

and the corresponding many-body eigenenergies are given by

${\cal E}_{\bf
m}=\sum_{k}(n^{+}_{k}-n^{-}_{k})E_{fk}=\sum_{k}(1-2n^{-}_{k})E_{fk},$ (193)

where to derive the first equality we have used the reflection symmetry (160)
and to derive the second we have used Eq. (192). Taking these into account we
simplify Eq. (191) to

$\displaystyle\langle\Psi(t)|{\cal
A}|\Psi(t)\rangle=\sum_{\\{{n_{k}^{-}}{}^{\prime}\\}}\sum_{\\{{n^{-}_{k}}\\}}e^{it\sum_{k}({n^{-}_{k}}-{n_{k}^{-}}{}^{\prime}){2E_{fk}}}\,w_{{\bf
m}^{\prime}}^{*}\,w_{\bf m}\,\langle\Psi_{{\bf m}^{\prime}}|{\cal A}|\Psi_{\bf
m}\rangle.$ (194)

Note that due to Eq. (192) the many-body eigenstate ${\bf m}$ (or ${\bf
m}^{\prime}$) here is completely fixed by the occupation number configuration
of the hole band $\\{n_{k}^{-}\\}\equiv{\bf m}^{-}$ (or
$\\{{n_{k}^{-}}{}^{\prime}\\}\equiv{\bf m}^{\prime}{}^{-}$). In Eq. (194), the
dynamical phases are governed by the frequencies associated with the energy
gaps $2E_{fk}$ of the single-particle band structure and the difference in
hole numbers ${n^{-}_{k}}-{n_{k}^{-}}{}^{\prime}$ between two states ${\bf
m},\,{\bf m}^{\prime}$.

Then we apply Eq. (194) to one-body operators, e.g. the correlation function
$C(t)$ which belongs to this class. Such operators take the general form as

${\cal A}=\sum_{\zeta\zeta^{\prime}}\sum_{kk^{\prime}}{\cal
A}_{kk^{\prime}}^{\zeta\zeta^{\prime}}e^{i(kj-k^{\prime}j^{\prime})}d_{k\zeta}^{\dagger}d_{k^{\prime}\zeta^{\prime}},$
(195)

where $d^{\dagger}_{k\zeta}\,,d_{k\zeta}$ are, respectively, fermionic
creation and annihilation operators at the Bloch momentum $k$ and at the band
$\zeta$, the coefficients ${\cal A}_{kk^{\prime}}^{\zeta\zeta^{\prime}}$ are
complex in general, and $j,j^{\prime}=1,\ldots,L$ label the unit cells. It is
easy to see that

$\displaystyle\langle\Psi_{{\bf
m}^{\prime}}|d_{k\zeta}^{\dagger}d_{k^{\prime}\zeta^{\prime}}|\Psi_{\bf
m}\rangle$ $\displaystyle=$
$\displaystyle\delta_{kk^{\prime}}\langle\Psi_{{\bf
m}^{\prime}}|d_{k\zeta}^{\dagger}d_{k\zeta^{\prime}}|\Psi_{\bf m}\rangle$
(198)
$\displaystyle\times\left\\{\begin{array}[]{ll}n_{k}^{\zeta}\prod_{k^{\prime}}\delta_{{n_{k^{\prime}}^{-}}{}^{\prime}{n_{k^{\prime}}^{-}}},&{\rm
for}\,\,\zeta=\zeta^{\prime}\\\
n_{k}^{\zeta^{\prime}}\delta_{(1-{n_{k}^{-}}{}^{\prime})n_{k}^{-}}\prod_{k^{\prime}\neq
k}\delta_{{n_{k^{\prime}}^{-}}{}^{\prime}{n_{k^{\prime}}^{-}}},&{\rm
for}\,\,\zeta=-\zeta^{\prime}\end{array}\right.,$

due to the constraint (192). The product of Kronecker symbols in the second
line shows that ${\bf m}^{\prime}$ is uniquely determined by ${\bf m}$ (more
precisely, its occupation number configuration ${\bf m}^{-}$ at the hole band)
for given $\zeta,\,\zeta^{\prime}$, and thus we denote it as ${\bf
m}^{\prime}({\bf m}^{-})$. With the substitution of Eqs. (195) and (198) into
Eq. (194), we obtain

$\langle\Psi(t)|{\cal A}|\Psi(t)\rangle={\cal A}_{0}+\delta{\cal A}(t),$ (199)

where ${\cal A}_{0}$ has no time dependence,

$\displaystyle{\cal A}_{0}=\sum_{k}e^{ik(j-j^{\prime})}\sum_{\zeta}{\cal
A}_{kk}^{\zeta\zeta}\sum_{{\bf m}^{-}}n_{k}^{\zeta}|w_{\bf m}|^{2},$ (200)

and $\delta{\cal A}$ depends on time,

$\displaystyle\delta{\cal A}(t)=\sum_{k}e^{ik(j-j^{\prime})}\sum_{\zeta}{\cal
A}_{kk}^{\zeta(-\zeta)}\sum_{{\bf m}^{-}}n_{k}^{-\zeta}\,w_{{\bf
m}^{\prime}({\bf m}^{-})}^{*}\,w_{\bf m}\,e^{it(2n_{k}^{-}-1){2E_{fk}}}.$
(201)

The first factors in Eqs. (200) and (201), $e^{ik(j-j^{\prime})}$, imply that
the expectation value is translationally invariant in the unit cell index.
Most importantly, we can rewrite Eq. (201) as

$\displaystyle\delta{\cal A}(t)$ $\displaystyle=$
$\displaystyle\sum_{k}e^{ik(j-j^{\prime})+i{2E_{fk}t}}\,{\cal
A}_{kk}^{+-}\,\sum_{{\bf m}^{-}}\delta_{1n_{k}^{-}}\,w_{{\bf m}^{\prime}({\bf
m}^{-})}^{*}\,w_{\bf m}$ (202) $\displaystyle+$
$\displaystyle\sum_{k}e^{ik(j-j^{\prime})-i{2E_{fk}t}}\,{\cal
A}_{kk}^{-+}\,\sum_{{\bf m}^{-}}\delta_{0n_{k}^{-}}w_{{\bf m}^{\prime}({\bf
m}^{-})}^{*}\,w_{\bf m}.$

Taking into account the symmetries (159) and (160), we thus arrive at the same
physical picture as in Sec. II.1. That is, for a generic one-body operator
${\cal A}$,

$\displaystyle\langle\Psi(t)|{\cal
A}|\Psi(t)\rangle\leftrightsquigarrow\boldsymbol{\varphi}=\boldsymbol{\omega}t\in\mathbb{T}^{N},$
(203)

i.e., the quantum expectation $\langle\Psi(t)|{\cal A}|\Psi(t)\rangle$
corresponds one-to-one to a point $\boldsymbol{\omega}t$ along the classical
trajectory on $\mathbb{T}^{N}$.

##  Appendix D Proof of the concentration inequalities (53) and (54)

We shall prove inequality (53) only. Inequality (54) can be proved in the same
way. Since $L$ and $N$ — the number of independent random variables
$\varphi_{m}$ — are finite, we adopt the theory of concentration of measure
[15-17].

First of all, because $\mathfrak{R}_{l,\sigma\sigma^{\prime}}(t)$ is a
quasiperiodic in $t$, applying Eq. (43) gives

$\displaystyle\lim_{T\rightarrow\infty}\int_{0}^{T}\frac{dt}{T}\mathfrak{R}_{l,\sigma\sigma^{\prime}}(t)=\int{d\boldsymbol{\varphi}\over(2\pi)^{N}}\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})=0.$
(204)

Then, by Eq. (46) we have

$\displaystyle\lim_{T\rightarrow\infty}\int_{|\mathfrak{R}_{l,\sigma\sigma^{\prime}}(t)|\geq\epsilon}\frac{dt}{T}=\mathbf{P}\left(|\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})|\geq\epsilon\right).$
(205)

The rest of the proof is to bound the probability on the right-hand side.
Using Eq. (52) we find that

$\displaystyle\partial_{\varphi_{m}}\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})=-\frac{1}{L}R_{l,\sigma\sigma^{\prime}}\left(k_{m}\right)\,\sin\varphi_{m}.$
(206)

Thus

$\displaystyle|\partial_{\boldsymbol{\varphi}}\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})|^{2}$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{N-1}(\partial_{\varphi_{m}}\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi}))^{2}=\frac{1}{L^{2}}\sum_{m=0}^{N-1}\left(R_{l,\sigma\sigma^{\prime}}\left(k_{m}\right)\right)^{2}\,\sin^{2}\varphi_{m}$
(207) $\displaystyle\leq$
$\displaystyle\frac{1}{L^{2}}\sum_{m=0}^{N-1}\left(R_{l,\sigma\sigma^{\prime}}\left(k_{m}\right)\right)^{2}.$

For $L\gg 1$ the sum converges well to the integral, giving

$\displaystyle\frac{1}{L^{2}}\sum_{m=0}^{N-1}\left(R_{l,\sigma\sigma^{\prime}}\left(k_{m}\right)\right)^{2}\approx\frac{1}{2L}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\left(e^{ikl}\alpha_{\sigma\sigma^{\prime}}(kl)+e^{-ikl}\alpha_{\sigma\sigma^{\prime}}(-kl)\right)^{2}\equiv\|\mathfrak{R}_{l,\sigma\sigma^{\prime}}\|_{\rm
Lip}^{2}.$ (208)

Since $\mathbb{T}^{N}=\mathbb{R}^{N}/\mathbb{Z}^{N}$ inherits the Riemann
structure from the flat Euclidean space $\mathbb{R}^{N}$, we have

$|\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})-\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi}^{\prime})|\leq\|\mathfrak{R}_{l,\sigma\sigma^{\prime}}\|_{\rm
Lip}\|\boldsymbol{\varphi}-\boldsymbol{\varphi}^{\prime}\|.$ (209)

That is, the function
$\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})$ is Lipschitz
continuous with respect to the Euclidean norm $\|\cdot\|$, and
$\|\mathfrak{R}_{l,\sigma\sigma^{\prime}}\|_{\rm Lip}$ is called the Lipschitz
constant. Combining this continuity property with a well-known theorem in the
theory of concentration of measure [17], we find that

$\displaystyle\mathbf{P}\left(|\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})|\geq\epsilon\right)\leq
2e^{-\frac{\delta^{\prime}}{\|\mathfrak{R}_{l,\sigma\sigma^{\prime}}\|_{\rm
Lip}^{2}}\,\epsilon^{2}},$ (210)

where $\delta^{\prime}$ is some absolute constant of order unity. With the
substitution of Eq. (208) into the right-hand side, we justify inequality (53)
and Eq. (55).

##  Appendix E Statistical independence of correlation matrix elements

In this Appendix we show that for large $L$ any random variable
$\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}},\,\mathfrak{I}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$
in the block $(C(t))_{l_{1}}$ is statistically independent of
$\mathfrak{R}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}},\,\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$
in a distinct block $(C(t))_{l_{2}}$ ($l_{1}\neq l_{2}$), or equivalently, any
random variable
$\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}},\,\mathfrak{I}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$
in the block $(\tilde{C}(\boldsymbol{\varphi}))_{l_{1}}$ is statistically
independent of
$\mathfrak{R}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}},\,\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$
in a distinct block $(\tilde{C}(\boldsymbol{\varphi}))_{l_{2}}$. Without loss
of generality we prove the statistical independence of
$\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$ and
$\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$ below. The statistical
independence of $\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$ and
$\mathfrak{R}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$ and of
$\mathfrak{I}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$ and
$\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$ can be proven in the same
way.

Consider the characteristic function of those two random variables, defined as

$\Phi_{2}(z_{1},z_{2})\equiv\lim_{T\rightarrow\infty}\int_{0}^{T}{dt\over
T}\exp\left(i\left(z_{1}\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(t)+z_{2}\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(t)\right)\right).$
(211)

Applying Eq. (43) we obtain

$\Phi_{2}(z_{1},z_{2})=\int{d\boldsymbol{\varphi}\over(2\pi)^{N}}\exp\left(i\left(z_{1}\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(\boldsymbol{\varphi})+z_{2}\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(\boldsymbol{\varphi})\right)\right).$
(212)

With the substitution of
$\mathfrak{R}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})$ and
$\mathfrak{I}_{l,\sigma\sigma^{\prime}}(\boldsymbol{\varphi})$ given by Eq.
(50), we obtain

$\displaystyle\Phi_{2}(z_{1},z_{2})$ $\displaystyle=$
$\displaystyle\int{d\boldsymbol{\varphi}\over(2\pi)^{N}}\exp\left({i\over
L}\sum_{m=0}^{N-1}\left(z_{1}R_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(k_{m})\cos\varphi_{m}-iz_{2}I_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(k_{m})\sin\varphi_{m}\right)\right)$
(213) $\displaystyle=$ $\displaystyle\prod_{m=0}^{N-1}J_{0}\left({1\over
L}\sqrt{\left(z_{1}R_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(k_{m})\right)^{2}-\left(z_{2}I_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(k_{m})\right)^{2}}\right),$

where $J_{0}(x)$ is the zero-order Bessel function. Recall that
$I_{l_{2},\sigma_{2}\sigma^{\prime}_{2}}$ is purely imaginary. Because of
$J_{0}(x\ll 1)\approx 1-{x^{2}\over 4}\approx e^{-{x^{2}\over 4}}$,

$\displaystyle\Phi_{2}(z_{1},z_{2})\stackrel{{\scriptstyle L\gg
1}}{{\longrightarrow}}\prod_{m=0}^{N-1}\exp\left(-{1\over
4L^{2}}\left(z_{1}R_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(k_{m})\right)^{2}-\left(z_{2}I_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(k_{m})\right)^{2}\right).$
(214)

Next, we consider the characteristic function of the single random variable
$\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}$, defined as

$\Phi_{1}(z_{1})\equiv\lim_{T\rightarrow\infty}\int_{0}^{T}{dt\over
T}\exp\left(iz_{1}\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(t)\right),$
(215)

which is

$\Phi(z_{1})=\int{d\boldsymbol{\varphi}\over(2\pi)^{N}}\exp\left(iz_{1}\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(\boldsymbol{\varphi})\right).$
(216)

By the same token we find that

$\displaystyle\Phi_{1}(z_{1})\stackrel{{\scriptstyle L\gg
1}}{{\longrightarrow}}\prod_{m=0}^{N-1}\exp\left(-{1\over
4L^{2}}\left(z_{1}R_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(k_{m})\right)^{2}\right).$
(217)

Similarly, we have

$\displaystyle\Phi_{1}(z_{2})$ $\displaystyle\equiv$
$\displaystyle\lim_{T\rightarrow\infty}\int_{0}^{T}{dt\over
T}\exp\left(iz_{2}\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(t)\right)=\int{d\boldsymbol{\varphi}\over(2\pi)^{N}}\exp\left(iz_{2}\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(\boldsymbol{\varphi})\right)$
(218) $\displaystyle\stackrel{{\scriptstyle L\gg 1}}{{\longrightarrow}}$
$\displaystyle\prod_{m=0}^{N-1}\exp\left({1\over
4L^{2}}\left(z_{2}I_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(k_{m})\right)^{2}\right).$
(219)

for the characteristic function of
$\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}$. From Eqs. (214), (217)
and (219) we obtain

$\displaystyle\Phi_{2}(z_{1},z_{2})\stackrel{{\scriptstyle L\gg
1}}{{\longrightarrow}}\Phi_{1}(z_{1})\Phi_{1}(z_{2})$ (220)

That is, $\Phi_{2}(z_{1},z_{2})$ factorizes into
$\Phi_{1}(z_{1})\Phi_{1}(z_{2})$ for large $L$. This justifies the statistical
independence of $\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(t)$ and
$\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(t)$. By Eqs. (212), (216)
and (218) it is equivalent to the statistical independence of
$\mathfrak{R}_{l_{1},\sigma_{1}\sigma_{1}^{\prime}}(\boldsymbol{\varphi})$ and
$\mathfrak{I}_{l_{2},\sigma_{2}\sigma_{2}^{\prime}}(\boldsymbol{\varphi})$.

##  Appendix F Some properties of $\boldsymbol{\check{\alpha}_{m}}$,
$\boldsymbol{\check{\beta}_{m}}$ and $\boldsymbol{\check{\gamma}_{m}}$

We recall that ${\check{\alpha}_{m}}$, ${\check{\beta}_{m}}$ and
${\check{\gamma}_{m}}$ are defined as

$\displaystyle\check{\alpha}_{m}\equiv\left.\left(\begin{array}[]{cc}\alpha_{{\bar{A}}{\bar{A}}}(k_{m})&\alpha_{{\bar{A}}{\bar{B}}}(k_{m})\\\
\alpha_{{\bar{B}}{\bar{A}}}(k_{m})&\alpha_{{\bar{B}}{\bar{B}}}(k_{m})\end{array}\right)\right.,\,\check{\beta}_{m}\equiv\left.\left(\begin{array}[]{cc}\beta_{{\bar{A}}{\bar{A}}}(k_{m})&\beta_{{\bar{A}}{\bar{B}}}(k_{m})\\\
\beta_{{\bar{B}}{\bar{A}}}(k_{m})&\beta_{{\bar{B}}{\bar{B}}}(k_{m})\end{array}\right)\right.,\,\check{\gamma}_{m}\equiv\left.\left(\begin{array}[]{cc}\gamma_{{\bar{A}}{\bar{A}}}(k_{m})&\gamma_{{\bar{A}}{\bar{B}}}(k_{m})\\\
\gamma_{{\bar{B}}{\bar{A}}}(k_{m})&\gamma_{{\bar{B}}{\bar{B}}}(k_{m})\end{array}\right).\right.$
(227)

Combining Eqs. (182)-(188) with Eqs. (165)-(171) and Eq. (180), we can show
straightforwardly the following properties for the elements of these three
matrices:

  * •

For $m>0$,

$\displaystyle\textrm{(i)\
}\gamma_{{\bar{A}}{\bar{A}}}(k_{m})=\gamma_{{\bar{A}}{\bar{A}}}(-k_{m})\in\mathbb{R};$
$\displaystyle\textrm{(ii)\
}\gamma_{{\bar{A}}{\bar{B}}}(k_{m})=e^{i\delta(k_{m})}\,\gamma^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m}),\,\,\,\,\,\,\,\,\textrm{with
}\delta(k_{m})=-\delta(-k_{m}),\,\,\,\gamma^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})=\gamma^{\prime}_{{\bar{A}}{\bar{B}}}(-k_{m}),\,\,\,\textrm{and}\,\gamma^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})\in\mathbb{R};$
$\displaystyle\textrm{(iii)\
}\alpha_{{\bar{A}}{\bar{A}}}(k_{m})=\alpha_{{\bar{A}}{\bar{A}}}(-k_{m})\in\mathbb{R};$
$\displaystyle\textrm{(iv)\
}\alpha_{{\bar{A}}{\bar{B}}}(k_{m})=e^{i\delta(k_{m})}\alpha^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m}),\,\,\,\textrm{with\
}\alpha^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})=\alpha^{\prime\,*}_{{\bar{A}}{\bar{B}}}(-k_{m});$
$\displaystyle\textrm{(v)\
}\beta_{{\bar{A}}{\bar{A}}}(k_{m})=-\beta_{{\bar{A}}{\bar{A}}}(-k_{m})\in\mathbb{R};$
$\displaystyle\textrm{(vi)\
}\beta_{{\bar{A}}{\bar{B}}}(k_{m})=e^{i\delta(k_{m})}\,\beta_{{\bar{A}}{\bar{B}}}^{\prime}(k_{m}),\,\,\,\textrm{with\
}\beta^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})=-\beta^{\prime\,*}_{{\bar{A}}{\bar{B}}}(-k_{m}).$
(228)

Recall that $\gamma^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})$,
$\alpha^{\prime}_{{\bar{A}}{\bar{B}}}(k_{m})$ and
$\beta_{{\bar{A}}{\bar{B}}}^{\prime}(k_{m})$ are respectively the off-diagonal
(${\bar{A}}{\bar{B}}$) component of $\check{\Gamma}^{\prime}_{m}\,$,
$\check{\alpha}^{\prime}_{m}$ and $\check{\beta}^{\prime}_{m}$.

  * •

For $m=0$, we have $k=0$ and

$\displaystyle\textrm{(vii)\
}\gamma_{{\bar{A}}{\bar{A}}}(0)\in\mathbb{R},\textrm{\ \ \ (viii)\
}\gamma_{{\bar{A}}{\bar{B}}}(0)\in\mathbb{R},\textrm{\ \ \ (ix)\
}\alpha_{{\bar{A}}{\bar{A}}}(0)\in\mathbb{R},$ $\displaystyle\textrm{\ (x)\
}\alpha_{{\bar{A}}{\bar{B}}}(0)\in\mathbb{R},\textrm{\ \ \ (xi)\
}\beta_{{\bar{A}}{\bar{A}}}(0)=0,\textrm{\ \ \ \ \ (xii)\
}\beta_{{\bar{A}}{\bar{B}}}(0)\,\,\,\textrm{ purely\, imaginary}.$ (229)

##  Appendix G Proof of the identity Eq. (126)

In this Appendix we prove the identity Eq. (126). Throughout this appendix we
suppress all the subscripts $m$ and the corresponding Bloch momenta $k_{m}$ to
make the formulae compact.

For $m=0$, by using the definitions of $\gamma$, $\tilde{\gamma}$, $\alpha$,
$\tilde{\alpha}$ and Eqs. (182), (183), (186) and (187), we obtain

$\displaystyle\gamma\alpha+\tilde{\gamma}\tilde{\alpha}=4ab\,\biggl{(}(au^{-})^{2}+(bu^{+})^{2}-\frac{1}{2}\biggr{)}u^{+}u^{-}+2ab\,(a^{2}u^{-}v^{-}+b^{2}u^{+}v^{+})(u^{-}v^{+}+u^{+}v^{-})=0.$
(230)

The last expression is obtained with the help of the identities:
$u^{-}u^{+}=-v^{-}v^{+}$, $(u^{\pm})^{2}+(v^{\pm})^{2}=1$ and $a^{2}+b^{2}=1$.

For $m>0$, similarly, by writing out the expressions explicitly we have

$\displaystyle\gamma=\Bigl{(}|a|^{2}(u^{-})^{2}+|b|^{2}(u^{+})^{2}\Bigr{)}-\frac{1}{2},\,\,\,\,\,\,\,\tilde{\gamma}=|a|^{2}u^{-}{v^{\prime}}^{-}+|b|^{2}u^{+}{v^{\prime}}^{+},$
(231)

and

$\displaystyle\alpha=2u^{+}u^{-}(ab^{*}+a^{*}b),\,\,\,\,\,\,\,\tilde{\alpha}=(u^{-}{v^{\prime}}^{+}+u^{+}{v^{\prime}}^{-})(ab^{*}+a^{*}b).$
(232)

So one can check that

$\displaystyle\gamma\alpha+\tilde{\gamma}\tilde{\alpha}$ $\displaystyle=$
$\displaystyle(ab^{*}+a^{*}b)\biggl{(}|a|^{2}\Bigl{(}2(u^{-})^{3}u^{+}+(u^{-})^{2}{v^{\prime}}^{-}{v^{\prime}}^{+}+u^{-}u^{+}({v^{\prime}}^{-})^{2}\Bigr{)}$
(233)
$\displaystyle+|b|^{2}\Bigl{(}2(u^{+})^{3}u^{-}+(u^{+})^{2}{v^{\prime}}^{-}{v^{\prime}}^{+}+u^{-}u^{+}({v^{\prime}}^{+})^{2}\Bigr{)}-u^{-}u^{+}\biggr{)}=0$

with the help of the identities:
$u^{-}u^{+}=-{v^{\prime}}^{-}{v^{\prime}}^{+}$,
$(u^{\pm})^{2}+({v^{\prime}}^{\pm})^{2}=1$ and $|a|^{2}+|b|^{2}=1$.

##  Appendix H Time evolution of generic entanglement probes for
$\boldsymbol{L\rightarrow\infty}$

For a generic entanglement probe $O$ such as $S$ and $S_{n}$, as shown in
Secs. VI and VII its variance obeys the same scaling law as Eq. (148). (Note
that the unimportant coefficients $a,b$ in Eq. (148) change, when $S$ is
replaced by $O$.) As a result, at fixed $L_{A}$ a vanishing variance results
in the limit $L\rightarrow\infty$, i.e.,

$\lim_{T\rightarrow\infty}\int_{0}^{T}{dt\over T}\left(O(t)-\langle
O\rangle\right)^{2}=0.$ (234)

Because $O(t)$ is continuous in time, $O(t)\neq\langle O\rangle$ cannot occur
in a zero measure set. Therefore, $(O(t)-\langle O\rangle)^{2}$ must decay for
sufficiently large $t$ (otherwise the limit does not exist), but this
contradicts the quasiperiodic behavior of $O(t)$. So the only possibility for
Eq. (234) to hold is that there exists a critical time $t_{c}$, such that

$O(t\geq t_{c})=\langle O\rangle.$ (235)

In words, $O(t)$ is strictly a constant $\langle O\rangle$ beyond $t_{c}$.
This kind of phenomena were first found numerically in Ref. [1] and proven
analytically in Ref. [3], but only for the entanglement entropy $S$.

References:

[1] P. Calabrese and J. L. Cardy, Evolution of entanglement entropy in one-
dimensional systems, J. Stat. Mech. 2005, P04010 (2005).
[2] B.-Q. Jin and V. E. Korepin, Quantum spin chain, Toeplitz determinants and
the Fisher-Hartwig conjecture, J. Stat. Phys. 116, 79 (2004).
[3] M. Fagotti and P. Calabrese, Evolution of entanglement entropy following a
quantum quench: Analytic results for the XY chain in a transverse magnetic
field, Phys. Rev. A 78, 010306 (2008).
[4] I. Peschel and V. Eisler, Reduced density matrices and entanglement
entropy in free lattice models, J. Phys. A: Math. Theor. 42, 504003 (2009).
[5] B. L. Altshuler, Fluctuations in the extrinsic conductivity of disordered
conductors, Pis’ma Zh. E’ksp. Teor. Fiz. 41, 530 (1985) [JETP Lett. 41, 648
(1985)].
[6] P. A. Lee and A. D. Stone, Universal conductance fluctuations in metals,
Phys. Rev. Lett. 55, 1622 (1985).
[7] P. Sheng, Introduction to wave scattering, localization, and mesoscopic
phenomena (2nd Ed., Springer, Berlin, Germany, 2006).
[8] E. Akkermans and G. Montambaux, Mesoscopic physics of electrons and
photons (Cambridge University Press, UK, 2007).
[9] K. B. Efetov, Supersymmetry in Disorder and Chaos (Cambridge University,
Cambridge, UK, 1997).
[10] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod.
Phys. 69, 731-808 (1997).
[11] M. Kac, Statistical independence in probability, analysis and number
theory (The Mathematical Association of America, 1959).
[12] V. I. Arnold, Mathematical methods of classical mechanics 2nd Ed.
(Springer, New York, U.S., 1989).
[13] A. M. Samoilenko, Quasiperiodic oscillations, Scholarpedia 2, 1783
(2007).
[14] R. Modak, V. Alba, and P. Calabrese, Entanglement revivals as a probe of
scrambling in finite quantum systems, J. Stat. Mech. 2020, 083110 (2020).
[15] M. Talagrand, A new look at independence, Ann. Probab. 24, 1 (1996).
[16] V. D. Milman and G. Schechtman, Asymptotic theory of finite-dimensional
normed spaces (Springer, Berlin, Germany, 1986).
[17] M. Ledoux, The concentration of measure phenomenon (AMS, Providence,
2001).
[18] E. Bogolmony, Spectral statistics of random Toeplitz matrices, Phys. Rev.
E 102, 040101(R) (2020).
[19] E. Bogolmony and O. Giraud, Statistical properties of structured random
matrices, Phys. Rev. E 103, 042213 (2021).
[20] F. Haake, Quantum signatures of chaos 2nd ed. (Springer, Berlin, 2001).
[21] F. Evers and A. D. Mirlin, Anderson transitions, Rev. Mod. Phys. 80, 1355
(2008).
[22] B. L. Altshuler, I. Kh. Zharekeshev, S. A. Kotochigava, and B. I.
Shklovskii, Repulsion between levels and the metal-insulator transition, Sov.
Phys. JETP 67, 625 (1988).
[23] C. Hammond and S. J. Miller, Distribution of eigenvalues for the ensemble
of real symmetric Toeplitz matrices, J. Theor. Probab. 18, 537 (2005).
[24] E. B. Bogomolny, U. Gerland, and C. Schmit, Models of intermediate
spectral statistics, Phys. Rev. E 59, 1315(R), (2001).
[25] M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, Universal structure
of transmission eigenchannels inside opaque media, Nature Commun. 6, 6893
(2015).
[26] A. Altland, A. Kamenev, and C. Tian, Anderson localization from the
replica formalism, Phys. Rev. Lett. 95, 206601 (2005).
[27] P. Fang, L. Y. Zhao, and C. Tian, Concentration-of-measure theory for
structures and fluctuations of waves, Phys. Rev. Lett. 121, 140603 (2018).
[28] S. Boucheron, G. Lugosi, and P. Massart, Concentration inequalities - A
nonasymptotic theory of independence (Oxford University Press, UK, 2013).
[29] L. Vidmar, L. Hackl, E. Bianchi, and M. Rigol, Entanglement entropy of
eigenstates of quadratic fermionic Hamiltonians, Phys. Rev. Lett. 119, 020601
(2017).
[30] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M.
Preiss, and M. Greiner, Quantum thermalisation through entanglement in an
isolated many-body system, Science 353, 794 (2016).
[31] J. von Neumann, Beweis des Ergodensatzes und des H-Theorems in der neuen
Mechanik, Z. Phys. 57, 30 (1929).
[32] S. Goldstein, J.L. Lebowitz, R. Tumulka, and N. Zanghi, Long-time
behavior of macroscopic quantum systems, Eur. Phys. J. H 35, 173 (2010).
[33] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos
and eigenstate thermalization to statistical mechanics and thermodynamics,
Adv. Phys. 65, 239 (2016).
[34] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the
emergence of statistical mechanics in closed quantum systems, Rep. Prog. Phys.
79, 056001 (2016).
[35] F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky, Quantum
chaos and thermalization in isolated systems of interacting particles, Phys.
Rep. 626, 1 (2016).
[36] C. Tian and K. Yang, Breakdown of quantum-classical correspondence and
dynamical generation of entanglement, Phys. Rev. B 104, 174302 (2021).