Local pairing of Feynman histories in many-body Floquet models

We study many-body quantum dynamics using Floquet quantum circuits in one space dimension as simple examples of systems with local interactions that support ergodic phases. Physical properties can be expressed in terms of multiple sums over Feynman histories, which for these models are paths or many-body orbits in Fock space. A natural simplification of such sums is the diagonal approximation, where the only terms that are retained are ones in which each path is paired with a partner that carries the complex conjugate weight. We identify the regime in which the diagonal approximation holds, and the nature of the leading corrections to it. We focus on the behaviour of the spectral form factor (SFF) and of matrix elements of local operators, averaged over an ensemble of random circuits, making comparisons with the predictions of random matrix theory (RMT) and the eigenstate thermalisation hypothesis (ETH). We show that properties are dominated at long times by contributions to orbit sums in which each orbit is paired locally with a conjugate, as in the diagonal approximation, but that in large systems these contributions consist of many spatial domains, with distinct local pairings in neighbouring domains. The existence of these domains is reflected in deviations of the SFF from RMT predictions, and of matrix element correlations from ETH predictions; deviations of both kinds diverge with system size. We demonstrate that our physical picture of orbit-pairing domains has a precise correspondence in the spectral properties of a transfer matrix that acts in the space direction to generate the ensemble-averaged SFF. In addition, we find that domains of a second type control non-Gaussian fluctuations of the SFF. These domains are separated by walls which are related to the entanglement membrane, known to characterise the scrambling of quantum information.


