Solution of a minimal model for many-body quantum chaos

We solve a minimal model for quantum chaos in a spatially extended many-body system. It consists of a chain of sites with nearest-neighbour coupling under Floquet time evolution. Quantum states at each site span a $q$-dimensional Hilbert space and time evolution for a pair of sites is generated by a $q^2\times q^2$ random unitary matrix. The Floquet operator is specified by a quantum circuit of depth two, in which each site is coupled to its neighbour on one side during the first half of the evolution period, and to its neighbour on the other side during the second half of the period. We show how dynamical behaviour averaged over realisations of the random matrices can be evaluated using diagrammatic techniques, and how this approach leads to exact expressions in the large-$q$ limit. We give results for the spectral form factor, relaxation of local observables, bipartite entanglement growth and operator spreading.


I. INTRODUCTION
Random matrix theory plays a central role in the understanding of chaotic quantum systems. 1 It is founded on the idea that for many systems of interest there is no privileged basis in Hilbert space. Some important phenomena in this field, however, arise only if there is spatial structure and a notion of locality. Diffusive transport in weakly disordered conductors is such an example for single-particle systems, and spreading of quantum information is a counterpart for many-body systems. It is natural to attempt to combine the simplifying features of random matrix theory with extended spatial structure. For diffusive conductors this is achieved in Wegner's n-orbital model, 2 in which hopping between sites of a tight-binding system is governed by n × n random matrices and disorder-averaged properties can be calculated exactly in the limit n → ∞. Our aim in this paper is to establish a comparable simplification for spatially extended many-body systems.
Chaotic many-body quantum systems lie at the focus of efforts to understand the foundations of quantum statistical mechanics. [3][4][5][6] Generic features expected in the dynamics of such systems include rapid equilibration of local observables 5 and ballistic propagation of quantum information, 7 as well as ballistic growth of bipartite entanglement. 8,9 Conservation laws play a central part in dynamics, and one expects that systems with a given set of conservation laws will form a distinct class. Our focus in the following is on evolution arising from a time-dependent Hamiltonian: without even energy as a conserved quantity, it constitutes a particularly simple example.
Random matrix approaches offer natural routes to constructing models with minimal structure, and unitary quantum circuits provide an attractive way to formulate the evolution operator for time-dependent quantum systems. Unitary circuits that are random in both space and time have recently yielded valuable insights into chaotic quantum dynamics. [10][11][12][13] Here we initiate an analytic study of unitary circuits that are random in space but periodic in time.
We study a Floquet operator acting on a onedimensional system consisting of q-state 'spins' at each site. The Floquet operator is constructed from unitary matrices that couple adjacent sites. These q 2 × q 2 matrices are drawn independently from the circular unitary ensemble (CUE) and we compute physical properties averaged over the ensemble. Our key simplification is to treat the large-q limit. We show that quantum dynamics in this system exhibits a range of features that are expected to be characteristic of ergodic many-body quantum systems: correlators of local observables decay rapidly in time and quantum information spreads ballistically, in the sense that the bipartite entanglement of an initial product state grows linearly in time and the out-of-timeorder correlator 14,15 (OTOC) shows the 'butterfly' effect. Our approach also provides access to the spectral properties of the Floquet operator and we show (on scales much larger than the level spacing) that its eigenvalue correlations are those of the CUE. At a technical level, our calculations are based on the application of diagrammatic techniques developed previously for single-particle problems in mesoscopic physics involving random matrices from the CUE. 16 Some features of the large-q limit are non-generic: as for random unitary circuits 11,12 at large q, the distinction between the velocities associated with entanglement and operator spreading vanishes; and sublinear growth of entanglement, expected at long times in one dimensional disordered systems because of weak links, 17 is absent.
The balance of this paper is organised as follows. In Sec. II, we define the model, observables and results. In Sec. III, we develop the diagrammatic approach for taking the ensemble average of a given observable. In Sec. IV, we sketch the proof of the results presented in Sec. II using the tools developed in Sec. III. Lastly, in Sec. V, we summarise. Technical details are described in a series of appendices.

II. MODEL, OBSERVABLES AND RESULTS
We seek a minimal model for quantum chaos in a spatially extended many-body system with local interactions. We formulate the model directly in terms of the time-evolution operator, rather than a Hamiltonian, and for simplicity we consider a Floquet problem. Taking a one-dimensional lattice, evolution over one period is separated into two half-steps. Each even site is coupled to its neighbour on the left in the first-half-step, and to its neighbour on the right in the second half-step. This quantum circuit is illustrated in Fig. 1.
To write the evolution operator explicitly, let U 2n−1,2n denote the unitary matrix that couples sites 2n − 1 and 2n in the first half-step, and let U 2n,2n+1 be the matrix that couples sites 2n and 2n + 1 in the second half-step. Each U i,i+1 is independently distributed with the Haar measure and we calculate physical properties as an average over this ensemble, denoted by . . . . L sites the full Hilbert space has dimension q L . Taking L even for definiteness, the first half-step is represented by an evolution operator acting in this space with the form The evolution operator for the second half-step is similarly where 1 1 q denotes the q × q unit matrix. The Floquet operator describing evolution over a single complete period is We denote the evolution operator for an integer number t of periods (the t-th power of W ) by W (t).
To demonstrate that this model has the features expected in a chaotic quantum many-body system, we examine a range of physical properties, as detailed below. Our results for entanglement spreading and for the OTOC are the same as those for random unitary circuits with the same structure as our model, if one takes the large q limit in previous results. 11,12 In this context, one of our significant conclusions is that random unitary circuits do indeed share important physical features with deterministic systems. At the same time, our model opens up discussion of the spectral properties of the evolution operator. This has no equivalent for random unitary circuits since they lack any such fixed operator.

