# Operator growth and eigenstate entanglement in an interacting integrable Floquet system

###### Abstract

We analyze a simple model of quantum dynamics, which is a discrete-time deterministic version of the Frederickson-Andersen model. We argue that this model is integrable, with a quasiparticle description related to the classical hard-rod gas. Despite the integrability of the model, commutators of physical operators grow as in generic chaotic models, with a diffusively broadening front, and local operators obey the eigenstate thermalization hypothesis (ETH). However, large subsystems violate ETH; as a function of subsystem size, eigenstate entanglement first increases linearly and then saturates at a scale that is parametrically smaller than half the system size.

A central theme in many-body physics is the mechanics of thermalization, decoherence, and information scrambling in isolated, interacting quantum systems Cazalilla and Rigol (2010); Polkovnikov et al. (2011); Basko et al. (2006); Nandkishore and Huse (2015); Langen et al. (2013); Kaufman et al. (2016). Generic systems (except for the many-body localized phase Basko et al. (2006); Nandkishore and Huse (2015)) are believed to obey the eigenstate thermalization hypothesis (ETH), which posits that local observables in a many-body eigenstate behave as they would in an appropriately chosen thermal state Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Kim et al. (2014); Dymarsky and Liu (2017). However, in one dimension, a number of physically important systems (e.g., the Heisenberg spin chain) are integrable—i.e., they have extensively many conservation laws and thus do not thermalize Rigol et al. (2008); Ilievski et al. (2016). The absence of thermalization in integrable systems, and the anomalously slow “prethermal” behavior of nearly integrable systems, have been experimentally observed Kinoshita et al. (2006); Gring et al. (2012); Langen et al. (2013); Tang et al. (2018); Erne et al. (2018); Li et al. (2018). Some integrable systems can be solved by mapping to free fermions Schultz et al. (1964), but many others, including the Heisenberg chain, cannot Faddeev (2016). In the latter class of “interacting” (i.e., Bethe-ansatz-solvable) integrable systems, it is challenging to compute the dynamics of physical observables, because physical operators have complicated representations in terms of the quasiparticles Caux and Calabrese (2006), although coarse-grained approaches to some dynamical questions have recently been developed Castro-Alvaredo et al. (2016); Bertini et al. (2016); Bulchandani et al. (2018); Doyon et al. (2018a). In the quantum context, most work on integrable systems has focused on Hamiltonian dynamics, though there has been some recent work on time-periodic, driven integrable systems Gritsev and Polkovnikov (2017); Vanicat et al. (2017).

This work presents and analyzes a simple integrable Floquet model for which many of these questions can be explicitly addressed. This model is a deterministic discrete-time version of the Frederickson-Andersen model Ritort and Sollich (2003); Garrahan (2017); Hickey et al. (2016), which is a standard model of kinetically constrained dynamics; we call it the Floquet-Frederickson-Andersen (FFA) model. The dynamics of the FFA model is in some ways analogous to the classical hard-rod gas Spohn (2012); Doyon et al. (2018a), a canonical interacting integrable system; as we discuss, there is a natural description in terms of ballistically propagating quasiparticles. Beyond its integrability, what renders the model tractable is that its dynamics maps each computational-basis product state to a unique computational-basis product state; this allows for efficient classical simulations of dynamics. Despite its simplicity the FFA model retains two key features of generic integrable systems: first, the relation between physical observables and quasiparticles is nontrivial, so physical observables are complicated in the quasiparticle language; and second, each quasiparticle’s motion (and even the geometry it feels Doyon et al. (2018b)) is modified by the distribution of other quasiparticles.

Our main results are as follows (Fig. 1). For physical observables the OTOC displays scrambling, with a front that broadens diffusively as expected for generic chaotic systems Xu and Swingle (2018a, b); Khemani et al. (2018); however, this front is anomalous on timescales exceeding the system size . In addition, small subsystems satisfy ETH, with a fraction of outlier states that appears to vanish (slowly) in the thermodynamic limit. However, large subsystems are strikingly nonthermal: for subsystem size , the eigenstate entanglement crosses over from a thermal volume law to a constant. This crossover sharpens for larger systems.