I. INTRODUCTION
This paper is concerned with understanding generic features of spectral and eigenstate correlations for a class of ergodic many-body systems. Random matrix theory (RMT) and the eigenstate thermalisation hypothesis (ETH) provide a baseline description of these features [1][2][3][4][5][6][7][8][9][10]. Our objective is to identify when these approximations are accurate, and to characterise the behaviour when they break down, focussing on the regime of large times and distances for which one can expect a degree of universal behaviour. As a route to doing so, we start from expressions for physical properties in terms of sums over Feynman histories, or paths in Fock space, and aim to identify the dominant contributions to these sums.
Analogous approaches involving sums over histories have been applied very successfully to single-particle quantum systems in a variety of settings. In particular, the Gutzwiller trace formula provides a connection between the periodic orbits of a classical system and the spectrum of its quantum counterpart [3,11]. Contributions from pairs of orbits that carry opposite phases are dominant in the semiclassical limit, and the restriction to these pairs is known as the diagonal approximation [12]. RMT spectral correlations are then a consequence of the nature of periodic orbits in classical systems that are chaotic, and the diagonal approximation is the starting point for a systematic semiclassical expansion [13][14][15]. As a second example, in the theory of mesoscopic conductors the spectral and transport properties are expressed in terms of paths of electrons undergoing multiple scattering by impurities [16][17][18]. Pairs of paths with opposite phases, known as diffusons and Cooperons, survive a disorder average and determine long distance and low energy properties. They are also the basis for an expansion in inverse powers of the mean free path that captures weak localisation effects in disordered conductors. Here we will refer quite generally to the approximation of retaining only those pairs of paths with opposite phases as the diagonal approximation.
The many-body systems we discuss are quantum circuits. These offer a minimal description of time evolution in lattice models: each lattice site carries a 'spin' with Hilbert space dimension q, and pairs of sites are coupled by unitary gates during discrete time steps of the evolution. The examples we treat are one-dimensional brickwork models, in which every site is coupled alternately to its nearest neighbours on either side in successive time steps. In an obvious way, matrix elements of the evolution operator for a single time step define amplitudes for the steps of a path in Fock space. The transition amplitude between particular initial and final states under evolution over multiple time steps is a sum of contributions from multiple paths, labelled by states at intermediate times and weighted by products of amplitudes for the individual steps. These sums over Feynman histories will be at the centre of our discussion.
Following the spirit of RMT, many useful insights have been obtained recently by considering quantum circuits in which gates are drawn from a random distribution and properties of a circuit are averaged over an ensemble. In this way, two classes of model arise: random unitary circuits (RUCs) [19][20][21][22][23][24], in which gates are chosen independently at different time steps, and random Floquet circuits (RFCs) [25][26][27][28][29][30]. RFCs may equally be viewed as examples of kicked spin chains, which have been stud-ied extensively in their own right [31][32][33][34][35]. In RFCs the evolution operator is the same over each whole period, and is known as a Floquet operator. RUCs with gates chosen from the Haar distribution are particularly simple to analyse because the diagonal approximation is exact following an ensemble average. By contrast, RFCs, even with Haar-distributed gates, are not amenable to an exact analysis except in the large-q limit [26]. There is nevertheless a strong interest in understanding their properties, because time evolution in a RFC provides information on the spectral and eigenstate correlations of the evolution operator for a simple form of local manybody dynamics. This motivates the study of RFCs that we present here.
The spectral statistics of a Floquet operator are characterised by the spectral form factor (SFF). Denoting the Floquet operator by W and its t-th power for integer time t by W (t), the SFF is K(t) = |TrW (t)| 2 . Similarly the matrix elements of a local operator τ in the basis of eigenstates of W are characterised, in the off-diagonal case, by the autocorrelation function Tr[τ W (t)τ W † (t)] and in the diagonal case by |Tr[τ W (t)]| 2 . A key technical fact is that the ensemble average of the SFF can be generated using a transfer matrix that acts in the spatial direction. One of the central ideas that we present in this paper is that questions about the pairings of Feynman histories that contribute to the ensemble-averaged SFF and matrix element correlators can be rephrased as questions about the eigenvectors associated with the leading eigenvalues of this average transfer matrix. Similarly, deviations of the SFF from RMT predictions are controlled by the behaviour of these eigenvalues as a function of t. Although the dimension of this matrix grows very rapidly with t (for a brickwork RFC it is of dimension q 4t × q 4t ), we are able to probe its properties numerically for sufficiently large t that we believe we have established its asymptotic behaviour. On this basis we argue that the character of the relevant pairs of histories is determined by the relative size of t and of the system length L. In large systems the diagonal approximation becomes exact for times which are large but much smaller than the inverse level spacing. On the other hand, at sufficiently large L for fixed large t, the contributing Feynman histories are locally paired, with multiple domains in the orbit pairing, and distinct pairings in neighbouring domains.
This work builds on and complements other recent research in a number of ways. Most directly, the possibility that domains arise in the pairing of Feynman histories was shown for a model solvable in the large-q limit by one of the present authors and others [27]. The current paper shows how to formulate and test this idea in a generic setting, demonstrates that it has consequences beyond the behaviour of the SFF, and links it to the notion of the entanglement membrane, which characterises the scrambling of quantum information in chaotic quantum systems [36][37][38]. The idea of using a transfer matrix to generate the SFF and to calculate other quantities has been applied previously in several settings [28,29,34,35,[39][40][41][42][43][44]. This approach arose from discussions of periodic orbits in many-body systems [45], was developed as a method for treating kicked spin chains [34], and has been elaborated further as a way of accessing the semiclassical limit and making connections with periodic orbit theory [35,40]. Such a transfer matrix also forms the basis for the analysis of self-dual kicked spin chains [28,46], which display exact RMT behaviour of the SFF. Away from the self-dual point, recent work has investigated the evolution of this transfer matrix from ergodic to manybody localised behaviour [42]. The important distinction between those works and ours is that we are concerned with generic behaviour unrelated to self-duality, and at late times.
Both the RFCs we consider and the kicked spin chains investigated by other groups lead naturally to a description in terms of Feynman histories in Fock space. This perspective has also been adopted in studies of thermalisation, through the introduction of an influence matrix [43,44]. We note that there have been complementary efforts to study many-body quantum dynamics in the semiclassical limit by relating the orbits of classical models to properties of their quantum counterparts. These include studies of the spectral properties, many-body versions of coherent backscattering, and the out-of-time-order correlator [47][48][49][50].
The models we study here differ in two ways from standard random-matrix systems, such as the Gaussian and circular ensembles [1]. In our case the non-zero matrix elements are local in space. There is a long history of work on systems with two-body random interactions [51,52], but without the restriction of locality. The behaviour of the SFF in ensembles of that type has been examined from the current perspective in Ref. [53] and we comment on this further in Sec. VI.
The remainder of this paper is organised as follows. In Sec. II we first gather definitions of the models and physical quantities we consider and describe how sums over histories appear in our calculations. We then give a brief overview of our results. In Sec. III we show how to construct the transfer matrices generating the SFF of a RFC, Haar-average these transfer matrices, and then study the spectral properties of the average. Our analysis reveals a picture of local pairings of paths, in this case closed orbits in Fock space, and correspondingly deviations of the SFF from RMT which diverge in the thermodynamic limit. Through the transfer matrices we will also highlight a connection between the spectral statistics of chaotic many-body systems and their behaviour under local coupling to a bath. Then, in Sec. IV, we show that local orbit-pairing implies strong correlations between the diagonal matrix elements of local observables. Relative to the predictions of ETH these correlations grow without bound in the thermodynamic limit, and we verify this growth numerically. Following this in Sec. V we consider the statistical fluctuations of the SFF. By writing the higher moments of the SFF in terms of multiple sums over histories, we show that a distinct freedom in their local pairing gives rise to non-Gaussian statistics. We also highlight a connection with the entanglement membrane. In Sec. VI we extend aspects of our discussion to a class of models with gates whose distribution can be tuned continuously from Haar to the identity, and discuss how their behaviour differs away from the Haar case. In Sec. VII we provide a summary, discussion and outlook. Various technical details are described in a series of appendices.

II. OVERVIEW
Before discussing our work in detail, we first introduce the Floquet models and the physical quantities that we consider throughout the paper. Following this we discuss sums over histories and the diagonal approximation in the context of quantum circuits. We then give an overview of our main results.

A. Models and correlators
We focus on one-dimensional Floquet circuits with brickwork structure [26], illustrated in Fig. 1. The unitary evolution operator over integer time t is W (t) ≡ W t , where W = W 2 W 1 is the Floquet operator. With local Hilbert space dimension q, each of W 1 and W 2 are tensor products of two-site (q 2 × q 2 ) unitary matrices U x,x+1 coupling alternate pairs of neighbouring sites (x, x + 1). We label sites x = 0, 1 . . . (L − 1) for a system of length L, and so with Fock-space dimension q L . For periodic boundary conditions, which necessitates L even, the two half-steps W 1,2 are given by With open boundary conditions and L even we replace U L−1,0 with the q 2 × q 2 identity matrix. With open boundary conditions and odd L we instead have where now q × q identity matrices act on site x = (L − 1) in W 1 and x = 0 in W 2 . With the exception of Sec. VI we are concerned with circuits constructed from Haarrandom gates, or Haar-RFCs. For all numerical investigations we focus on a local Hilbert space dimension q = 2 (a form of kicked spin-1 2 chain). To demonstrate that our choice of Haar-random gates is not crucial to the results we obtain, in Appendix G we compare with behaviour for a kicked Heisenberg model.
The central quantity considered in this work is the spectral form factor (SFF), which probes correlations in the level density. The SFF is defined, for integer t, by Here θ n is the quasienergy of the Floquet operator W associated to eigenstate |n , W |n = e iθn |n . RMT, in the standard Wigner-Dyson sense [1], will be a key reference point. For W a Haar-random N × N unitary matrix, TrW (t) is normally and isotropically distributed in the complex plane in the limit of large N [54]. In RMT for all N we have the average SFF for integer t An important timescale is the Heisenberg time t H = N , set by the mean level spacing 2π/N . For circuit models N = q L , and except where stated explicitly we consider times below t H .
For many examples of chaotic quantum systems it is known that RMT behaviour of the SFF sets in only beyond a characteristic timescale. In some circumstances this scale reflects specific microscopic features of the system, but for spatially extended models it may diverge with system size and arise from physical processes that are common to a class of systems. Diffusive mesoscopic conductors provide a key example, for which this timescale is known as the Thouless time, and is given in terms of the diffusion constant D and the linear system size L by t Th = L 2 /D. For the chaotic many-body systems we consider, we will also refer to the timescale at which RMT behaviour of the SFF sets in as the Thouless time [27,[55][56][57][58] (see also Ref. [53]). We note however that there are a number of alternative definitions, based on the behaviour of local observables [59], the sensitivity to boundary conditions [60], and on many-body return probabilities [61].
The characterisation of the behaviour of a system involves the matrix elements of local observables, and here Two of the diagonal orbit pairings that contribute to the spectral form factor K(t) = |TrW (t)| 2 . The outer circles represent the forward orbit, and the inner the backward, corresponding to terms in the sum-over-histories representations of TrW (t) and TrW * (t), respectively. The dashed lines represent the pairing of forward and backward orbits. Diagonal pairings have ar = a * r+s , with addition defined modulo t. On the left and right we show the s = 0 and s = 1 pairings, respectively.
the predictions of the ETH provide a useful reference. In Floquet systems the ETH is as follows: for a set of eigenstates drawn from a sufficiently narrow window of quasienergies, the statistical properties of the matrix elements of local observables are as for random vectors.
The diagonal matrix elements of the operator τ at site x can be written in terms of the reduced density matrices at this site, ρ x (n) = Tr x |n n| for the eigenstate |n , in the form Tr[τ ρ x (n)]. Here and throughout this paper, we use Tr x to denote a trace over site x and Tr x to denote a trace over its complement, the other (L − 1) sites. A correlator between the eigenstates |n and |m can then be defined as Tr x [τ ρ x (n)]Tr x [τ ρ x (m)]. To avoid referencing a specific observable it is useful to average over a complete orthonormal set of operators (see Sec. IV). This leads us to consider the correlator of reduced density matrices Tr x [ρ x (n)ρ x (m)]. It is convenient to work in the time domain, and so we define the reduced form factor (RFF) at site x, We will see also that a straightforward generalisation of this, describes spectral structure in the off-diagonal matrix elements of operators local to site x.

B. Histories in circuits
We now discuss sums over histories in quantum circuits, and how they appear in the quantities of interest.
Consider the probability for a Floquet system to evolve from an initial state |a 0 to final state |a t over time t, This is a discrete sum over all pairs of paths between state a 0 and a t . The forward (a 0 a 1 . . . a t ) and backward (a 0 a * 1 . . . a t ) paths are labelled by integers a r , a * r = 0 . . . (N − 1). The transition amplitudes are simply the matrix elements of the Floquet operator.
The weights of the N t−1 diagonal pairs of paths, with a r = a * r for all r, give real non-negative contributions to Eq. (7). If the weights of the N t−1 (N t−1 −1) off-diagonal pairs behave like a set of independent random variables, their sum gives a real contribution of the same order as the sum of the diagonal pairs. The off-diagonal contributions may, however, vanish after a suitable average. We then arrive at the diagonal approximation. The dynamics within the diagonal approximation is described by the diagonal propagator, which could also be thought of as the analogue of a diffuson in Fock space, For example, provided the ensemble average can be performed independently for each step, the transition probability Eq. (7) becomes P t ata0 (see Sec. VI for a discussion of the distinction between P t and P t ). For definiteness in the following we use the term diagonal approximation to refer to calculations based on the average diagonal propagator P. Unitarity constrains P to be a doubly stochastic matrix (the sum along any row or column is unity) and so it has a leading eigenvalue of unity.
To construct the SFF in the diagonal approximation we write TrW (t) as a sum over all closed paths of t steps, We will refer to these closed paths as orbits. The SFF is a sum over all pairs of forward and backward orbits, The amplitude of an orbit (a 0 , a 1 . . . a t ) is invariant under a cyclic permutation. Consequently, in a system without time-reversal symmetry, for a typical forward orbit in Eq. (10) there are t backward orbits with the complex conjugate amplitude. In Fig. 2 we illustrate two of these diagonal orbit pairings. The diagonal approximation to the average SFF is For a propagator with just one unit-modulus eigenvalue, TrP t → 1 at late times. In the regime t < t H , the diagonal approximation to the SFF therefore approaches the RMT result of Eq. (4). As a concrete example consider again the SFF within RMT. For W Haar-random, and in the limit of large N , the ensemble average of Eq. (10) is [62,63] for nonzero integer t N . The sum over s = 0 . . . (t − 1) is the sum over t diagonal orbit pairings. Here the diagonal approximation coincides with the exact result.
In circuit models we are concerned with sums over paths in Fock space. Choosing our Fock-space basis states |a r to be product states, the sum in for example Eq. (9) can be recast as a multiple sum over the local orbits of individual sites. TrW (t) and the SFF can then be written in terms of transfer matrices acting on these orbits (see Sec. III for details). For the Haar-RFCs we study, all choices of site basis are statistically equivalent.
These transfer matrices are associated with individual gates U x,x+1 . Writing the transfer matrices generating K(t) as T x,x+1 (t), with open boundary conditions we find for example Here B L | and |B R are vectors encoding the boundary conditions at the left and right ends of the system, respectively. For independently and identically distributed gates the average SFF is determined by powers of the average transfer matrix T (t), Fixing t and taking the limit of large L, the average SFF is dominated by contributions from the leading eigenvalues of T (t). We show, for a system without time-reversal symmetry, that there are exactly t leading eigenvalues, and we denote these by λ(ω, t), where ω is a t-valued symmetry label associated with time-periodicity.

C. Results
In this work we study behaviour beyond RMT in chaotic many-body Floquet systems with local interactions. For Haar-RFCs the deviations of the SFF from RMT arise from a particular class of off-diagonal pairings of paths.
Simplicity, in the form of the diagonal approximation or something beyond, is only to be expected after a degree of averaging. Moreover, the SFF is not self-averaging [64], and typically exhibits large system-dependent fluctuations around the RMT result [54]. Averaging can in principle be approached in various ways. For a specific model or a specific realisation drawn from an ensemble, one can average the SFF over a time window. This window should be narrow on the scale set by its mid-point, to avoid distortion, but wide enough to contain many Floquet periods, to ensure efficient averaging [64]. Alternatively, one can average over a small region in the space of possible models. In this case the average may be representative of all systems in this region, but systemdependent fluctuations will only be washed out at late times. A final alternative, which is the one we follow, is to average over a wide variety of systems: in our case the Haar distribution for gates. In practice however quite limited averaging is sufficient. In Fig. 3 we show that an average over just one gate is enough to dramatically suppress fluctuations of the SFF.
In Sec. III we show that the spectral structure in Haar-RFCs, with local Hilbert space dimension q = 2, is not captured by the diagonal approximation. To understand the form of the SFF we develop the picture of local orbit pairing which first appeared in the large-q limit in Ref. [27], demonstrating here its appearance in a generic setting. To do this we show how to express the SFF of a Haar-RFC in terms of transfer matrices T (t) acting in the space direction. For independently and identically distributed gates the average SFF K(t) is then expressible in terms of powers of the average transfer matrix T (t). Through this average transfer matrix the accuracy of the diagonal approximation, as well as the corrections arising in large systems, acquire a sharply defined meaning in the language of local orbit pairing.
The transfer matrices are too large to study directly in the regimes of interest. To probe the spectral properties of T (t) we instead impose a variety of boundary conditions on the model, thereby coupling to the eigenvectors in controlled ways. By analysing the length-scaling of the SFF and related objects, we separate the leading eigenvalues of T (t) according to their symmetry sector and determine their magnitudes. In practice our approach requires only very small systems of L ≤ 8 sites for the times of interest, and our results allow us to reconstruct the SFF for arbitrarily large L. By studying the corre-sponding eigenvectors we then directly probe the local orbit pairing.
In Fig. 4 we illustrate the different regimes of the many-body spectrum for the Haar-RFC and relate these regimes to the transfer matrix spectrum. At fixed L, increasing t brings us into the diagonal regime (for t < t H = q L ). Here K(t) is dominated by the sum over the t diagonal orbit pairings, each contributing unity, so we find the RMT result K(t) = t. We will refer to these contributions as the global diagonal orbit pairings. By contrast, taking the limit of large L at fixed t the SFF is dominated by orbits which are diagonally paired only locally, with distinct diagonal pairings in neighbouring domains. These two regimes can be understood by considering only the leading eigenvalues λ(ω, t) of the average transfer matrix, of which there are exactly t. The deviations of these eigenvalues from unity control the contributions of domain walls to the SFF. For example, increasing t these eigenvalues approach unity and we move from a picture of local orbit pairing to one of global orbit pairing. Conversely, taking the limit of large L at fixed t, the largest of these close-lying eigenvalues dominates. Small deviations of the eigenvalues from unity are then responsible for large deviations of the SFF from RMT. We show that under certain assumptions the timescale for the crossover between these two regimes, which we refer to as the Thouless time, scales logarithmically with system size L. A third regime is entered on increasing t beyond the Heisenberg time. There the subleading eigenvalues of T (t) control the SFF.
A natural question is whether this picture of local orbit pairing has implications for local observables and eigenstate correlations, and we investigate this in Sec. IV. The trace structure of the RFF R x (t) in Eq. (5) is, away from the site x, identical to that of the SFF, and we show that this implies exponential growth of R x (t) with increasing L. Remarkably, on transforming this behaviour from the time domain to that of quasienergies we find that the deviations from ETH are most prominent on small quasienergy-scales, where the ETH is conventionally expected to be most accurate. We find enhancements of the correlations relative to RMT by two orders of magnitude even in systems of L ≤ 14 sites, so our results represent a substantial correction to the ETH for one-dimensional chaotic Floquet systems.
We investigate the statistical fluctuations of the SFF in Sec. V. Focussing on its second moment we highlight the presence of another kind of domain wall, in the pairing of the multiple copies of orbits which appear. While these are distinct from the domain walls in the orbit pairing considered in earlier sections, their role also becomes more prominent with increasing system size. These domain walls in the pairing of copies of orbits are closely related to entanglement membranes. Our results here also provide an understanding of the striking effect of the single-gate average in Fig. 3(b).
In Sec. VI we discuss deviations of the SFF from RMT which arise within the diagonal approximation. These deviations result from the subleading eigenvalues of the diagonal propagator [25,56,58,65]. We consider this effect both for Haar-RFCs and more generally.

III. LOCAL ORBIT PAIRING
The diagonal approximation to sums over pairs of orbits is blind to locality. In this section, which contains the central arguments of the paper, we demonstrate how the picture of paired orbits must be extended in onedimensional many-body systems with local interactions. First in Sec. III A we show that the diagonal approximation fails to describe the spectral correlations in Haar-RFCs. We then discuss the model of Ref. [27], where the breakdown of the diagonal approximation is evident even in the large-q limit. There an exact treatment is possible, and this reveals the picture of local orbit pairing. Motivated by this, in Sec. III B we express the SFF of a brickwork model (with arbitrary q) in terms of transfer matrices for the local orbits which act in the spatial direction. On averaging these transfer matrices the picture of local orbit pairing emerges.
Through Secs. III C and III D we extract information on the average transfer matrix by imposing various boundary conditions on the orbits. This allows us to determine the leading eigenvalues, as well as the relationship between the corresponding eigenvectors and the local orbit pairing. As a test of our results, we show that this information is sufficient to reconstruct the average SFF accurately in large systems (see Fig. 16). In Sec. III E we further explore the connection between the eigenvectors of the average transfer matrix and the local orbit pairing, and then in Sec. III F discuss the behaviour of individual circuit realisations.

A. Breakdown of the diagonal approximation
First we demonstrate the necessity of moving beyond the diagonal approximation. In particular, we will show that for Haar-RFCs K(t) = t within the diagonal approximation. This is also the large-q result [26], and should be compared with numerical results for q = 2 in Fig. 5. The diagonal propagator for brickwork models takes the form where all subscripts are many-body indices (taking q L values). Each independent gate appears once in P, as does its conjugate. P is given by averaging over the gates. For example, the average over the gate U 0,1 acting on sites 0 and 1 in the first half-step W 1 is where we have made the single-site indices explicit. The result is the same for gates acting in the second half-step.
Multiplying the expressions Eq. (16) for each independent gate, and summing over the internal c, c * indices, we find the average propagator a matrix with constant entries. The leading eigenvalue of P is unity and all others are zero. As a result, the diagonal approximation to the SFF K(t) = t as in RMT. By contrast in Fig. 5 there are obvious deviations of the SFF from RMT. These deviations grow with increasing system size, and furthermore are significantly larger with open boundary conditions than with periodic. Through this section we set out to understand this behaviour. The diagonal approximation to the SFF involves an overall factor of t, the number of diagonal orbit pairings. The breakdown of the diagonal approximation in a chaotic Floquet system with local interactions was demonstrated in Ref. [27]; there a calculation of the SFF in the large-q limit revealed deviations from RMT which grow exponentially with increasing L. Here we briefly review certain aspects of that work.
The Floquet operator W = W 2 W 1 of Ref. [27] is made up of two half-steps. The first, W 1 , consists of independent q × q Haar-random gates u x acting on individual sites of a chain, W 1 = u 0 ⊗ . . . ⊗ u L−1 , for L sites. The second, W 2 , consists of weak interactions between adjacent sites with coupling strength ε.
In K(t) each of the q × q gates u x appears t times in TrW (t) and its conjugate appears t times in [TrW (t)] * . To calculate K(t) we must average independently over each u x , and this gives a sum over orbit pairings s x = 0 . . . (t − 1) as in Eq. (12). This immediately promotes the orbit pairing to a local degree of freedom. K(t) is then given by a sum over the orbit pairings s x at each FIG. 5. Average SFF K(t) in the q = 2 brickwork model with (a) periodic and (b) open boundary conditions. The system size 3 ≤ L ≤ 12 is shown on the legend. L is necessarily even with periodic boundary conditions. In the diagonal regime K(t) t, and beyond the Heisenberg time tH = q L we have , and takes the form of a partition function Here and throughout this paper, we will refer to a configuration s x = s x+1 as a domain wall. The statistical weights of domain walls are suppressed by factors e −εt ; the coupling strength ε appears here as the domain wall line tension. In this way K(t) has been expressed in terms of the transfer matrix of a ferromagnetic t-state Potts model. At late times the domain walls are suppressed, and we recover the RMT result K(t) = t, corresponding to a sum over t global diagonal orbit pairings.
The first corrections to RMT with open boundary conditions come from configurations with one domain wall, and with periodic boundary conditions they come from configurations with two domain walls, The Although these results were derived in the limit of large q, in this work we show that the picture of local orbit pairing is more general. For example, it describes the spectral statistics of Haar-RFCs at small q. Indeed, the deviations from RMT displayed in Fig. 5, including their dependence on both the system size and on the boundary conditions, illustrate the phenomenology of domain walls in the orbit pairing. Specifically, for fixed t the deviations from RMT behaviour grow with L and are larger with open than with periodic conditions.

B. Transfer matrix for orbits
Here we show how to construct the transfer matrices generating the SFF of a brickwork model, discuss their spectral properties, and connect the behaviour of their average to the behaviour observed in Fig. 5. Similar approaches have appeared in the study of kicked Ising models [34]. By writing the average spectral form factor of the brickwork model in terms of a transfer matrix, we make the local orbit pairing degrees of freedom explicit.
Throughout this section we will work with single-site orbits, and so here we introduce some notation. Two unitary gates act on each site during a time step, so at time t each single-site orbit is a string of 2t integers, each taking a value 0 . . . (q − 1). We denote the forward orbits appearing in TrW (t) by (a 0 b 0 . . . a t−1 b t−1 ). For integer r, a r represents the state of the site at time r, and b r represents the state of the site at time (r + 1/2), in the middle of the step. We denote the backward orbits appearing in [TrW (t)] * by strings (a * 0 b * 0 . . . a * t−1 b * t−1 ). In Fig. 6 we illustrate the SFF of a brickwork model, as well as a pair of forward and backward single-site orbits. Dashed lines exiting at the top of the figure are connected to those entering from below, giving independent traces TrW (t) and TrW * (t). Right: segment of the doubled time evolution, highlighting notation for forward (a0b0 . . .) and backward (a * 0 b * 0 . . .) orbits

Construction of the transfer matrix
Consider a unitary gate U x,x+1 acting on sites (x, x+1) in the first half of the Floquet step and so with x even (see Fig. 1). This matrix evolves the state of the two sites from a time r to (r + 1 2 ), with r integer. We can also think of U x,x+1 as a non-unitary matrix acting on the orbit of site (x + 1) at the times (r, r + 1 2 ), and refer to this matrix asŨ x,x+1 . Writing the components of U x,x+1 as [U x,x+1 ] brb r ,ara r , where unprimed indices a r , b r correspond to site x and primed a r , b r to site (x + 1), the components ofŨ x,x+1 are [Ũ x,x+1 ] arbr,a r b r = [U x,x+1 ] brb r ,ara r , as illustrated in Fig. 7.
In TrW (t) each gate U x,x+1 appears t times. By taking a tensor product of the t copies ofŨ x,x+1 , we create a q 2t × q 2t matrixŨ ⊗t x,x+1 acting on the entire forward orbit of site (x + 1). To construct TrW (t) from the standard matrix multiplication of these operators, we introduce an orthogonal matrix S which acts on orbits as a translation of one half-step in time, The matrices S andŨ ⊗t are illustrated in Fig. 8. With periodic boundary conditions where L is the number of sites, which is necessarily even with these boundary conditions. Here tr is a trace over Components U br b r ,ar a r of the q 2 × q 2 unitary matrix U , and componentsŨ ar br ,a r b r of the q 2 × q 2 non-unitary matrixŨ . Time runs vertically and space horizontally. Here, U acts in the first half-step of the Floquet operator, so it describes evolution from time step r to (r + 1 2 ). Where U acts on the state of a pair of sites (x, x + 1) at time r,Ũ acts on the state of site (x + 1) at times (r, r + 1 2 ).
single-site orbits, to be distinguished from Tr, a trace over many-body states. Since the evolution operator is periodic we have [S 2 ,Ũ ⊗t ] = 0, where S 2 translates an orbit one full-step in time. Clearly S 2t = 1, where 1 is the q 2t × q 2t identity acting in the space of single-site orbits. Using these properties we can write TrW (t) in terms of one type of transfer matrix, SŨ ⊗t , A similar expression can be derived for open boundary conditions, where the boundary state |B has components and is invariant under time translation by a full step, The transfer matrices which generate the SFF K(t) = |TrW (t)| 2 are straightforwardly expressed as tensor products of the transfer matrices which generate TrW (t) and [TrW (t)] * , and are also illustrated in Fig. 8. The transfer matrix T x,x+1 (t) acts on the product space of forward and backward orbits at site (x + 1), which has dimension q 4t . These matrices have some obvious symmetries associated with translation in the time direction, which we now discuss (see also [42] for a related treatment in the context of the kicked Ising model).
The transfer matrices T x,x+1 (t) commute with the translation operators for a full time step of either of the forward and backward orbits, Since (S 2 ) t = 1, T x,x+1 can be block-diagonalised into t 2 sectors, having eigenvalues e iω+ under S 2 ⊗ 1 and e iω− under 1 ⊗ S 2 . The frequencies are ω ± = 2πn/t for n = 0 . . . (t − 1) defined modulo t. The transfer matrices for each bond (x, x + 1) are different. Our focus, however, is on the full ensemble average. Following such an average all the transfer matrices are the same.

Average transfer matrix
The notion of orbit pairing takes on a concrete meaning in the ensemble average of the SFF. For identically and independently distributed gates, K(t) is determined by powers of a single averaged transfer matrix, T (t). From here on, T (t) will be at the centre of our discussion.
To calculate T (t) we must Haar-average the tensor productŨ ⊗t ⊗ [Ũ ⊗t ] * . Writing the first indices ofŨ ⊗t and [Ũ ⊗t ] * as orbits , respectively, the only non-vanishing matrix elements of their average tensor product have , where σ and τ denote permutations of t objects [62,63]. The same is true of the second indices. To arrive at a succinct expression for the average transfer matrix, it is convenient to introduce vectors |σ, τ in the product space of forward and backward orbits which have nonzero entries only where indices are paired in this way: Note that the vectors |σ, τ are neither normalised nor orthogonal. The average ofŨ ⊗t ⊗ [Ũ ⊗t ] * can be expressed in terms of the vectors |σ, τ and the (q 2 -dependent) Weingarten functions Wg(στ −1 ), here taking the composed permutation στ −1 as an argument. From the Haar average over the gate we find the average transfer matrix where the sum is over all pairs of permutations of t objects. We have also introduced the doubled half-step time-translation operator S = S ⊗ S. Only t terms in the sum in Eq. (29) contribute in the large-q limit, and these are essentially the local orbit pairings discussed in Sec. III A. We elaborate on this in Appendix C.
As discussed above, the transfer matrices have a block structure associated with symmetry under time translation. Focussing on the average T (t), we write the left and right eigenvectors in the block ω + , ω − as ω + ω − α L ; t| and |ω + ω − α R ; t , respectively, and the corresponding eigenvalues as λ(ω + ω − α; t). Here α = 0, 1, . . . label eigenstates in descending order of magnitude, so λ(ω + ω − , 0; t) is the leading eigenvalue in the block ω + , ω − . The spectral decomposition of the average transfer matrix T (t) is then We are now in a position to write down expressions for the average SFF in terms of the average transfer matrix. With periodic boundary conditions, where L is necessarily even, For the case of open boundary conditions we define the doubled boundary vectors |B R ≡ |B ⊗ |B and B L | ≡ B R | S T , for the right and left ends of the chain, respectively. These vectors are in the ω + = ω − = 0 sector. From Eq. (24) we then have and this is illustrated in Fig. 9. An immediate question is how, for general L, RMT level statistics emerge from the average transfer matrix T (t) at late times. We anticipate that in the diagonal regime the result K(t) = t is expressible as a sum over the t global diagonal orbit pairings, each contributing unity. In the language of the average transfer matrix T (t) this corresponds to having t leading eigenvalues which at late times approach unity. Different behaviour sets in beyond t H = q L , where the average SFF plateaus at K(t) = q L . To see how this behaviour could arise, note that the number of nonzero subleading eigenvalues of T (t) grows very rapidly with t, but that the contributions of these eigenvalues to the SFF are suppressed for large L. In Appendix D we show, based on exact diagonalisation of T (t) for t ≤ 5, that the proliferation of subleading eigenvalues with increasing t is responsible for the plateau.
Our analysis is focused on the regime 1 t t H for large L, and therefore on the late-time behaviour of the t leading eigenvalues of T (t). General features then emerge because the only microscopic timescale in T (t) is q 2 , the Heisenberg time for a single gate. As we show in Sec. III D, each of the sectors ω = ω + = −ω − contains one of the t leading eigenvalues, and since our attention is effectively limited to these and the corresponding eigenvectors, it is convenient to introduce the shorthand notation In the regime of interest the transfer matrix is too large to compute directly, so in the next two sections we study the length-scaling of the average SFF and related objects with a variety of boundary conditions. Through this we will determine λ(ω, t) as well as certain properties of the corresponding eigenvectors.

C. Open boundary conditions
In this section we study the length-scaling of the average SFF K(t) with open boundary conditions. From Eq. (32) this probes the ω + = ω − = 0 block of the average transfer matrix. Based on the discussion in Sec. III A we anticipate that the leading correction to the SFF relative to RMT can be understood in terms of a domain wall, and so here we aim also to infer an effective domain wall tension.
If the leading eigenvalue in the ω + = ω − = 0 sector dominates the SFF in Eq. (32), we have (34) where the ellipses indicate contributions from subleading eigenvalues. By following the L-dependence of K(t) for each time t we extract the leading eigenvalue λ(0, t) and the overlap of the corresponding eigenvector with the boundary states, B L |0, t; R 0, t; L|B R . The results are shown in Fig. 10, and details of the analysis are presented in Appendix F. We find that λ(0, t) approaches unity from above at large t. We find also that the overlap, shown in the inset, approaches t at late times. This is exactly the value expected if the leading eigenvector of T (t) represents a locally diagonal orbit pairing, as we discuss in Sec. III E.
Expanding the average SFF around the RMT result, at large t we then have where we have defined δλ(0, t) = λ(0, t) − 1. We can interpret the deviation from RMT as a domain wall contribution in analogy with Eq. (19). Here (L − 1) is the translational entropy, and δλ(0, t) in Eq. (35) plays the role of (t − 1)e −εt , the sum of statistical weights of the (t − 1) different domain walls. This interpretation motivates the definition of an effective domain wall tension FIG. 11. Evolution of the effective tension ε eff (t) defined in Eq. (36), calculated using numerical results for the leading eigenvalue λ(0, t) in Fig. 10.
In the model of Ref. [27] all domain walls have the same tension ε, so there ε eff = ε, which is also timeindependent. This is not the case in general [27,66]. We show in Fig. 11 that the effective domain wall tension in the q = 2 Haar-RFC decreases monotonically with time.
To understand this decrease, suppose that a domain wall can be located on a particular bond with corresponding gate U , and that we can associate a tension ε(U ) with this gate. The average of the statistical weight e −ε(U )t over U is at late times increasingly dominated by those gates with small ε(U ). Consequently, it is natural for the effective tension derived from this average to decrease with time.
We can relate this behaviour to the Thouless time t Th as follows. From Eqs. (35) and (36) we see that the crossover to RMT occurs when the decreasing statistical weight e −ε eff (t)t associated with a domain wall overwhelms the translational entropy (L − 1). If the effective domain wall tension approaches a nonzero constant ε eff (∞) at late times, the timescale for the crossover of the SFF to RMT behaviour is t Th = ln L/ε eff (∞).
In this section we have determined one of the t leading eigenvalues of the average transfer matrix. Next we will determine the other (t − 1) and provide further evidence for domain walls in the orbit pairing.

D. Twisted boundary conditions
In order to access t different symmetry blocks of T (t) here we impose local diagonal orbit pairings at the two ends of the system. This approach allows us to force domain walls into the many-body orbit pairing, and gives us information on λ(ω, t), ω, t; L| and |ω, t; R for all ω. By determining all of these leading eigenvalues, we will show that they alone are sufficient to calculate the SFF with periodic boundary conditions to a remarkable degree of accuracy.

Domain walls
The local diagonal orbit pairings are represented by vectors |s in the product space of forward and backward orbits. The pairing |s is a normalised member of the set of vectors |σ, τ introduced in Eq. (28), with σ = τ the permutation mapping r → (r + s) modulo t, for s = 0, 1 . . . (t − 1). The components are Imposing the orbit pairing s L on the left and s R on the right of an L-site system we find For s L = 0 and s R = 1 these boundary conditions correspond to the two pairings shown in Fig. 2. We illustrate Z(s, t) in Fig. 12. For s L = s R we force equal-time pairings, and for s L = s R we force a domain wall into the orbit pairing. In practice we impose these boundary conditions using a Monte-Carlo method, and give details on the implementation in Appendix B. For the particular case of Z(0, t) these boundary conditions can also be thought of as local couplings to Markovian baths. We discuss this further in Sec. III F. Via the ensemble average, we probe the average transfer matrix T (t). Suppose that the RMT behaviour K(t) = t at late times arises from the sum over the t possible global diagonal orbit pairings. The boundary conditions imposed for Z(0, t) are compatible with one of these pairings, whereas those for Z(s = 0, t) are not compatible with any of them. We therefore expect that Z(0, t) → 1 and Z(s = 0, t) → 0 with increasing t, on the grounds that contributions from pairs of many-body orbits with domain walls should vanish at late times. We show that this is indeed the case in Fig. 13 for s = 0 and s = 1 for various L, and in  Fig. 13 is then understood as a consequence of this entropy, and furthermore demonstrates the failure of the diagonal approximation. In principle, and in contrast to Ref. [27], we expect the decay rate of Z(s = 0, t) to depend on s (see Appendix C), and behaviour of this kind is evident in Fig. 14(a).

Leading eigenvalues
Using our numerics on Z(s, t) we can extract the leading eigenvalues λ(ω, t). The local diagonal pairings |s defined in Eq. (37) are eigenvectors of S 2 with unit eigenvalue, so they are linear combinations of vectors in sectors with ω + = −ω − . We define their Fourier transform in FIG. 14. Z(s, t) and Z(ω, t) for q = 2 and L = 8. (a) Z(s, t), with s on the legend. Z(0, t) tends to unity and Z(s = 0, t) decays approximately exponentially as in Fig. 13. (b) Z(ω, t), defined in Eq. (41) as the Fourier components of Z(s, t) with respect to s; the legend shows n = ωt/2π. sector Using Eq. (40) and the spectral decomposition of T (t) in Eq. (30) we have the Fourier transform of Z(s, t), For fixed t and large L, Z(ω, t) = λ L−1 (ω, t) ω|ω, t; R ω, t; L|ω + . . . , (42) where the ellipses represent the contributions of subleading eigenvalues. We show Z(ω, t) for L = 8 and various ω as a function of t in Fig. 14(b). The qualitative behaviour of the leading eigenvalues and corresponding eigenvectors at large t can be deduced as follows. First note that Z(0, t) → 1 and Z(s = 0, t) → 0 at large t for all L. Substituting this behaviour into the first line of Eq. (41), we find Z(ω, t) 1 for large t, for all L and ω. From Eq. (42), this implies that the leading eigenvalues λ(ω, t) 1, and the overlaps ω|ω, t; R ω, t; L|ω 1, at large t. Moving beyond this qualitative analysis, the leading eigenvalues λ(ω, t) and the overlaps ω|ω, t; R ω, t; L|ω can be extracted from the length-scaling of Z(ω, t), and the results are shown in Figs. 15(a) and (b), respectively. We give details of the analysis in Appendix F, and repeat it for another model in the same symmetry class in Appendix G. For comparison with length-scaling of Z(ω, t) we show (i) the exact results for t ≤ 5, determined from the exact diagonalisation of T (t) (see Appendix D), and (ii) the calculation of λ(0, t) from Fig. 10, and find excellent agreement between the different approaches. We see that λ(0, t) ≥ 1 and λ(ω = 0, t) ≤ 1, with all t of the leading eigenvalues approaching unity at late times. In Fig. 15(b), ω|ω, t; R ω, t; L|ω also approaches unity, and we discuss this in Sec. III E.
We can conduct a stringent test of the accuracy of the eigenvalues that we have extracted in Fig. 15(a), and of the adequacy of focussing on only the t leading eigenvalues. To do so we use our results to reconstruct the SFF with periodic boundary conditions, and the results are shown in Fig. 16. We find remarkable agreement with exact diagonalisation from Fig. 5(a). To demonstrate the power of the transfer matrix approach, and the exponential growth of deviations from RMT, we also extrapolate to larger systems than are directly accessible. For fixed t, our restriction to only the leading eigenvalues of T (t) is exact in the large-L limit.
Within the framework of the average transfer matrix the mechanisms giving rise to the RMT result K(t) = t in the diagonal regime are different with periodic and open boundary conditions. In the periodic case the t leading eigenvalues each contribute unity. In the open case just one of these eigenvalues, in the ω = 0 sector, contributes to K(t), but this contribution is enhanced by a factor B L |0, t; R 0, t; L|B R t.

Domain wall tensions
Having established the approximate equality Fig. 16, we will now relate the average SFF with periodic boundary conditions to properties of the domain walls studied in for example Fig. 13. At late times and for a nonzero twist s, the Fourier transform of Eq. (41) gives Here we have written δλ(ω, t) = λ(ω, t) − 1 and used ω|ω, t; R ω, t; L|ω 1. The ellipses represent terms higher order in δλ(ω, t) and the contributions from subleading eigenvalues. With no twist we instead have Z(0, t) 1 from Fig. 13. Eq. (43) and the observed decay of Z(s, t) in Fig. 13 motivate the definition of the domain wall tensions ε(s, t) through Since Z(s = 0, t) = 1 for large t, independently of L, we see that ω δλ(ω, t) = 0. The statistical weight e −ε eff t defining the effective tension ε eff in Eq. (36) is then simply the average over s = 0 of the statistical weights e −ε(s,t)t .
To write an expression for the SFF with periodic boundary conditions in terms of the domain wall tensions we invert Eq. (44) and use K(t) ω λ L (ω, t). Expanding around the RMT result we find to second order in δλ(ω, t). We see then that the first correction to RMT with periodic boundary conditions can be interpreted as arising from pairs of many-body orbits with two domain walls, as opposed to pairs with just one domain wall in the case of open boundary conditions.

E. Pairing domains
In this section we investigate the character of the orbitpairing domains. We probe the space-time-local pairing of orbits across domain walls and within the different domains. Our tools are the twisted boundary conditions from Sec. III D, and a suitable correlator defined below.
Our results on the leading eigenvectors of T (t) in Figs. 10 and 15(b) highlight a close connection with the vectors |ω , and therefore with the local diagonal orbit pairings |s . In particular, at late times we find B L |0, t; R 0, t; L|B R t, to be compared with B L |ω = 0 ω = 0|B R = t. Additionally, we have ω|ω, t; R ω, t; L|ω 1. Since |ω, t; R and ω, t; L| are vectors in a space of very high dimension for large t, it is quite remarkable to find that they have such large overlap with |ω , and even more striking that this overlap approaches unity at large t. However, because the left and right eigenvectors of T (t) form a biorthogonal set as opposed to being orthonormal, this does not imply equality of the leading eigenvectors and the |ω vectors.
To probe pairing beyond these overlaps, we construct a correlator as follows. Writing the forward orbit of site x as (a 0 b 0 . . . a t−1 b t−1 ) and the backward orbit as (a * 0 b * 0 . . . a * t−1 b * t−1 ), our space-time-local pairing correlator for the step r will be chosen to take the value 1 if a r = a * r and −1/q otherwise. In this way the correlator vanishes if a r and a * r are uncorrelated q-valued random numbers. To calculate this correlator we introduce the operator C(r), acting in the product space of forward and backward single-site orbits, Imposing the local diagonal orbit pairing s L on the site x = 0, and s R on site (L − 1), the average correlator on a site with x odd is given by inserting C(r) into a product of transfer matrices in Eq. (39). If we probe the orbit pairing of a site with x even we redefine C → SCS T . For an individual realisation our correlator for the site x is Here we have inserted C between the transfer matrices T x−1,x (t) and T x,x+1 (t) in the numerator, and have omitted the argument of C which is in this case arbitrary. The denominator is Z(s R − s L , t). We are interested in average properties of the transfer matrix, and although both the numerator and denominator are nonzero and positive after averaging, this is not the case in all individual realisations. Consequently it is not possible to average the right-hand side of Eq. (47) directly. We therefore define our average correlator C(s L , s R , t, x) as the ratio of the average numerator to the average denominator, Consider imposing a particular pairing domain using the boundary conditions s L = s R = s, as in Z(0, t). At late times if the forward and backward many-body orbits are diagonally paired we expect C(s, s, t, x) to be equal to unity if s = 0 and equal to zero if s = 0. At early times, based on the results of Fig. 13, we anticipate domain wall contributions. The corresponding multidomain orbit pairing configurations suppress C(0, 0, t, x) and enhance C(s, s, t, x) for s = 0. In Fig. 17(a) we show C(s, s, t, x) calculated for various x and for s = 0 and 1 in a system of L = 8 sites, and find exactly the behaviour expected.
To study the correlations in the pairing across a domain wall, we impose s L = 0 and s R = s = 0 as in Z(s, t). In this case, as we increase x and thereby sweep across the domain wall, we expect our correlator C(0, s, t, x) to decrease toward zero, and in Fig. 17(b) we confirm that this is indeed the case. The results of Fig. 17 demonstrate quite directly the existence of domains in the orbit pairing.

F. Bath interpretation of Z(0, t)
It turns out that the boundary conditions used to define Z(0, t) also arise for a system coupled at its ends to Markovian baths. This fact highlights a connection between thermalisation and, through the transfer matrices, the spectral statistics.
Consider first evolution under only the Floquet operator W .
The SFF can be written K(t) = Tr[W (t) ⊗ W * (t)], where W ⊗ W * is the doubled Floquet operator, and here the trace is over the q 2L -dimensional doubled Fock space. To couple the system to a bath we act on each of the sites x = 0 and (L − 1) with independent q × q Haar-random unitary matrices at each time step. Writing the composite index for the end sites 0 and (L − 1) using Greek letters α, β, and the composite index for the (L − 2) central sites using Roman letters a, b, the components of the doubled Floquet operator W ⊗ W * are W aα,bβ W * a * α * ,b * β * . Averaging over the bath couplings we find that the only nonzero matrix elements have α = α * and β = β * . This average therefore forces local diagonal pairings of the orbits of each of the sites x = 0 and (L − 1). For a more detailed discussion, see Appendix B.
Defining the Kraus operators Q αβ via their components Q αβ ab = (1/q)W aα,bβ , the quantum channel describing the evolution of the central (L−2) sites is represented by the operator Q acts on the q 2(L−2) -dimensional doubled space of the central (L − 2) sites. From this a sum over orbits with local diagonal pairing at the ends of the system. The matrix Q has leading eigenvalue unity, and the corresponding eigenvector describes a density matrix proportional to the identity. In generic circuit realisations the other eigenvalues of Q lie within the unit circle, and so at late times Z(0, t) approaches unity. We show this behaviour for individual circuit realisations in Fig. 18, and find also that for s nonzero, Z(s, t) decays to zero with increasing time. In the language of Sec. III D this implies an approach toward a global orbit pairing in individual circuit realisations. From the perspective we have just described on the other hand, the statement that Z(0, t) approaches unity at late times is a statement about thermalisation.

IV. EIGENSTATE CORRELATORS
Our attention has so far been limited to spectral properties. However, the picture of local orbit pairing, and in particular the possibility of domain walls in the pairing of many-body orbits, has implications for local correlations between eigenstates. We now show that there are correlations between the diagonal matrix elements of local operators, or equivalently the reduced density matrices of eigenstates. These correlations, relative to the predictions of the ETH, grow without bound in the thermodynamic limit. We note that apparently related correlations have been observed in a recent numerical study [67].
We first introduce probes for these correlations. Restricting ourselves to a single site x, let τ j,x (j = 0 . . . (q 2 − 1)) be a set of orthonormal q × q Hermitian operators acting only on this site, which for brevity we will refer to as observables. Then Tr[τ i,x τ j,x ] = δ ij , and we choose τ 0,x = 1/ √ q so τ j =0,x are traceless. The reduced density matrix of the eigenstate |n on site x is A correlator between diagonal matrix elements of the observables can then be defined as where Tr x denotes a trace over the site x. In a similar way, for the off-diagonal matrix elements, consider summing | n|τ j,x |m | 2 over τ j . The result is Here ρ x (n) = Tr x |n n| is the reduced density matrix of the eigenstate |n over the (L − 1)-site complement of the site x, and Tr x is a trace over that region. In the time domain the correlations between diagonal matrix elements are characterised by the reduced form factor (RFF) where in the summand we have the correlator defined in Eq. (52). Spectral structure in the off-diagonal matrix elements is instead encoded in where now the object introduced in Eq. (53) appears. Note that R x (t) is the autocorrelation function of the operator τ j,x , averaged over choices of τ j,x , and is therefore accessible in principle to experimental measurements. In Fig. 19 we show diagrams of the correlators R x (t) and R x (t) in circuit notation, as well as the relevant orbit illustrations (as in Fig. 2). In light of the picture of local orbit pairing developed through Sec. III we see that R x (t) and R x (t), and therefore the diagonal and off-diagonal matrix elements, behave very differently.
In R x (t) the evolution operator appears as the product of Tr x W (t) and its Hermitian conjugate. Away from site x, R x (t) has the same trace structure as K(t). This trace structure is associated with a freedom in the local orbit pairing, or more formally with the transfer matrices T x,x+1 discussed in Sec. III. Therefore, just as K(t) grows exponentially with L, so does R x (t). Exact calculations are straightforward for the model of Ref. [27] in the largeq limit, yielding R x (t) = (q/t)K(t), to be compared with R x (t) = q within RMT. Therefore at times t < t Th , for which K(t) > t, the average RFF R x (t) exceeds its RMT value. The freedom in the local orbit pairing gives rise to correlations between the diagonal matrix elements, and this effect is stronger in larger systems.
In Appendix H we show how the correlator R x (t) can be evaluated using an extension of the transfer matrix method of Sec. III. Here we present numerical results, obtained in the quasienergy domain using exact diagonalisation and Lanczos methods (see Appendix A for details), that demonstrate divergent departures from ETH.
Note first that the completeness of eigenstates im-   It is useful to parametrise the correlator as follows Comparing with Eq. (56) we see that the ETH prediction is r x,nm = −1. More generally, the sum rule satisfied by the correlator ensures that the average of r x,nm over n or m is −1. It is reasonable to expect that the ensemble average is determined only by the quasienergy difference, so r x,nm = r x (ω = |θ n − θ m |). We define also the average of r x (ω) over x, r(ω), and for open boundary conditions exclude the two sites at each of the two ends of the chain from this average. r(ω) is a correlation function for the diagonal matrix elements of local observables, averaged over position and summed over an orthonormal set of observables at each site.
In Fig. 21 we show r(ω) for both open and periodic boundary conditions. The deviations from the ETH are striking, and they are most prominent at the smallest quasienergy separations, where the ETH is usually thought to be applicable. The observed low-frequency peak has its origin in the late-time decay of the leading eigenvalues of the transfer matrix, λ(ω, t), toward unity (see Figs. 10 and 15(a)). From this we conclude that at large L the height of the peak grows exponentially with L while its width is given by the inverse of the Thouless time, which decreases with increasing L. With open boundary conditions and for only L = 14 sites we find already a relative enhancement of the correlator r(ω) by over two orders of magnitude. This result constitutes a substantial correction to the ETH for chaotic one-dimensional Floquet systems.
The picture of local orbit pairing does not suggest a significant modification to the ETH for the off-diagonal matrix elements of local observables, at least at the level of R x (t). To see this, note that R x (t) is constructed from the product of Tr x W (t) and its Hermitian conjugate. Based on the trace structure, there is only a freedom in the local orbit pairing at the site x. This fact is demonstrated by a large-q calculation in the model of Ref. [27] where we find R x (t) = q L−1 (1 + (t − 1)e −2εt ) for 1 ≤ x ≤ (L − 2), to be compared with the RMT result R x (t) = q L−1 . The deviations of R x (t) from RMT do not grow with L.
Objects analogous to R x (t) and R x (t) can also be de-fined for operators with arbitrary spatial support. From Eqs. (54) and (55), as well as Fig. 19, we see that R x (t) and R x (t) are complementary. Considering operators with support on sites (where = 1 in the above), there is a freedom in the local orbit pairing over (L − ) sites in the analogue of R x (t) and over sites in the analogue of R x (t). The substitution ↔ (L − ) converts between the analogues of R x (t) and R x (t). Whereas the spectral structure in the diagonal matrix elements grows with (L − ), the spectral structure in the off-diagonal matrix elements grows with .

V. SPECTRAL FLUCTUATIONS
We now turn to an investigation of the sample-tosample fluctuations of spectral correlations within the ensemble of Haar-RFCs. In Fig. 3 we have already shown that a local average dramatically suppresses statistical fluctuations of the SFF, a global quantity. Here we set out to understand this and the distribution of spectra more generally.
Information on the distribution of spectra is buried in the moments of K(t). The n th moment, K n (t), is an average over the product of n copies of the forward orbits TrW (t) and n copies of the backward TrW * (t). Focusing on only the second moment, in Sec. V A we make some first steps toward adapting our theory of local orbit pairing to these higher order objects, and demonstrate non-Gaussian statistics of the SFF. In Sec. V B we discuss the connection with the entanglement membrane.

A. Non-Gaussian statistics
To ground the discussion, consider first the case where W is a single random matrix drawn from the Haar distribution. In the limit of large matrix dimension the object TrW (t) is then normally distributed, with mean zero, in the complex plane [54]. In calculating the second moment of the spectral form factor, K 2 (t) = |TrW (t)| 4 , Wick's theorem gives a sum over the two possible pairings of copies of TrW (t) and its conjugate, where we have denoted the two copy pairings '+' and '−'. Each of the '+' and '−' copy pairings contributes K 2 (t), so K 2 (t) = 2K 2 (t). Here the copy pairing is necessarily global since our evolution operator is just one Haar-random matrix.
The global pairing of copies in Eq. (58) is to be contrasted with the situation in a circuit model. We have seen already that the local freedom in the pairing of orbits gives rise to deviations of the average SFF from the RMT prediction. In K 2 (t) there is an additional local freedom in the pairings of copies of these orbits. We show in Fig. 22 that this gives rise to non-Gaussian statistics of the spectral form factor at early times. Specifically we show that K 2 /(2K 2 ) is signficantly larger than unity even in small systems, and grows with L at fixed t.
As a point of comparison, we now sketch the calculation of K 2 for the model of Ref. [27] in the large-q limit; see also [68]. In K 2 (t) = |TrW (t)| 4 each one-site Haarrandom gate appears t times in each of the two copies of TrW (t), and similarly its conjugate appears t times in each of the two copies of TrW * (t). Haar-averaging a one-site gate at large q we find a sum over pairings of orbits analogous to Eq. (12), and additionally a sum over the two pairings of copies of these orbits, as in Eq. (58). K 2 (t) then involves a sum over all local orbit pairings, as in Eq. (18), as well as a sum over all local copy pairings. Whereas a domain wall in the orbit pairing has statistical weight e −εt , it can be shown that a domain wall in the copy pairing here has statistical weight e −2εt . With open boundary conditions, at late times The first term on the right-hand side of Eq. (59) involves contributions from domain walls in the orbit pairing, as in Eq. (19), whereas the second term arises from a domain wall in the pairing of copies of orbits. The factor (L − 1) is the translational entropy of this domain wall, and the factor t 4 arises from sums over the orbit pairing at each of the two ends of the two copies of the chain. With periodic boundary conditions the leading term in K 2 (t) − 2K 2 (t) at late times comes from configurations with two domain walls in the copy pairing. Fig. 22 shows exactly this type of behaviour. A useful analogy with Sec. III is now evident. There a great deal was learnt about K(t) via calculations in which boundary conditions were imposed to force diagonal pairings of the local orbits. In this approach the emergence of RMT level statistics on small energy scales in the Haar-RFC can be understood as arising from a global pairing of orbits at late times. Looking at Eq. (58), we see that in order to recover RMT statistics of the spectral form factor beyond only the first moment, we also require a global pairing of copies of orbits.
Proceeding in a similar way then to Sec. III D, here we investigate how, in the calculation of the second moment of the spectral form factor, a local freedom in the pairing of copies of orbits gives way to a global pairing. To do this we again impose local diagonal pairings at opposite ends of the system, but now to force a domain wall in the pairing of copies of orbits. This setup is ilustrated in Fig. 23. In principle we could also impose a relative twist on the orbit pairing, as in Z(s = 0, t) (see Eq. (39) and Fig. 12), but we choose to restrict ourselves to equal-time pairing, as in Z(0, t).
In the notation of Eq. (58), Fig. 23 illustrates a '+' pairing on one end of the chain, and a '−' pairing on the other. The resulting object, which we denote Z 2 (+, −, t), is a multiple sum over two copies of the forward orbits and two copies of the backward orbits, with a do- main wall in the pairing of copies running along the time direction. Z 2 (+, −, t) is to be contrasted with Z 2 (+, +, t) = Z 2 (0, t), in which we impose equal-time pairings of the same copies of orbits at both ends of the system.
We show the averaged quantities Z 2 (+, −, t) and Z 2 (+, +, t) in Fig. 24. We see that Z 2 (+, −, t) grows with L at fixed t, and therefore so do the deviations from RMT. At late times, t 10 for L ≤ 8, Z 2 (+, −, t) goes to zero. In other words the contributions from those sets of many-body orbits with domain walls in the copy-pairing vanish. As expected, Z 2 (+, +, t) approaches unity in this regime. The appearance of domain walls in the copy pairing implies non-Gaussian statistics of the SFF at times t 10 for these system sizes, as revealed by direct calculation in Fig. 22. As in Sec. III, the deviations from RMT are more prominent with open boundary conditions than with periodic because domain walls can only appear singly in the first case.
We now return to the effect of the single-gate average observed in Fig. 3. To understand why an average of the spectral form factor over only one gate of the circuit is enough to suppress statistical fluctuations, consider performing this average and subsequently calculating the variance of the spectral form factor over the ensemble of the remaining gates. The tendency toward global copy pairing observed in Fig. 24 implies that the second moment can be approximated by a sum over the two global pairings ('+' and '−') at late times, as in Eq. (58). However, in performing the one-gate average beforehand, we have locally selected for just one of these pairings. Consequently the second moment is approximately halved, and so the variance is suppressed to near zero.

B. Entanglement membrane
The domain walls in the copy pairing which we have discussed through this section are related to the entanglement membrane [38]. To see this we recall the growth of entanglement under a unitary circuit, with or without time periodicity.
Dividing space into subregions A and B, the purity of subregion A is defined as e −S2 = Trρ 2 A (t), where ρ A (t) is the reduced density matrix in subsystem A and S 2 is the second Renyi entropy. Writing the initial density matrix as ρ, we have Here a boundary condition at time 0 is set by ρ. The boundary condition at time t is determined by the trace structure in Eq. (60), and is illustrated in Fig. 25. The boundary condition at time t imposes a domain wall, the entanglement membrane, in the pairing of the indices of the evolution operator. A similar domain wall in the pairing of copies appears in Z 2 (+, −, t). In that case, as shown in Fig. 23, it is imposed by boundary conditions at the left and right ends of the chain, while the boundary conditions in time are periodic.

VI. DIAGONAL APPROXIMATION
As discussed in Sec. III A, in Haar-RFCs the behaviour of the SFF beyond RMT comes from the contributions of off-diagonal pairings of many-body orbits. However, deviations from RMT can also arise within the diagonal approximation through a different mechanism. This involves the subleading eigenvalues of the diagonal propagator P. These contribute in two ways: through having nonzero average values, and through their fluctuations. In this section we discuss both types of contribution.
First we consider a general class of RFCs in which the gates can be tuned from Haar-random to the identity. Away from Haar-randomness the subleading eigenvalues of the diagonal propagator have nonzero averages. We parametrise the two-site gates U of the circuit via their spectral decompositions, U = V e −iϕ V † . Here ϕ is a diagonal matrix, and we take the unitary matrix V to be Haar random. If we choose ϕ to be distributed as the eigenphases of a Haar-random unitary then U is Haar random. On the other hand if ϕ is proportional to the identity so is U . By tuning the distribution of ϕ we can interpolate between these two extremes.
In a brickwork model the diagonal propagator is of the form Eq. (15). The gate acting on sites 0 and 1 in the first half-step appears in P a0...a L−1 ,b0...b L−1 as U c0c1,b0b1 U * c * 0 c * 1 ,b0b1 , where the indices c, c * are to be summed over. The average over U now involves an average over V and an independent average over the ϕ-distribution. The result is where Clearly for ν = 0 we recover the Haar-random result Eq. (16). On the other hand for U proportional to the identity we have ν = 1. The average propagator, with matrix elements P a0...a L−1 ,b0...b L−1 , is given by multiplying expressions of the form Eq. (61) and summing over indices c, c * . The result is a matrix-product operator (MPO), and with periodic boundary conditions this takes the form where M ab is a 2 × 2 matrix with components The eigenvectors of the MPO Eq. (63) are product states, which we label σ This is achieved by taking v 0 to have constant entries, and v σ =0 to make up the rest of an orthogonal set (the sum of entries of any v σ =0 is therefore zero). The resulting 2 × 2 matrices µ σ , with components µ ij σ , are outer products and the eigenvalues of P are Tr[µ σ0 . . . µ σ L−1 ]. The leading eigenvalue is unity as required by unitarity of W , and the next-to-leading eigenvalues are ν 2 with multiplicity L(q − 1). Within the diagonal approximation the average SFF K(t) = tTrP t is then where in the second line we have expanded around the late-time result. For 0 < ν < 1, K(t)/t decreases monotonically with t. Writing ν 2t = e −2 ln(1/ν)t we see that the Thouless time is t Th = ln L/(2 ln(1/ν)) within the diagonal approximation. As indicated below Eq. (8) a full exploration of the diagonal approximation requires a discussion not only of P but also of P t . We now discuss this aspect in the context of Haar-RFCs. In this case, as we have shown following Eq. (17), the average P has a single eigenvalue of unity, with all subleading eigenvalues zero. Individual realisations, however, typically have non-zero subleading eigenvalues which contribute to P t . We probe the consequences of this by evaluating the SFF using the approximation K(t) ≈ tTrP t . In Fig. 26 we compare TrP t to K(t)/t for Haar-RFCs, and our results demonstrate that corrections to RMT arising from the subleading eigenvalues of P cannot account for the behaviour of Haar-RFCs.
In the general case the SFF has contributions both from subleading eigenvalues of the diagonal propagator and from domain walls. Determining which effect is dominant requires a comparison of the associated timescales. For Haar-RFCs, with ν = 0, we have seen in Sec. III that the domain wall tensions ε are finite. Since ν sets the proximity of the gates to the identity it is natural to expect these tensions to decrease with increasing ν. On the other hand the subleading eigenvalues of the diagonal propagator vanish for ν = 0. The behaviour we have studied in Haar-RFCs, in which domain walls control the SFF, must therefore extend to at least a finite range of ν, and so is not restricted to the Haar case.
A contrasting instance is provided by quantum circuits which are built from two-site gates but without any spatial structure imposed. It has recently been argued that models of this type have a 'ramp' (or Thouless) time that is logarithmic in system size [53]. In the absence of spatial structure a domain wall interpretation cannot be appropriate, and it would be interesting to investigate whether this behaviour can be attributed to subleading eigenvalues of P.

VII. DISCUSSION
The spectral statistics of many-body quantum systems in the ergodic phase coincide with the predictions of RMT below a certain (quasi)energy scale, the inverse of our Thouless time t Th . On timescales much greater than t Th , but much shorter than the Heisenberg time t H , this result can be understood through the diagonal approximation to the SFF. In this work we have determined the regime of validity of the diagonal approximation, and the form of the corrections to it, in systems having local interactions and no conserved densities. One of our main results is to reveal a generic mechanism setting the Thouless time in this context. Our approach has been to develop a theory of local orbit pairing, centred on the properties of a transfer matrix which acts on pairs of local orbits and which generates the SFF. In large systems the average SFF is controlled by the leading eigenvalues of the average transfer matrix; we have calculated these eigenvalues and shown that there is a connection between the corresponding eigenvectors and the diagonal approximation to sums over pairs of manybody orbits. At fixed time and in the limit of large systems the dominant contributions to the SFF come from orbits which are locally paired as in the diagonal approximation, but with distinct diagonal pairings in neighbouring spatial domains. This domain structure is associated with the exponential growth of the SFF with system size at fixed time, and the corresponding divergence of the Thouless time. Conversely, in a system of fixed large size there is a wide window of time t Th t t H in which the pairing for the whole system forms a single domain, so the diagonal approximation is accurate.
We believe the structures we have uncovered are universal, in the sense that our results should apply quite generally to lattice Floquet models with local interactions, and with the form of the deviations of the SFF from RMT governed by the cost of domain walls in the orbit pairing. On timescales much greater than the Floquet period we expect that systems in the same symmetry class that share a domain-wall cost, but are different in microscopic detail, have similar spectral statistics. This kind of universality is familiar from studies of single-particle disordered conductors. In that case the form of the departure of the spectral statistics from RMT on timescales t < t Th is set by the diffusion constant, a coarse-grained quantity [69].
The picture of local orbit pairing has also revealed strong correlations between the diagonal matrix elements of local observables. We have shown that, relative to the predictions of the ETH, these correlations grow exponentially with increasing system size. Our results represent a prominent correction to the ETH for Floquet models with local interactions. In the same way, we have argued that deviations from the ETH arise in the off-diagonal matrix elements of non-local operators, and that these features grow exponentially with the support of the operator as opposed to the system size. This result has implications for operator spreading, and furthermore it has been noted that the SFF can be related to the autocorrelation functions of operator strings [53]. It would be interesting to understand the connection between these two approaches.
Our results should be compared against calculations within the diagonal approximation, in which the con-tributing many-body orbits are globally paired. In that approximation Haar-RFCs behave like Haar-random matrices, and the non-RMT behaviour we have observed in Haar-RFCs cannot be understood. There is, however, a separate mechanism, additional to the one involving pairing domains, that leads to non-RMT behaviour even from global pairing. To illustrate this mechanism we have introduced a generalised class of RFCs with gates that can be tuned between Haar randomness and the identity. Away from the Haar case, these circuits generate within the diagonal approximation a Thouless time growing logarithmically with system size L. The mechanism behind this, observed in Ref. [25], is distinct from the domain wall contributions we have investigated here, and the effects are sub-dominant for ensembles in the vicinity of Haar randomness.
In interacting self-dual models, however, the diagonal approximation to the average SFF is exact. This is because the transfer matrix which generates it has leading eigenvalues equal to unity [30,42,46]. Within our framework this implies an infinite domain wall tension, and therefore a vanishing contribution to the SFF from pairs of many-body orbits involving domain walls. Additionally, the subleading eigenvalues of the diagonal propagator in the Ising basis vanish at the self-dual point [25]. The vanishing of both the domain wall tension and the subleading eigenvalues of the diagonal propagator conspire to give an SFF equal to the RMT result. Perturbing away from the self-dual point causes the leading eigenvalues of the transfer matrix to deviate from unity [42], giving rise to behaviour similar to that observed here in Haar-RFCs, for example the exponential growth of the SFF with system size at fixed time.
It is natural to ask what other types of correction to the diagonal approximation can arise in this setting. Important candidates are provided by the class of interference effects studied by Sieber and Richter [15]. In systems with time-reversal symmetry these are crucial to recover the exact RMT behaviour, while they cancel in systems without time-reversal symmetry. In time-reversal symmetric Floquet quantum circuits these interference effects were studied in Ref. [25], where the authors recover both the leading and dominant subleading terms in the expansion K(t) = 2t − 2t 2 /t H + . . . that is expected from RMT for t t Th . The physical mechanism responsible for Sieber-Richter corrections is distinct from the pairing domains which have been our focus. In large systems these effects appear in different time regimes, t t Th and t < t Th , respectively.
Moving beyond the average SFF, we have initiated a similar investigation into the higher moments (see also Refs. [68] and [46]). Our discussion highlights another source of deviations from RMT, associated with a freedom in the local pairing of the multiple copies of orbits which appear. The presence of domain walls in the pairing of copies of orbits enhances these higher moments, and consequently gives rise to non-Gaussian statistics of the SFF. We have discussed the connection between these domain walls and the entanglement membrane [36][37][38]. Although our transfer matrix approach was here restricted to the first moment of the SFF, it is clear that one can construct analogous transfer matrices to generate the higher moments in generic Floquet models. Such an investigation would connect with recent work on rare region effects. Weak links for example have been shown to play an important role in entanglement growth [70]. For the average transfer matrix generating the first moment of the SFF, the possibility of weak links would be captured by the decrease of the effective domain wall tension at late times.
Although our work has focused on systems without any locally conserved densities, their introduction is known to lead to prominent deviations of the SFF from RMT even within the diagonal approximation [56,58,65]. With a conserved scalar charge, for example, these deviations persist up to a Thouless time which scales diffusively t Th ∼ L 2 . This is far greater than the Thouless time which we have shown to arise from pairing domains, and it would be interesting to understand the interplay between these effects. There is of course an outstanding question of how to apply our ideas to Hamiltonian systems. In that case the local conservation of energy surely plays a role.
Our approach to studying the average transfer matrix involves imposing local orbit pairings on the doubled time-evolution of the system. Fixing different pairings at the two ends -our twisted boundary conditionswe forced domains walls into the many-body orbit pairing. We have also shown that the untwisted case, where the same pairing is imposed at each end, can be thought of as introducing local couplings to a bath. Applying these same boundary conditions to a many-body localised phase we are faced with its destabilisation [71], and this behaviour must be reflected in the transfer matrices. A central question in that context, which we address elsewhere [72], concerns the behaviour of the transfer matrix eigenvalues on crossing a phase boundary between an ergodic and a many-body localised phase.
tary Haar distribution we use Monte Carlo. We sample this distribution as follows [73]. We first generate complex matrices A with independent, unit-normally distributed entries. Performing the QR decomposition A = QR we then define the diagonal matrix D with entries D ii = R ii /|R ii |. The unitary matrices U = QD are then Haar-random.
In Fig. 3(b) we have averaged the SFF over 10 4 realisations of one gate, and in (c) over 10 6 realisations of the Floquet operator. Our calculations of the average SFF K(t) in Fig. 5 use exact diagonalisation (ED) with 10 6 realisations for L ≤ 8 and with 10 5 for L ≥ 9. These calculations were then used for the analysis in Figs. 10 and 11 in conjunction with the methods in Appendix F. Our approach for the calculations in Figs. 13 and 14 is detailed in Appendix B, and the analysis leading to Fig. 15 is detailed in Appendix F.
In Figs. 20 and 21 we have carried out calculations in the quasienergy domain using ED in system sizes 8 ≤ L ≤ 12 and using a Lanczos method for L = 14. The Lanczos method is particularly efficient since the action of the Floquet operator W on many-body states is specified by its definition in terms of local unitary gates. Standard algorithms require Hermitian operators, and so instead of working with W we determine first a set of eigenvectors |n of 1 2 (W + W † ), with eigenvalues cos θ n . There are no symmetries or degeneracies, so this approach unambiguously determines the eigenstates of W . The sign of θ n is then determined by acting on the eigenvectors with 1 2i (W − W † ). In all system sizes we sample eigenstates of W whose quasienergies reside in bins of fixed width 20(2π)/(2 14 ) centred on θ = 2πk/50 for k = 0 . . . 49. For each eigenvector |n we only store its eigenphase θ n and the diagonal matrix elements of a complete set of local Hermitian operators τ x,j (with x = 0 . . . (L − 1) and j = 0 . . . (q 2 − 1)), n|τ x,j |n . Square bins for the sampling of quasienergies implies triangular bins for the sampling of quasienergy differences ω = |θ n − θ m |. Having fixed the bin widths, the number of circuit realisations we use varies with L in such a way that we have over 10 7 contributions to each Tr[ρ x (n)ρ x (m)] data point in Fig. 21. This number is based on using 10 3 realisations of the Floquet operator for L = 14.

Appendix B: Monte-Carlo pairing
In Sec. III we have imposed local pairings between forward and backward histories, and through this we have calculated Z(s, t), defined in Eq. (39). A similar approach was used in the calculation of Z(+, +, t) and Z(+, −, t) in Sec. V. In this Appendix we describe the numerical method used to impose these local pairings of histories.
The local diagonal pairing s is represented in the product space of forward and backward histories by the state |s , defined in Eq. (37). Each of the forward and back-ward histories is a vector in a space of dimension q 2t , and so the product space has dimension q 4t . Since we are interested in late-time behaviour, it is not feasible to perform calculations in this space directly. One alternative, which we do not follow here, is to study the doubled Floquet operator W ⊗ W * , which acts on a space of dimension q 2L . Our Monte-Carlo approach instead requires us only to work with operators acting on a space of dimension q L . We calculate properties of the forward and backward histories independently at first, and impose pairing through averaging.
Focussing on a fixed Floquet operator W for a chain with open boundary conditions, we discuss this construction for Z(s, t). At the time step r, where r runs from 0 to (t − 1), we act on each of the left-hand and right-hand sites with independent q × q Haar-random matrices, u (r) L and u (r) R , respectively. These matrices are also independent for different steps r. The forward evolution operator for step r is then where we act with the identity on each of the central sites. This gives the sum over forward orbits For the backward evolution, at step r we instead act with u (r+s L ) L on the left-hand site, and u (r+s R ) R on the righthand site, where s L and s R are integers 0 . . . (t − 1) and addition is defined modulo t. The backward evolution operator for step r is the conjugate of and the sum over backward orbits is The double sum over forward and backward orbits is now TrW f (t)TrW * b (t). If we average this double sum over the matrices u (r) L for example, using we impose the local diagonal orbit pairing s L at the left end of the system. Similarly an average over u (r) R imposes the pairing s R at the right end. This procedure gives Z(s R −s L , t). The configurations of single-site gates u (r) L,R used to impose local diagonal pairings are illustrated in Fig. 27 for Z(0, 2) and Z(1, 2).
We use a very similar method to calculate Z 2 (+, +, t) and Z 2 (+, −, t) in Sec. V. These objects each involve two copies of the sum over forward many-body orbits and two copies of the sum over backward many-body orbits. The different single-site Haar-random unitary matrices are then used to fix local pairings of the multiple copies L,R , drawn as open circles, used to impose the local diagonal pairings at the ends of the system which define Z(s, t) (see Eq. (38)), here for t = 2. Grey and black circles represent independent matrices, and furthermore the matrices at the left-hand sites and the right-hand sites are independent. Other conventions as in Fig. 6. Averaging over the single-site matrices in the left-hand diagram fixes sL = 0 and sR = 0, so we find Z(0, 2). In the right-hand diagram this average fixes sL = 0 and sR = 1, so we find Z(1, 2). of the orbits, and their arrangements are illustrated in Fig. 28.
In the calculations of Z(s, t) in Fig. 18, where we consider individual realisations of the Floquet operator, we use 10 4 realisations of the single-site unitary matrices u In this Appendix we discuss the large-q limit of the average transfer matrix T (t), defined in Eq. (29) in terms of Weingarten functions Wg(στ −1 ) and unnormalised permutation states |σ, τ .
The Weingarten functions are maximised for σ = τ , taking the value Wg(1) = q −2t in the large-q limit [62,63]. Furthermore, the states |σ, τ that maximise σ, τ |S|σ, τ are those in which σ = τ are permutations mapping r → r + s modulo t. Normalising these states, we find |s defined in Eq. (28). This motivates the large-q approximation to the average transfer matrix |s s| . (C1) FIG. 28. As in Fig. 27, but here to calculate Z2(+, +, t) (lefthand diagram) and Z2(+, −, t) (right-hand diagram). Here t = 1. The two copies of the forward evolution operator are shown light, and the two copies of the backward evolution operator are shown dark. Averaging over the single-site matrices fixes a '+' pairing between the copies of the many-body orbits on the left-hand site in both of the diagrams, and also on the right-hand site in the left-hand diagram. In the right-hand diagram, the average imposes a '−' pairing on the right-hand site.
The different states |s correspond to the different local diagonal orbit pairings, and so in the large-q limit the notion of local orbit pairing has a very sharply defined meaning. Note however that the states |s are not orthogonal. Approximating T (t) by Eq. (C1) at general values of q, it is straightforward to calculate the SFF, and this approach reveals a number of interesting features. Note that Eq. (C1) is the transfer matrix for a clock model: s|s depends only on |s−s |. We first diagonalise T (t) by writing it in terms of the orthogonal (but not normalised) states |ω = (1/ We will see below that the inner products s|s for s = s decay at least as quickly as q −t , and so ω|ω approaches unity at late times. Also, since s|s > 0, the largest eigenvalue of T (t) is in the ω = 0 sector, as at small q (see for example Fig. 15). From the eigenvalues ω|ω we can calculate the average SFF with periodic boundary conditions With open boundary conditions on the other hand Because the eigenvalues ω|ω approach unity at late times, we recover the RMT result K(t) = t. The eigenvalues ω|ω are given by the inner products s|s through Eq. (C3), and these are Evaluating the sums on the right-hand side we see that s|s is determined by the number of cycles in the permutation mapping r → r + |s − s | mod t. Denoting this cycle number N t,|s−s | , There is clear s dependence in this expression, and therefore in the domain wall tension and Z(s, t) (defined in Eq. (39)), as observed for q = 2 in Fig. 14(a). This is to be contrasted with the s-independent domain wall tension found in the large-q limit in the model of Ref. [27].
There is an interesting difference between the SFF at odd and even times. For an even value of t, the largest inner product 0|s for s = 0 is 0|t/2 = q −t . On the other hand, for t an odd multiple of 3, the largest are 0|t/3 = 0|2t/3 = q −4t/3 . The deviations of ω|ω from unity can be significantly larger for t even than for t odd. From Eqs. (C4) and (C5) we see that this implies larger deviations of the SFF from RMT at even times. We have observed the same behaviour at small q in Fig. 5.
The q 4t × q 4t matrix T (t) has cokernel spanned by the vectors |σ, τ , of which there are (t!) 2 . For small t it is then feasible to calculate T (t) exactly using Eq. (29). In this Appendix we set out the technical details required for this calculation.
With this structure in place we can separate the sum in Eq. (29) into separate sums over roots |σ, τ ; 0, 0 and the states within the corresponding classes. The average transfer matrix Eq. (29) involves the outer products S |σ, τ σ, τ |, and having chosen the roots we can write where the primed root |σ , τ ; 0, 0 is defined relative to |σ, τ ; 0, 0 . The integers a + (σ, τ ), a − (σ, τ ) depend on our initial choice of roots. We are now in a position to write down T (t) in terms of roots and classes. We find where in the summand r ± ≡ r ± + a ± (σ, τ ) and σ , τ are defined via Eq. (D3). Here we have used the fact that the Weingarten function depends only on the conjugacy class of its argument.
In terms of the orthonormal basis states |i we have The transfer matrix block ω + , ω − is then The different ω + , ω − blocks of T can be diagonalised numerically, giving exact results for the eigenvalues and eigenvectors. The leading eigenvalues computed through this method, as well as the overlaps of the corresponding eigenvectors with the local diagonal states, are shown in Fig. 15 for t ≤ 5, and we find excellent agreement with our other approaches.
Appendix E: Subleading eigenvalues of T (t) In the main text we have focused on the leading eigenvalues and eigenvectors of the average transfer matrix T (t) for q = 2, which control the domain structure in the orbit pairing. On the other hand, the subleading eigenvalues of T (t) control the behaviour beyond the Heisenberg time t H . Here, using the results of Appendix D, we examine some aspects of their distribution.
As we have shown in Fig. 15, the t leading eigenvalues tend to unity at late times. Therefore their contribution to K(t), with periodic boundary conditions for example, is simply t at late times. However, for t > t H = q L the SFF plateaus at K(t) = q L . For t ≤ q L the subleading eigenvalues appear to do nothing, whereas for t = q L + δt their contribution to K(t) is −δt.
In practice we have access to all of the eigenvalues of T (t) for t ≤ 5, so for q = 2 we can verify the role of the subleading eigenvalues across the Heisenberg time (t H = 4 for L = 2). We indeed find K(4) = K(5) = 4. On the other hand for q ≥ 3, t H ≥ 9, so we cannot probe behaviour beyond the Heisenberg time. The eigenvalue distributions for times t = 4 and t = 5 are shown in Fig. 29, for local Hilbert space dimensions q = 2, 3, 4. At these times there is a clear gap between the magnitudes of the t leading eigenvalues and the magnitudes of the subleading ones, which becomes more pronounced with increasing q.

Appendix F: L-scaling
In this Appendix we describe the methods used to determine the leading eigenvalues of the transfer matrix, as well as the overlaps of the associated eigenvectors with the boundary states B L |, |B R and the local diagonal states |ω .
In Sec. III C we have shown that with open boundary conditions and large L the behaviour of the leading eigenvalue of T (t) in the ω = 0 sector, λ(0, t), controls the behaviour of the average SFF. Taking the logarithm of Eq. (34) we have ln K(t) (L − 1) ln λ(0, t) + ln B L |ω, t; R ω, t; L|B R , where we have neglected contributions from the subleading eigenvalues. Fixing t and varying L we have extracted ln λ(0, t) and B L |ω, t; R ω, t; L|B R from linear fits to ln K(t) versus (L−1), and the results are shown in Fig. 10. In practice, for a given t we fit only to data with sufficiently large L that t < 1 2 t H ≡ 1 2 q L . By restricting to L such that t is well below the Heisenberg time, we avoid contributions to K(t) from subleading eigenvalues. In Fig. 30 we show ln K(t) versus L for various t, and find excellent linear scaling.
In Figs. 31 and 32 we show ln Z(ω, t) for the sectors ω = 0 and ω = 2π/t, respectively, as a function of L and for times 1 ≤ t ≤ 20. The linear relationship suggested by Eq. (F2) holds very well even for systems of only L = 3 sites (and so only two gates), and from fits we can read off λ(0, t) and ω|ω, t; R ω, t; L|ω . The results are shown in Fig. 15. Fits for other values of ω are of similar quality. As we discussed through Secs. III D and III E the local diagonal pairings |ω are closely related to the leading eigenvectors of T (t). Consequently we expect that the contributions of subleading eigenvalues to Eq. (F2) are suppressed. Moreover, while there is an abrupt change in the behaviour of K(t) at t = t H , this is not the case for Z(ω, t). For these reasons, in contrast to Fig. 30, we include all available data in the fits in Figs. 31 and 32. The quality of these fits is evidence that the subleading eigenvalues are playing only a minor role, if any.

Appendix G: Heisenberg interactions
To investigate the generality of our results, in this section we consider another Floquet model in its er-  godic phase. We again use a local Hilbert space dimension q = 2, but whereas in the main text nearestneighbour interactions were implemented via 4 × 4 Haarrandom unitary matrices U , here we instead write U = [B ⊗ B ]e iπJΣ [A ⊗ A ], where A, A , B and B are independent 2 × 2 Haar-random unitary matrices acting on the individual sites, Σ is the two-site swap operator (essentially a Heisenberg coupling), and J is a fixed interaction strength. The model is therefore a kicked spin-half Heisenberg chain with random fields acting at each site, and as in the Haar-random case the Floquet operator does not have time-reversal symmetry. This model is many-body localised for small J < J c , as we have discussed elsewhere [72], but here we restrict ourselves to J = 1/4 deep in the ergodic phase.
In Fig. 33 we calculate Z(s, t) and Z(ω, t) for L = 8 sites and different values of s and ω. As in Fig. 14(a) of the main text, Z(0, t) approaches unity at late times whereas Z(s = 0, t) decays to zero. This implies the observed behaviour of Z(ω, t). From the scaling of Z(ω, t) with L we then extract the leading eigenvalues of the average transfer matrix T (t) for this model, as well as the overlaps of the corresponding eigenvectors with the paired states |ω . The results, shown in Fig. 34, are similar to those for the original model in Fig. 15. In this Appendix we show how to construct the form factor giving information on correlations between the diagonal matrix elements of a local observable τ , |Tr[τ W (t)]| 2 , in terms of the transfer matrices T x,x+1 which also generate the SFF. The time-domain approach here is complementary to the calculations in the quasienergy domain discussed in Sec. IV.
For a q × q operator τ acting only on site x, we definẽ τ acting on the q 2t -dimensional space of forward orbits at the site x via its matrix elements We therefore introduceτ and it can be straightfowardly verified that the righthand side of Eq. (H2) is unchanged by the replacement ofτ byτ (t) . From Eq. (H2) and its conjugate, the form factor |Tr[τ W (t)]| 2 can be constructed as × ω + , ω − , α L ; t|(τ (t) ⊗ [τ (t) ] * )|ω + , ω − , α R ; t , where in the second line we have used the spectral decomposition of the average transfer matrix. This result highlights the exponential growth with system length of correlations between the diagonal matrix elements of local observables.