A. Spectral form factor
Our results for the spectral form factor show that the Floquet operator for the model has exactly the same eigenvalue correlations in the large q limit as an ensemble of Haar-distributed unitary matrices. This insensitivity to the spatial structure of the system is a striking emergent feature.
The spectral form factor K(t) is the Fourier transform of the two-point correlation function of the eigenvalue density. Denote the eigenphases of W by {θ m } for m = 1, . . . q L . Then For N × N matrices from the CUE one has 1 The behaviour of K(t) on scales |t| N reflects level correlations at separations much larger than the mean spacing, and the linear dependence of K(t) on |t| in this regime is a consequence of Coulombic suppression of long wavelength fluctuations in the eigenvalue density. Conversely, the form of K(t) for |t| ∼ N encodes spectral correlations on the scale of the level spacing.
We obtain exactly the CUE form, K(t) = |t| for t = 0, in the large q limit at fixed L and t. We stress that this result is a consequence of coupling between all sites, and is not simply a reflection of the properties of individual random matrices U i,i+1 . To illustrate the point, consider an alternative model in which the couplings U 2n,2n+1 on even bonds are replaced with unit matrices 1 1 q 2 . For this toy system of L/2 independent pairs of sites one has much larger fluctuations in eigenvalue density, with Returning to the original model, it is interesting to speculate on behaviour in regimes other than the one we are able to analyse exactly. At large but fixed q and fixed t, it is natural to expect K(t) to increase with L if L is sufficiently large, since distant parts of the system should be only weakly correlated; results for this regime will be presented elsewhere. 18 In the opposite limit, fixing q and L and varying t, the form of K(t) at |t| ∼ q L probes correlations on the scale of the level spacing. We expect these to be of CUE form for all q and L, provided only that q L 1. They should therefore show a transition from K(t) = |t| for |t| q L to K(t) = q L for |t| ≥ q L . The transition to constant K(t) at large t is, however, well-known to be inaccessible in a perturbation expansion in inverse powers of q. 19,20

B. Dynamics
We next examine physical properties chosen to reveal both the dynamics of local degrees of freedom and the spreading of quantum information.

Local relaxation
A characteristic feature of chaotic dynamics is that local observables relax rapidly towards equilibrium. To probe relaxation in our model we introduce an operator O x representing an observable at site x and obeying and In our Floquet model the infinite-temperature density matrix ρ = q −L 1 1 q L is the natural choice and we denote the normalised many-body trace by tr ≡ q −L Tr. One anticipates that [O(x, t)O(x)] av will relax to the value [O(x)] 2 av (which is zero with this form for O(x)) on a microscopic timescale of order the Floquet period. At large q we in fact find complete relaxation within a single period, obtaining for q → ∞ the result We discuss the leading corrections to this result for finite q and small t in Sec. IV B.

Entanglement spreading
We probe entanglement spreading via the time dependence of the ensemble-averaged bipartite entanglement purity (and, more generally, ensemble averages of exponentials of Rényi entropies) for an initial state |ψ that is a direct product over sites. Specifically, let A be the left half of the system (with site labels 1 ≤ n ≤ L/2) and B the right half, and consider the reduced density matrix The Rényi entropies S α (t) are given by and the entanglement purity is P(t) = e −S2(t) . We compute e (1−α)Sα(t) for α = 2 and α = 3. The purity shows an exponential decay in time until it falls to a value P(t) ∼ q −L/2 typical of random states in the many-body Hilbert space. In detail, we obtain at large q The function f α (t) grows exponentially with t, and in particular, we find f 2 (t) = 4 t and f 3 (t) ((4+3 √ 2)/2) 2t for t 1. At times t > L/4, the Renyi entropies saturate and we conjecture that K α = Cat(α), where Cat is the α-th catalan number. The latter result is proven for α up to 5 (see Appendix F). This is consistent with the generalized Page formula for bipartite entanglement in random states 2122 and with the results from random unitary circuit 23 .
An interpretation of these results is that the reduced density matrix for the pure initial state spreads at time t over the Hilbert space spanned by basis states at 2t sites. As an illustration, suppose that ρ A (t) has q 2t nonzero eigenvalues that are all equal. Then Tr A ρ A (t) = q −2(α−1)t , which is consistent with the leading order behaviour of (10) for t ≤ L/4. This demonstrates that the entanglement spreads ballistically with a velocity at large q that is the same as the naive light-cone velocity, v = 2, introduced below in Fig. 2. The values of f 2 (t) and f 3 (t) give information on the distribution of non-zero eigenvalues of ρ A (t), and it is noteworthy that these two quantities are distinct.
The growth of entanglement found here is similar to that obtained for integrable systems using conformal field theories 24 . In integrable systems this behaviour is associated with the presence of quasiparticles that travel ballistically; a quite different physical picture is required for ergodic systems.

