Temporal fluctuations of correlators in integrable and chaotic quantum systems

We provide bounds on temporal fluctuations around the infinite-time average of out-of-time-ordered and time-ordered correlators of many-body quantum systems without energy gap degeneracies. For physical initial states, our bounds predict the exponential decay of the temporal fluctuations as a function of the system size. We numerically verify this prediction for chaotic and interacting integrable spin-1/2 chains, which satisfy the assumption of our bounds. On the other hand, we show analytically and numerically that for the XX model, which is a noninteracting system with gap degeneracies, the temporal fluctuations decay polynomially with system size for operators that are local in the fermion representation and decrease exponentially in the system size for non-local operators. Our results demonstrate that the decay of the temporal fluctuations of correlators cannot be used as a reliable metric of chaos or lack thereof.


I. INTRODUCTION
The spreading of local observables under unitary time evolution has been used as a measure of information scrambling in quantum systems out-of-equilibrium and has spurred a lot of interest across various areas of physics.A central quantity used to assess this spreading is the out-of-time ordered commutator (OTOC) between two operators, Ŵ (t) and V (0), where one is fixed at time t = 0 and the other evolves in time, that is, where [.] indicates the real part of the out-of-timeordered correlation function Ŵ † (t)V † (0)W (t)V (0) .Initially, the commutator is small or zero, and the spreading of the operator in time is manifested in the growth of the commutator.
For chaotic systems with a well-defined semiclassical limit, the initial growth of the OTOC with time is exponential, with a rate determined by the positive Lyapunov exponents of the corresponding classical system [1].As such, OTOCs are natural candidates to explore the chaotic features of a given system.However, this initial exponential growth happens also in integrable models due to instability [49][50][51][52][53][54], as experimentally confirmed in [67].There have also been examples of interacting-integrable systems without a well-defined semiclassical limit, where the OTOC exhibits a diffusive front broadening as in nonintegrable models [16,17], and of chaotic models with local conserved quantities, where the OTOC grows algebraically [11,15,[68][69][70].Recently, it was also shown that for a class of many-body local circuits, exponential OTOC decay is not a good indicator of chaos [71].
While the short-time behavior of the OTOC does not categorically distinguish between integrable and chaotic quantum systems, one may wonder whether its long-time behavior could.In [36], the authors used the size of the temporal fluctuations after the saturation of the OTOC as a way to differentiate between chaos and integrability.
Motivated by [36] and various other studies on the temporal fluctuations of observables [72][73][74][75][76][77][78][79] and OTOCs [40,41], we investigate how the magnitude of the temporal fluctuations of time-ordered correlation functions and out-of-time-ordered correlation functions depend on the system size L. We obtain analytical bounds that show that for systems without energy gap degeneracies (chaotic or not), the fluctuations decay at least exponentially with the system size.We confirm this result numerically by considering three spin-1/2 models in onedimensional (1D) lattices: a chaotic model with firstand second-neighboring couplings, the integrable interacting XXZ model, and the integrable noninteracting XX model.In the first two cases, energy gap degeneracies are absent, so the scaling of the fluctuations with L cannot set them apart and the fluctuations decay exponentially with L. For the XX model, where energy gap degeneracies are present, the decay can be polynomial or exponential depending on the local operators used in the correlators.
The paper is structured as follows.In Sec.II, we introduce the correlators that we study and the general bounds for the decay of their temporal fluctuations as a function of system size.In Sec.III, we present the models and initial states that we consider.In Sec.IV, we study numerically the decay of the fluctuations with L and verify that our bounds are tight.In Sec.V, we present analytical results for the XX model, for which the general bounds do not apply.In Sec.VI, we summarize our main results and outline possible research directions.
Derivations and supporting material are provided in the appendices.

II. CORRELATORS AND BOUNDS ON THEIR TEMPORAL FLUCTUATIONS
We consider the time-ordered correlation function, and the out-of-time-ordered correlation function, where |Ψ 0 is an initial state, Â is a local operator, which in our case is Hermitian and unitary, Â (t) = e i Ĥt Âe −i Ĥt , and Ĥ is the Hamiltonian of the system.We investigate the infinite-time average of the correlators, and the magnitude of their temporal fluctuations around this value,

A. Bounds on Fluctuations
To obtain general bounds on the temporal fluctuations of F A 2,4 , we generalize the results of Refs.[74][75][76] for the fluctuations in an observable Â (t) .Our results for F A 2 are complementary to those in Ref. [41], where a bound on the fluctuations of F A 2 was obtained for thermal initial states and for systems exhibiting weak ETH, as well as to those in Refs.[80,81], where the fluctuations of ktime-ordered correlation functions were bounded on average (see also [82] for studies on temporal fluctuations of nonequilibrium currents).
The Hamiltonian associated with the evolution of the correlators has eigenvalues E n and eigenstates |E n .It is a Hermitian operator that can be written as Ĥ = n E n Pn , where Pn = Kn q=1 |E nq E nq | is a projector onto the degenerate subspace with K n equal eigenvalues E n , and all E n 's in the sum for Ĥ are distinct by construction.

Fluctuations of the time-ordered correlation function
Using the projectors, the time-ordered correlation function in Eq. (2) becomes Since all E n 's are distinct, the infinite-time average is and the fluctuations around F A 2 are given by ∆ 2 We now assume that the energy gaps are nondegenerate, which means that if the gaps E n − E m = E k − E l for any given n, m, k, l, then either n = m and k = l or n = k and k = l.Since n = m and k = l are excluded from the sums in Eq. ( 9), we obtain the following expression for the infinite-time averaged fluctuations in Eq. ( 5), where we defined, Using the Cauchy-Schwarz inequality, Eq.( 10) gets bounded as (see Appendix A), where A is the matrix norm corresponding to the largest eigenvalue of Â.The quantity tr(ω 2 ) corresponds to the inverse participation ratio (IPR) of the initial state in the basis of the eigenstates of the Hamiltonian, where C (0) nq = E nq |Ψ 0 .The result in Eq. (13) implies that for Hamiltonians with non-degenerate energy gaps and initial conditions that have weight on exponentially many eigenstates of the Hamiltonian (which may happen when the Hamiltonian describes a many-body quantum system), the fluctuations decay exponentially with the system size.It is important to stress that, similarly to Refs.[74,76], the obtained bound does not require the system to be chaotic or the spectrum to be non-degenerate.

Fluctuations of the out-of-time-ordered correlation function
To bound the fluctuations of F A 4 (t), we first introduce the following simplified notation, The infinite-time average of F A 4 (t) is given by and the fluctuations around the infinite-time average are ∆ 2 where δ(x) is a Kronecker delta and the prime indicates that all terms with S nmkl , S n m k l = 0 are not included in the sums.The Kronecker delta implies that S nmkl = S n m k l = 0. Similarly to the fluctuations of F A 2 (t), we assume that the nonzero S nmkl 's are unique up to trivial permutations n ←→ k and m ←→ l that leave S nmkl invariant.In other words, for nonzero S nmkl , only the following S nmkl are equal, Using this assumption, we can reduce the constraint into four sums, ∆ 2 where the prime over the sum includes all constraints on (nmkl) that ensure no double counting between the permutations.
To obtain the bound for ∆ 2 , we show (see Appendix A) that the first term on the right-hand-side of Eq. ( 18) is dominating, and write it in the form n,m,k,l where ω is defined in Eq. ( 11), and (20) As shown in Appendix A, we can bound the fluctuations as This is similar to the result in Eq. ( 13) and implies again that if IPR 0 is proportional to the dimension of the Hilbert space of a many-body quantum system, then the fluctuations decay exponentially with the system size.

III. MODELS AND INITIAL STATES
In this section, we describe the models and initial states that we use to numerically confirm in Sec.IV that the bounds derived above are tight.We consider three spin-1/2 chains as described next.
(i) The XX model, where Ŝx,y,z i = σ x,y,z i /2 are spin-1/2 operators operating on site i and σ x,y,z i are Pauli matrices, L is the size of the system, J is the coupling strength, and h b is a border defect, which we introduce to break parity and spinreversal symmetries [83], without breaking the integrability of the model.Throughout the work, we set J = 1 and h b = 0.1.This model can be exactly mapped to noninteracting fermions using the Jordan-Wigner transformation.
(ii) The XXZ model, which is integrable via the Bethe ansatz [84].We set the anisotropy parameter to J z = 0.48 to avoid additional symmetries that appear for J = J z and at the roots of unit, such as J z = J/2 [85].
(iii) The XXZ model in Eq. ( 23) with additional nextnearest neighbor (NNN) couplings, We choose λ = 1, which guarantees that the system is chaotic.
We use open-boundary conditions for the three models to break translational symmetry.All of them conserve the total magnetization in the z-direction, Ŝz tot = i Ŝz i .We work within the Ŝz tot = 0 subspace, which has the Hilbert space dimension D = L L/2 .Interacting Hamiltonians without symmetries, like those defined in Eqs. ( 23) and ( 24), typically satisfy the condition of non-degenerate energy gaps, which is needed for the bounds in Eq. ( 13) and Eq. ( 21).This contrasts with noninteracting models, or systems that can be mapped to free fermions, such as the XX model in Eq. ( 22), which have a large number of gap degeneracies even when the spectrum is non-degenerate.

A. Initial States
In this work, we use three initial states: a random state drawn from the Haar measure, which we refer to as the

IV. NUMERICAL RESULTS
In this section, we show that when the bound that we derived in Sec.II applies, it is tight.We numerically investigate the evolution and temporal fluctuations of the time-ordered and out-of-time-ordered correlation functions for The corresponding correlators are denoted by F X 2,4 and F Z 2,4 , respectively.We consider the three models presented above and take a random state from the Haar measure as our initial state.The results do not change qualitatively for the other initial states listed in Sec.III A, but they differ in details (see Appendix B).
Since we are interested in the long-time dynamics, we use exact diagonalization, which limits the accessible system's size to L ≤ 18.

A. Results
Before analyzing the temporal fluctuations, we present in Fig. 1, the evolution of F Z 2,4 (t) [(a)-(c)] and F X 2,4 (t) [(d)-(f)] from time t = 0 up to their saturation; for the XX, XXZ, and NNN models; for various system sizes.With the exception of F Z 2,4 (t) for the XX model [Fig.1(a)], the correlation functions saturate to a very small value and exhibit small temporal fluctuations that decrease as the system size increases.The fluctuations of F X 2,4 (t) are noticeably smaller than those of F Z 2,4 (t).This is likely caused by the conservation of the total z-magnetization in the studied models, which is also known to slow down the decay in time of F Z 4 (t) [70].The behavior of F Z 4 (t) after saturation for the XX model [Fig.1(a)] stands out.While for the other models, F Z 4 (t) fluctuates around zero at long times, for the XX model, it approaches 1, as discussed also in Ref. [26].In contrast, F Z 2 (t) for the XX model [Fig.1(a)] does decay to zero, even though it exhibits larger fluctuations than for the XXZ [Fig.1(b)] and NNN [Fig.1(c)] models.In Sec.V, we elucidate the behavior of F Z 2,4 (t) for the XX .The models are (a) the XX defined in Eq. ( 22), (b) the XXZ from Eq. ( 23), and (c) the NNN in Eq. (24).
model by analytical arguments.
The saturation of F X 4 (t) for the XX model [Fig.1(d)] is preceded by quasi-periodic oscillations, which are absent for the chaotic NNN model [Fig.1(f)].Smaller oscillations at intermediate times are also visible for the XXZ model, although they happen at a plateau [Fig.1(e)].Interestingly, a plateau that gets longer as L increases also appears for F Z 4 (t) in the XXZ model [Fig.1(b)].Contrary to F X 4 (t), the time-ordered correlation function F X 2 (t) decays fast to a small value without oscillations for all three models.
With the general picture of the behavior of F X,Z 2,4 (t) provided in Fig. 1, we study the dependence of the temporal fluctuations, ∆ F X,Z 2,4 as a function of system size.Fig. 2 shows this dependence for the Haar state (for a comparison with the results for the Néel and domain wall states, see Appendix B), and confirms that the bound derived in Sec.II is tight.Namely, that ∆ F X,Z 2,4 decays exponentially with system size for the interacting-integrable [Fig.2(b)] and chaotic [Fig.2(c)] systems.For a given chain length, even the magnitude of the fluctuations is comparable for both models.
The behavior for the XX model [Fig. 2 decreases slower than exponentially.As stated in Sec.II, this model has gap degeneracies, so the proof of Sec.II does not apply to it.Nevertheless, since it can be mapped to free fermions, in the next section, we provide numerical and analytical results for the dependence of ∆ 2 on L, as well as for the infinite-time averages F Z 2,4 .

V. XX MODEL
Since the XX model maps exactly to noninteracting fermions, its infinite-time average, F Z 2,4 , and temporal fluctuations, ∆ 2 , can be computed analytically.The numerical analysis of the previous section can also be extended to much larger system sizes, because the complexity of the calculations increases only polynomially with system size.
We compute F Z 2 (t) and F Z 4 (t) by writing them in terms of fermionic creation, ĉ † L/2 , and annihilation, ĉL/2 , operators acting on site L/2, that is, For quadratic Hamiltonians, the time-evolution of the creation and annihilation operators is given by where u ik (t) = i e −ihst k is the single-particle propagator and h s the single-particle Hamiltonian.
After expressing F Z 2,4 (t) in terms of fermionic operators, we calculate the correlations using Wick's theorem (see Appendix C).For generic initial states, the final expression is cumbersome, however, it considerably simplifies at infinite temperature, yielding and We can then obtain the infinite time-average, where |α and ε α are the eigenstates and eigenvalues of h s .In the last equality, we assumed that the singleparticle spectrum is non-degenerate, which is true for the XX model with the border impurity.We see that F Z 2 is the IPR of the initial state corresponding to a single particle at the center of the chain, L/2, and projected in the  , as a function of system size; for the XX model and an initial state at infinite temperature.Solid lines represent the numerical data and dashed lines, the analytical results.
basis of the single-particle eigenstates of h s .For delocalized single-particle initial states, we have that F Z 2 ∝ L −1 .Using F Z 2 , we write the temporal fluctuations of F Z 2 (t) as and obtain the variance A similar derivation follows for the temporal fluctuations of F Z 4 (t).The infinite-time average of the out-oftime ordered correlation function is given by a result that was first discussed in Ref. [26].Since F Z 4 (t) is dominated by the second term in the parenthesis of Eq. ( 29), its temporal fluctuations decay with system size similarly to what happens for F Z 2 (t) [cf.Eq. ( 28)].
To confirm these analytical estimates numerically, we compute F Z 2 (t) and F Z 4 (t) using Eq. ( 28) and Eq. ( 29) for a number of system sizes.Figure 3(a) shows that the infinite-time average decays as L −1 , and Fig. 3(b) shows that the temporal fluctuations around that average decay as L −2 , in perfect alignment with the analytical estimates.
The calculation of the temporal fluctuations of F X 2,4 (t) is considerably harder, because both correlators expressed in terms of fermionic operators are non-local, so one has to perform Wick's contraction of the order of L operators.While this can be done numerically, we were not able to obtain analytical estimates, which would help to explain the exponential decay of the fluctuations with the system size.

VI. DISCUSSION
We rigorously showed, that for any quantum system with non-degenerate energy gaps, the temporal fluctuations around the saturation value of time-ordered and out-of-time-ordered correlation functions are bounded by the square root of the inverse participation ratio of the initial state.Since most physical initial states are composed of exponentially (in the system size L) many eigenstates of the Hamiltonian, this implies that for such initial states, the temporal fluctuations decay at least exponentially with L. We verified numerically that the bounds are tight; the fluctuations decay exponentially for the interacting integrable XXZ spin-1/2 chain and for its chaotic version with next-nearest-neighbor couplings.Our results demonstrate that the decay of the temporal fluctuations of correlators as a function of the system size cannot serve as a reliable metric of chaoticity.This has to be contrasted with the decay of temporal fluctuations when the initial state is one eigenstate of the Hamiltonian.For such initial states, the fluctuations are related to the off-diagonal matrix elements, which do show different decay with system size for chaotic and integrable systems [86].
The only distinguishing behavior that we identified was for the 1D XX model, which is exactly mappable to noninteracting fermions.Since this model has degenerate energy gaps, the bounds that we obtained on the decay of the fluctuations do not apply; however, they can be calculated both analytically and numerically.We find that for this system, the temporal fluctuations of the timeordered and out-of-time-ordered correlation functions in the z-direction, F Z 2,4 (t), decay as L −2 .On the other hand, we observed numerically that the temporal fluctuations of F X 2,4 (t) decay exponentially with system size.We argue that, in this case, the exponential decay stems from the non-locality of F X 2,4 (t), when written in terms of fermionic annihilation and creation operators.It remains to put this claim on more solid ground.It would also be interesting to investigate if there are scenarios in which different power-law decays of the fluctuations are possible.

ACKNOWLEDGMENTS
This research was supported by a grant from the United States-Israel Binational Foundation (BSF, Grant No. 2019644), Jerusalem, Israel, and the United States National Science Foundation (NSF, Grant No. DMR-1936006).LFS thanks Vinitha Balachandran, Marcos Rigol, and Dario Poletti for valuable discussions.
Appendix A: Proofs of the general bounds In this section, we provide detailed proofs of the bounds in Eqs. ( 13) and ( 21) of the main text.We showed in Eq. ( 10) that ∆ 2 where ω and ω A are defined in Eq. (11).Using the Cauchy-Schwarz inequality, setting V = Âω A and W = Â † ω, and using the cyclic property of the trace, we further bound Eq. ( 10) as follows, where we used that for two positive matrices A and B, tr(AB) ≤ A tr(B).We can bound tr(ω 2 A ) by This gives For the fluctuations of F A 4 , we start with Eq. ( 18).We designate a given permutation of (nmkl) by σ (nmkl), and notice that it is bijective.Then, using the Cauchy-Schwarz inequality, we have nmkl where in the last inequality we renamed the indexes of the second term in the square root.We also removed the constraints on the sum, using the fact that all elements of the sum are now positive.This allows us to write Using the definitions of ω in Eq. ( 11) and ω AAA in Eq. ( 20), which are positive operators, we have n,m,k,l Using the Cauchy-Schwarz inequality combined with the cyclic property of the trace, We now proceed by bounding tr(ω 2 AAA ), which is given by Since Â |Ψ 0 Ψ 0 | Â † and the matrix to the left of it are positive, we can write and eliminate the sum over l using the cyclic property of the trace and the fact that m P m = I.Therefore, We can now proceed along the same lines as before until we get Combining the expressions, we obtain Appendix B: Dependence on the initial state Here, we compare the results for the variance of the temporal fluctuations obtained for the Haar state as the initial state with those for the Néel and the domain-wall states.The expectation of the Hamiltonian, E 0 = Ψ 0 | H |Ψ 0 of these two states depends on the model according to [87], Néel state: Domain wall state: For the XX model, where J z , λ = 0, the energies of all three states are in the middle of the spectrum, E 0 = 0.For the XXZ model with large L and for our choices of parameters, |E 0 | ∼ JL/8, where E 0 is negative (positive) for the Néel (domain wall) state.For the NNN model, the energy of the Néel state is close to the middle of the spectrum, E 0 ∼ 0, where chaos is strong and the eigenstates are very delocalized, while the energy of the domain wall state for large L is closer to the positive edge of the spectrum, E 0 ∼ JL/4.In Fig. 4, we plot ∆ 2  In the following, we present the detailed derivation of F Z 4 (t) defined in Eq. ( 29) of the main text, where Ŵ (t) = σ z i (t) and V = σ z j .We rewrite Ŵ and V in terms of fermionic operators as where ni (t) = ĉ † i (t) ĉi (t) and nj = ĉ † j ĉj , so that Expanding the previous expression, some terms trivially cancel and it gets reduced to where we used n2 i (t) = ni (t) and n2 j = nj .In the following, we expand the expectation value and work on each of the terms above.But before doing that, let us express the time evolution of ĉ † i (t) and ĉi (t) as where u ik (t) = i e −ihst k is the single-particle propagator and h s the single-particle Hamiltonian.Using Eq.(C5), we then express the first term of Eq.(C4) as It is convenient to rewrite all the expressions above using the Green's function While we can obtain a general equation for any initial state for which the Wick's theorem holds, it is rather cumbersome.In particular, for the infinite-temperature state considered in this work, it considerably simplifies, since nj = ni (t) = 1/2, and the remaining terms can be readily written in terms of Green where we used Eq.(C5).Finally, substituting Eq. (C15) in Eq. (C14), we obtain the expression for F Z 4 (t) at infinite temperature, which we used in the main text, where u ij (t) is the single-particle propagator.In the text, we set i = j = L/2.

Time-ordered correlation function
From the previous subsection, one can check that in terms of fermionic operators, the time-ordered correlation function F Z 2 (t) defined in Eq. ( 28) of the main text is given by

FIG. 1 .
FIG. 1. Dynamics of the correlators F Z 2,4 (t) in the top row (a)-(c) and F X 2,4 (t) in bottom row (d)-(f) evaluated numerically for the Haar state initial condition and three spin-1/2 models.Solid lines correspond to F X,Z 4 (t) and dashed lines to F X,Z 2 (t).For F Z 4 (t), a number of system sizes are considered, where the larger sizes are indicated with darker colors [legend in (b) is for (a)-(c) and legend in (e) for (d)-(f)].For F Z 2 (t), the system size is L = 18 and for F X 2 (t), it is L = 16.The XX model from Eq. (22) is depicted in the left column [(a), (d)], the XXZ model from Eq. (23) in the middle column [(b), (e)], and the chaotic NNN model from Eq. (24) in the right column [(c), (f)].

FIG. 2 . 4 ( 2 ( 2 F
FIG. 2. Variance of the temporal fluctuations of F X,Z 4 (t) (circles) and F X,Z 2 (t) (diamonds) for the Haar state initial condition as a function of system size, L. Filled symbols represent ∆ 2 F Z 2,4 and empty symbols correspond to ∆ 2 F X 2,4

, 4 for 2 F X,Z 4 for 2 F Z 4 2 F X 4 decays
the three different initial states.Their results are qualitatively similar, but there are subtle differences associated with the position of the energy of the initial state in the spectrum.While for the Haar state [Fig.4(a),(d)], there is practically no difference in the values of ∆ the interacting-integrable and the chaotic model, the fluctuations for the Néel state [Fig.4(b),(e)] are smaller for the chaotic model, since in this case, |Ψ 0 is in the middle of the spectrum, being thus more delocalized than for the XXZ model.A similar explanation can be given for the smaller fluctuations associated with the Néel state for the NNN model [Fig.4(b),(e)] when compared with those for the domain wall state under the same model [Fig.4(c),(f)].With respect to the XX model, the magnitude of the fluctuations is analogous for all three initial states.The difference lies in the operator considered, as already discussed in Fig.2(a), that is, ∆ decreases slower than exponentially [Fig.4(a)-(c)] and ∆ exponentially with L [Fig.4(d)-(f)].