Model.—The model we consider was recently introduced in Ref. Gopalakrishnan and Zakirov (2018); the system is a spin- chain, subject to repeated application of the unitary

(1) |

where consists of the following rule, applied to each odd spin : apply the Pauli operator unless the two neighboring even sites, and , are both in the state. This rule can be composed from the gate sequence , in which controlled NOT and Toffoli gates are applied to the target site from its two neighbors. repeats this process with even and odd sites exchanged. The full-cycle unitary has a strict causal light cone that expands at two sites per period (dashed red lines in Fig. 1). In what follows we treat periodic boundary conditions; the effects of other boundary conditions are deferred to future work. We also refer to odd (even) sites as the () sublattices respectively.

Eq. (1) can be regarded as a reversible block cellular automaton Schumacher and Werner (2004); Lesanovsky et al. (2018). A number of works have addressed cellular automata and quantum circuits involving Clifford gates Gütschow et al. (2010); Werner et al. (2010); Chandran and Laumann (2015); Nahum et al. (2017); Gopalakrishnan and Zakirov (2018). Clifford gates induce simple operator dynamics, mapping Pauli strings to other Pauli strings, and therefore do not give rise to operator entanglement Hosur et al. (2016); Jonay et al. (2018). A related but Clifford-only model, without the Toffoli gate, was analyzed in Ref. Gopalakrishnan and Zakirov (2018) and shown to map to free particles. In contrast, the Toffoli gate in Eq. (1) induces operator entanglement for local operators, even after a single step of time evolution: for example, the operator evolves to , which cannot be factored into on-site operators. Our methods here cannot be used to compute operator entanglement, however, so we defer a full discussion to future work.

Quasiparticle picture.—We first present a simple quasiparticle picture Yang and Yang (1969) of the dynamics of this model; this picture is not derived microscopically, but is inferred from simulating the dynamics, and justified by its explanation of various exact numerical results. The structure of quasiparticles is easiest to describe when the density of up spins is low. In this limit, the elementary quasiparticles are pairs of two adjacent up spins surrounded by down spins. There are two inequivalent quasiparticles, respectively right- and left-moving, based on whether their sublattice structure is or . (The dynamics is symmetric under exchanging sublattices and time-translating by half a period, but not under each action separately.) All right-movers and all left-movers have the same speed, so collisions necessarily involve a right-mover and a left-mover; also, each three-body collision can only occur in one sequence, precluding diffractive processes. Each collision induces a time delay of one time step, which is the same for all collisions [Fig. 2(a)]; thus, this model resembles hard rods with negative rod length. The quantization condition for right (left) movers simply depends on the total number of left (right) movers. Since all quasiparticles have the same velocity, a state is characterized by the spacings between adjacent left-movers and between adjacent right-movers. Adjacent right-movers must have at least one empty site between them, since the configuration in Fig. 2(b) is not just a pair of right-movers; thus the state space for right-movers is constrained, with the same constraints as in Ref. Chandran et al. (2016).

At finite density [Fig. 2(c)], the quasiparticle velocity is decreased by an amount proportional to the density of other quasiparticles, which in turn is (approximately) proportional to the density of occupied sites; however, the dynamics still consists of ballistically moving quasiparticles.
At high density, the model remains integrable in terms of these quasiparticles. However, one can identify the quasiparticle content of a product state by recasting the model in terms of bonds, as follows ^{1}^{1}1D.A. Huse, private communication: assign a quasiparticle to each bond where both spins are up, and assign two quasiparticles to any spin configuration that has the sequence . Note that the number of physical up spins fluctuates, though the number of quasiparticles remains conserved, because quasiparticles transiently “merge” during a collision [Fig. 2(a-b)]. The quasiparticle structure is explored further in sup (), by simulating the “free expansion” Rigol and Muramatsu (2005) of a general initial state.