Out-of-time-order correlator
The spreading with time of local operators is characterised by the behaviour of the (ensemble-averaged) commutator which measures how a local perturbation at site x affects measurements at a later time at the site y. With our normalisation for local operators, one has and the second term on the right-hand side is the OTOC. For short times t and large separations x − y, the operators O(x, t) and O(y) commute and C(x − y, t) is zero, while for large times the OTOC is small and C(x − y, t) approaches unity. We obtain at large q a sharp light-cone behaviour Hence operator spreading occurs with a butterfly velocity that, like the entanglement velocity, is equal to the naive light-cone velocity v = 2 in the large q limit. Results on random unitary circuits 11,12 suggest that all three velocities should be distinct for q finite. From the same comparison, we also expect at finite q that the step function of (13) will broaden.

III. ENSEMBLE AVERAGING
We now set out in several steps the general approach that we use to obtain these results.
First, it is useful to extend the notation of Fig. 1 in various ways, so as to represent pictorially the quantities defined in Sec. II. An example for Tr[O(x, t)O(x, 0)] is shown in Fig. 2. The vertical timelines of sites join rectangles representing factors of U i,i+1 and indicate matrix multiplication. Repeated copies of W denote multiple time steps; W † is shown as a differently shaded version of W ; and local operators appear as squares. The matrix trace is shown by joining the timelines of sites to form closed loops. Finally, product states in the site basis (which do not appear in this example) are shown using circles at the ends of the timelines.
A straightforward simplification in many instances is that some factors of U and U † cancel, as illustrated in Fig. 2b, and the naive light-cone velocity v = 2 clearly emerges.
As a second step, it is helpful to fold the pictorial representations, so that while the timelines in W run upwards, those in W † run downwards, as shown in Fig. 2c.
These folded pictures provide a direct depiction of physical quantities but are cumbersome. A simpler representation is possible if we focus on the time evolution of a single site. To this end, and in anticipation of the disorder average, we switch to an alternative notation for U i,i+1 and U † i,i+1 in which individual sites appear separately, as shown in Fig. 3.
With these preliminaries in place, we can set out a diagrammatic representation for ensemble averages. It is a many-body extension of the one introduced for averages over the CUE by Brouwer and Beenakker. 16 It can be applied to an arbitrary observable f (W ) but to be definite we illustrate it for K(t) with t = 2 and L = 2, so that W (2) = [U 1,2 ] 2 . In general three steps are involved.
(i) The observable f (W ) is represented as a collection of single-site diagrams using the notation introduced in Fig. 3. For the example of K(t=2), one starts from the many-body diagram shown in Fig. 4 (right) and obtains the single-site representation shown in Fig. 5 (left).
(ii) The ensemble average f (W ) is computed by generating a collection of single-site 'contracted' diagrams G = {G i } L i=1 as follows. On each site, filled U -dots are connected to filled U † -dots of the same color on the same site with a dashed line (a contraction) in all possible ways, and likewise for the empty dots.
Since U s and U † s act on neighbouring pairs of sites, these contractions must be matched: for even i the choice of contractions between blue dots must be the same in the diagrams G i−1 and G i ; similarly for the red dots and the diagrams G i , G i+1 . We refer to this as the bond constraint.
(iii) Each contracted site diagram G i gives rise to an algebraic expression A(G i ), obtained as the product of two factors These factors are associated with loops of two kinds, called T -loops and U -loops in Ref. 16. The T -loops simply record pairings of matrix labels, and give rise to powers of q in the contribution of a diagram. The U -loops distinguish different contributions from the average of factors of U and U † . Examples of these two types of loop are illustrated in Fig. 6.
A T -loop is a closed sequence of alternating single and dashed lines. It carries an index corresponding to one of the q basis states at a site. We generate the contribution associated with each T -loop by summing over its index. This leads to a factor of q for a T -loop of single lines if it does not pass through any operator insertion, or of k } bi k=1 be the sets of lengths of red and blue U -loops in G i . Then from the theory of CUE averages 16 where the explicit form of the coefficients V Ri ≡ V c1...cr i and V Bi ≡ V c1...c b i (also known as the Weingarten functions) is given in App. A. Note that the exponent 1/2 in (15) arises because we distribute the contribution arising from unitary operators equally over the two sites on which the operators act. The final value of the ensemble average reads where the sum runs over all contracted diagrams G satisfying the bond constraint. For the example of K(t=2) , , which gives K(t=2) after ensembleaveraging, using the notation of the left and right sides of The factors of q within the bracket are contributions from the T -loops. We have taken the large-q limit on the far right, using At large q the first and fourth diagrams are leading order, while the second and third are sub-leading. The above procedure is exact for any q, but the sum in (16) is problematic in general as it may involve an extremely large number of terms. For example, a total of (t!) 2(L−1) diagrams contribute to K(t). It is therefore useful to establish which terms dominate in the largeq limit. For large q one has A U (G i ) ∼ q ui−ni where u i = r i +b i is the total number of U -loops and n i = k c k is the number of contractions on the i-th site. Let τ i be the number of T -loops in G i that contribute a factor q. Then we have the large-q expansion where we have introduced the order O(G) and omitted a proportionality constant independent of q. In Appendix B, we discuss ways of enumerating the leading order diagrams. As a general rule, since the total number of contractions i n i is fixed, we should retain in (16) only the diagrams that maximize the total number of T -and U -loops. A natural approximation at large q is to treat the elements of U as independent Gaussian random variables, so that a standard Wick theorem applies. This approximation corresponds to including all diagrams where all U -loops are of unit length, and omitting all others. We refer to this set as the Gaussian diagrams. As an example, in Fig. 5, the first and fourth diagrams on the right hand side are Gaussian. For these diagrams, Eq. (18) holds with proportionality constant one; therefore counting the overall number of U -and T -loops is sufficient to obtain the leading contribution of the diagram.
In this paper, the leading contributions to all quantities calculated (except for the autocorrelation function of a local observable for t > 0) are Gaussian. A procedure that goes beyond the Gaussian approximation is necessary, firstly, for the proof of this statement, and secondly, for future discussion of the sub-leading contributions. Since the number of diagrams contributing to each of the quantities we consider is finite at fixed L and t, our results are exact in the limit q → ∞ with L, t fixed.

