Mesoscopic fluctuations in entanglement dynamics

Understanding fluctuation phenomena plays a dominant role in the development of many-body physics. The time evolution of entanglement is essential to a broad range of subjects in many-body physics, ranging from exotic quantum matter to quantum thermalization. Stemming from various dynamical processes of information, fluctuations in entanglement evolution differ conceptually from out-of-equilibrium fluctuations of traditional physical quantities. Their studies remain elusive. Here we uncover an emergent random structure in the evolution of the many-body wavefunction in two classes of integrable—either interacting or noninteracting—lattice models. It gives rise to out-of-equilibrium entanglement fluctuations which fall into the paradigm of mesoscopic fluctuations of wave interference origin. Specifically, the entanglement entropy variance obeys a universal scaling law in each class, and the full distribution displays a sub-Gaussian upper and a sub-Gamma lower tail. These statistics are independent of both the system’s microscopic details and the choice of entanglement probes, and broaden the class of mesoscopic universalities. They have practical implications for controlling entanglement in mesoscopic devices.

Understanding fluctuation phenomena plays a dominant role in the development of many-body physics.The time evolution of entanglement is essential to a broad range of subjects in many-body physics, ranging from exotic quantum matter to quantum thermalization.Stemming from various dynamical processes of information, fluctuations in entanglement evolution differ conceptually from out-of-equilibrium fluctuations of traditional physical quantities.Their studies remain elusive.Here we uncover an emergent random structure in the evolution of the wavefunction in a class of integrable models.It gives rise to out-of-equilibrium entanglement fluctuations which, strikingly, fall into the paradigm of mesoscopic fluctuations of wave interference origin.Specifically, the entanglement entropy variance obeys a universal scaling law, and the distribution displays a sub-Gaussian upper and a sub-Gamma lower tail.These statistics are independent of both the system's microscopic details and the choice of entanglement probes, and broaden the class of mesoscopic universalities.They have practical implications for controlling entanglement in mesoscopic devices.
When an isolated many-body system evolves, entanglement tends to spread.Owing to the diversity of the fate of the wavefunction evolution (e.g., localized or delocalized, thermalized or not thermalized), a wealth of entanglement patterns develop [1][2][3][4][5][6][7].These patterns are the building blocks of the physics of recently discovered exotic phases of matter [4,7,8], and are central to the foundations of statistical mechanics [6,7].Understanding the long-time evolution of entanglement, and especially its universal aspects, is indispensable in the study of pattern formation.
To address this issue, one often investigates mesoscopic rather than macroscopic systems.Recent advancements in quantum simulation platforms, ranging from cold atoms, trapped ions to superconducting qubits, have made possible the measurement of information-theoretic observables and the experimental study of entanglement evolution [6,7,9].In these investigations, quantum coherence is maintained across the entire sample, as required also for mesoscopic electronic and photonic devices [10,11].At the same time, the relationship between the evolution of entanglement and quantum thermalization in isolated systems is currently under investigations [6,7].Since various scenarios for the latter [12][13][14][15][16][17] are built upon a basis of wavefunctions with finite spatial extent, emphasis naturally has to be placed on the dynamics of entanglement on the mesoscopic scale.
A prominent feature of mesoscopic systems is the occurrence of unique fluctuation phenomena, when randomness due to quenched disorders [10,11] or chaos [18,19] is present.Notably, the conductance -a basic probe of mesoscopic transport -fluctuations have a universal variance, independent of sample size and the strength of randomness [20,21].Mesoscopic fluctuations * ct@mail.itp.ac.cn are of wave interference origin and conceptually different from thermodynamic fluctuations.They are related to various entanglement properties [22,23].Understanding their universalities is of fundamental importance to mesoscopic physics.Here we uncover a 'random' structure emergent from the dynamical phases in the wavefunction evolution.Treating the information-theoretic observable as an unconventional 'mesoscopic' probe, we explore out-of-equilibrium fluctuation phenomena in entanglement evolution, whose origin are similar to that of mesoscopic fluctuations in genuine disordered samples.
In fact, there is a rapid increase of interests in entanglement fluctuations.In particular, understanding out-ofequilibrium entanglement fluctuation properties is a key to the statistical physics of isolated systems [24,25].So far studies have focused on the kinematic case [16,[26][27][28][29], where fluctuations arise from random sampling of some pure state ensemble, initiated by Page [26].Since kinematic theories cannot describe wave effects and dynamical properties of the Schrödinger evolution [30], outof-equilibrium entanglement fluctuations are beyond the framework of those theories.
Here we develop an analytical theory for fluctuations in long-time dynamics of entanglement in a class of integrable lattice systems, including the Rice-Mele model and the transverse field Ising chain.We find that the wavefunction evolution endows the correlation matrix a random structure, even though the system is neither chaotic nor disordered.Specifically, the time dependence enters through N ≈ L 2 dynamical phases (ω 1 t, . . ., ω N t) ≡ ωt, with L being the number of unit cells, so that the instantaneous correlation matrix C(t) is given by some Nvariable (matrix-valued) function C(ϕ) for ϕ = ωt; due to the incommensurality of ω an ensemble of random matrices C(ϕ) then results.Each C(ϕ) is determined by ϕ, the virtual disorder realization uniformly distributed over a N -dimensional torus (Fig. 1).It describes a virtual disordered sample, and determines entanglement properties FIG. 1. Emergence of mesoscopic fluctuations in entanglement evolution.a.We simulate entanglement entropy evolution of Rice-Mele model up to t = 10 4 in unit of /J, J the amplitude of hopping between two nearest sites (inset, LA = 25, L = 124).Its fluctuation statistics (histograms) is shown to be equivalent to the statistics of entanglement entropy fluctuations in an ensemble of virtual disordered samples (dashed line, for 5 × 10 5 disorder realizations ϕ).b.These long-time fluctuations differ from growth and damped oscillations appearing in early entanglement evolution (inset).c.Simulations show that the nearest-neighbor spacing distribution characterizing spectral fluctuations of the correlation matrix C(t) (inset) is indistinguishable from that for an ensemble of truly random matrices C(ϕ), and is semi-Poissonian.d.Physically, as system's wavefunction evolves, the dynamical phases ϕ = (ω1t, . . ., ωN t) sweeps out an ensemble of 'mesoscopic samples' C(ϕ).
of that sample in the same fashion as C(t) determines system's instantaneous entanglement properties.Consequently, when system's wavefunction evolves, the trajectory ϕ = ωt sweeps out the entire disorder ensemble, trading the temporal fluctuations of various informationtheoretic observables, such as the entanglement entropy and the Rényi entropy, to mesoscopic sample-to-sample fluctuations [20,21].In particular, we find that these out-of-equilibrium entanglement fluctuations arise from wave interference, similar to mesoscopic fluctuations.Interestingly, this kind of trajectories play important roles in Chirikov's studies of the relations between mesoscopic physics and quantum chaos [31].
However, there are important differences between ordinary quenched disorders and the randomness emergent from entanglement evolution.As shown below, the latter has a strength ∼ 1/ √ L and diminishes for L → ∞.This situation renders canonical mesoscopic theories based on diagrammatical [10,11] and field-theoretical [32] methods inapplicable, since they require the disorder strength to be independent of the sample size.In addition, because C(t) is a (block-)Toeplitz matrix and very little [33] is known about the spectral statistics of random Toeplitz matrices, mesoscopic theories based on random matrix methods [34] are inapplicable either.Here we develop a different approach based on the modern nonasymptotic probability theory [35], that relies merely on the statistical independence of the components of ϕ and applies to any L. A related approach has recently been used to find novel universalities in mesoscopic transport [36].
Uncovering the random structure, we show that fluctuations in entanglement evolution exhibit intriguing universal behaviors, independent of microscopic details.First, when the variance Var(S) of the entanglement entropy S as well as L and L A (the subsystem size) are rescaled by appropriate microscopic quantities, the universal scaling law: follows.Second, the statistics of S is universal and the distribution is asymmetric with respect to its mean S , displaying a sub-Gaussian upper and a sub-Gamma lower tail.In particular, the probability for large deviation is where b ± ∝ Var(S) and c > 0 depends on the ratio L A /L. Third, Eqs. ( 1) and ( 2) hold for other probes, e.g., the Rényi entropy.These universal fluctuation behaviors are irrespective of the location of S in Page's curve [26].By Eq. ( 1) at fixed L A the variance vanishes in the limit L → ∞ (cf.Fig. 2a), implying the full suppression of temporal fluctuations beyond some critical time, in agreement with a benchmark result of entanglement evolution [1].In contrast, at fixed L, as L A increases the variance grows as ∼ L 3 A eventually (cf.Fig. 2b), which is much faster than ∼ L A displayed by typical extensive quantities.We shall see below this enhanced growth results from quantum interference.
Having summarized the main results, we outline the derivations and present numerical verifications.A complete description is given in Supplemental Information (SI).We focus on the Rice-Mele model subjected to the periodic boundary condition.Generalizations to other models mentioned above are straightforward.Let the system be at the half-filling ground state Ψ(0).At t = 0 we suddenly change parameters of the Hamiltonian.So the pre-quench state Ψ(0) evolves unitarily under the new Hamiltonian to state Ψ(t) at later time t.Because Ψ(0) is a Gaussian state and the system is fermionic, the instantaneous entanglement entropy can be expressed as using the method in [37][38][39][40].Here e(λ) = −λ ln λ − (1 − λ) ln(1−λ) is the binary entropy function.Tr A δ(...) gives the spectral density of the correlation matrix C(t) defined in the unit cell and sublattice sector, labelled by i and σ, respectively; its element iσ ) being the fermion annihilation (creation) operator.The trace is restricted to the subsystem A.
When replacing e(λ) by an appropriate function of λ, we obtain other entanglement probes such as the Rényi entropy.This kind of expressions indicate that the evolving spectral density underlies out-of-equilibrium behaviors of different entanglement probes.They are analogous to the expressions for probes of mesoscopic transport.Indeed, if we replace C(t) by the product of transmission matrix and its hermitian conjugate, we transform Eq. (3) to the Landauer formula for conductance with e(λ) changed to λ, and to formulaes for other transport probes with e(λ) changed to appropriate functions of λ [34].
Because the eigenenergy spectrum displays a reflection and a particle-hole symmetry, when particle eigenenergies ωm 2 (Planck's constant set to unity) at Bloch momenta , are given, all other particle and all hole eigenenergies are known.Due to the translational invariance of the system, the time parameter enters the correlation matrix through the dynamical phases ωt associated with ω ≡ (ω 1 , ..., ω N ).Specifically, we can define a function C(ϕ) = C 0 + C 1 (ϕ) on the N -dimensional torus.Leaving its detailed form for SI, here we only expose its key properties.First, C 0,1 are block-Toeplitz matrices, with elements (C 0,1 ) ii being 2 × 2 blocks defined in the sublattice sector and depending on the unit cell indexes i, i via (i − i ), i.e., (C 0,1 ) ii ≡ (C 0,1 ) i−i .Second, C 0 is ϕ-independent, whereas C 1 is not and its elements take the form of where the elements of blocks, R's, I's, are complex and depend on k m (as well as post-quench Hamiltonian parameters).Then C(t) is given by C(ϕ) at ϕ = ωt.Similarly, with the introduction of S(ϕ) ≡ dλ e(λ) Tr A δ(λ− C(ϕ)) in parallel to Eq. ( 3) (for notational simplicity we use the same symbol S despite differences in the arguments.),S(t) is given by S(ϕ) at ϕ = ωt.This implies that, like C(t), an evolving entanglement probe depends on t through the dynamical phases ωt.Such dependence has an immediate consequence.That is, because in general the components of ω are incommensurate, after initial growth [1] and damped oscillations [41] due to the traversal of quasiparticle pairs or the incomplete revival of wavefunction (Fig. 1b), an entanglement probe displays quasiperiodic oscillations (Fig. 1a, inset), which are reproducible under the same initial conditions.To understand fluctuation properties of quasiperiodic oscillations we note that the trajectory ϕ = ωt generates an ensemble of random matrices C(ϕ), each of which is determined by the 'disorder realization', ϕ, and thus is separated into two parts: nonrandom C 0 and random C 1 (ϕ).The probability measure of this ensemble is induced by the uniform distribution of ϕ via Eq.( 4).This ensemble has some prominent features (see SI for detailed discussions): First, since ϕ m 's are statistically independent, Eq. ( 4) implies that each element randomly fluctuates around its mean, with a magnitude ∼ 1/ √ L. Thus for fixed ratio L A /L the randomness diminishes in the limit of large matrix size.Second, the elements of two distinct blocks are statistically independent.Third, the average elements decay rapidly with their distance to the main diagonal.These features lead to a semi-Poissonian nearest-neighbor spacing distribution [42], as shown in simulations (Fig. 1c).Strikingly, despite that the Rice-Mele model is integrable and has no extrinsic randomness, the evolving correlation matrix can exhibit level repulsion: P 0 (s → 0) ∼ s, which is a distinctive property of quantum chaos [18,19].We can demonstrate that the statistical equivalence of the ensemble of C(ϕ) and the time series C(t) (Fig. 1c) hinges only on the incommensurabilty of ω (see SI when this condition is not met).Furthermore, much like that a transmission matrix determines transport properties of a mesoscopic sample, a matrix C(ϕ) determines S(ϕ) and other entanglement probes of a virtual mesoscopic sample at the disorder realization ϕ; consequently, the statistical equivalence between C(t) and C(ϕ) leads to the statistical equivalence between out-of-equilibrium and sample-to-sample fluctuations of various entanglement probes, in agreement with simulation results (Fig. 1a).
Exploiting this equivalence, we proceed to study the statistics of entanglement entropy fluctuations.To overcome the difficulties discussed in the introduction with the unusual disorder structure, below we combine the continuity properties of the N -variable function S(ϕ) with the nonasymptotic probabilistic method, so-called concentration inequality [35].This allows us to work out a statistical theory for mesoscopic sample-to-sample fluctuations of S(ϕ) at total system size L, which is finite so that the disorder strength does not vanish.
In order to study the distribution of S(ϕ), we introduce the logarithmic moment-generating function G(u) ≡ ln e u(S− S ) , with u being real and • denoting the average over ϕ.Consider the downward fluctuations (i.e., S − S < 0) first.Because the N components of ϕ are statistically independent, we can apply the so-called modified logarithmic Sobolev inequality [35] to obtain with φ(x) = e x − x − 1 and u ≤ 0.Here S − m is the maximal values of S(ϕ), when ϕ m varies and other arguments are fixed.Observing that the leading u-expansion of the we separate the right-hand side of the inequality into two terms, b− 2 and the remainder.Then, we show that the latter is bounded by c − dG du with c − being a negative constant.So we cast the inequality (6) to d du which can be readily integrated to give 1+|c−|u .Such bound holds also for Gamma random variables.It generalizes the tail behaviors of the Gamma distribution, giving the so-called sub-Gamma tail [35].Specifically, following standard procedures, we can use Markov's inequality to turn this bound for G into a bound for the probability of downward fluctuations.The result is for any >0.This gives a sub-Gamma lower tail, which crosses over from a Gaussian to an exponential form at Similarly, we can study the upward fluctuations (i.e., S − S > 0).We replace S − m in the inequality ( 6) by S + m , which is the minimal S(ϕ) when ϕ m varies and other arguments are fixed, and consider u > 0. Upon separating for any > 0, which is a sub-Gaussian upper tail.The inequalities (9) and (10) show that S(ϕ) concentrates around S albeit with different bounds for upward and downward fluctuations.Simulations further show that the exact deviation probability for large downward (upward) fluctuations agrees with the form given by the right-hand side of the corresponding concentration inequality, with b ± and c − as fitting parameters (Fig. 2c).Therefore, for large deviation, the upper (lower) tail distribution has the universal form given by the first (second) line in Eq. ( 2), and the parameters b ± and c in Eq. ( 2) are proportional to b ± and c − , respectively.So for large the upper tail is always Gaussian e − 2 /(2b+) while the lower is always exponential e − /(2c) , different from the distribution tails of thermodynamic fluctuations which are symmetric and Gaussian.
With Eq. ( 2) we find that the variance Var(S) is given by b ± .To calculate the latter, note that by the mean value theorem there exists φm between ϕ m and ϕ ± m (at which S ± m is reached), so that (S − S ± m ) 2 is given by (ϕ m − ϕ ± m ) 2 (∂ φm S) 2 .Then, for large L the Fourier series of ∂ ϕm S with respect to ϕ m are truncated at the second harmonics, giving (∂ φm S) 2 ∼ dϕm 2π (∂ ϕm S) 2 .Applying these analyses to the definitions of b ± , we obtain This relation is confirmed numerically (Fig. 3a), and the proportionality coefficient is found to be ≈ 1/8.Equation (11) uncovers a relation between entanglement entropy fluctuations and continuity properties of the N -variable function S(ϕ).It resembles the so-called concentrationof-measure phenomenon, a modern perspective of probability theory [43,44], where fluctuations of an observable are controlled by its Lipschitz continuity.This continuity is a key ingredient of universal wave-to-wave fluctuations in mesoscopic transport [36].By definition of S(ϕ), we have , we expand the logarithm in C 1 up to the first order.Taking into account that C 0 is short-ranged, we obtain where H 0 = ln(C −1 0 − I) is the entanglement Hamiltonian in the absence of disorder.Substituting Eq. ( 12) into Eq.( 11), we find that the two terms in Eq. ( 12) contribute to the variance separately.The contribution by the first term is a/L and that by the second is bL 3  A /L 2 , and the former (latter) is found to be a subsystem's edge (bulk) effect.Here the coefficient a is proportional to the square of the size of subsystem's edge, and both a and b have no dependence on L, L A .Upon rescaling: L, L A by a/b and Var(S) by √ ab, we obtain the scaling law (1), which is confirmed by simulations (Fig. 3bd).By Eq. (1), one enters the regime Var(S) = L −1 for L A L 1/3 (b) and the regime Var(S Let us consider other entanglement probes such as the second-order Rényi entropy S 2 .As said above, in this case we have an expression similar to Eq. ( 3), with e(λ) changed (see SI). Repeating the analysis above, we find for S 2 the same relation as (11).Furthermore, we can calculate |∂ ϕ S 2 | 2 in the same way as |∂ ϕ S| 2 .As a result, we find that Var(S 2 ) obeys the same scaling law as Eq.(1).These statistics of S 2 are confirmed numerically (Fig. 3).In SI we further show that Eqs.(1), ( 2) and (11) hold for more general probes.
To understand physically the scaling behavior we use the concept of coherent entangled quasiparticle pair [1].
Consider a quasiparticle inside the subsystem A. When pairing with another outside, it contributes to the bipartite entanglement.Due to the Heisen- Wave origin of entanglement fluctuations.The pairing amplitude of the coherent entangled quasiparticle pair (solid lines) fluctuates with time.Constructive interference between two paths due to virtual hopping (blue and red dashed lines), that underlies such fluctuations, leads the variance of a generic entanglement probe to exhibiting the universal scaling behavior ∼ L 3 A /L 2 .berg uncertainty, this particle's position fluctuates with time, leading to the temporal fluctuation Φ(t) of the pairing amplitude.In the simplest case, the particle hops virtually from a site i to j (in A as well) and back to i. Since the entangled pair is a correlation effect, Φ(t)∼ ij (C 1 (ωt)) ij (C 1 (ωt)) ji and thus by Eq. ( 4) Φ(t)∼ 1 L 2 ij mn e i(km−kn)(i−j) e i(ωm+ωn)t/2 , where k m , ωm 2 are respectively the Bloch momentum and the particle eigenenergy associated with the hopping i→j, and k n , ωn 2 with j→i.The variance of a generic entanglement probe is given by where (k m − k n )(i − j) and (k m − k n )(i − j ) are the phases of the paths: i→j→i and i →j →i , re-spectively.Because ω's are incommensurate, we obtain (m, n)=(m , n ) or (n , m ).So the first sum is dominated by those terms with two phases being identical.As a result, dt|Φ(t)| 2 ∼ L 3 A /L 2 , with the numerator (denominator) given by the first (second) sum: This is the second term in Eq. (1).We see that it arises from the constructive interference between the two hopping paths (Fig. 4).
Our theory essentially hinges upon the relation between the wavefunction evolution and the trajectory ϕ = ωt on a high-dimensional torus, and the informationtheoretic observable as a function on that torus.So it can be extended to more general contexts.First, it applies to truly disordered systems and to other characteristics of entanglement, e.g., the multipartite entanglement entropy and the largest eigenvalue of the entanglement Hamiltonian.Second, when the initial state is not a Slater determinant state or has spontaneous symmetry breaking, or when the interaction is present, the evolving state is not Gaussian.In this case, Eq. ( 3) does not apply.However, it is conceivable that using the reduced density matrix one can still establish the statistical equivalence between the fluctuations with time and with disordered samples, and study the fluctuation statistics by the same token.Finally, because each virtual disordered sample corresponds to a pure state, our work suggests a simple way of producing a random pure-state ensemble to which great experimental efforts [45] are made.That is, we evolve an initial pure state by a single Hamiltonian and collect states at distinct sufficiently long time.

SUPPLEMENTAL INFORMATION: MESOSCOPIC FLUCTUATIONS IN ENTANGLEMENT DYNAMICS
Lih-King Lim, 1 Cunzhong Lou, 1 and Chushun Tian In this Supplemental Information (SI) we present a complete description of the theory and details of numerical simulations, written in the self-contained manner.It is organized in the following way: • In Sec.I as the preliminary, we introduce the Rice-Mele model and give the overall structure of its correlation matrix.
• In Sec.II we develop for several lattice fermionic models a formalism for long-time dynamics of entanglement, and demonstrate its connection to classical dynamics on a high-dimensional torus.
• In Sec.III we use the formalism to establish a statistical principle for out-of-equilibrium fluctuations of various information-theoretic observables.
• In Sec.IV, armed with that principle we show that, despite the models considered are all integrable and not disordered, a random structure emerges from the unitary pure-state evolution, and as a result mesoscopic sample-to-sample fluctuations emerge from entanglement evolution.
• In Sec.V we develop a general theory for the statistics of emergent mesoscopic fluctuations, that allows us to study the full distribution and the variance of a generic information-theoretic observable.
• In Sec.VI we apply the general theory to study entanglement entropy fluctuations.In particular, we uncover a universal scaling law for the variance and derive the full distribution.
• In Sec.VII we describe the methods of numerical simulations in details.We also show numerically that the nth-order Rényi entropy exhibits the same universal behaviors as the entanglement entropy.
• In Sec.VIII we further show that the mesoscopic fluctuations emergent from entanglement evolution disappear for commensurate ω.
• In Appendix A we describe the dynamics of entanglement in the transverse field Ising chain.
• In Appendices B-H we give some additional technical details.
Except Secs.I, II and VIII, we consider incommensurate ω throughout this SI.

I. PRELIMINARY: CORRELATION MATRIX OF THE RICE-MELE MODEL
The Rice-Mele model describes the motion of many fermions on a one-dimensional discrete lattice with two sublattices: Ā, B. In the lattice space, its Hamiltonian is given by Here J, J are the hopping amplitudes, M is the staggered onsite mass, c † iσ , c iσ (σ = Ā , B) are, respectively, fermionic creation and annihilation operators at the σ-sublattice sites belonging to the i-th unit cell with a total of L unit cells.The periodic boundary condition is imposed.We set the lattice constant to be unity, and set the distance between adjacent Ā -and B -sites to be one half.
In this work, we consider generic sudden global quenches of Hamiltonian from H 0 to H f , corresponding to the sudden change of the Hamiltonian parameters in Eq. ( 1): at time t = 0, where the subscript 0 (f ) denotes the parameter set in the pre-quench (post-quench) Hamiltonian H 0 (H f ).Moreover, the initial state Ψ(0) is taken to be the ground state of H 0 at half-filling, which corresponds to a fermion number of L, and the pure-state evolution results in Ψ(t) = e −iH f t/ Ψ(0).Hereafter the Planck constant is set to unity.
To study the evolution of general entanglement probes, namely information-theoretic observables, after quench, we will see that all information is encoded in the time evolution of the correlation matrix C(t) ≡ {C iσ,i σ (t)} restricted to the subsystem A with L A contiguous unit cells.Its matrix elements C iσ,i σ (t) = Ψ(t)|c † iσ c i σ |Ψ(t) , with the indexes i, i = 1, . . ., L A and σ, σ = Ā, B, can be computed to give where the coefficients α's, β's and γ's have no time dependence, but depend only on the Bloch momentum k, the sublattice indexes σ, σ as well as the parameters of Hamiltonians H 0 and H f .E f k is the positive branch of the energy spectrum of H f : It shows that the correlation matrix is a block-Toeplitz matrix, sharing the same form as the correlation matrix for the transverse field Ising chain (TFIC) [1,2].For TFIC, we refer to the Appendix A for the description.To be specific, its matrix element C iσ,i σ depends on the unit cell indexes i, i only through the displacement So we write the coefficients in the block-matrix form in the sublattice sector as follows with six independent complex coefficients.Their explicit forms are given in the Appendix B.
A closely related model is the Su-Schrieffer-Heeger (SSH) model, by setting M → 0 in Eq. ( 1).All results pertaining to the Rice-Mele model hold also for the SSH model.

II. FORMALISM FOR LONG-TIME DYNAMICS OF ENTANGLEMENT
In this section we develop a formalism that yields a connection between dynamics of entanglement and classical dynamics on a high-dimensional torus.It is important that throughout this work we consider finite L, unlike canonical analytical studies of entanglement evolution [1,3] that address infinite L. Of course, the statistical theory developed here can be extended to the case of infinite L. We focus on the bipartite entanglement and the condition: 1 L A ≤ L/2 is assumed, so that the subsystem is smaller than its complement but large enough.
A. Relating evolving correlation matrix to classical trajectory on T N Due to a particle-hole symmetry and a reflection symmetry of the energy spectrum in Eq. ( 4) (see also Eqs. (B5) and (B6)), we see that the time parameter enters in the instantaneous correlation matrix C(t) in Eq. ( 3) only through associated with the energy gaps at the Bloch momenta k m .We actually have for the correlation matrix with the notation: (m → −m) denoting the term obtained from that in front of the + sign by making the replacement: m → −m.Mathematically, we can first introduce a N -variable (matrix-valued) function on the N -dimensional torus T N defined as with ϕ ≡ (ϕ 0 , . . ., ϕ N −1 ) ∈ T N .Then, we obtain the instantaneous correlation matrix C(t) from C(ϕ) via We shall call C(ϕ) the correlation matrix as well.Its physical meanings will become clear later.Note that by definition C(ϕ) is periodic in each argument ϕ m , Thus, these phases ϕ = ωt in C(ϕ) completely determine C(t) and entail a classical motion that is a rotation in the N -dimensional torus T N with constant angular velocity ω (cf.Fig. 1(d) in the main text): In other words, given the initial state Ψ(0), the correlation matrix corresponds one-to-one to a point ωt along the classical trajectory.This relation plays important roles in investigating fluctuations in long-time entanglement evolution.We remark that this relation makes no reference to the number-theoretical properties of ω, e.g., commensurate or incommensurate, but, as we will demonstrate later, the behaviors of long-time entanglement evolution are very sensitive to these properties.We note that similar correspondence also holds generally for the evolution of quantum expectation value Ψ(t)|A|Ψ(t) of one-body operator A upon a global quench, see Appendix C.

Reduced density of matrix and entanglement entropy
The time evolution of the correlation matrix completely determines the entanglement evolution.Indeed, the instantaneous reduced density of matrix provides full information on entanglement evolution.It is given by with the trace restricted to the subsystem A, where H A (t) is the instantaneous entanglement Hamiltonian.Because the systems considered are noninteracting and fermionic and the evolving state Ψ(t) is Gaussian, H A (t) takes a quadratic form [4]: and the 2L A × 2L A matrix H A (t) is determined by the correlation matrix C(t) via So we see that the instantaneous reduced density of matrix is a functional of the instantaneous correlation matrix.Moreover, all 2L A instantaneous eigenvalues of C(t), denoted as p ν (t) (ν = 1, . . ., 2L A ), belong to the interval [0, 1].Physically, p ν (t) gives the occupation probability of the instantaneous eigenmode ν of the correlation matrix.
In most of this work, we study the entanglement entropy, which is the von Neumann entropy associated with the reduced density matrix and defined as

S(t) ≡ −Tr
Combining it with Eqs. ( 13)-( 15) we find that or equivalently, So, it is again a functional of C(t), like ρ A (t). Furthermore, with the introduction of the binary entropy function one can readily express Eq. ( 17) as For the reduced density matrix, we also have a relation similar to Eq. (10).Indeed, due to Eqs. ( 10) and (15) we have where HA (ϕ) may be regarded as the entanglement Hamiltonian associated with C(ϕ).Then, upon introducing the following reduced density of matrix associated with C(ϕ): Tr A e − HA (ϕ) , we obtain By the same token, we can introduce the following entanglement entropy associated with C(ϕ): then by Eq. ( 10) the relation: follows.Since all entanglement probes associated with C(ϕ) such as HA , ρA and S are obtained from their counterparts associated with C(t) by the replacement: C(t) → C(ϕ), hereafter we use the same symbols: H A , ρ A , S, etc. as corresponding probes associated with C(t), keeping in mind the difference in the arguments.

More general entanglement probes
Let us consider a larger class of entanglement probes or information-theoretic observables O(t), which are required only to be a functional of C(t), i.e., like Eqs. ( 15) and (18).In words, O(t) depends on the time parameter only through C(t), albeit its dependence on C(t) is highly nonlinear in general.Therefore, O(t) and C(t) can be diagonalized simultaneously at any t.Besides the entanglement entropy and the reduced density of matrix, such entanglement probes include the n-th order Rényi entropy and the entanglement spectrum, namely, the eigenvalues of the reduced density matrix labelled by m = {n ν }, which is the configuration of occupation numbers n ν (= 0, 1) at the eigenmode ν.Owing to the relation (10), O(t) must depend on t via ϕ = ωt, i.e., in the same fashion as Eqs.( 23) and (25).That is, with the introduction of the N -variable function: O(ϕ) ≡ O[ C(ϕ)], the relation: then follows.Note that by definition O(ϕ) is 2π-periodic in each ϕ m also.
In order to appreciate more the importance of C(t) to the time evolution of generic entanglement probe more, we proceed to find explicitly the functional dependence of S n (t) on C(t).First of all, because of − ln Moreover, from Eqs. ( 14) and ( 15) we obtain Then, by substituting Eqs. ( 30) and ( 31) into Eq.( 27), we obtain where In the limit: n → 1, e n (λ) → e(λ).Therefore, we recover the expression for the entanglement entropy (20).
Observing the right-hand sides of Eqs. ( 20) and ( 30)- (32), we find that they have the common structure: and differ only in the function O(λ).Remarkably, this structure relates the time evolution of different entanglement probes to the same quantity, the instantaneous spectral density of the correlation matrix.Likewise, all O(ϕ) can be related to the same quantity, the spectral density of C(ϕ), in the way as This expression bears a firm analogy to the expression for canonical quantities in mesoscopic physics [5][6][7][8][9][10], as we will demonstrate in Sec.IV C. So the evolving correlation matrix C(t) is the building block of our out-of-equilibrium fluctuation theory.Moreover, as we shall see in the next section, the relation ( 29) allows us to trade out-of-equilibrium entanglement fluctuations to mesoscopic fluctuations from one virtual disordered sample to another, which are seemingly unrelated to the dynamics of entanglement at first glance.Equations ( 10), ( 12), ( 29), ( 34) and ( 35) constitute a formalism connecting various information-theoretic observables, that characterize entanglement evolution, to classical rotation on T N .

III. STATISTICAL EQUIVALENCE
It is well known that dynamical properties of classical rotation on the torus T N are extremely sensitive to the number-theoretic properties of the angular velocity ω.Notably, if the equation: then ω 0 , . . ., ω N −1 are incommensurate or said to bear the arithmetic linear independence [11], which should be distinguished from the concept of statistical independence to be introduced shortly; if the condition (36) is not satisfied, the N frequencies are commensurate.For incommensurate ω the rotation is ergodic, while for commensurate ω it is periodic [12].By the dispersion relation (4) the angular velocity ω governing the evolution of C(t) (cf.Eq. ( 3)) is incommensurate.So we will focus on the incommensurate case from now on.In Sec.VIII we will study the consequences of commensurate ω.By definition O(ϕ) introduced in Sec.II B is 2π-periodic in each argument ϕ m .As a result, O(t) is a quasiperiodic function of t with frequency basis ω; in other words, after the initial process the time evolution O(t) transits to quasiperiodic oscillations [13], as exemplified by the entanglement entropy (Fig. 1 of the main text) and the secondorder Rényi entropy (see Figs. 1(a) and 1(b) below): These quasiperiodic oscillations are conceptually different from the entanglement entropy oscillations due to the traversal of quasiparticle pairs or the incomplete revival of system's wavefunction; the latter oscillations rapidly damp in the course of time [14].The quasiperiodic oscillations do not exclude the occurrence of a big bump or dip at certain times, and the big dip may signal the recurrence.It should be emphasized that the quasiperiodic oscillations of different probes, all of which can be attributed to those of the correlation matrix, are quantum oscillations (cf.Eq. ( 8)).
For quasiperiodic O(t) we have lim by the ergodic theorem [12].Here stands for the average with respect to the probability measure P over T N , which has a uniform density of 1/(2π) N .Therefore, O(t) quasiperiodically oscillates around O .The ergodic theorem further gives that for arbitrary interval ∆, lim Here the left-hand side gives the frequency for the time series O(t) to appear in ∆, while the right-hand side gives the probability for O(ϕ) to be in the same interval, i.e., P(O(ϕ) ∈ ∆).Consequently, out-of-equilibrium fluctuations displayed by the quasiperiodic oscillations O(t) must have the same statistics as fluctuations of O(ϕ) with ϕ, provided ϕ is drawn randomly from P. This statistical equivalence holds for any entanglement probes introduced in Sec.II B, and is confirmed numerically (Fig. 1a in the main text).More generally, for any set A ⊂ T N , we have lim with Eq. ( 39) as a special case.The right-hand sides of Eqs. ( 39) and ( 40) entail a uniform joint probability density, 1/(2π) N , for the N components of ϕ.Therefore, (with respect to this probability) these components are statistically independent, and each of them is uniformly distributed over the 1D torus T. That is, the above probability measure P is a product of N uniform probability measures over T. Such product structure of probability measures allows the applications of the theory of concentration of measure [15][16][17], and leads to important consequences later on.

IV. EMERGENCE OF MESOSCOPIC FLUCTUATIONS
The Rice-Mele model, the Su-Schriefer-Heeger model, and TFIC are all integrable.Besides, they are deterministic, namely, free of any extrinsic randomness or stochasticity.Despite of these, armed with the general principle established in the last section we show in this section that some canonical quantum chaotic phenomena can still emerge from the time evolution of the correlation matrix.Most importantly, we show that out-of-equilibrium fluctuations of various entanglement probes can be traded to their fluctuations from one virtual, mesoscopic disordered sample to another.

A. Emergent of ensemble of random correlation matrices
According to Eq. ( 9), the correlation matrix C(ϕ) can be separated into two parts, i.e., Each part inherits the block-Toeplitz matrix structure from C, i.e., (C 0,1 ) ii = (C 0,1 ) i−i =l with (C 0,1 ) l being 2 × 2 blocks defined in the sublattice sector.The first part, C 0 , has no dependence on ϕ and is well defined for L → ∞, whereas the second part, C 1 , depends on ϕ, and has no limiting behaviors for L → ∞ because of the incommensurality of ω.Using Eqs.(B17), (B18), (B21) and (B22), we find after lengthy but straightforward calculations that the coefficient of cos ϕ m (sin ϕ m ) in the summand is real (purely imaginary), which is a function of k m and denoted as R l,σσ (k m ) (I l,σσ (k m )).So Eq. ( 43) is rewritten as where R l,σσ (ϕ) and I l,σσ (ϕ) are the real and imaginary part of (C 1 (ϕ)) l,σσ , respectively.It is easy to show that for l = 0 and σ = σ the imaginary part vanishes, i.e., the diagonal of C 1 is real.By Eqs. ( 10) and ( 41), we have with Therefore, C(t) displays quasiperiodic oscillations around C 0 .
In Appendix D we show that, thanks to Eq. ( 40), both R l,σσ (t) and I l,σσ (t) or correspondingly R l,σσ (ϕ) and I l,σσ (ϕ) are sub-Gaussian centered random variables.More precisely, they have zero mean and their deviation probabilities satisfy the following concentration inequalities: lim for large L and for any positive .Here whose numerical coefficients have no dependence on L and thus are not important.Both η l,σσ and η l,σσ are positive.We see that C 1 (ϕ) varies around zero with variance ∼ 1/L, i.e., Physically, this scaling can be understood as follows.When ω is incommensurate, cos(ω 0 t), cos(ω 1 t), ..., cos(ω N −1 t) are N statistically independent random numbers of zero mean [11].As a result, Eq. ( 46), that corresponds to R l,σσ (respectively I l,σσ ), is a sum of N ∼ L random numbers.For large N one may expect the central limiting theorem to apply, which then gives R l,σσ ∼ √ N /L ∼ 1/ √ L and likewise for I l,σσ .In Appendix E we further show that any random variable R or I in the block (C(t)) l (respectively ( C(ϕ)) l ) is statistically independent of that in any distinct block (C(t)) l =l (respectively ( C(ϕ)) l =l ).In contrast, from Eq. ( 46) we see that random variables R, I in the same block (C(t)) l (respectively ( C(ϕ)) l ) are not statistically independent.
So, when the classical trajectory: ϕ = ωt sweeps out T N uniformly, the infinite time series C(t) generates a specific ensemble of random block-Toeplitz matrices C(ϕ), denoted as E. For each member of E, its ϕ-dependent part, C 1 (ϕ), is composed of L A statistically independent 2 × 2 blocks (C 1 (ϕ)) l (l = 0, 1, . . ., L A − 1), the real or imaginary part of whose matrix elements are sub-Gaussian centered random variables, with the variance ∝ η l,σσ /L, η l,σσ /L.The probability measure over E is induced by the uniform probability measure P over T N in the way as follows: Given arbitrary intervals ∆ r, i l,σσ for every l, σ, σ , the probability for C to satisfy: Re(C 1 ) l,σσ ∈ ∆ r l,σσ and Im(C 1 ) l,σσ ∈ ∆ i l,σσ is P(ϕ ∈ B), where the set B is defined as With this induced probability measure, the ensemble E is equivalent to the time series C(t) statistically.This ensemble differs from the ensemble of random hermitian Toeplitz matrices studied very recently [18,19] in two aspects.First, each member in E, though being hermitian and Toeplitz-type, carries a block structure; that is, the element of Toeplitz matrix is now a 2 × 2 block, instead of a complex number.Second, according to Eq. ( 51) at fixed ratio L A /L the variance of its matrix elements decays with the matrix dimension and vanishes in the large matrix dimension limit, whereas in canonical studies of random matrices that variance is independent of the matrix dimension.

B. Spectral statistics of correlation matrix
We proceed to study the statistical properties of the spectrum of the correlation matrices from either the time series C(t) or C(ϕ) drawn randomly from the ensemble E. Since they are statistically equivalent, we will mainly focus on the properties derived from the time series C(t).In this work it suffices to consider the nearest-neighbor spacing distribution, P 0 (s), with the eigenvalues rescaled by the mean eigenvalue spacing.
In order to spell out the block-Toeplitz structure in a transparent way, the correlation matrix C(t) can be rewritten in the following form: where with defined in the sublattice sector and with expressions of ũfk (t), ṽfk (t) given in the Appendix B. The statistical property C(t) is therefore determined by Γ(t).In parallel, we can generalize the above expressions Eqs. ( 54) and (55) to define the corresponding C(ϕ) = Owing to the structure described by Eqs. ( 54) and ( 55), the eigenvalues of Γ(t) come in pair: one positive and the other negative, with the same absolute value.Therefore, its spectrum splits into two subspectra symmetric with respect to zero (cf.Fig. 3(b) below): One consists of negative eigenvalues and the other positive.We focus on the former below.
First, we follow the method of Ref. [19] to estimate small-s behaviors of P 0 .For this purpose we need to consider only block-Toeplitz matrices of minimal dimension.Such matrices carry the structure of Γ but in lower matrix dimension, and are perturbed from the block-diagonal matrix, with each block being diag(ϑ, −ϑ) and ϑ < 0 being some fixed constant.To satisfy the requirements we set L A in Eq. ( 54) to be two and the diagonal elements of Π0 to be ±(ϑ + δϑ), with δϑ being real.Because of L A = 2 the subspectra consist of two eigenvalues.By analysis in Sec.IV A, these matrices have only two independent random variables, one from the block Π0 and the other from Π1 .So the joint probability of the two eigenvalues can be expressed as an integral of these two variables.Thus the probability for both eigenvalues to be within a small distance of s from ϑ is ∼ s 2 , with the power 2 accounting for the number of independent random variables.Since the joint probability is the cumulative spacing distribution s 0 P 0 (s )ds , we find that Thus P 0 vanishes at s = 0.This is the so-called level repulsion phenomenon -a defining property of quantum chaos [20].Strikingly, here the phenomenon occurs in a peculiar system, such that it is integrable and free of any extrinsic randomness and thus has nothing to do with chaos, at first glance.
Next, we estimate large-s behaviors of P 0 .For this purpose we return to the original matrix, and note that the matrix elements in Πl and Πl decay rapidly with l, namely, the distance to the main diagonal for which l = 0, σ = σ .So Γ and Γ are random band matrices, and their matrix elements are short-ranged correlated since as shown in Appendix E the matrix elements of distinct blocks are statistically independent.Recall that for random band matrices with all their elements being statistically independent, it is generally believed that the eigenvectors may exhibit localization [21].But it is not clear to us whether this might occur here, because the random Toeplitz band matrix is structured, i.e., all its elements of the same distance to the main diagonal are identical, and additionally the disorder strength has a size dependence shown in Eq. ( 51).From the Toeplitz structure we can expect only that the level repulsion exists up to the scale of mean eigenvalue spacing.So we can follow the arguments of Ref. [22] to estimate large-s behaviors.Supposing that a band of large width s includes N (s) eigenvalues on the average, by the central limiting theorem we find that the variance of the number of eigenvalues within this band Var(N (s)) ∼ N (s) .Since the probability for the band to include N (s) eigenvalues is ∼ e −(N (s)− N (s) ) 2 /Var(N (s)) , the probability for the band to be empty is ∼ e − N (s) 2 /Var(N (s)) ∼ e −const.N (s) , with const.being some positive constant.Since P 0 (s) is the probability for the band to be empty and N (s) = s, a Poissonian tail, P 0 (s 1) ∼ e −const.N (s) = e −const.s(58) then follows.The full Poissonian distribution was conjectured for random real symmetric Toeplitz matrices in Ref. [23].Finally, the limiting behaviors (57) and (58) can be unified via the following simple form: with the numerical coefficients fixed by the conditions: The distribution (59), that holds for arbitrary s, is called semi-Poissonian distribution [24].It was seen for the first time in Ref. [18] for truly random hermitian Toeplitz matrices, but without the block structure.This distribution is confirmed in simulations of the spectral statistics of the evolving correlation matrix C(t) (see Fig. 1(c) in the main text).

C. A new class of mesoscopic fluctuations
We have shown that for long time various entanglement probes O(t) quasiperiodically oscillate around their means O .This implies that when the entanglement evolution enters into the stationary state at long time, various probes are not time independent but display fluctuations around the stationary value.As shown above the statistics of these out-of-equilibrium fluctuations is equivalent to that of the fluctuations of the same probe with ϕ, provided ϕ is drawn from T N with a uniform probability.Now we show that these fluctuations with ϕ connect the entanglement evolution and mesoscopic physics.
To this end we observe that the fluctuations of entanglement probes with ϕ have following prominent features.First, because a generic probe O(ϕ) is a functional of the correlation matrix C(ϕ), the fluctuation statistics of different probes can all be attributed to that of C(ϕ) or more precisely its spectral density.Second, because the random, namely, ϕ-dependent part of C(ϕ) vanishes in the limit L → ∞, these fluctuation phenomena can occur only for finite L A at fixed ratio L A /L (> 0).Third, these phenomena inherit the quantum coherence nature from quasiperiodic oscillations.By these three features the fluctuations of various entanglement probes with ϕ fall into the paradigm of mesoscopic sample-to-sample fluctuations [5][6][7][8][9][10]; see Table I for comparisons.Indeed, in mesoscopic electronics or photonics various transport characteristics (e.g. the conductance or transmittance, the shot-noise power, etc.) can be expressed in the same form as the right-hand side of Eq. (35).Specifically, upon passing to various mesoscopic transport probes, the hermitian matrix: C(ϕ) is replaced by the hermitian matrix: T ≡ tt † , whose dimension gives the number of channels for wave transport in disordered samples and is finite also [25,26].Here t is the transmission matrix that encodes all information of wave propagation from the input to the output end.Similar to that C(ϕ) is determined by ϕ randomly distributed over T N , the matrix t and thereby T are determined by the disordered potential or the dielectric configuration, denoted as V ; furthermore, since typically the potential and dielectric fluctuations at distinct spatial points are independent and Gaussian, each disorder realization can be viewed as a vector in a highdimensional Euclidean space, V ≡ (V r1 , V r2 , . ..) with r 1 , r 2 , . . .labelling distinct spatial points, Gaussian distributed in that space [27].So, similar to that varying ϕ gives rise to an ensemble of C(ϕ), varying V gives rise to an ensemble of random matrices T (V ), and the distribution of V induces the distribution of T (V ), albeit in a highly complicated manner.With the replacement described above, the right-hand side of Eq. ( 35) gives various transport characteristics, e.g. the conductance or transmittance for O(λ) = λ and the shot-noise power for O(λ) = λ(1 − λ) [10].As such, T (V ) represents a finite disordered sample.In the same fashion, as C(ϕ) determines various entanglement probes O(ϕ), it mimics a finite disordered sample, with ϕ being the 'disorder realization'.So the fluctuations of various entanglement probes with ϕ are traded to mesoscopic sample-to-sample fluctuations.Furthermore, since the former fluctuations are statistically equivalent to the out-of-equilibrium fluctuations displayed by quasiperiodic oscillations, those oscillations have the same nature as the mesoscopic sample-to-sample fluctuations.
However, standard random matrix theories require the variance of matrix elements to be independent of the matrix dimension.This is totally opposed to the present situations: Equation (51) shows that at fixed ratio L A /L > 0 the variance of the matrix elements of C(ϕ) decays with the matrix dimension 2L A as ∼ 1/L A , so that nonrandom C(ϕ) results in the limit L A → ∞.Moreover, analytical or mathematical studies of random Toeplitz matrices are highly challenging, as put explicitly by Bogolmony and Giraud: "seem to be inaccessible to known analytic methods" [19]; even for the simplest random Toeplitz matrices, only recently have the studies of their spectral statistics been initiated and very few results been reported [18,19].For these reasons, we are prohibited to invoke standard mesoscopic theories to study the mesoscopic sample-to-sample fluctuations of entanglement probes.Therefore the long-time dynamics of entanglement brings us a new class of mesoscopic fluctuations.

V. STATISTICS OF EMERGENT MESOSCOPIC FLUCTUATIONS: GENERAL THEORY
To overcome the difficulties discussed in the end of last section with the unusual disorder structure, below we develop an approach to study the statistics of mesoscopic fluctuations emergent from entanglement evolution.The idea starts from that a generic evolving entanglement probe O(t) is given by a N -variable function, O(ϕ), over T N according to Eq. ( 29).So, when T N is equipped with uniform probability, as imposed by Eq. (40), O(ϕ) becomes a random variable depending on N -which is finite but large -statistically independent angular variables ϕ 0 , ϕ 1 , . . ., ϕ N −1 .Note that, however, the dependence of O on those angular variables is highly nonlinear.Then, taking the advantage of the very fact that N (i.e., L) is finite, we can combine the continuity properties of the function O(ϕ) with the nonasymptotic probability theory, so-called concentration inequality [28], to study the statistical properties of fluctuations of O with ϕ.A similar idea has recently been developed to find the universal wave-to-wave fluctuations in mesoscopic transport [27].

A. Concentration inequalities
To implement the idea above, we introduce the logarithmic moment-generating function defined as for a generic entanglement probe O.Because ϕ 0 , ϕ 1 , . . ., ϕ N −1 are independent and identically distributed random variables, we then have the so-called modified logarithmic Sobolev inequality [28], read Here O m is an arbitrary function of ϕ m ≡ (ϕ 0 , ϕ 1 , . . ., ϕ m−1 , ϕ m+1 , . . ., ϕ N −1 ), independent of ϕ m .By definition it is also a random variable, but depending on (N − 1) instead of N independent ϕ's.It is important that inequality (61) holds for any N .This is fundamentally different from the situations in traditional probability theory, where the number of independent random variables has to be sent to infinity eventually.In other words, the traditional probability theory is a limiting theory, whereas the concentration inequality is not.In this section we need two special classes of O m 's, defined as the infimum (supremum) of O(ϕ 0 , . . ., ϕ m , . . ., ϕ N −1 ) over ϕ m with the rest of (N − 1) unprimed variables held fixed.The function is denoted as O + m (O − m ) correspondingly.For the convenience below we also introduce the following quantity: Roughly speaking, it is the total variance of the probe O with respect to each argument.

The upper tail of the distribution
Following Ref. [28], we set u > 0 and O m = O + m for inequality (61) to study the distribution of the upward fluctuations (i.e., O − O > 0).However, our subsequent treatments of that inequality differ substantially from those in Ref. [28].Observing that φ −u(O − O + m ) can be Taylor expanded as we separate from φ the contribution: to the first term, and rewrite inequality (64) as 2 from its mean.Substituting Eq. ( 65) into the right-hand side of inequality (61), we arrive at where It is important to note that although the right-hand side of inequality (66) bears some similarities to familiar ones in mathematical literatures [28], there are essential differences.In particular, we do not resort to the bounding property of the function: , and consequently the first term makes no references to the bounding properties of that function.Whereas in mathematical literatures b + is determined by the bounding properties of that function.
Next, noting that O − O + m ≥ 0 and the equal sign is taken only at a set of zero (Lebesgue) measure, we find that the right-hand side of inequality (61) As a result, there must exist some constant const.and some sufficiently large u * , so that This implies that the right-hand side of the modified logarithmic Sobolev inequality is bounded by the curve ∼ 1 u in the large u regime.On the other hand, it is easy to show that for u > 0, dG(u)  du > 0 and grows unboundedly with u.Taking this and inequality (69) into account, we find that the function δF/ dG du must be bounded for u > 0. Let its supremum be c + .We obtain As we shall see immediately, that dG/du is positive definite (for u > 0) entails the sign of c + important consequences.
Combining inequalities (66) and (70), we arrive at the following differential inequality: Moreover, by Eq. ( 60) this differential inequality is implemented by the boundary condition: G(u) u = 0 at u = 0. Now, let us study the cases of c + < 0 and c + > 0 separately: It can be solved readily.The result is Then it is a standard exercise to use the Markov inequality [28] to obtain the concentration inequality for the upper tail of the distribution.Specifically, by the Markov inequality we have for arbitrary positive and u, Combining it with inequality (73) we obtain Since u is arbitrarily positive, upon minimizing the exponent over positive u it gives This concentration inequality defines a sub-Gaussian upper tail distribution.
• c + > 0: We solve inequality (71).The result is This kind of bounds hold also for Gamma random variables, and thus generalize the tail behaviors of the Gamma distribution, giving the so-called sub-Gamma tail distribution [28].Specifically, by the same method of deriving inequality (76), we obtain The constants b + and c + are called the variance factor and scale parameter, respectively [28].Of course, it does not mean that b + gives the exact variance, but gives a bound for the exact variance in general.In the present work, by the choice of b ± defined by Eq. ( 63) (b − is for the lower tail; see concentration inequalities ( 80) and (81) below), which are substantially different from canonical choices in mathematical literatures, we further find that b ± are proportional to the exact variance, with the help of numerical simulations (see Sec. VI E for detailed discussions).Finally, we remark that as increases, the sub-Gamma bound displays a crossover from the Gaussian form e − 2 /(2b+) to the exponential form e − /(2c+) form at ∼ b + /c + .

The lower tail of the distribution
The distribution of downward fluctuations (i.e., O − O > 0) can be studied in the same way.We set u < 0 and O m = O − m for the modified logarithmic Sobolev inequality (61).With these two replacements, we have with dG/du < 0, similar to the bounding (70).The subsequent analysis can be carried out in exactly the same way.But the results are reversed: i.e., a sub-Gamma lower tail.

B. Calculations of the variance factors b±
In this part we study the variance factors b ± defined by Eq. ( 63) in details.First of all, we keep in mind that O ± m (ϕ m ) in Eq. ( 63) actually belongs to a special class of (real-valued) functions of ϕ, such that they are independent of Note that by definition ϕ ± m depends on ϕ m as well as O.Then, by the mean value theorem there exists φ± m between ϕ m and ϕ ± m which depends on ϕ, so that Here the subscript: ϕ m = φ± m stands for that the arguments of ∂ ϕm O are (ϕ 0 , . . ., φ± m , . . ., ϕ N −1 ), and to make the formula compact we suppress all the arguments on the right-hand side except φ± m .From Eq. ( 82) follows immediately.
Then we perform the Fourier expansion of ∂ ϕm O with respect to ϕ m .Taking into account the general expression of O (cf.Table I) and C 1 = O(1/ √ L), we can expand O in C 1 for large total system size L 1. Keeping the expansion up to the second order (see Sec. VI A for explicit calculations for the entanglement entropy) gives where we have made use of Eq. ( 44) and the coefficients a's, b's are ϕ m independent.We can rewrite Eq. ( 84) as where A's and φ's are determined by a's and b's.With the help of this expression it can be readily shown that and Inequality (86) and Eq. ( 87) give Actually the factor 6 can be improved but this is not important.This bound suggests that (with the components ϕ n =m treated as parameters) (∂ ϕm S) 2 at given ϕ m is the same order as its ϕ m -average, i.e., Substituting it into Eq.( 83) we obtain where ϕ ± * m ≡ (ϕ 0 , . . ., φ± m , . . ., ϕ N −1 ).Because the first factor characterizes the distance between ϕ m and ϕ ± m and the second describes the ϕ m average of (∂ ϕm O) 2 , we assume that they are statistically uncorrelated and obtain Furthermore, by symmetry (ϕ m − ϕ ± m ) 2 r m (ϕ ± * m ) must be the same for distinct m, at least approximately.As a result, we reduce Eq. ( 91) to where is an overall numerical coefficient.As mentioned above, in the present work owing to the special choice of b ± , which are given by Eq. ( 63), b ± are proportional to the variance Var(O).As a result, where This formula relates fluctuations of the entanglement probe O to the continuity of the N -variable function O(ϕ) defined by (35).As such, the mesoscopic fluctuations emergent from entanglement evolution fall into the class of the so-called concentration-of-measure phenomena, where fluctuations of an observable are controlled by its Lipschitz continuity.The concentration-of-measure phenomena are a new perspective of probability theory [15], and are first reported for mesoscopic electronic and photonic transport in Ref. [27].In Fig. 3(a) of the main text, we provide the numerical verification of Eq. ( 94) for the entanglement entropy.In Fig. 2(a) we provide further numerical verifications for the n-th Rényi entropy S n .

B. The contribution Var1(S)
To calculate Eq. ( 101) we note that the trace operation can be cast to Tr A (•) = L A i=1 σ (•).That is, we perform that operation in two steps: In the first step, we sum over the sublattice index σ = Ā, B; in the second, we sum over the cell index.Now, because the matrix block (C 0 ) ij decays rapidly with |i − j|, we can divide the sum over i into two contributions, one from the bulk of the subsystem A and the other from its edge.The size of the edge, L e , is order of the decay length of C 0 , and thus is determined merely by the parameters of the Hamiltonian; while the size of the bulk is ≈ L A .The bulk and the edge of the subsystem are denoted as A b and A e , respectively.The cells in A b are labelled by i b and in A e by i e .So, Below we calculate the two terms separately.
For the bulk term, owing to the rapid decay of C 0 we can extend all the intermediate cell indexes involved in the matrix product to the total system.Taking this into account, we find that where the trace tr is restricted to the sublattice sector and I is the unit matrix in that sector.Moreover, we introduce the matrix defined in the sublattice sector, αm with matrix elements α σσ (k m ), and likewise βm and γm .
With the introduction of Γm ≡ γm − I/2, we rewrite Eq. ( 104) as Then we make use of the properties (i)-(xii) of γm , αm and βm given in Appendix F to calculate Eq. ( 105).
When we expand the logarithm in Γm in Eq. ( 105) and perform the trace, we find that each γ Ā B factor is paired with a T B Ā (T = γ, α (or β)) factor.Therefore, all the phase factors e iδ (k) , that are carried by the off-diagonal components of γm , αm and βm (cf. the properties (ii), (iv) and (vi)), cancel out.Therefore, the expansion corresponds, term by term, to the Γ m -expansion of tr − ln where Γ m , α m and β m differ from Γm , αm and βm only in that their off-diagonal components have no phase factor.Note that Γ m is a real symmetric matrix.Taking this and the properties (v) and (vi) into account, we find that the coefficient of cos ϕ m vanishes for m > 0; and for m = 0 taking this and the properties (xi) and (xii) into account, we find that the coefficient also vanishes.So the second term in Eq. ( 106) vanishes, and thus Eq. ( 105) is simplified as with αm ≡ α m + α −m , where we have used the properties (i) and (ii) to derive the second line.
To calculate Eq. ( 107) we expand the logarithm in Γ m , and consider tr(( Γ m ) r αm ) for any r ∈ N. To this end we note that Γ m and αm have the following general structure: Here all matrix elements are real and their explicit expressions can be easily found by using Eqs.(B16), (B17), (B20) and (B21).In Appendix G we prove the following identity: With its help one can readily show that As a result, So, only the edge term contributes to Tr A ln C −1 0 − I ∂ ϕm C 1 .For the edge term, we have similar to Eq. ( 104).Note that the edge size L e does not scale with L A .Combining this result with Eq. ( 111) we obtain where Then we substitute Eq. (113) into Eq.(101).For L 1 the sum over m converges to the integral over the momentum.We finally obtain where the proportionality coefficient It is important that a depends neither on L nor on L A .

C. The contribution Var2(S)
For the moment let us ignore the second line in Eq. (102) in order to simplify discussions.Because (( Here κ ≡ ((I − C 0 )C 0 ) −1 ii which is independent of i.With the substitution of Eq. ( 44), it is written as where are matrices in the sublattice sector and are even in m, n.Taking the square of Eq. ( 118), we obtain 16 terms.They have the form as Here g 1,2 and h 1,2 stand for the symbols of functions: cos, sin.It is easy to see that only the four terms, with g i = h i (i = 1, 2) and n = n , dominate the ϕ-averaged square.The average of other terms either vanishes or is smaller by an order of 1/L.Taking these into account, we have where the last three terms in the bracket are obtained by replacing Jcs in the first term by respectively Jsc , Jcc and Jss , and we have used the fact that the leading contributions to the sum over m, n come from those terms with nm = 0. Next, we perform the sum over the indexes i, j, i , j .Let us observe each trace product in Eq. (120).By Eq. (119) the first trace includes four terms.Each term carries a phase factor, and the four factors are different, which are e ±i 2π(n−m)(i−j) L and e ±i 2π(n+m)(i−j) L .Similarly, each term in the second trace carries a phase factor, and the four factors are e ±i 2π(n−m)(i −j ) L and e ±i 2π(n+m)(i −j ) L .Therefore, the trace product has 4 × 4 = 16 terms, each of which carries a phase factor e iΞ , with the phase Due to n ± m = 0 the phase factor e iΞ rapidly oscillates with i, j, i , j .As a result, the sum over i, j, i , j is dominated by the 'stationary phase' configuration, for which Ξ vanishes, requiring η 2 = η 2 and η 1 = η 1 , (i − j) + (i − j ) = 0 (or where For L 1 the sum over m, n in Eq. ( 122) can be well approximated by the continuum limit.We obtain It is important that b has no dependence on L A and L, and is completely determined by system's microscopic parameters.Now let us retrieve the second line in Eq. ( 102).We repeat the derivations above.As a result, we find that only the coefficient b in Eq. ( 124) is modified and, similar to b , the modified coefficient b has no dependence on L A and L. Finally, we have Interestingly, the scaling behavior ∼ L 3 A /L 2 coincides with a previous result about the entanglement entropy variance in a completely different context [29].In that context fluctuations are due to randomly drawing a member from a free-fermion eigenstate ensemble; thus that result is a purely kinematic consideration.

D. Universal scaling law for Var(S)
Combining Eqs. ( 115) and (126), we obtain the variance of the entanglement entropy, Importantly, from the derivations above we have seen that the first and second term describe a subsystem's edge and bulk effect, respectively.Upon the rescaling: L ≡ L/ , LA ≡ L A / ( ≡ a/b), Eq. ( 127) is rewritten as with s 0 ≡ √ ab, which is Eq. ( 2) in the main text (where we used the same symbols L, L A for notational simplicity).This result is universal with respect to system's detailed constructions and initial states (required to be Gaussian, however), that enter only into the microscopic parameters and s 0 .It is thus suggested that the behaviors of entanglement entropy fluctuations are completely controlled by two dimensionless macroscopic lengths, L and LA .
Equation (128) implies that in the regime of LA L1/3 , the edge term dominates over the bulk term and Var(S) ∼ L−1 .Thus in this regime the fluctuations are very weak.In the special case of L → ∞ at fixed LA , the variance vanishes.Combined with the fact that S(t) oscillates quasiperiodically, this result implies that fluctuations are fully suppressed and the entanglement entropy is strictly a constant beyond some critical time, in agreement with a celebrated result in entanglement evolution [1,3]; see Appendix H for details.In the regime of LA L1/3 the fluctuation behaviors are totally opposite.In this case the bulk term dominates over the edge term and Var(S) ∼ L3 A / L2 .So when the ratio LA / L = L A /L is fixed, the variance increases linearly with L, and thus enlarging the entire system can drive very strong entanglement fluctuations.This is contrary to the belief [30] that temporal fluctuations observed in experiments on entanglement evolution of finite systems would eventually diminish by increasing the system size.Thus, we see that the behaviors of out-of-equilibrium entanglement entropy fluctuations can be completely different in approaching the limit L → ∞, depending on how L A scales with L.
The situations above resemble those in the conductance fluctuations of quasi-one-dimensional disordered wires in several aspects.First, the conductance fluctuations are also controlled by some dimensionless macroscopic parameter, namely, the sample length rescaled by the localization length, and the microscopic details of the wire, notably, the disorder strength, enter only into the localization length.Second, varying the rescaled sample length leads to distinct fluctuation behaviors: Increasing the sample length drives the wire from a metallic regime, where the universal conductance fluctuations [5,6] follow, to a localized regime, where the conductance distribution is very broad so that the variance of logarithmic conductance increases with the sample length linearly [10].

E. The distribution tail
By the general results obtained in Sec.V A, the upper tail of the distribution is controlled by two parameters (b + , c + ) and the lower by (b − , c − ).While to calculate c ± analytically is an intractable task, we resort to numerical analysis.We compute the function δF/ dG du with the use of ab initio data of S(ϕ) obtained from numerical simulations (see Sec. VII for details).Then, according to Eqs. ( 70) and (79) the maximal value of δF/ dG du over positive u gives c + while the minimal value over negative u gives c − .We also compute b ± numerically by directly using the definitions (63).Table II gives the numerical value of (b ± , c ± ) for different L and the ratio L A /L.
We find that for all L A /L, L considered c ± are negative.(The only exception is that for L A /L = 0.5, c + is slightly positive, and the reasons are under investigations.)Consequently, by the results obtained in Sec.V A the distribution displays a sub-Gaussian upper and a sub-Gamma lower tail, described by concentration inequalities (76) and (80), respectively.We also see that b + ≈ b − for all L A /L, L, consistent with the analysis in Sec.V B. Furthermore, Table II shows that b − decreases rapidly as L A /L decreases from 0.5, whereas c − does not change too much.This implies that as L A /L decreases the distribution becomes more and more asymmetric.Specifically, for L A /L = 0.5 it is near symmetric and near Gaussian; as L A /L decreases the upper tail is always sub-Gaussian, and with more and more weights transferred from the lower to the upper tail, a heavier and heavier sub-Gamma lower tail develops.
The tail behaviors described above can be understood in a simple picture based on the concept of the coherent entangled quasiparticle pair [1] as follows.Upon quench such pairs are locally created across the whole system.Because they have opposite group velocities and, moreover the system is finite, so that a quasiparticle reenters into the system after reaching system's boundary, after long time a stationary configuration of entangled pairs is form.In this configuration the pairs are randomly distributed in the system subjected to Pauli's exclusion principle.Let δ be the probability for a quasiparticle inside subsystem A to pair with another outside -an event contributing to the entanglement entropy S, and (1−δ) 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., δ = 1/2; it is well known that in this case the binomial distribution reduces to a Gaussian distribution.For small L A /L, δ ≈ 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 (76) and (80) can be promoted to equalities, i.e., for sufficiently large , with the parameters b ± ∝ b ± and c ∝ |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 B we provide further simulation results to show that Eq. ( 129) holds also for the second-order Rényi entropy S 2 .Let us further discuss semi-quantitatively how the rate of the exponential decay ∼ e − 2c ( b − /c) of the lower tail behaves, when L approaches the infinite and the ratio L A /L is fixed.First of all, Eq. (127) implies that Var(S) ∼ L for sufficiently large L and fixed L A /L.Then, expanding dG/du in u we obtain dG/du = uVar(S) + O(u 2 ).So the coefficient of the leading term is ∝ L, and thus the leading expansion of L −1 dG/du is well defined in the limit L → ∞.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. (67), δF (u) is a sum of N ∼ L terms.Thus we estimate δF (u) ∼ L. Finally, combining the estimations for dG/du and δF (u) with the inequality (79), we find that c = O(L 0 ), i.e., depends only on the ratio L A /L.This is in agreement with the data in Table II.

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

A. The entanglement entropy S
We first diagonalize numerically for each time t the correlation matrix C(t) of Eq. ( 3) to obtain its eigenvalues {p ν (t)} 2L A ν=1 .We then substitute them into Eq.( 17) 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 /J, where J is the amplitude of hopping between nearest Āand B-site.
When the eigenfrequencies ω 0 , ω 1 , . . ., ω N −1 are incommensurate, generic for the energy eigenspectrum, the statistics of the time series S(t) and the random function S(ϕ) are equivalent, as shown by Eq. ( 39).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 ω 0 , ω 1 , . . ., ω N −1 are commensurate, this statistical equivalence does not hold; see Sec.VIII.

B. The nth-order Rényi entropy Sn
As an example for more general entanglement probes discussed in Sec.II B, we numerically compute the nth-order Rényi entropy according to Eqs. (32) and (33).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 (94) for S n , for n = 2, 3, 10.Even more surprisingly, simulations show that the scaling law (127) holds also for the variance 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.
C. Numerical computation of (b ± , c ± ) In Sec.V A, we introduce the variance factor and the scale parameter (b ± , c ± ) in solving the modified logarithmic Sobolev inequality (61) for the logarithmic moment-generating function G(u) of S(ϕ).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 II.For b ± , following Eq.(63): • We realize a disorder realization ϕ and compute S(ϕ), in precisely the same way as what is described above.
• With the same ϕ, we repeat the second step for each component and obtain N values of S ± m , and then determine • We repeat the previous three steps for 2 × 10 3 disorder realizations ϕ and average the outcome.
To obtain the scale parameter c ± , we numerically evaluate the two quantities δF (u) and dG(u)/du for a range of u around the origin.c + is determined by the maximum of δF/(dG/du) in the u > 0 interval, whereas c − is determined by the corresponding minimum in the u < 0 interval.

D. 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.

VIII. ENTANGLEMENT EVOLUTION WITH COMMENSURATE ω
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., ω is incommensurate.In this section we study the commensurate case.In this case, by Eq. ( 8), C(t) is periodic in time.Thus all entanglement probes defined by Eq. ( 26) 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 ω m is set to ω 0 + m r, with r being a suitably chosen rational number and ω 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.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., With the help of Ψ(t) = e −iH f t/ Ψ(0) as well as Eqs.(B7) and (B8), the matrix elements of the correlation matrix C iσ,i σ (t) can be computed to give with the summation over all Bloch momenta, since we take the initial state as the half-filling ground state of the pre-quench Hamiltonian where the subscript 0 denotes the eigenstate of H 0 .For σ = σ , with the help of Eq. (B12) we find that Inserting it into Eq.(B10) we get For σ = σ , with the help of Eqs.(B12) and (B13) we find that where v ± f k = e iφ f k v ± f k , and φ f k is defined in the same way as φ k but with post-quench Hamiltonian parameters.Inserting it into Eq.(B11) we get and δ(k) ≡ k/2 + φ f k .We note that u ± f k , v 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 Ψ(0) to a state Ψ(t) at time t, given by Ψ(t) = e −iH f t Ψ(0). (C1) In the basis of many-body eigenstates Ψ m (labelled by m), Eq. (C1) is written as where E m 's are the many-body eigenenergies and w m 's are the superposition coefficients of Ψ(0).Under this evolution, the expectation value of a generic operator A at time t is So the time parameter enters the evolution Ψ(t)|A|Ψ(t) through a large number ∼ L 2 of dynamical phases: (E m − E m )t with E m = E m .Here 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][33][34][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. (C3), 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 m corresponds to a configuration, {n ζ k }, of the occupation number n ζ k = 0, 1 at the single-particle eigenstate (u ζ f k , v ζ f k ) T of H f , with ζ = + (−) 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 Ψ(0) is a half-filling ground state, the configurations corresponding to the many-body eigenstates of H f superposing Ψ(t) -not all the many-body eigenstates -must satisfy n + k + n − k = 1, for every Bloch momentum k, (C4) and the corresponding many-body eigenenergies are given by where to derive the first equality we have used the reflection symmetry (B6) and to derive the second we have used Eq.(C4).Taking these into account we simplify Eq. (C3) to Note that due to Eq. (C4) the many-body eigenstate m (or m ) here is completely fixed by the occupation number configuration of the hole band {n − k } ≡ m − (or {n − k } ≡ m − ).In Eq. (C6), the dynamical phases are governed by the frequencies associated with the energy gaps 2E f k of the single-particle band structure and the difference in hole numbers n − k − n − k between two states m, m .Then we apply Eq. (C6) to one-body operators, e.g. the correlation function C(t) which belongs to this class.Such operators take the general form as Taking into account the symmetries (B5) and (B6), we thus arrive at the same physical picture as in Sec.II A. That is, for a generic one-body operator A, i.e., the quantum expectation Ψ(t)|A|Ψ(t) corresponds one-to-one to a point ωt along the classical trajectory on T N .
Appendix D: Proof of the concentration inequalities (47) and (48) We shall prove inequality (47) only.Inequality (48) can be proved in the same way.Since L and N -the number of independent random variables ϕ m -are finite, we adopt the theory of concentration of measure [15][16][17].
First of all, because R l,σσ (t) is a quasiperiodic in t, applying Eq. ( 37) gives

FIG. 2 .
FIG.2.Entanglement entropy distribution.We perform statistical analysis of the temporal fluctuations in the simulated entanglement entropy evolution.a. Variation of the distribution with increasing L at fixed LA. b.Same as a., but with increasing LA at fixed L. c.The large deviation probability P(|S − S | ≥ ), with upper and lower tail respectively, is well fitted by Eq. (B1) (dashed lines), implying that the upper (squares) tail distribution is sub-Gaussian and the lower (circles) is sub-Gamma.The ratio LA/L is 0.1 (yellow), 0.2 (green) and 0.5 (blue).

2 FIG. 3 .
FIG. 3. Universal scaling behaviors of the variance.We perform simulations for both Rice-Mele (RM) model and transverse field Ising chain (TFIC) with different sizes and different quench protocols (I-III) to study the variance of two entanglement probes O = S, S2. a.For both O the data confirm the relation (B8).b-c.They also confirm the limiting scaling behavior described by the first (second) term of Eq. (1) for sufficiently small (large) L 3 A /L. d.After rescaling Var(O), L, LA all data collapse to the universal curve described by Eq. (1).All theoretical predictions are presented by dashed lines.

FIG. 1 .
FIG. 1. Numerical simulations of the time evolution of the second-order Rényi entropy S2(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 LA is shown in (c).For different ratios LA/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. (129) (dashed lines), respectively (d).The subsystem size LA = 25 and for (a,b) the total system size L = 124.

FIG. 2 .
FIG. 2. Simulation results for the variance of nth-order Rényi entropy Sn, n = 2, 3, 10.For each Sn, we vary either LA 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 (94).(b)-(d) They confirm the first (b) and the second (c) terms, and the overall (d) of the scaling law (127).All theoretical predictions are presented by dashed lines.In (d) Var(Sn), L, LA are all rescaled in the manner similar to what is performed in the scaling analysis of Var(S).

FIG. 3 .
FIG. 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.

FIG. 4 .
FIG. 4. (a) Numerical simulations confirm that, when ω is commensurate, the pattern of the entanglement entropy is periodic and displays complete revival at multiple periods (yellow curves).The pattern in the incommensurate ω 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 ωm = ω0 + m r, ω0 = 3.162, r = 0.043, m = 0, . . ., 62.The dash lines are the physical energy spectrum (B5).The total system size L = 124 and the subsystem size LA = 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.
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.

)θ k 2 −e iφ k cos θ k 2 (θ k 2 e iφ k sin θ k 2 ,
That is, the spectrum bears the reflection symmetry.The eigenstates are SU (2) spinors.The eigenstate corresponding to E k is (B8)where cos θ k = M/E k , tan φ k = B y (k)/B x (k) = − J−J J+J tan( k 2 ).It is obvious that . The time-dependent form factors in Eqs.(B10) and (B11) are given byũfk (t) = a k u − f k e iE f k t + b k u + f k e −iE f k t (B12)andṽfk (t) = a k v − f k e iE f k t + b k v + f k e −iE f k t .(B13)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

± f k
are real, whereas a k , b k , and v ± f k are generally complex except for k = 0 when they are also real.Furthermore,u ± f k , v ± f k are even in k, δ(k) is odd in k, and a * k = a −k , b * k = b −k .

TABLE I .
Comparisons of different fluctuations.