OTOCs.—The clearest evidence that the model remains integrable at high densities comes from studying OTOCs. In a classical system, the OTOC corresponds to the local overlap between two histories with identical initial conditions except for a disturbance at the origin Das et al. (2017). To establish integrability we perform the following numerical experiment [Fig. 2(c-d)]: we create an initial state with a region where the density of up spins in part of the system is low (so we can reliably create a single quasiparticle there) but the rest of the system is at high density. We then move this quasiparticle, and overlap histories with different initial positions of the quasiparticle. We find that translating a single quasiparticle does not have a “butterfly” effect: the OTOC is simply the time trace of the quasiparticle that was moved [Fig. 2(d)]. This is direct evidence that the model remains integrable at high densities: in a chaotic system, “firing” a quasiparticle into the system at time time rather than at would cause a butterfly effect, which is clearly absent here.

Although moving quasiparticles does not cause a spreading disturbance, adding or removing quasiparticles does, as the presence of a new quasiparticle modifies the phase shifts of all the others. Physical spin operators typically create and/or move quasiparticles, depending on the underlying product state. In contrast with, e.g., the Ising model, all simple operators have an amplitude for both creating and translating quasiparticles. For instance, is a state with four quasiparticles, while is a state with only two quasiparticles. Therefore all simple operators spread.

For concreteness we focus on the simplest OTOC, . (We have also checked the OTOCs of more complicated operators sup () but they do not behave appreciably differently, in contrast with the Ising case Lin and Motrunich (2018).) The OTOC averaged over many initial product states grows with a light-cone typical of chaotic systems (Fig. 1). The magnitude of the OTOC near its growth “front” precisely matches the recent prediction for chaotic systems Xu and Swingle (2018a, b); Khemani et al. (2018), behaving to leading order as , where [Fig. 3]. The velocity of the OTOC front is half the causal light-cone speed; this is expected, since a random state is at half-filling. The FFA model is thus distinct from other large- or “classical” limits, in which the OTOC front is sharp Nahum et al. (2017); Xu and Swingle (2018a, b) (note that in the FFA model the front is sharp if the initial state is a computational-basis product state but not otherwise). Starting from a random initial distribution of quasiparticles, the diffusive broadening of the front is a natural consequence of the random time delays due to collisions.

The late-time behavior of has some unusual features, which can be seen in the lower two panels of Fig. 3. Even at times before the operator has wrapped around the system, its behavior inside the front is unexpected. Rather than decorrelating completely, the two histories remain weakly anticorrelated within the front (so the value of the OTOC inside the front is rather than ). After a timescale , this behavior changes and the two histories compared by the OTOC become weakly correlated, giving rise to the diamond-like shape seen in the space-time plot [Fig. 3, lower left]. The OTOC then saturates at a value until the much longer revival timescale . These saturation values vary depending on the operator. These effects stem from high-density initial configurations sup ().

Our quantitative analysis of the OTOC involved averaging over eigenstates. The OTOC in a single randomly chosen eigenstate also spreads out with a light-cone, but there is less broadening, and the front “refocuses” on the timescale sup (). This can be understood within the quasiparticle picture. In a particular eigenstate, the distribution of left and right moving quasiparticles, and the spacings among the left- and right-movers, are fixed, but one averages over the point in the classical trajectory at which the new quasiparticle is introduced. Thus, the sequence of time delays experienced by (say) a right-mover is randomized, but the total time delay is fixed by the total number of left-movers. Nevertheless, we expect the broadening of the front to be Gaussian at times , since in this temporal regime the average runs over randomly timed collisions.

Eigenstates.—We now turn to the eigenstates of the FFA model. Since takes each product state to a product state, the dynamics of an initial product state consists of chains of transitions . A random eigenstate can therefore be constructed Gopalakrishnan and Zakirov (2018) by picking a random initial product state and summing over its orbit, with appropriate phases, i.e.,