IV. DIAGRAMMATIC EVALUATION OF RESULTS
We now show how this diagrammatic approach can be used to generate the results given in Sec II. We sketch the main ideas here, deferring formal proofs to appendices.

A. Spectral form factor
Consider the spectral form factor K(t) introduced in (5). Applying the procedure described in Sec. III (and generalising the example given for L = 2 and t = 2 in Fig. 5), we obtain for L > 2 a many-body diagram consisting of single-site diagrams as shown in Fig. 7a. Here each timestep in W (t) contributes with two unitaries, leading to 4t dots divided into four types -blue or red, and empty or filled. Three examples from a total of t leading order diagrams for K(t) at t > 0.
Next we consider ensemble-averaging this diagram, making all possible contractions. On a single site, one can easily check that one of the diagrams with the maximum number of loops is, for example, the site diagram in Fig. 7b. Since this diagram is Gaussian and τ i = u i = 2t and n i = 4t, its contribution for large q is simply 1.
There are t equivalent diagrams of this kind, obtained by cyclically shifting the right dots with respect to the left ones, as shown in Fig. 7c and Fig. 7d. Additionally, once one of these configurations has been chosen on the site i, the bond constraint forces all other sites to be in the same configuration in order to maximise the number of loops. All other diagrams are smaller by powers of q for q → ∞. In consequence we get the result K(t) = |t|. A proof of this statement is given in Appendix C.

B. Relaxation of local observables
Contributions to the autocorrelation function are generated by contractions of the site diagrams shown in Fig. 8a and 8b. The leading contribution G i at a site i = x is from a contraction of the form shown in Fig. 8c. However, if this contraction is made at every site i = x, then because of the bond constraint it also applies at site x and yields Fig. 8d. An example of an alternative contraction, for which A(G x ) = 0, is shown in Fig. 8e. With such a choice at i = x, the bond constraint imposes contractions at nearby sites i = x that are sub-leading in q. We have also evaluated the leading non-zero contributions to autocorrelation function for large q at small values of t. We find These results are interesting for two reasons. First, the equivalent quantity for time evolution with a random unitary circuit is identically zero at all t = 0. The results hence expose a difference between our Floquet model and random unitary circuits. The finite relaxation rate at finite q in the Floquet model is consistent with expected generic behaviour, whereas complete relaxation for any t = 0 is likely to be a special feature of random unitary circuits. Second, and quite separately, the dominant contributions arise from non-Gaussian diagrams: for example, at t = 3 the largest Gaussian term is 2q −9 , but is cancelled by non-Gaussian contributions.

C. Purity and Rényi entropies
We now show that the large-q calculation of the purity can be reduced to counting domain wall (DW) configurations with positive weights, in a problem analogous to one from classical statistical mechanics. Related ideas apply (albeit more elaborately) to the evaluation of averages of the exponential of the Rényi entropy S α (t) for general positive integer α and we discuss the case α = 3 in Appendix F. Domain structure similar to what we derive here also appears in recent treatments of random unitary circuits 10,11 and it is striking to find it as an emergent property of our Floquet model.
The main steps in our procedure are summarised in Fig. 9. We wish to average Tr A [ρ A (t) α ] for α = 2. Using the conventions of Fig. 2, Eq. (9) has the pictorial representation shown in Fig. 9a. Since each ρ A (t) = W † (t) |ψ ψ| W (t), we have four sectors, containing W (t), W † (t), W (t) and W † (t) respectively. It is convenient to fold this diagram, as discussed in Sec. III and Fig. 2. This procedure leads to a folded representation containing four layers, as shown in Fig. 9b top. The folded site diagrams in regions A and B are shown in Fig. 9b bottom left. Note in particular how the timelines of sites in the regions A and B connect differently at the top of these diagrams, because of the structure of traces in (9). After averaging, it turns out that the leading contributions come from contractions of the U 's and U † 's that are brought to lie in a stack on top of each other by the folding operation. We call a contraction of this type local. Such a stack of unitaries is indicated with a large blue box in Fig. 9b top. In the evaluation of the purity (where two U 's and two U † 's are involved) there are two possible local contractions (Fig. 9b bottom right).
An intuitive explanation of the fact that only local contractions contribute for q → ∞ is that any contraction between two dots in distant blocks would necessarily lead to a longer loop; as the total loop length is fixed, this implies a smaller number of loops and a lower order according to (18). This statement is proved in Appendix E.
Given the restriction to local contractions, a further simplification of the diagrammatic representation is possible. It is no longer necessary to represent the unitaries within a stack individually. Instead we can simply depict the stack, using a label to indicate the type of local contraction. We label the two types of local contraction that appear in a calculation of the purity A and B, because they involve the same types of pairing as are induced by the trace structure in regions A and B of the system. The diagrammatic notation applied to one stack is shown in Fig. 9b bottom right. The same notation is used for a larger system in Fig. 9c.
We prove in Appendix E that the order of a diagram with only local contractions is given by q −h AB , where h AB is the number of segments of horizontal wall of length one lattice unit between an A and a B block. Such walls are shown with horizontal red lines in Fig. 9c. In this framework, the problem of finding the leading contributions admits a simple geometrical interpretation: they are represented by the set of minimal-length DW diagrams that separate the A and B contractions. Given the fixed boundary conditions for the blocks on the top of the diagram, a DW must connect the top centre to either the side or the bottom of the diagram (Fig. 9c).
For t ≤ L/4, a minimal-length DW can only connect the top centre to the bottom edge of the diagram. Since there are a total of 2t rows, the minimal h AB is 2t, and the order is q −2t . Moreover, there are two choices (for the DW to go left or right) below every row, so the number of leading order diagrams is 2 2t . For t > L/4, the lengths of the DW-s are minimized when they connect the top centre and the side of the diagram with the shortest possible paths. This implies that there are only two leading order diagrams: those where the DW is directed exclu- sively towards either the left or the right. For these cases, h AB = L/2 and the order is q −L/2 . In summary, Gaussian diagrams give the leading behaviour at large q: Subleading terms eliminate the discontinuity in (20) at t = L/4. For larger α, one can repeat the construction of Fig. 9. Each block has now α! possible states, corresponding to the possible local contraction inside a block of α U 's and α U † 's. So, generically the problem reduces to counting, but for large α the procedure quickly becomes problematic. The case α = 3 is discussed in Appendix F.

D. Out-of-time-order correlator
Consider the diagrammatic representation of the OTOC defined in Eq. (12). After the cancellation of U s with U † s wherever possible, we can distinguish two alternative scenarios.  (11), we see that the commutator vanishes. This is consistent with the fact that the light-cone spreading from x with velocity v = 2 has not yet reached the operator in y. Conversely, for t ≥ |x − y|/2 there is a region of unitary operators that cannot be cancelled.
Calculation of the OTOC for t ≥ |x − y|/2 is shown in Fig. 10 left using a folded diagram. The diagram contains four layers of W s, as for the purity, and we employ the block representation introduced in Sec. IV C. Similarly to Sec. IV C, it can be shown that the leading order diagrams necessarily involve only local contractions. Moreover, the trace structure of OTOC sets boundary conditions along the top of the diagram that act as Ablocks, and along the bottom that act as B-blocks. From this one would expect that the leading order diagrams are minimal-length DWs that cross the diagram horizontally, separating an upper domain of As from a lower domain of Bs (Fig. 10 right). These diagrams however vanish, because they involve factors of either tr O(x) ≡ 0 or tr O(y) ≡ 0. In Appendix G, we prove that the leading contributions to the OTOC for t ≥ |x − y|/2 go to zero for q → ∞, and hence Eq. (13) is recovered in the large-q limit.

V. SUMMARY AND OUTLOOK
In summary, we have introduced a minimal randommatrix model with extended spatial structure to study the chaotic Floquet dynamics of a many-body quantum system. We have presented a diagrammatic technique to compute several quantities using a systematic and controlled expansion in the inverse of the local Hilbert space dimension q.
Our study of a minimal model and the techniques we have developed are complementary to a variety of other recent works. In particular, the semiclassical approach to quantum chaos in few-body systems 25 has been extended to bosonic systems at high density by considering interfering paths in Fock space that arise from solutions to the Gross-Pitaevskii equation, 26,27 and to periodically driven spin systems in the large-spin limit. 28 These semiclassical techniques have the attraction of applying directly to specific physical systems that are of wide interest, rather than simply to a minimal model. They been used so far to identify particular many-body interference phenomena, but have not been developed to allow general computation of the dynamics of quantum information. In a separate advance, the Keldysh technique has been generalised to permit calculation of out-of-time order correlators, with applications to a variety of microscopic models. 29 These calculations are well-controlled in a quasiclassical regime, analogous to our large-q limit, while the augmented Keldysh contour of Ref. 29 is (unsurprisingly) mirrored quite closely by the structure of the diagrammatic calculations we describe here. Further topical research, 30 addressing a spatially extended version of the SYK model, 31 is set apart from the results we present by the fact that the zero-dimensional SYK model exhibits much greater structure than the individual random matrices of our model.
There are some obvious and interesting directions for additional investigations using the techniques we have set out. First, in contrast with random unitary circuits, our model has a well-defined Floquet operator, opening the possibility for the study of its spectral properties. It seems likely that further work in this direction will be useful, going beyond our evaluation of the spectral form factor in the random matrix regime. Second, the results we have presented are averages over an ensemble of systems. It would be useful to understand the magnitude of sample-to-sample fluctuations, by evaluating quantities such as (e −Sα(t) ) 2 − e −Sα(t) 2 . More ambitiously, it would be appealing to use higher order terms in the 1/q expansion to search for the expected differences between the naive light-cone velocity, the entanglement-spreading velocity and the butterfly velocity, and to investigate broadening of the step in the OTOC given in (13) for the large-q limit.
There are also several generalisations. Our model and the techniques we have developed can be naturally extended to higher dimensions. Analogous models could also be developed for the other symmetry classes in random matrix theory. In the context of quantum transport, it would be interesting to incorporate the presence of conserved quantities, by modifying the local structure of the Floquet operator, as has been done recently for random quantum circuits. 32,33 Finally, and speculatively, it is possible that an understanding of the ergodic phase in our model for chaotic many-body systems may play a role analogous to the treatment of a diffusive metal in the theory of disordered conductors, and provide a starting point for a theory 34 of the many-body localisation transition. 35 After completing this work, we learned of recent related calculations of the spectral form factor for a manybody system, in which the equivalent of Eq. (5) is derived for the orthogonal symmetry class 36 . There are however important differences between the two approaches. In particular, the result of [36] develops from a specific long-range Hamiltonian using ideas of periodic orbit theory. Instead, the present work utilises a limit of large onsite Hilbert space and disorder averaging still retaining a local short-range structure. Therefore, one important potential use of our method is to address ergodicity in disordered systems 37 .
where V 0 ≡ 1 and N = q 2 is the dimension of the unitary group from which the Haar-distributed unitary operators are drawn.