(2) |

Such states are evidently eigenstates of so long as for some .

A central quantity in our analysis is the recurrence time for a given initial state, which measures how much of configuration space is accessible under unitary evolution from a given initial state. In an “ergodic” system, essentially every configuration would be visited, and would grow exponentially with system size. This is not what happens in the FFA model (Fig. 4); instead . This follows from the quasiparticle picture: a left-mover traverses the system on a timescale set by the number of right-movers, and vice versa. Since both numbers are of order , their least common multiple is , although there are many configurations with a much smaller least common multiple, and for these is much smaller (Fig. 8). Note that is a recurrence time specific to a particular initial computational-basis product state. The recurrence time of a random initial vector is the least common multiple of each orbit’s recurrence time. By sampling many initial states and computing their , we find that this global recurrence time grows at least as fast as sup (), so it is in general exponentially larger than the Hilbert space dimension. The quasiparticle picture also suggests that the scaling is exponential with system size: according to this picture, is the least common multiple of all possible quasiparticle periods , which scales as the product of all primes ; asymptotically this product grows exponentially in by the prime number theorem. Thus, scales exponentially rather than double-exponentially with system size, in contrast with generic chaotic models Hosur et al. (2016).

We previously showed Gopalakrishnan and Zakirov (2018) that upper-bounds the second Renyi entanglement entropy Nielsen and Chuang (2002) through any bipartite cut; an implication is that for any subsystem size and system size . This is what we find, by directly computing for various subsystem sizes (Fig. 1). As Fig. 4 shows, there are strong state-to-state fluctuations in the saturation value of , but the entanglement of every eigenstate saturates for subsystems well below half the system size. Small subsystems, on the other hand, appear to satisfy ETH: both the Renyi entropy for small subsystems of length and the expectation values of on-site operators are narrowly distributed, with state-to-state spread that narrows with system size (Fig. 4). This narrowing is slower than for generic ETH systems, however. More details about the outlying states are discussed in sup ().

Eigenvalues and level statistics.—From the construction of eigenstates, it follows that each eigenvalue is of the form , where labels inequivalent orbits. There is an -fold degeneracy at quasienergy zero (i.e., eigenvalue unity for the unitary ), and also other large degeneracies at special frequencies that divide many orbit periods. Owing to these degeneracies the level statistics are neither Poisson nor random-matrix.

Off-diagonal matrix elements.—Finally, we consider the off-diagonal matrix elements of local operators (specifically, the two-spin-flip operator ) between eigenstates. For convenience we restrict ourselves to matrix elements between states in the quasi-energy zero sector; as noted above, the spectrum is highly degenerate and each distinct “orbit” contributes one state to this sector. Most pairs of eigenstates have strictly zero matrix element; the fraction of nonzero matrix elements decreases exponentially with system size, approximately as (Fig. 5). This behavior has also been seen in other integrable models ^{2}^{2}2V. Khemani, private communication.. The distribution of the nonzero matrix elements broadens with system size but does not seem to approach a Gaussian at the accessible sizes.

Discussion.—This work analyzed the FFA model, a simple interacting integrable Floquet system. Because this model has a special basis with classical dynamics, various quantities can be computed here that are harder to compute in more realistic interacting integrable systems such as the Heisenberg chain. In its integrable dynamics the FFA model is analogous to a discretized hard-rod gas (except that collisions here delay the propagation of quasiparticles instead of speeding it up). Despite the integrability of the model, physical operators exhibit a chaotic butterfly effect, with a diffusively broadening front that quantitatively matches the predictions of Ref. Xu and Swingle (2018a). We believe this feature should be generic for interacting integrable systems, as it arises from the medium-dependence of the velocities, which is a generic feature. Moreover, small subsystems (relative to the system size) obey the eigenstate thermalization hypothesis. Thus, integrability is “hidden” from these diagnostics. However, large subsystems clearly violate the ETH: the eigenstate entanglement has a sharpening crossover from volume-law growth to saturation, at a subsystem size that grows logarithmically with the full system size. Metrics such as the level statistics also diagnose the non-thermal character of this model.