Appendix B: Enumeration of leading order diagrams
In the following, we describe two methods for efficiently enumerating the leading order diagrams for the physical quantities we have computed. The contraction addition method described in App. B 1 can be used to eliminate sub-leading single-site diagrams efficiently, which is particularly useful if the many-body diagram of interest has the same site-diagrams across all sites (e.g. the spectral form factor). The domain wall (DW) approach described in App. B 2 allows us to obtain an upper bound to the order, a global property of a diagram, by making only local calculations between neighbouring domains. We apply the ideas from App. B 1 in App. C, D and G, and those from App. B 2 in App. E and G.
The foundation for both these methods is Eq. (18), which gives the order O(G) of a diagram G as a product of factors q τi+ui−ni from each site. As the numbers τ i of T -loops and u i of U -loops contribute to the order in the same way, we do not need to distinguish between T -and U -loops in the calculation of order discussed below.

Method of contraction-addition
An ensemble-averaged many-body diagram consists of L site diagrams labelled by i, each with a number n i of contractions. The idea of this approach is that, given a site diagram, we can first remove all the contractions, and then re-construct the diagram by adding contractions one at a time. The order can be evaluated by considering the effect of each contraction-addition. The crucial point is that the procedure is independent of the sequence of contraction-additions. Each addition either leaves the order in q unchanged (if it increases both τ i + u i and n i by one), or reduces the order (if the addition increases only n i , or reduces τ i +u i and increases n i ). A convenient sequence (one that contains order-reducing contractionadditions in the first few steps) can therefore be used to eliminate sub-leading diagrams efficiently.
To be more precise, given a diagram G, we consider the single-site configuration G i on site i. We choose a sequence of diagrams {G Here, ∆ ) is the change in the order when one contraction is added.
As we have established that we do not need to distinguish between T -and U -loops here, each contracted site diagram can be represented as a collection of loops, using the same notation for both types of loop. The same representation can also be used for each intermediate G . The addition of a contraction has the effect either of merging two loops into a single longer one, or of breaking one loop into two. To see this, recall the prescription for construction of loops in a single-site configuration G i that is implied by Fig. 6. It can be expressed as follows. Start at any point on a single (for T -loops) or double (for U -loops) line; follow the direction along this line; when a dashed line is encountered, go to the other end of the dashed line and continue along the single or double line in the same sense; repeat until the starting point is reached: the path traced is a loop. Addition of an extra contraction to this prescription simply re-routes the paths of the two loop segments that meet the contraction from either side. Moreover, since the two ends of the dashed line depicting a contraction are attached either to two separate loops or to two different portions of the same loop, we can evaluate ∆ (m) i omitting all information except what concerns the one or two loops involved in the contraction. Once one focuses on the relevant details in this way, there are only eight distinct contraction-addition scenarios. We list them all in Fig. 11, using a single unbroken line to denote generic loops. A minor complication is that it is necessary to distinguish loops that pass through one or two operators representing local observables from ones that do not, since these operators may also affect the weight of a diagram.
We provide a simple example of how this approach to determining the order of a site diagram works in Fig. 12.

Domain wall approach
In Sec. IV C we set out a way of considering contributions to some quantities of interest, in terms of spacetime domains within which contractions are all of the same type. The power of this approach lies in the fact that one can associate a bound on the diagram's order with each DW between two neighbouring domains. This allows us to calculate a bound on the diagram's order,  a global property, by making only local calculations between neighbouring domains. We explain details here, using for concreteness the 2nd Rényi entropy (see Sec. IV C) and referring to the block representation introduced in App. E 1 and Fig. 13.
First, we rewrite the order of a diagram G in terms of loops. We define the length c k of a loop k as half of the number of non-dashed lines it contains. We make the convention that single lines at the top edge (where W and W † are multiplied) are counted twice, so that lengths are always integer numbers (e.g. at the top edge of Fig. 5 there are four single lines in each single-site diagram). This generalises the definition given in Sec. III to both U -loops and T -loops. We also define the total length l t = k c k where the sum runs over all loops. Then the order of a diagram in (18) can be re-written where Γ c = q 1−c and n = i n i is the total number of contractions. The values of l t and n are fixed by the observable of interest. In particular, for the 2nd Rényi entropy l t = n. The intuition behind Eq. (B2) is as follows. Since a given diagram has a fixed number of length segments, its order is high (low) when there are many short (a few long) loops. The pre-factor q lt−n represents the theoretical maximum order when all length segments are used to form 1-loops. Γ c is the cost associated with a loop of length c, and its value decreases as c increases. For example, Γ 1 has value 1 since 1-loop is the shortest possible loop. Γ 2 has value q −1 because one could have formed two 1-loops (giving a factor of q 2 in Eq. (18)) instead of one 2-loop (giving a factor of q).
Next, we rewrite the order of a diagram in terms of boundaries that intersect the loops. We define a boundary as a horizontal segment of unit lattice spacing that separates two blocks (for example, the red lines in Fig. 15 top), so that each block has four boundaries (except for the blocks representing the trace structure at the upper edge, e.g. zero-th row in Fig. 15 top). Every loop instersects two such boundaries per unit length. To each boundary w, we can therefore associate a cost where k labels the loops that the boundary w intersects, and c k is the corresponding loop length. In this way, the order of a diagram G is expressed in terms of the cost of each boundary, so that Since Γ 1 = 1, loops of unit length do not contribute to O(G). We obtain an upper limit on the cost of a boundary by assuming longer loops to have length not more than c k = 2. This leads to where n w is the number of loops crossing the boundary w with length longer than 1.
Referring to examples in Fig. 15, C(w) for the top sub-figure is bounded by (and equal to) q −1 with n w = 4, and C(w) for the bottom sub-figure is bounded by q −3/4 with n w = 3. As we will see, a costly boundary (one with n w > 0) is always sandwiched between blocks with different contractions. It is therefore natural to call such a boundary a domain wall (DW). With minor modifications, this method can also be applied to diagrams that contain local operators.

Appendix C: Spectral form factor
Here we use the method of contraction-addition (App. B 1) to enumerate the leading order diagrams for K(t) . In Eq. (B1) we have n i = 4t and O(G There are only two relevant contraction-addition scenarios illustrated in Fig. 11: (i) An intra-loop addition (item 1 in figure), where the two legs of the new contraction line land on the same loop. Then ∆ = q 0 according to Eq. (18) because τ + u and n both increase by one. (ii) An inter-loop addition (item 2 in figure), where the two contraction legs land on two different bare loops. Since the two loops merge into a single bigger loop due to the new contraction, τ + u decreases by one, and n increases by one. So ∆ = q −2 .
The many-body diagram for K(t) comprises site diagrams before contraction as in Fig. 7a. For t = 0 we have K(t) = q 2L since W (0) = 1 1 q L . To compute the large-q limit of K(t) for t > 0, we first note that there are diagrams with multiple contractions that are O(1), such as Fig. 8c. This is in fact a highest order diagram for t > 0 according to Eq. (B1), because there must be at least one inter-loop contraction-addition which is associated with ∆ = q −2 per site, and because, from Fig. 11, the later contraction-additions can be made without increasing the order of the diagram. So any diagrams with order smaller than O(1) are sub-leading.
There are t leading order site diagrams on site i, by the following argument. Using the fact that the order determination is independent of the sequence of contractionaddition, we choose to contract on site i the (filled and un-filled) blue dots on the top layer from left to right. The first filled blue dot can be contracted with any one of the t filled blue dots on the bottom layer. This interloop contraction costs q −2 . In order to obtain a leading diagram, the later contractions must all be of O(1), i.e. intra-loop contractions. A general feature of an O(1) contraction is that it must partition a bigger loop into two smaller loops, such that the q factor associated with the extra loop (since τ or u increases by one) cancels the q −1 factor of the contraction (since n increases by one) (see Fig. 11). Furthermore, on each of the two smaller loops, there must be equal numbers of un-contracted blue dots on the top and bottom layer. Otherwise, there will be an inter-loop contraction which will render the diagram subleading. It is straightforward to see that there is a unique choice of contraction of the second blue dot that satisfies this requirement. Similarly, there are unique choices for the rest of the blue dots on this site. In order to not incur further cost on site i, each red dot on the top layer must be contracted with the only other red dot on the same loop on the bottom layer. So there are unique choices for the red dots as well.
Due to the bond constraint, the (i + 1)-th site diagrams inherit the choice of either the blue or red dot contractions on site i. We can repeat the above analysis site-by-site and conclude that there are t leading order diagrams of O(1) for t > 0 (Fig. 7b,c and d).
Lastly, each of these t diagrams is translated algebraically to the factor 1, simply observing that it is Gaussian and for each site τ i +u i = n i in (18). In other words, for t > 0, K(t) = t, and we have arrived at Eq. 5.
Here we prove Eq. (7). For t = 0, the calculation is straightforward since W (0) = 1 1 q L . For t = 1, there is only one configuration at each site because there is only one pair of dots of each kind. Since we have chosen Tr O x = 0, the algebraic term associated with this diagram vanishes. For t > 1, we use the same argument as in App. C to show that every possible diagram has an order that vanishes as q → ∞. In the formulation of contraction-addition (App. B 1), we have n i = 4t and O(G (0) i ) = 1. On site x, we choose to first contract the blue dot on the bottom left of the diagram (Fig. 8b). This bottom left blue dot cannot be contracted with the bottom right blue dot non-trivially since Tr O x is zero. However, if this bottom left blue dot is contracted with any other blue dots, it will partition a bigger loop into two smaller loops, each of which has an unequal number of dots from U s and U † s. This implies that there will be at least one addition of an inter-loop contraction with ∆ ∝ q −1 (item 6 in Fig. 11). Since all other contractionadditions can only further reduce the order according to

Appendix E: 2nd Rényi entropy
The diagramatic representation used in Fig. 9c requires a straightforward extension when it is empolyed in the presentation of detailed proofs, and we start by describing this extension. It arises for the following reason. The rectangles of Fig. 9c Fig. 9c is divided horizontally into two in this Appendix.

Block representation of exp[−S2(t)]
In detail, we introduce the block representation as follows. Beginning from the folded representation introduced in Sec. IV C, we consider a stack of four unitary operators (blue box in Fig. 9). We use Fig. 3 to represent each unitary operator in terms of dots and double lines. A block is defined as a region within this box that encloses only filled (or only empty) dots. In this way, to each stack of four unitaries, we associate two blocks containing respectively filled and empty dots (so that, the four-legged symbol in the bottom right of Fig. 9b corresponds to two blocks in Fig. 13 bottom). For exp[−S α (t)] , there are 4t rows of such blocks, and drawing each of them as an empty rectangle leads to the representation in Fig. 13.
We categorize the possible contractions of the dots within a block into seven types T = {a, b, a 1 , a 2 , b 1 , b 2 , x}.  a a a a a a a a a b b b b b b Table I). Bottom: The correspondence between blocks and the symbols used in the bottom right of Fig. 9b. As explained in the main text, a and b only involve local contractions (within the block), while x involves only non-local contractions. The other types are of mixed local and non-local character, and are defined in Fig. 14.
An upper bound ω(c, c ) for the cost of the boundary between two neighbouring block configurations c, c ∈ T is shown in Table E 1. Two examples of the evaluation of ω(c, c ) are illustrated in Fig. 15 using the method introduced in Appendix B 2.