Many questions remain for future work, including perturbations of the model that restore chaos and/or quantum fluctuations; identifying the local conserved operators and formally demonstrating the integrability of the model, e.g., through methods for integrable cellular automata Bruschi et al. (1992); Tokihiro et al. (1996); Joshi and Lafortune (2005); and exploring operator entanglement.

Note added.—While this manuscript was being prepared, a preprint appeared Rowlands and Lamacraft (2018), mapping operator spreading in noisy quantum circuits to the (stochastic) Frederickson-Andersen model. The results seem conceptually unrelated to ours, however.

Acknowledgments.—S.G. thanks David Huse, Vedika Khemani, Brian Swingle, and Romain Vasseur for helpful discussions and comments on a draft of this paper, and Lincoln Carr, Andrew Potter, and Bahti Zakirov for helpful discussions on related topics. This work was supported by NSF Grant No. DMR-1653271.

## References

- Cazalilla and Rigol (2010) M. Cazalilla and M. Rigol, New Journal of Physics 12, 055006 (2010).
- Polkovnikov et al. (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
- Basko et al. (2006) D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
- Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
- Langen et al. (2013) T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J. Schmiedmayer, Nature Physics 9, 640 (2013).
- Kaufman et al. (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
- Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
- Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
- Rigol et al. (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
- Kim et al. (2014) H. Kim, T. N. Ikeda, and D. A. Huse, Phys. Rev. E 90, 052105 (2014).
- Dymarsky and Liu (2017) A. Dymarsky and H. Liu, arXiv preprint arXiv:1702.07722 (2017).
- Ilievski et al. (2016) E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, Journal of Statistical Mechanics: Theory and Experiment 2016, 064008 (2016).
- Kinoshita et al. (2006) T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
- Gring et al. (2012) M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, Science , 1224953 (2012).
- Tang et al. (2018) Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Physical Review X 8, 021030 (2018).
- Erne et al. (2018) S. Erne, R. Buecker, T. Gasenzer, J. Berges, and J. Schmiedmayer, arXiv preprint arXiv:1805.12310 (2018).
- Li et al. (2018) C. Li, T. Zhou, I. Mazets, H.-P. Stimming, Z. Zhu, Y. Zhai, W. Xiong, X. Zhou, X. Chen, and J. Schmiedmayer, arXiv preprint arXiv:1804.01969 (2018).
- Schultz et al. (1964) T. D. Schultz, D. C. Mattis, and E. H. Lieb, Reviews of Modern Physics 36, 856 (1964).
- Faddeev (2016) L. Faddeev, in Fifty Years of Mathematical Physics: Selected Works of Ludwig Faddeev (World Scientific, 2016) pp. 370–439.
- Caux and Calabrese (2006) J.-S. Caux and P. Calabrese, Phys. Rev. A 74, 031605 (2006).
- Castro-Alvaredo et al. (2016) O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
- Bertini et al. (2016) B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016).
- Bulchandani et al. (2018) V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore, Phys. Rev. B 97, 045407 (2018).
- Doyon et al. (2018a) B. Doyon, T. Yoshimura, and J.-S. Caux, Phys. Rev. Lett. 120, 045301 (2018a).
- Gritsev and Polkovnikov (2017) V. Gritsev and A. Polkovnikov, SciPost Physics 2, 021 (2017).
- Vanicat et al. (2017) M. Vanicat, L. Zadnik, and T. Prosen, arXiv preprint arXiv:1712.00431 (2017).
- Ritort and Sollich (2003) F. Ritort and P. Sollich, Advances in Physics 52, 219 (2003).
- Garrahan (2017) J. P. Garrahan, arXiv preprint arXiv:1709.09208 (2017).
- Hickey et al. (2016) J. M. Hickey, S. Genway, and J. P. Garrahan, Journal of Statistical Mechanics: Theory and Experiment 2016, 054047 (2016).
- Spohn (2012) H. Spohn, Large scale dynamics of interacting particles (Springer Science & Business Media, 2012).
- Doyon et al. (2018b) B. Doyon, H. Spohn, and T. Yoshimura, Nuclear Physics B 926, 570 (2018b).
- Xu and Swingle (2018a) S. Xu and B. Swingle, arXiv preprint arXiv:1802.00801 (2018a).
- Xu and Swingle (2018b) S. Xu and B. Swingle, arXiv preprint arXiv:1805.05376 (2018b).
- Khemani et al. (2018) V. Khemani, D. A. Huse, and A. Nahum, arXiv preprint arXiv:1803.05902 (2018).
- Gopalakrishnan and Zakirov (2018) S. Gopalakrishnan and B. Zakirov, arXiv preprint arXiv:1802.07729 (2018).
- Schumacher and Werner (2004) B. Schumacher and R. F. Werner, arXiv preprint quant-ph/0405174 (2004).
- Lesanovsky et al. (2018) I. Lesanovsky, K. Macieszczak, and J. P. Garrahan, arXiv preprint arXiv:1804.09794 (2018).
- Gütschow et al. (2010) J. Gütschow, S. Uphoff, R. F. Werner, and Z. Zimborás, Journal of Mathematical Physics 51, 015203 (2010).
- Werner et al. (2010) R. F. Werner, V. Nesme, and J. Gütschow, Discrete Mathematics & Theoretical Computer Science (2010).
- Chandran and Laumann (2015) A. Chandran and C. R. Laumann, Phys. Rev. B 92, 024301 (2015).
- Nahum et al. (2017) A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys. Rev. X 7, 031016 (2017).
- Hosur et al. (2016) P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, Journal of High Energy Physics 2016, 4 (2016).
- Jonay et al. (2018) C. Jonay, D. A. Huse, and A. Nahum, arXiv preprint arXiv:1803.00089 (2018).
- Yang and Yang (1969) C.-N. Yang and C. P. Yang, Journal of Mathematical Physics 10, 1115 (1969).
- Chandran et al. (2016) A. Chandran, M. D. Schulz, and F. J. Burnell, Phys. Rev. B 94, 235122 (2016).
- (46) D.A. Huse, private communication.
- (47) See online supplemental material. .
- Rigol and Muramatsu (2005) M. Rigol and A. Muramatsu, Modern Physics Letters B 19, 861 (2005).
- Das et al. (2017) A. Das, S. Chakrabarty, A. Dhar, A. Kundu, R. Moessner, S. S. Ray, and S. Bhattacharjee, arXiv preprint arXiv:1711.07505 (2017).
- Lin and Motrunich (2018) C.-J. Lin and O. I. Motrunich, Phys. Rev. B 97, 144304 (2018).
- Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, Quantum computation and quantum information (AAPT, 2002).
- (52) V. Khemani, private communication.
- Bruschi et al. (1992) M. Bruschi, P. Santini, and O. Ragnisco, Physics Letters A 169, 151 (1992).
- Tokihiro et al. (1996) T. Tokihiro, D. Takahashi, J. Matsukidaira, and J. Satsuma, Phys. Rev. Lett. 76, 3247 (1996).
- Joshi and Lafortune (2005) N. Joshi and S. Lafortune, Journal of Physics A: Mathematical and General 38, L499 (2005).
- Rowlands and Lamacraft (2018) D. Rowlands and A. Lamacraft, arXiv:arXiv:1806.01723 (2018).

## Supplemental Material

In what follows, we discuss the quasiparticle structure of the FFA model at general densities, quantify states that are outliers in their entanglement properties, estimate of the global recurrence time, and present numerical results on the OTOCs of more complicated operators as well as OTOCs measured within single product states or single eigenstates.

Structure and counting of quasiparticles.—To extract the quasiparticle content of a general initial state, we perform the following numerical experiment. We create an initial state with no up spins outside of a region of size , and observe its time dynamics. The initial state is composed of ballistic quasiparticles, and under time evolution these fly apart; at timescales , the left- and right-movers will have completely separated, allowing us to resolve the quasiparticle content of the state. Fig. 6 shows this expansion dynamics for a random initial state, as well as a very high-density initial state in which all spins are up.

These results suggest the following picture for the quasiparticle content of a state, and for entropy vs. number of quasiparticles. (For specificity we consider right-movers; all statements in the rest of this paragraph hold under exchanging right- and left-movers.) First, in general, quasiparticle states can be specified by the asymptotic spacings between adjacent right-movers. Second, two right-movers cannot occupy adjacent states; there is a hardcore constraint, as noted in the main text. Third, the number of right-movers that can fit in a system of size is a function of the number of left-movers, because left-movers effectively “shrink” the right movers and allow them to fit closer together. In a system of size , the number of available right-moving states is where is the number of left-movers. If we guess that the count of left-movers is , the average number of right-movers in the state is (using results from the Fibonacci chain Chandran et al. (2016)), so this guess is indeed self-consistent.

OTOC at very high density.—The diamond-shaped pattern in the OTOC is in large part a result of initial states that are at very high densities. As Fig. 7 shows, a state with all spins initially up undergoes a “blinking” pattern in its time evolution; flipping a single spin creates two propagating domain walls in the blinking pattern, which lead to anticorrelation (rather than just decorrelation) in the OTOC. For a state with a high density of up spins, the OTOC remains strongly anticorrelated until the domain wall wraps around the system, and then decorrelates.

Outlier states.—We now turn to outlier states that are maximally nonthermal Kim et al. (2014); our proxy for non-thermality is an anomalously low . Fig. 8 histograms the distribution of across eigenstates at a given system size; note the “outlier” peak at . Exponentially many such outliers can be constructed, if one begins with states that have an atypical quasiparticle distribution. What is less certain is whether these outliers constitute a finite fraction of all states in the thermodynamic limit. Fig. 8 histograms the distribution of across eigenstates at a given system size; there is an “outlier” peak at . Since our numerical approach is restricted by the growth of , we can track the weight of this outlier peak for quite large systems (up to sites). As system size is increased, the outliers rapidly drop to about of the states. At large system sizes, the outlier “density” apparently decreases further with increasing system size (Fig. 8), but we have not been able to identify the functional form for this decrease.

Global recurrence times.—Our approach to estimating global recurrence times is simply to compute for many random initial configurations (typically ) and then find the least common multiple of all the . This is guaranteed to be a lower bound for the global recurrence timescale. This numerical lower bound grows with system size as (Fig. 9). For the smallest system sizes we have checked explicitly against exact diagonalization that our procedure gives the correct global recurrence timescale. We expect it to be reasonably accurate for the system sizes we have studied, as sampling over 2500 initial conditions instead of 5000 never misses more than one or two frequencies (so our estimate is probably off by only a few orders of magnitude).

OTOC for “pair-hopping” operator.—Given the quasiparticle structure of this system, it seems natural that the operator should scramble, as it usually changes the quasiparticle number. Correspondingly we have compared our results for this against the results for the OTOC with and with . The results are shown in Fig. 10; qualitatively these behave much like the results presented in the main text.

OTOC in single product states and in eigenstates.—We now present results for the OTOC in single product states and for the case where the quantum average is performed over a single random eigenstate. In the former case, the OTOC is always binary, but its pattern shows clear signs of the quasiparticle “tracks” going through the system (Fig. 11). When the OTOC is computed in a single eigenstate, as discussed in the main text, it spreads to some extent and then refocuses on a timescale of order , as seen clearly in Fig. 11. While this single-eigenstate OTOC fills in the light-cone, the spread of the front is narrower than for the infinite-temperature-averaged OTOC, although our results are too noisy to draw any definite numerical conclusions about the shape of the front in this case.