Evaluation of exp[−S2(t)]
We discuss diagrams for exp[−S 2 (t)] in terms of different block pairings. We first fix our convention as follows. We label the rows of blocks from the top as in Fig. 13. We refer to the boundaries immediately above the p-th row of blocks as the p-th row of boundaries.
Due to the simplification illustrated in Fig. 2, pairs of U and U † can be cancelled outside of the light cone originating from the sub-system boundary (shown with dashed lines in Fig. 13). Alternatively, instead of cancelling them, it is equivalent to assume that the blocks on the far left and far right have respectively a-and b-contractions, as long as we remain in the large-q limit. Consistently, this choice implies costless boundaries according to Table I. Therefore, in particular, a configuration for any given odd row of boundaries can be parametrised as in Fig. 16, and we denote this as (ad 1 d 2 d 3 . . . d k b) wall with variables d 1 = a, d k = b, and d i = d i+1 for i = 1, . . . , k − 1. a, b and d i represent (connected) domains of blocks with the same type of contractions. Note that for every change of domain, there is a DW associated with a factor of q −m with m > 0. contractions. We classify the leading diagrams at t by the width r of the j i -region on the bottom row of blocks. To count the number of leading diagrams, we use an inductive approach from time t to t + 1. Every leading diagram at t + 1 can be generated from a leading diagram at t by adding four rows of blocks at the bottom of the diagram at t. A newly-added even row must have a configuration identical to the one above. A newly-added odd row has configuration depending on the one above according to the following rules (Fig. 18). For each even row in which the width of the j i -region is zero, there are five possible configurations for the odd row below: two with width 0 and three with width 1. For each even row with width r > 0, there are four possible configurations for the odd row below, two with width r, one with width (r − 1) and one width (r + 1). An example of a leading diagram is given in Fig. 19.  18. The recursive rules for enumerating the p+1-th row of the leading order diagrams given the p-th row of the leading order diagram. These rules are derived from the fact that, in a leading diagram, the difference between the DW positions at even row p and odd row p + 1 is exactly one lattice spacing.
We can write this recurrence relation in matrix form, by introducing a vector v(p) = (v 0 (p), v 1 (p), . . . ) T where v r (p) denotes the number of diagrams with total number of rows p that have a j i -region of width r in the last row. The degeneracy D(t) can be obtained by summing up the components of v(t). These in turn can be found by acting 2t times with a transfer matrix representing the recursive rules on the vector v(t = 0) = (1, 0, 0, . . . (F1) On the right side of this expression we have used the fact that, for 1 t < L/4, D(t) is dominated by the largest eigenvalue of the transfer matrix.
For t > L/4, there are five leading order diagrams. The first two have ab DW directed only towards the left or the right. The other three diagrams have only j icontractions within the light-cone region (dashed lines in Fig. 19), with the aj i and j i b DW-s directed only towards the left and the right respectively After translating these diagrams into algebraic terms, we obtain (10) for α = 3.
For higher Rényi entropies, we can prove for α = 4 that leading diagrams consist of local contractions only, and that K 4 = Cat(4) in Eq. 10. If we assume the analogous statements about local contractions for α > 4, we can apply the above approach and show that K α = Cat(α) up to α = 10.  a a a a a a a a a b b b b b b b  i ) = 1. We choose to contract blocks of U s and U † s bondby-bond from the left to the right in Fig. 10 . For the left-most bond, there are always two blocks of filled and empty dots for all t ≥ x/2. For each of these blocks there are two choices of contraction: A or B. A choice of contraction-A in one of the two blocks can be followed by other contractions in two alternative ways. Either all other blocks within its upward light cone (in this case a stripe of blocks) have A-contractions: this results in a factor ∝ (Tr O) 2 = 0 at site x, due to contractionaddition of type seven in Fig. 11. Alternatively there is a q −1 cost due to at least one inter-loop contractionaddition of type two, three or six in Fig. 11. Similarly, the choice of contraction-B would give a vanishing contribution because of the operator in y. These arguments are illustrated in Fig. 20. So we have shown that all diagrams have an order bounded by q −1 for t ≥ |x − y|/2.