Spectral decoupling in many-body quantum chaos

We argue that in a large class of disordered quantum many-body systems the late time dynamics of time-dependent correlation functions is captured by random matrix theory, specifically the energy eigenvalue statistics of the corresponding ensemble of disordered Hamiltonians. We find that late time correlation functions approximately factorize into a time-dependence piece which only depends on spectral statistics of the Hamiltonian ensemble, and a time-independent piece which only depends on the data of the constituent operators of the correlation function. We call this phenomenon"spectral decoupling,"which signifies a dynamical onset of random matrix theory in correlation functions. A key diagnostic of spectral decoupling is $k$-invariance, which we refine and study in detail. Particular emphasis is placed on the role of symmetries, and connections between $k$-invariance, scrambling, and OTOCs. Disordered Pauli spin systems, as well as the SYK model and its variants, provide a rich source of disordered quantum many-body systems with varied symmetries, and we study $k$-invariance in these models with a combination of analytics and numerics.


Introduction
New connections between many-body quantum chaos and random matrix theory (RMT) have emerged over the last several years [1][2][3][4][5][6][7][8][9][10][11][12]. Whereas traditional treatments of singlebody or few-body quantum chaos have emphasized characteristic energy level statistics of chaotic Hamiltonians [13][14][15], recent research has forged connections between RMT and the dynamics of quantum many-body systems, as probed by correlation functions. Particular emphasis has been placed on disordered systems, which naturally give rise to ensembles of Hamiltonians. Such ensembles will be the setting for much of the analysis in this paper.
There is an important structural distinction between RMT energy level statistics and correlation functions of local quantum many-body systems. Energy level statistics of an ensemble of local Hamiltonians are mostly agnostic to local physics, since energy eigenvalues are basis-independent data of a Hamiltonian. By contrast, the behavior of time-dependent correlation functions crucially depends on whether the constituent operators are local. Said simply, local correlation functions behave differently than non-local correlation functions, if the Hamiltonian dynamics is local.
Interestingly, the dynamics of time-dependent correlation functions of highly non-local operators is not sensitive to the specific choice of constituent operators [4]. Accordingly, the time-dependence of such correlation functions is completely described by basis-independent data of the Hamiltonian dynamics, i.e. energy eigenvalue correlation functions of the Hamiltonian ensemble which are the core objects of study in RMT. For instance, if A, B are highly non-local operators and E H is a Hamiltonian ensemble, we schematically have approximate equalities like where the integral over E H performs a disorder average with respect to the Hamiltonian ensemble. Note that the function [spectral quantity](t) will only depend on the joint eigenvalue statistics of the Hamiltonian ensemble E H , and not on the highly non-local operators A, B.
In fact, it often suffices that only some of the operators in the correlation function be highly non-local. Eqn. (1.1) exemplifies "spectral decoupling", namely that after a certain timescale, the dynamics of correlation functions (in a certain class of systems) only depend on spectral quantities (i.e., the joint eigenvalue statistics of the Hamiltonian ensemble), schematically denoted by [spectral quantity](t). This dynamical data is decoupled from the data of the operators A, B, via an approximate factorization as per Eqn. (1.1). We will elaborate precisely on many details in the sections below. Thus, we can predict the time dependence of disorder-averaged, highly non-local correlation functions simply by calculating energy eigenvalue correlations of the corresponding Hamiltonian ensemble. Conversely, we can calculate energy eigenvalue correlations of a Hamiltonian ensemble by studying the dynamics of disorder-averaged, highly non-local correlation functions.
This bridge between non-local correlation functions and RMT energy eigenvalue statistics is attractive, but it is often more natural to consider correlation functions of local operators. However, consider the following: take a local operator O, and evolve it with a Hamiltonian in the Heisenberg picture out to a time t as O(t) = e iHt O e −iHt . After a sufficiently long time, O(t) will be a highly non-local operator. Then it is possible that correlation functions containing O(t) will satisfy an equation like Eqn. (1.1) after sufficiently long times. In this way, we can imagine relating RMT energy eigenvalue statistics to the late time dynamics of correlation functions of initially local operators.
In this paper we make the above idea precise, and study how late time dynamics of correlation functions of disordered systems are captured by RMT energy eigenvalue statistics via spectral decoupling. A quantitative diagnostic of the dynamical onset of spectral decoupling is k-invariance, introduced in [4], which is central to our analysis. We find an intimate relationship between k-invariance, scrambling, and out-of-time-order correlators (OTOCs) [16][17][18], thus drawing on recent diagnostics of many-body quantum chaos. A subtle issue is the role of symmetry, which we study in great detail. We generalize k-invariance to systems with symmetries, including time reversal symmetries, particle-hole symmetries, etc. Disordered Pauli spin systems and the SYK model [17,19,20] will be our analytical and numerical testing ground for these ideas. These models have quenched disorder, and come with a wide variety of symmetries. We also analyze k-invariance in random circuits, leveraging recent technical results to perform explicit computations [21].
An outline of the paper is as follows: • In Section 2, we review the essentials of random matrix theory, and discuss timedependent correlation functions with respect to random matrix ensembles with different symmetries.
• In Section 3, we define k-invariance, explain its relation to the dynamical onset of chaos in time-dependent correlation functions, and make connections to scrambling and operator growth. We also give a generalization of k-invariance to finite temperature, and explore k-invariance in random quantum circuits.
• In Section 4, we introduce and study a refinement of k-invariance that accounts for symmetries of a Hamiltonian ensemble. As examples, we treat in explicit detail manifestations of time-reversal symmetry and describe a general construction.
• In Section 5, we apply our techniques to various disordered Pauli spin models and several versions of the SYK model, numerically confirming our results, and explore the role of approximate symmetries. Connections to quantum gravity and Jackiw-Teitelboim matrix models in particular are also explained.
• In Section 6, we discuss our results and future directions.
• The Appendices contain various formulas and derivations, including the first two moments of the Haar unitary, Haar orthogonal, and Haar symplectic ensembles.

Correlation functions in RMT
Consider a Hilbert H space of dimension d. We also consider a Hamiltonian ensemble E H which may be viewed as subset of the space of Hamiltonians on H, equipped with a probability measure dH such that Physically, this ensemble E H may be viewed as probabilistic instances of some disordered system. Before discussing the behavior of correlation functions in disordered quantum manybody systems, we will first review the behavior of correlation functions in conventional random matrix ensembles. The standard random matrix theory ensembles are the Gaussian Unitary Ensemble (GUE), Gaussian Orthogonal Ensemble (GOE), and Gaussian Symplectic Ensemble (GSE). We will not require the precise measures of these ensembles for our analysis here, and so refer the reader to reviews in [14,22] for details. The GUE, GOE, and GSE each belong to different symmetry classes; for instance, Hamiltonians in the GOE or GSE class represent different realizations of time-reversal symmetry. The GUE has no symmetries, meaning that there is no symmetry possessed by all Hamiltonians in the GUE. We start by focusing on the Gaussian Unitary Ensemble (GUE), an ensemble of d × d random Hermitian matrices, defined as having a Gaussian probability density P (H) ∝ e − d 2 tr(H 2 ) and a unitarily invariant measure dH = d(U HU † ). Therefore, if U is any unitary on the Hilbert space H, then for any function f (H) we have P eigval (λ 1 , ..., λ d ) (2.5) where in the second line we have used P eigvec (e 1 , ..., e d ) = 1/vol(U (d)) since all eigenbases are equally likely. We see that the probability density over the GUE depends only on eigenvalues. This feature is also present in both the GOE and GSE. 1 Accordingly, random samples from the GUE are rather non-physical. Suppose that d = 2 N for some integer N . If we fix a tensor factor decomposition of our Hilbert space into qubits as H = N i=1 C 2 , then a random sample of the GUE will have N -body interactions. Indeed, these interactions are as non-local as possible. 2 The simple reason is that locality of interactions with respect to our subsystems of qubits is basis-dependent. However, the GUE ensemble treats all bases equitably, and so random samples are agnostic to locality. As we will see shortly, this will have interesting implications for time-dependent correlation functions with respect to evolution by the GUE, namely that the locality of constituent operators will not affect the dynamics in any way.
To analyze time-dependent correlation functions with respect to the GUE, we leverage that a key corollary of unitary invariance is Haar-unitary invariance. That is, Eqn. (2.2) implies where on the left-hand side, the outer integral is over the Haar measure on U (d). Now consider two operators A and B acting on H, and their infinite-temperature time-dependent 2-point function By Eqn. (2.6) this is equivalent to where in the first line we have used the change of variables H → U HU † . We can compute the Haar-unitary integral to obtain [4] A(t)B(0) E GUE (2.9) To clean up the answer somewhat, suppose that A and B are traceless, and define · := 1 d tr( · ). We further define the 2k-spectral form factor with respect to a Hamiltonian ensemble E H by ..,i k =1 j 1 ,...,j k =1 e −i (λ i 1 + ··· +λ i k −λ j 1 − ··· −λ j k ) t (2.10) 2 For a more detailed analysis along these lines, see [23]. which depends solely on eigenvalue statistics for any ensemble E H . When the context is clear, we will often write R 2k (t) without a superscript. Spectral form factors have been studied extensively in the context of random matrix theory (for instance, see [22,24,25]). The GUE spectral form factor, and its features and time scales, are depicted in Figure 1. For the GUE, we rewrite Eqn. (2.9) as Notice that the dynamics of A(t)B(0) E GUE decouple from the details of the operators A and B. In particular, the only time dependence is due to the spectral form factor R E GUE 2 (t) which is an RMT quantity that diagnoses eigenvalue statistics. Indeed, Eqn. (2.9) has precisely the form that was advertised in Eqn. (1.1) in the introduction.
We can generalize Eqn. (2.11) in many ways, including to higher point functions at multiple times, to finite temperature, and so on. As an example, letting A and B be distinct Pauli operators, the infinite-temperature 2k-OTOC (i.e., out-of-time-order correlator) is, to leading order in large d [4] given by the 2k-spectral form factor as The time-dependence of the lower-order terms suppressed above similarly factorize as [spectral form factor](t) × [trace of operators] .
As emphasized in [4], the only essential feature in the vast simplication of correlation functions seen above is the Haar-unitary invariance of the Hamiltonian ensemble E H controlling time evolution. The GUE possesses an exact Haar-unitary invariance, and accordingly its time-dependent correlation functions solely depend on basis-invariant data of consistuent operators, and the joint eigenvalue distribution of the GUE. Unfortunately, the GUE is not a physical ensemble, since generic samples have N -body interactions and so are completely non-local. In fact, it appears that Haar-unitary invariance is at odds with the locality of interactions of a Hamiltonian ensemble. Indeed, if we consider a Hamiltonian ensemble comprising solely of Hamiltonians with, say, k-body interactions for k < N , then this ensemble cannot possibly be exactly Haar-unitary invariant. Thus it seems difficult to make contact with more physical ensembles Hamiltonians which possess few-body interactions.
To make progress, we need to change our perspective slightly. Instead of consider ensembles of Hamiltonians, we instead consider ensembles of unitaries generated by Hamiltonian time evolution. More precisely, consider an ensemble of Hamiltonians E H with measure dH. Now we construct an associated 1-parameter family of unitary ensembles E t = {e −iHt , H ∈ E H } (parameterized by t) with measure dU , such that for any function f (U ) we have Equivalently, the probability density of a unitary e −iHt in E t is the same as the probability density of H in E. At first, this appears to be merely a change of notation. For instance, we can equivalently rewrite the 2k-spectral form factor R E H 2k (t) in Eqn. (2.10) as and equivalently rewrite A(t)B(0) E GUE in Eqn. (2.11) as The crucial point is that the simplification of correlation functions explained above is only contingent on the (approximate) Haar-unitary invariance of E t . We will find that ensembles E H of local Hamiltonians which themselves are not (approximately) Haar-unitary invariant can still have their corresponding E t become (approximately) Haar-invariant at sufficiently late times, leading to a dynamical simplification of late time correlation functions. Spectral decoupling is due in part to the physics of scrambling.
To make the above ideas precise, we must: In the next section, we provide an answer to the first of these points, by reviewing and refining k-invariance [4]. The remaining points will be discussed throughout the remainder of the paper.

k-invariance and many-body chaos
In this section we introduce and discuss k-invariance as a measure of spectral decoupling in disordered quantum many-body systems. The goal will be to understand and quantify the ergodic nature of chaotic Hamiltonian evolution. Consider the ensemble of unitary time evolutions at a fixed time t, generated by an ensemble of Hamiltonians E H as with measure dU , as defined by Eqn. (2.13). The ensemble of Hamiltonians E H might be a spin system with quenched disorder, the SYK model, a random matrix ensemble, etc. One question that we can ask about the ensemble of unitaries E t is how rapidly and uniformly it spreads out over the unitary group U (d). A precise measure of the distance between E t and U (d) is provided by unitary k-designs. In words, if a unitary ensemble forms a k-design, its first k moments agree with the corresponding moments of the Haar unitary ensemble U (d). But as discussed in [4], k-designs fail to capture aspects of late time ergodicity even for random matrix ensembles (i.e., the GUE, GOE and GSE). An alternative to k-designs which encapsulates properties of random matrices is k-invariance. Next we will define k-designs and k-invariance and discuss their connection to correlation functions of a disordered theory.

k-designs
We begin with unitary k-designs, which are subsets of U If we instead consider a weighted subset of the unitary group, i.e. an ensemble of unitaries Such an ensemble E U ⊂ U (d) might be discrete or continuous. In the latter case, Eqn. (3.3) is upgraded to an integral with an appropriate probability density P (U ) as One might ask when averages over the ensemble E U look like averages over the full unitary group. An exact unitary k-design is an ensemble E U for which the k-fold channels are equal for all operators O acting on H ⊗k For the first moment, Pauli operators (or any basis of operators of B(H)) form an exact 1-design.
There are examples of exact 2 and 3-designs [26][27][28][29], but for higher k, exact constructions of unitary k-designs are not known. One might instead ask when an ensemble is merely close to replicating moments of the Haar measure. An approximate k-design is an ensemble of unitaries E U for which the distance between k-fold channels in the diamond norm is small, where the diamond norm is defined in Appendix A. We define approximate designs in terms of the diamond norm, but note that there are other definitions of an approximate design in the literature involving different norms. However, different norms bound each other up to factors of d ; see [30] for details in the context of k-designs.

Frame potential
A more tractable measure of approximate designs and the Haar randomness of an ensemble is the 2-norm on quantum channels, namely Usually, we consider the 2-norm on the space of operators B(H) acting on a Hilbert space H. However, in the present setting we are considering superoperators acting on B(H ⊗k ). These superoperators also have a natural 2-norm (more precisely, the 2-norm of the moment operators), which we discuss in Appendix A. A convenient way of representing Eqn. (3) is in terms of frame potentials. The k-th frame potential for an ensemble of unitaries E U is defined as [31,32] F (k) written here for a continuous ensemble equipped with a probability measure. The key fact is that we have Φ (k) As a result, we see that the frame potential F (k) E U is lower bounded by the corresponding frame potential computed for the Haar measure [32] F (k) where the Haar value is F For higher moments k > d, the Haar value is given by a known combinatorial expression. Closeness in the 2-norm, as seen by the frame potential, is a weaker notion of approximate k-design as defined in Eqn. (3.6). Nevertheless, the difference in frame potentials bounds the difference in k-fold channels up to factors of the dimension of the Hilbert space, as shown in Appendix A.
The frame potential has more recently been understood as a diagnostic of chaotic dynamics [4,[33][34][35][36]. Specifically, the k-th frame potential was related to operator-averaged OTOCs, making precise the connection between the chaotic decay of correlation functions and the approach to Haar randomness. We define the 2k-OTOCs averaged over an ensemble E t in Eqn. (3.1) by where dH is the measure over E H , B(t) = e iHt Be −iHt , and we take the expectation in the infinite temperature state. Note that these OTOCs can be naturally generalized to an ensemble of unitaries E U instead of just E t . The k-th frame potential is then related to operator-averaged 2k-OTOCs as [33] F (k) where we take the sum of each operator over an orthonormal basis of operators Since the frame potential is lower bounded by its Haar value F (k) E Haar = k!, the chaotic decay of OTOCs signifies the approach to randomness.

k-designs are not enough
In some models of stochastic, time-dependent evolution up to time t, such as random quantum circuits or Brownian circuits, it is known that for late enough times t the ensembles converge to approximate unitary k-designs [21,[37][38][39][40]. But for time-independent evolution by a disordered Hamiltonian this is not the case. Considering the ensemble E t = {e −iHt , H ∈ E H } of evolution by disordered Hamiltonians up to time t, there is a simple argument for why E t will never become Haar random at very late times [33]. The eigenvalues of a sample from the Haar unitary ensemble have the form e −iθ j for j = 1, ..., d, which experience level repulsion on the complex circle. However, sampling a Hamiltonian H from E H with eigenvalues λ j for j = 1, ..., d, the eigenvalues of e −iHt are e −iλ j t , which become Poisson distributed on the complex circle as t approaches infinity. Thus e −iHt for large t cannot look like a Haar random unitary.
However, at intermediate time scales, e −iHt can look approximately Haar random. Since the eigenvalues λ j of random Hamiltonians experience level repulsion, e −iλ j t still may experience level repulsion on the complex circle so long as t is not too large. Considering Hamiltonian evolution by a standard random matrix ensemble such as the GUE, we can analytically compute the frame potentials in terms of the spectral functions [4]. Hamiltonians drawn from the GUE become close to Haar random in 2-norm at an intermediate time scale called the "dip time" [4], see Figure 1. We believe that this phenomenon of becoming close to Haar random at an intermediate time holds more generally for chaotic quantum many-body systems, once the symmetries of the system are taken into account. Nevertheless, as evolution by random matrix Hamiltonians is not Haar random at late times, we require a new quantity to measure ergodicity and the onset of random matrix theory.

k-invariance
We now introduce k-invariance as a measure of the onset of random matrix behavior, and specifically spectral decoupling. We are interested in the ensemble of unitary time evolutions generated by disordered Hamiltonians, E t = {e −iHt , H ∈ E H }. For systems that break all symmetries, we expect that the ensemble will achieve (approximate) unitary invariance. This is motivated by a defining property of the GUE, namely that the measure on the space of Hermitian matrices is invariant under unitary conjugation, as discussed above. Given the ensemble E t , we define an ensemble which is unitarily invariant: E t = {e i(U HU † )t , H ∈ E H and U ∈ U (d)} where the U 's are Haar distributed and the H's are still distributed according to dH. We will often refer to E t as the "Haar'ed" version of E t . The distance between the ensembles E t and E t captures how basis-invariant the dynamics of E t are at the time t. We say that an ensemble E t is k-invariant if the k-fold channels of the two ensembles are equal (3.12) Given the time evolution of a chaotic quantum many-body system, we expect that the ensemble of unitary time-evolutions become approximately k-invariant at late times as information is scrambled and the dynamics become invariant under an arbitrary change of basis. Thus, we want to quantify an how close the two ensembles E t and E t are at a given time t. We might consider a strong notion of k-invariance in terms of the diamond distance between the ensembles which captures the distinguishability of the ensembles. In the interest of working with more tractable quantities, we will instead focus on the difference in frame potentials to understand Again, the two frame potentials quantify a 2-norm distance between the ensembles and the difference in frame potentials bounds the diamond distance up to factors of the dimension; see Appendix A. We will refer to F (k) Et as a Haar'ed frame potential. As we see, k-invariance as measured by the frame potentials defines a distance to unitary-invariance and will signify the onset of random matrix behavior in late time quantum many-body dynamics. In particular, an ensemble E t can become (approximately) k-invariant at a timescale t k-inv such that t 1-inv ≤ t 2-inv ≤ · · · , thus introducing a new hierarchy of timescales into quantum many-body chaos. Intuitively, t k-inv is the timescale when 2k-OTOCs with adjacent operators inserted t k-inv apart undergo spectral decoupling. More explicitly, a 2k-OTOC with t ≥ t k-inv will be approximately spectrally decoupled. Similarly, a 2k-point Keldyshordered correlation function of the form will be spectrally decoupled for t ≥ t k-inv . Here we have assumed that ρ 0 does not depend on the ensemble and so ρ 0 is not, for instance, a finite-temperature Gibbs state. It is possible to generalize the above discussion to finite-temperature OTOCs and finite-temperature Keldyshordered correlation functions. This utilizes a finite-temperature version of k-invariance, which we discuss later in this section. The definition of k-invariance in Eqn. (3.12) is weaker than the definition of k-designs in Eqn. (3.5). In fact, if a unitary ensemble forms a k-design, it is necessarily k-invariant. This is readily seen since if Φ (k)

1-invariant frame potentials
Given some ensemble of Hamiltonians time-evolutions E t , the question of whether the ensemble is k-invariant amounts to comparing the frame potential of E t to that of the invariant ensemble E t . By construction, the invariant frame potential can be computed in terms of spectral functions of the Hamiltonian ensemble.
Here we quickly review the calculation of the 1-invariant frame potential [4], which will lower bound F (1) Et with equality if and only if the ensemble is 1-invariant. We have (3.17) Integrating using the second Haar moment (Eqn. (A.16)), we compute the 1-invariant frame potential in terms of the 2-point form factor of E H as The early time decay of the function is ≈ R 2 (t) 2 /d 2 , and for a chaotic system we approach a late time value of F (1) Et ≈ 2 when t ∼ d. We derive the reproduce the 2-invariant frame potential in terms of spectral functions in Appendix B.

Relation to correlation functions
We have given a general definition of k-invariance in Eqn. (3.14), and in order to make the definition less abstract we now relate it to correlation functions of observables. Recall that we can exactly relate the frame potential of an ensemble to operator-averaged 2k-OTOCs as in Eqn. (3.11). From the definition of k-invariance, we have As such, approximate k-invariance can be interpreted in terms of the distance between correlation functions of the two ensembles E t and E t . We will now explore Eqn. (3.19) in detail in the case of 2-point functions and the onset of 1-invariance.

1-invariance and variance of correlators
The condition for 1-invariance may be written as the following difference of 2-point functions where P denotes the basis of generalized Pauli operators. Also, as in the previous section, · Et is the average of the correlator over the ensemble The brackets · Et are defined similarly, but with respect to E t . We can compute the averaged correlation functions A(t)B Et over the unitarily invariant ensemble explicitly. For Pauli operators A and B, A(t)B Et is only nonzero for A = B. For non-identity Paulis, A(t)A Et = (R 2 (t) − 1)/(d 2 − 1), and thus where have cancelled the identity contributions from F (k) Et and F (k) Et and summed over nonidentity Paulis P . As a sanity check, we note that which is our previously derived expression for the unitarily invariant frame potential in terms of spectral form factors.
Using the above expressions, we can understand the decay of the difference in frame potentials, and thus the onset of approximate 1-invariance, in terms of 2-point functions. First, note that the 2-point spectral form factor can be written as an average over nonidentity Paulis [4] (3.25) We define notation for the average over non-identity Pauli operators as (3.26) If we assume that the 2-point functions with A = B are approximately zero at all times can be written as the variance of 2-point functions Therefore, if the difference in frame potentials becomes small and A(t)B Et ≈ 0, then the variance of 2-point functions Var A(t)A Et is small.

Chaos from approximate 2-invariance
Here, we explain how the onset of 2-invariance gives rise to scrambling and late time operator growth. We will characterize scrambling from the decay of the 4-point OTOCs as well as quantum mutual information of subsystems of the time evolution operator. Additionally, we show that 2-invariance implies universal spectral behavior of late time operator dynamics. We note that thermalization by random Hamiltonian evolution was studied in [41][42][43][44][45], and these results hold for 2-invariant ensembles.

OTOC decay with 2-invariant Hamiltonians
To explore why 2-invariance gives rise to scrambling, we first look at OTOCs as a probe of chaotic dynamics. At infinite temperature, we consider the OTOCs A(t)BC(t)D Et . We want to show that if the ensemble E t is approximately 2-invariant, then OTOCs necessarily decay. We compute the full OTOC for the invariant ensemble E t in Appendix B, and for traceless operators we find to leading order in 1/d The spectral functions R 4 (t) and R 4,1 (t) are defined by where we note that R 4,1 (t) is generically complex. For non-identity Pauli operators with A and C not equal to B or D, the decay of the 2-invariant OTOC is entirely governed by the spectral 4-point form factor R 4 (t), as shown in [4]. Moreover, for more conventional OTOCs of the form A(t)BA(t)B Et , again assuming non-identity Pauli operators, we find Simply assuming a non-degenerate spectrum, the long-time average of the 4-point form factor is R 4 (t) ≈ 2d 2 . Then for 2-invariant ensembles, these OTOCs decay to a late time value of 1/d 2 . Moreover, in many examples we can explicitly compute the early time decay of the spectral functions. As such, 2-invariance of a disordered system implies the decay of OTOCs.

Information scrambling with 2-invariant Hamiltonians
Now we turn to study the scrambling of quantum information under evolution for 2-invariant Hamiltonians. In the following, we assume that the ensemble E t has reached a time where it is approximately 2-invariant. The setup we consider will be to study the scrambling by looking at the mutual information of the unitary time-evolution operator [33,46,47]. More precisely, we bipartition the input and output of e −iHt intoĀB and CD respectively, and takeĀ andB to be entangled with reference systems A and B, respectively. We then look at the 4-partite state on regions ABCD. By studying the mutual information between various subregions of the state, we probe the delocalization and scrambling of quantum information under e −iHt . This setup is similar to studying the entanglement of the Choi state of the operator [46], and has an operational interpretation in terms of the Hayden-Preskill decoding protocol [47,48].
Consider an initial state ρ = |ψ ψ|, where |ψ AĀBB is a pure state on a 4-partite system. The state can be depicted by We evolve theĀB subsystems as ρĀB(t) = e −iHt ρĀB e iHt , and partition the output of the time-evolution into CD, such that the evolved state ρ(t) is To be clear, the unitary e −iHt is expressed as a map HĀ ⊗ HB where HĀ ⊗ HB and H C ⊗ H D are different bipartitions of the same Hilbert space (i.e., HĀ ⊗ HB H C ⊗ H D ).
To determine whether information has delocalized and scrambled over the entire system, we want to compute the mutual informations I(A, BD) and I(A, C) in the state ρ(t). If I(A, C) is small for any small input region A, then information has delocalized and we cannot learn about the state on A by performing measurements on C. When I(A, BD) becomes large, the information in A has scrambled and we can reconstruct the state by acting on the BD systems alone.
To make the computation more tractable, we instead consider the Rényi-2 mutual information where S (2) A = − log tr(ρ 2 A ) and the other terms are defined similarly. We compute the Rényi-2 entropies by calculating the purities of subsystems of ρ(t) averaged over 2-invariant ensembles. For instance, considering tr(ρ BD (t) 2 ), we use the invariance of the Hamiltonians to Haarconjugate the expression as where Λ is the unitary time evolution operator e −iHt written in the diagonal basis and U is a Haar random unitary. We then average over U using the 4-th moment of Haar random unitaries.
Computing the quantity, we find that in terms of the purities of the initial state, we have to leading order in 1/d and similarly mutual information and trivially we also have tr(ρ A (t) 2 ) Et = tr(ρ 2 A ). Combining these, the full expression for the mutual information I (2) (A, BD) at early times is Early times : which is simply the mutual information of the initial state. At late times we find Late times : . (3.36) Now assume that we have maximally entangled inputs, ρ AĀBB ≈ ρ AĀ ⊗ ρB B , where we take A to be maximally entangled withĀ and take B to be nearly maximally entangled withB. Then tr(ρ 2 A ) = 1/d A and tr(ρ 2 B ) ≈ 1/d B , and at late times we find the the mutual information is nearly equal to its maximal value Late times : where n A is the number of qubits in A, indicating that reconstruction of the state A is possible by only acting on BD. Here we have assumed that d 1 and that the B and C systems are larger than their complements, We note that if the system is 2-invariant, we only really need that t 1 so that the spectral function R 4 (t) has decayed. We do not necessarily need to be probing exponentially long times.
Similarly, we can compute the mutual information between A and C and at early times the quantity is large, but at late times becomes small I (2) (A, C) ≈ 0, indicating that A and C are no longer entangled. We have Early times : . (3.39) and at late times Late times : .

(3.40)
Again assuming A is maximally entangled withĀ and B is nearly maximally entangled with B, we find that at late times I (2) (A, C) ≈ 0. In summary, for ensembles E t which become 2-invariant, we have that at time scales greater than t ∼ 1, indicating that 2-invariance is sufficient for the delocalization and scrambling of information.

Late time operator growth
Chaotic many-body systems exhibit a ballistic growth in the support of an initially local operator [49][50][51][52][53][54]. Here we show that many-body systems which become approximately 2invariant exhibit a universal behavior in the late time growth of operators. In particular, suppose we have a traceless operator O in the Heisenberg picture which evolves by which is summed over an orthonormal basis of operators. Also, γ A (t) is the support of the growing operator on a particular basis element A, and A∈P |γ A (t)| 2 = 1. We further denote O(t = 0) = O 0 , which is taken to be a non-identity Pauli operator. Then we can express the coefficients γ A (t) by and Following the operator growth calculation in Appendix E of [55], suppose we disorder average over an ensemble E t where the constituent Hamiltonian ensemble E H has no symmetry. For t greater than the approximate 2-invariance time, we have This captures a type of ergodicity of late time operator growth -namely, how an operator spreads itself uniformly across operator space. For a chaotic Hamiltonian ensemble with no symmetries, we expect that at late times of t ∼ d, This is the timescale when the initial operator O 0 has evenly spread itself around operator space. So before times of t ∼ d but after the approximate 2-invariance time, the spectral statistics of the Hamiltonian ensemble evidently capture the approach to ergodicity. We note that in the approximation for |γ O 0 (t)| 2 in Eqn. (3.45), which captures the decay of the support of O 0 (t) on the initial operator O 0 , we have assumed that the approximate 2-invariance time is shorter than the time scale where the support on all operators becomes uniform (t ∼ √ d).

Finite temperature
We can generalize our discussion of k-invariance to finite temperature by using the thermal frame potential defined and computed in [4,33], as well as an alternative definition better suitable to "physical" correlation functions. We begin with the former definition first. Recalling that the k-th frame potential can be expressed as the operator-averaged 2k-point OTOCs, we can also consider a thermal 2k-point OTOC Here, we have equitably distributed 2k interstitial factors of e −βH/2k , and then normalized by tr(e −βH ). This convention is packaged in the notation for the correlator, · β,1 , where the "1" subscript denotes that this is a thermal correlator of the "first kind." (A thermal correlator of the "second kind" will follow shortly.) Taking the square of the correlator and averaging over orthonormal bases of the A's and B's, we find the finite-temperature frame potential of the first kind: This frame potential has a pleasing form, which is why it is so-defined. Note that as β → 0, we recover the infinite-temperature frame potential. Of course, thermal correlators that arise in the study of physical systems 3 do not have the form of Eqn. (3.46). Rather, there is a single Gibbs state in the correlator so that it takes the form (3.48) which we call a thermal correlator of the second kind. Similarly taking the square of the correlator and averaging over orthonormal bases of the A's and B's, we arrive at the finitetemperature frame potential of the second kind: Similar to before, as β → 0, we recover the infinite-temperature frame potential. The finite-temperature frame potentials can help characterize k-invariance for thermal correlators. In particular, for j = 1, 2 (i.e., for thermal correlators of either the first or second kind), we have which is the finite-temperature version of Eqn. (3.20). The proof of this identity is a trivial modification of the infinite-temperature analysis in Appendix B.
≈ 0 3 We are referring to the placement of the e −βH/2k factors as being unphysical, not the lack of time order (which is also unphysical).
for t ≥ t k-inv , then we expect 2k-point finite-temperature Keldysh-ordered correlation functions will be spectrally decoupled for t ≥ t k-inv . Naturally, the Haar'ed finite-temperature frame potentials can be expressed solely in terms of spectral statistics. For the Haar'ed finite-temperature frame potential of the first kind, we have In a similar vein, the Haar'ed finite-temperature frame potential of the second kind can be written as F where we now define One can obtain similar expressions for the higher-order finite-temperature frame potentials as well. We expect that most of the results in this paper generalize to the finite-temperature settings.

k-invariance in random circuits
We end the section with an example where the k-invariance time is essentially exactly computable. Random quantum circuits (RQCs) are solvable models of strongly-coupled local unitary dynamics, and as such are a valuable resource for understanding many-body chaos. It is known that random circuits are rapid information scramblers [47,56], decouplers [57], and generators of randomness [37,38]. Random circuits are known to form approximate k-designs in depth O(n poly(k)), although it is likely that the dependence on k is linear [21], reaching the optimal lower bound O(nk). It has also been shown that random circuits on D-dimensional lattices form approximate designs in O(n 1/D poly(k)) depth [58]. As explained earlier in this section, a unitary ensemble forming a k-design is a sufficient condition for that ensemble to be k-invariant. The convergence of random circuits to k-designs directly implies that they become k-invariant. Furthermore, we motivated k-invariance as a measure of chaotic dynamics in time-independent Hamiltonian systems, as opposed to the strongly time-dependent dynamics of random circuits. Nevertheless, it is still instructive to compute the k-invariance times in a solvable model with local dynamics. t Figure 2. A random quantum circuit on n qudits of local dimension q. Each layer is comprised of parallelized 2-site unitaries, where each gate is randomly sampled from U (q 2 ).
The random circuits we consider act on n qudits, each of local dimension q, arranged in a 1-dimensional chain. The circuits are built out of layers of 2-site random unitaries, acting in parallel on even links at even time steps and odd links at odd time steps. The 2-site unitaries are each independently samples from U (q 2 ). For a diagrammatic depiction of the circuit, see Figure 3.4. Evolution to time t is given by t layers of the circuit, and the total unitary evolution is denoted by U t . The ensemble of the accumulated random unitary circuit up to time t will be notated as E RQC t . See [51,52,59,60] for a recent treatment of entanglement and operator growth in these models.
We start by discussing the first moment of the random circuit ensemble E RQC t . The k = 1 frame potential for random circuits is exactly 1 after a single time step, i.e. F 1-invariant for any n, q, and t > 0. The second moment is nontrivial. The frame potential for random quantum circuits was recently computed in [21] using a mapping to a statistical mechanics model [51,60]. The k = 2 frame potential can be written as where w(t, q, n) is a sum of time-dependent contributions, which exhibits an exponential decay in time since w(t, q, n) ∼ n/q 2t . The function w(t, q, n) can be computed exactly, but the key fact is that the frame potential is exactly expressed as a constant term plus an exponentially decaying function in time. Recalling that F (2) E RQC t is lower bounded by its Haar value of 2, this proves that random circuits form 2-designs. 4 4 We note that although the frame potential decays in log(n) time, defining an approximate design in terms of the diamond norm and bounding it in terms of the difference in frame potentials we pick up an additional factor of n. So the 2-design time is O(n) for local random circuits.
We can also compute the invariant frame potential by computing the form factors for RQCs. First note that the Haar values for the higher-point spectral form factors that appear in the 2-invariance calculation are R 4 = tr(U ) 2 tr(U † ) 2 E Haar = 2, R 4,1 = tr(U 2 )tr(U † ) 2 E Haar = 0, and R 4,2 = tr(U 2 )tr(U †2 ) E Haar = 2. Since random circuits will become 2-designs and thus become 2-invariant, the previous Haar values of the spectral correlators will be achieved by the RQC. In summary, the exact values of the form factors for Haar random unitaries are Haar: 57) which are achieved by the RQC at the 2-design time. It turns out that R 4,1 = 0 for all times under random circuit evolution. This fact is transparent for t = 1 when we have a disjoint tensor product of unitaries. But for general times a proof of this involves expressing R 4,1 as the partition function for a statistical mechanics model and seeing that all spin configurations are disallowed. Furthermore, the form factors R 4 and R 4,2 can be computed as the partition function of a lattice model, and the resulting expressions involve the same time dependence as the frame potential. For t > 0, the exact values of the form factors for local random quantum circuits are RQCs: where again w(t, n, q) ≈ n/q 2t , which we emphasize is an exactly calculable quantity for random circuits [21]. Using the expressions for the RQC form factors in Eqn. (3.58), we can compute the 2-invariant frame potential (see Eqn. (B.7) in Appendix B) to find With these expressions it becomes clear that the decay to 2-invariance precedes the decay to a 2-design, with a very short timescale in between. The decay to 2-invariance and 2-design in random circuits can be summarized as Unsurprisingly, the time scale that both F is of order log(n). Asking when the diamond norm distance between the ensembles is small, the 2-invariance and 2-design times are both O(n). We emphasize that approximate 2-invariance is achieved before the system becomes an approximate 2-design. That is, for some small tolerance ε, we have F E Haar = ε, even though the difference in the two times is parametrically suppressed. This is in line with the fact that the 2-invariance time should come before the 2-design time, as being a 2-design implies being 2-invariant. Since 2-invariance implies scrambling, as previously discussed, we expect a hierarchy t scramble ≤ t 2-inv ≤ t 2-design , where we have shown the latter inequality is a strict inequality for random circuits.

Connection to ETH
A related framework for connecting chaotic quantum systems with random matrix theory is provided by the Eigenstate Thermalization Hypothesis (ETH) [61][62][63]. There are many refinements of the original conjectures, but here we will only sketch one of the original formulations, more or less following [64]. Our main aim is to argue that ETH and k-invariance are distinct but compatible notions.
Consider a Hermitian operator A. If we have a single instance of a chaotic Hamiltonian with eigenstates are smooth functions of E and ω, S(E) is an energy-smoothed thermodynamic entropy, and R mn are random numbers with mean zero and unit variance.
To illustrate how the ETH ansatz is applied, consider the infinite-temperature 2-point function ) for a fixed Hamiltonian H. Then expanding in an eigenbasis of H, we have where A mn = E m |A|E n . Taking an infinite time average, we find that only the diagonal terms in the sum contribute, giving 1 d n |A nn | 2 . This is consistent with the diagonal ansatz in Eqn. (3.61), as A(E n ) ≈ A nn . Moreover, for a traceless operator, we expect the ensemble averaged diagonal matrix elements to have A nn E H = 0 and fluctuations A 2 nn E H ≈ 1/d. This is also consistent with the late time form of the 1-invariant 2-point function, A(t)A Et ≈ R 2 (t)/d 2 , which goes to 1/d.
The first term in Eqn. (3.61), known as the diagonal ansatz, helps explain features of thermalization and the equilibrium values of observables. The second term in Eqn. (3.61), known as the off-diagonal ansatz, aims to describe the dynamical approach to equilibrium. In our present context, the off-diagonal ansatz in Eqn. (3.61) will give a prediction for the decay in time of the 2-point function in Eqn. (3.62). We see that this agrees with the 1-invariant form and the decoupling of the matrix elements from the sum over energies in Eqn. (3.62) if the function f A (E, ω) is constant. Indeed, in energy windows ω E Th where E Th is set by a scale called the Thouless energy, f A (E, ω) is constant (although this can be subtle in systems with diffusive transport [65]). Such behavior in the off-diagonal matrix elements has been observed numerically [63,66]. Then it is consistent with k-invariance that in small energy windows (and thus at late time scales) ETH predicts decoupled forms of correlation functions. Aside from this consistency check, it would be interesting to more fully understand the precise interplay between k-invariance and ETH, and how the timescales involved in each compare with one another.

Spectral decoupling and symmetry
Our discussion of many-body chaos and the onset of spectral decoupling has thus far been centered on systems which break all symmetries. However, many physical systems do possess symmetries, such as time reversal symmetry, particle-hole symmetry, and so on. Here, we refine our analysis to accommodate symmetries. In particular, we analyze the extent to which spectral decoupling can occur in systems with symmetry, and diagnose this with a symmetrized form of k-invariance.

Symmetries in random matrix theory
Random matrix theory is meant to capture the universal eigenvalue statistics of 'generic' Hamiltonians, constrained by symmetry. For a system exhibiting features of quantum chaos, and which is also constrained by symmetry, a heuristic intuition is that its Hamiltonian behaves as if it was randomly sampled from some universal ensemble of Hamiltonians with appropriate symmetries.
Here, we emphasize the role of symmetry. Suppose we have a system with a Hamiltonian H, which generates unitary evolution via U (t) = e −iHt . If the system is endowed with symmetries, this means that there is a group G = {g} with unitary and anti-unitary representations V g such that [H, V g ] = 0 for all g ∈ G. The unitary U (t) acts on states in a Hilbert space H, which is itself endowed with an inner product structure. Indeed, the inner product of two states is invariant under unitary time evolution. Conversely, this invariance structure is the defining property of unitary evolution.
Then according to Dyson's threefold way [67], one can block diagonalize the Hamiltonian generating the time evolution such that the blocks are each either (i) a complex Hermitian matrix, (ii) a real symmetric matrix, or (iii) a real quaternionic matrix. The latter two possess different manifestations of time-reversal symmetry. For a review, see [68]. If the Hamiltonian is sufficiently 'generic,' it is expected that each block, necessarily falling into the categories (i), (ii), or (iii), will behave like a random matrix sampled from a corresponding universal ensemble. Indeed, Wigner and Dyson constructed such universal random matrix ensembles by considering Gaussian random matrices consistent with (i), (ii) or (iii). The resulting ensembles are referred to as the Gaussian Unitary Ensemble (GUE), the Gaussian Orthogonal Ensemble (GOE), and the Gaussian Symplectic Ensemble (GSE), respectively.
There are other situations in which the invariance structure of the Hilbert space under Hamiltonian time evolution is enriched. Such is the case for systems with fermions. Then Dyson's threefold way is expanded to the tenfold way of Altland and Zirnbauer [69]. In this paper, we only consider the threefold classification, although our analysis could be generalized to the extended ensembles.

Symmetric k-invariance and time-reversal symmetry
So far, we have only defined k-invariance for ensembles E t generated by Hamiltonians with no symmetries. In the presence of symmetry, we need to modify our definition of k-invariance. The idea is to find a modified version of the Haar'ed ensemble E t which we will call E sym t .
The ensemble E sym t will be compatible with the symmetries of E t , and we will use to quantify the distance between the two ensembles. We begin with the specific examples of time-reversal symmetry and particle-hole symmetry, and then build up to a general construction.
In the spirit of Dyson's classification, we want to consider many-body systems invariant under time-reversal symmetry. For such Hamiltonians, there is an antiunitary operator T which commutes with the Hamiltonian and squares to ±1, i.e. T 2 = ±1. If T squares to unity, then the Hamiltonian can be written as a real symmetric matrix with respect to a class of bases constructed from T , and which do not depend on H. For chaotic quantum systems in the T 2 = 1 symmetry class, the spectral statistics are expected to be that of the GOE, and thus at late times we expect E t to become orthogonally invariant, i.e. the measure over E t becomes invariant under conjugation by an orthogonal matrix. (This is in contrast to the previous setting with no symmetry, in which we expect unitary invariance.) For such systems, "orthogonal k-invariance" is determined by the distance between E t and the orthogonally Haar'ed ensemble where the O's are Haar distributed on O(d) and the H's are distributed according to dH. Indeed, the measure on E O t has the desired orthogonal invariance. If instead T 2 = −1 then we expect GSE-type spectral statistics and at late times we expect E t to become symplectically invariant, i.e. the measure over E t becomes invariant under conjugation by a symplectic matrix. Then symplectic k-invariance quantifies the distance between E t and where the S's are Haar distributed on Sp(d) (where we assume d is even) and the H's are distributed according to dH. Also, S D is the symplectic transpose, given by S D := JS T J −1 such that J is the canonical symplectic form where 0 is the d/2 × d/2 zero matrix, and 1 is the d/2 × d/2 identity matrix. For a given k, we compute Eqn. (4.1) using the orthogonally and symplectically invariant frame potentials. Recall that in the case with no symmetries, we compute the k = 1 unitary invariant frame potential for E t in terms of spectral functions using the Haar invariance of the integration measure. We can proceed similarly to compute the invariant frame potentials for time-reversal symmetric systems.
First, consider symmetric k-invariance for systems with T 2 = 1, where we look at the distance to orthogonal invariance. The k = 1 invariant frame potential for E O t is written as In the T 2 = −1 case, the symplectic invariant frame potential for k = 1 can also be computed similarly, this time using the second moment of the Haar measure on Sp(d), written explicitly in Eqn. (A.24) in Appendix A. The result is (4.7)

1-invariance and time-reversal symmetry
Let us explore the relation between approximate 1-invariance and 2-point functions in systems with time-reversal symmetry. Here we assume a multiqubit system, d = 2 n , and use the basis of Pauli strings. If T 2 = 1, then we expect the late time dynamics to achieve approximate orthogonal invariance, such that F (k) The difference in frame potentials may be written as the average over 2-point functions depending on whether the Pauli string A is even or odd, i.e. A T = ±A. We will simply refer to the Pauli strings as Pauli operators. The 2-point functions A(t)B with A = B vanish identically.
As a consistency check, we can write out the full expression for the orthogonally invariant frame potential. The d 2 −1 non-identity Pauli operators contain (d+2)(d−1)/2 even operators and d(d − 1)/2 odd operators. Therefore, 10) which is the same as the expression we derived above in Eqn. (4.6).
We would like to understand orthogonal 1-invariance in terms of 2-point functions. Assuming that the contribution from A = B correlation functions are negligible, we can write Since the operator average of the 2-point function A(t)A Et is given by the spectral form factor after some algebra we can write the condition for orthogonal 1-invariance as . (4.14) To leading order in d this is simply Therefore, achieving orthogonal invariance implies that at late times, there will be a residual variance in the 2-point functions. We will see this explicitly in our numerics for Pauli spin systems in Section 5, specifically in Figure 8. Although our analysis of this variance is for systems with T 2 = 1 time-reversal symmetry, a similar analysis holds for T 2 = −1 timereversal symmetry. Also, we note that if one did not know the system had time-reversal symmetry and went ahead computing the unitary 1-invariance, the late time value would be off by this additive factor of one, which is indeed what we observe in our numerics for T -invariant spin systems.

Orthogonal invariance from unitary invariance
In deriving the expression for orthogonal 1-invariance above, we assumed that the Hamiltonian was real, and thus e −iHt is a symmetric unitary. This glosses over an issue of basis dependence when computing the symmetric frame potentials. In general, making no assumptions about the symmetry of the Hamiltonian ensemble, the expression for the orthogonally-invariant frame potential is where R T (t) := tr(e iHt e −iH * t ) . (4.17) We note that when the Hamiltonian is of GOE class and commutes with an antiunitary operator T , there exists a set of bases (constructed independently of H) in which H is real (see Appendix C). In such bases, R T = d and the expression reduces to the orthogonal 1-invariant frame potential previously derived in Eqn. (4.6).
As a consistency check, we can compute R T (t) for the GUE, and find Plugging this into the expression for F More generally, if we wanted to check orthogonal invariance in the absence of time-reversal symmetry, we should use Eqn. (4.16) where R T no longer reduces to d in a T -invariant basis.

k-invariance and a Z 2 global symmetry
Now we consider a system with a Z 2 global symmetry. This example will be a good warm-up for our discussion of SYK in the next section. Suppose we have a Hamiltonian ensemble such that each Hamiltonian is block diagonal with respect to Z 2 symmetry operator. In particular, there are two blocks of equal size: We further suppose that each block H + and H − has no further symmetry. Also, the subscripts + and − are purely for convenience, and do not mean that the eigenvalues of H + and H − have a definite sign. Then the correctly Haar-symmetrized version of the ensemble E t is We can compute the invariant frame potential by (4.21) As it will be convenient to refer to later, we write this out explicitly as We can proceed by computing the first two terms using the second Haar unitary moment and the second two terms using the first Haar unitary moment, which gives where d + = d − = d/2 are the dimensions of the two sectors and the form factors R 2,± are those evaluated in each respective sector, i.e. involving only the eigenvalues in that part of the spectrum. In the first two terms are two decoupled versions of the normal Haar-invariant frame potential for each sector. We have also defined an interaction term, which is a mixed moment depending on both the + and − sectors where we see that the eigenvalues in each sector are summed over separately.

General construction
Having treated time-reversal symmetry and particle-hole symmetry on a case-by-case basis, we want to generalize our definition of k-invariance to accommodate general symmetries. To appreciate the need for a modification of k-invariance in the presence of symmetry, it is worth revisiting the definitions of our main ensembles of study. for α, β ∈ R. More generally, we have a system of linear constraints given by [H, V g ] = 0 for all g ∈ G. Let us denote the subspace of Hermitian operators satisfying these symmetry constraints by S G . In particular, In the trivial case for which G = {1}, we have S G = Herm(H). We now have the following picture: the ensemble of Hamiltonians E H with symmetry group G is supported on the subspace S G of Herm(H). This is notated as (4.30) A diagram can be seen in Figure 3.
To make this explicit, recall that the ensemble E H is equivalent to a probability measure P (H) dH on Herm(H). Let π : Herm(H) → S G be the projection map, and let δ S G (H) be the delta function of the subspace S G , satisfying for any function f (H). In the above equation, d(π * H) is the pushforward of the volume form dH to the subspace S G . Then we can interpret Eqn.   which is clearly no longer supported on S G alone. A schematic diagram is shown in Figure 4. Note that U † S G U is the subspace obtained from S G by conjugating each constituent Hermitian matrix by U † . So far, we have seen that an ensemble E H with (non-trivial) symmetry group G is supported on a subspace S G of the space of all Hermitian operators, whereas the Haar'ed version of the ensemble E H is supported on the entire space of Hermitian operators, and in fact is not supported on any proper subspace thereof. Thus, for ensembles with symmetry, we do not expect E H ≈ E H .
To ameliorate this issue, we appropriately modify the definition of E H so that it is consistent with the symmetries of the original ensemble E H . In other words, if E H has symmetry group G so that supp(E H ) ⊆ S G , we want to find a "canonical" definition of E G H such that supp( E G H ) ⊆ S G . To do so, consider a symmetry group G and the associated invariant subspace of Herm(H) which we have denoted by S G . There is a natural action of the unitary group U (d) on Herm(H) by U ·H = U HU † . Then we can consider a subgroup of U (d) defined by  This is precisely the subgroup of U (d) which leaves S G invariant under the group action. The idea is to "symmetrize" the ensemble E H with respect to Invar U (d) (S G ). More precisely, we define E G H by the probability distribution where dW is the Haar measure on Invar U (d) (S G ). Plugging in P (H) = p(H) δ S G (H) as per Eqn. (4.32), we find is itself a group and thus contains both W and W † .) In less mathematical terms, E G H is the most Haar'ed version of E H , consistent with its symmetries. We sketch a schematic diagram in Figure (5). Now we return to the ensemble E t . This ensemble of unitaries has probability measure P t (U ) dU = dP t . Defining g t : we have dP t := d(g t * P ) , (4.40) where d(g t * P ) is the pushforward measure of P with respect to g t . More simply, we have for any function f (U ). Then if E H has symmetry group G, the measure dP t (U ) has the form (4.42) We can naturally define E G t by the probability distribution and this probability distribution satisfies for ε small. Some other norm may be used instead, if desired. We say that E H has G as an approximate symmetry. While E H may not strictly lie in S G , the ensemble will be concentrated near S G . Considering an ensemble E t as usual, we can ask if the ensemble becomes k-invariant at late times. We would expect that the ensemble does not become k-invariant, but rather may become approximately symmetric k-invariant with respect to G. The reason is that if we consider the Haar'ed version of E H , namely E H , then the new ensemble will not be concentrated near S G . The simple reason is that Haar-averaging with respect to the unitary ensemble does not 'know' about the approximate symmetry group G. On the other hand, we expect that the G-symmetrized version of E H , namely E G H , is concentrated near S G . In summary, if E H and thus E t has an approximate symmetry group G, we may anticipate approximate symmetric k-invariance with respect to G. In Section 5.2 we will explore this numerically.

k-invariance in spin systems and SYK
In this section we examine k-invariance numerically in a number of disordered many-body systems, including disordered Pauli spin systems and various versions of the SYK model. We proceed by numerically computing the frame potential for the ensemble E t of a given disordered system at a fixed time. In particular, we generate a set of Hamiltonians by randomly sampling from the ensemble of disordered Hamiltonians, and use these samples as an approximation to the full Hamiltonian ensemble to compute F (k) Et . We then compute the invariant frame potentials by numerically computing the spectral form factors of the ensemble. Some additional details regarding our numerics are given in Appendix C.
The quantity we compute is the difference of frame potentials which we sometimes refer to informally as the 'distance to 1-invariance.' In computing this for systems breaking all symmetries, we take E t to be the unitarily invariant ensemble and use the 1-invariant frame potential in Eqn. (3.18). For time-reversal symmetric systems we use the invariant expressions derived in Section 4. But for SYK, which nontrivially realizes a particlehole symmetry due to the interplay with charge-parity sectors, we derive the appropriate 1-invariant expressions for F (1) Et below.

Pauli spin models
The simplest examples of disordered spin systems are Pauli spins with disordered interactions. We consider several versions with nonlocal and geometrically local interactions as well as systems which realize or break time-reversal symmetry. Numerics studying k-invariance in Pauli spin systems were first reported in [70].

Random Nonlocal 3-body
The first spin model we consider is a system of N all-to-all randomly coupled spins with 3-local interactions where each J ijk,αβγ is an i.i.d. Gaussian random variable with zero mean and variance 1/N 2 .
In the Hamiltonian, we sum i, j, k over all possible 3-body interactions between all triplets of sites, and sum α, β, γ over {x, y, z}, i.e. local Paulis at each site. The main feature of this Hamiltonian is that it consists of nonlocal (but still 3-local) random interactions which break all symmetries. In Figure 6, we plot F the distance to 1-invariance is zero at t = 0 is because E t=0 = {1}. So then the distance to 1-invariance must increase and then decrease. At late times, the dynamics appear approximately unitarily invariant as the distance drops to zero. Note that on the log-log plot we take the absolute value of the difference in frame potentials, and so the late time fluctuations are actually around zero. The late time floor appears to be zero within the precision of our numerics for all values of N we checked. Therefore, the random 3-local system appears to achieve approximate invariance at late times.

Random Nonlocal 2-body
The next Pauli spin model we consider is a system of N spins with all-to-all interactions via random couplings, summed over all possible 2-body Pauli interactions Here, each J ij,αβ is an i.i.d. Gaussian random variable with zero mean and variance 1/N , and we sum i and j over all pairs of sites and sum α and β over local Paulis at each site. First, we note that H 2-local realizes time-reversal symmetry and commutes with the antiunitary operator where K is the antiunitary complex conjugation operator. As T reverse the sign of single-site Paulis under conjugation, i.e., T σ α i T −1 = −σ α i , the 2-body Hamiltonian above commutes with T . Note that T 2 = 1 if the number of spins N is even, and thus H 2-local belongs to the GOE symmetry class. Also, T 2 = −1 if the number of spins N is odd, meaning H 2-local belongs to the GSE symmetry class.
Therefore, when computing F for the random 2-local model with N = 5, 6, 7 and 8 spins and observe that when accounting for the symmetry the system achieves approximate 1-invariance at late times.
In Figure 7, we plot the difference in frame potentials to show approximate unitary 1invariance and symmetric 1-invariance at late times. As in the previous example, the difference in the frame potentials increases at early times, and with a magnitude that increases with N . We observe the presence of a late time floor in the unitary invariance at a small but non-zero value. Once the symmetry is accounted for using the appropriate symmetrically invariant ensemble, we find that the late time floor goes to zero and the ensembles display approximate orthogonal or symplectic invariance at late times.
We note that a 4-local analog of H 2-local is also time-reversal invariant, and appears to achieve orthogonal/symplectic 1-invariance at late times better than the 2-local case (see [70] for some numerics demonstrating this). We believe that this is because, in some appropriate sense, the 4-local Hamiltonian is more "chaotic" than the 2-local ensemble.

2-point functions and form factors
Recall that there is an exact relation between the average of 2-point functions over the basis of Pauli strings and the 2-point spectral form factor [4] 1 This expression is true for any ensemble of Hamiltonians (even if that ensemble is concentrated on a single Hamiltonian). We previously argued that approximate unitary 1-invariance implies that the variance of 2-point functions is small in Section 3.1, and that approximate 1-invariance for time-reversal symmetric systems implies a larger variance in 2-point correlators.
In Figure 8   correlators are tightly clustered around R E H 2 (t), but in the case of the 2-local system the correlators spread out around the form factor in bands of varying Pauli weights.
One of the more interesting properties of 1-invariance is it implies that 2-point functions A(t)A Et are close to their spectrally decoupled versions A(t)A Et ≈ R 2 (t)/d 2 AA . For the 3-local system, which becomes unitarily invariant, all 2-point functions should become close to A(t)A Et = R 2 (t)−1 d 2 −1 . For time-reversal symmetric systems, the invariant 2-point function depends on if the operator is even or odd, but again are expressions in terms of R 2 (given in Eqn. (4.9)). In Figure 9 slower than correlation functions of multi-site Paulis. In the left panel of Figure 9, for the 3-local system the band above the approximately 1-invariant 2-point functions consists of 2point functions of single site Paulis. In the right panel of Figure 9, for the 2-local system the top band consists of correlation functions of 2-site Pauli operators. It would be interesting to better understand these features.

Geometrically local spin systems
Above we considered disordered spin systems with random few-body interactions that are not geometrically local. Here we consider more physical, geometrically local systems, namely a 1d system of N Pauli spins with all allowed nearest-neighbor and next-to-nearest-neighbor interactions We have permitted 3-body interactions in order to break time-reversal symmetry, opting to do so with 3-body terms instead of 1-body terms to avoid regimes where the system might manybody localize. As usual, we take each J i, αβγ and K i, αβ to be Gaussian random variables. In Figure 10, we plot F (1) for N = 5, 6, 7, 8 spins, observing a decay to approximate 1-invariance at late times as the difference in frame potentials becomes small. We note that the floor value is small, but does not appear to approach zero up to numerical precision in contrast with the case of the geometrically local 3-body spin system discussed previously.

Approximate symmetry in Pauli spin models
Now we consider a system with "approximate" symmetry. We return to Hamiltonians with non-local interactions for convenience. Consider a family of Hamiltonians of the form In Figure 11, we plot F (1) 1d with N = 6 spins, tuning λ from 0 to 1. We find that as we perturb away from the time-reversal symmetric point λ = 1, the system becomes increasingly unitarily 1-invariant at late times. Approximate orthogonal 1-invariance is achieved for all values of λ. This is because for λ = 1 the Hamiltonian is GOE class, but as we deform the Hamiltonian to break time-reversal symmetry the Hamiltonian becomes unitarily invariant, and unitary invariance is a sufficient condition for orthogonal invariance. Here orthogonal invariance in the time-reversal broken regime was computed using the full expression in Eqn. (4.16), evaluated in a T -invariant basis. For details on how to construct such a basis, see Appendix C.

SYK models
We now turn to studying k-invariance in versions of the Sachdev-Ye-Kitaev model (SYK), a system of N strongly-interacting all-to-all randomly coupled Majorana fermions χ i . The SYK model [17,19,20] has garnered much interest as a tractable model of quantum manybody chaos and captures features of low-dimensional black holes, exhibiting an emergent reparametrization invariance and maximal chaos [17,20]. The Hamiltonian for the q-local SYK model is written as where the Majorana fermions are Hermitian, χ i = χ † i , and obey the anticommutation relations {χ i , χ j } = δ ij . We disorder average the Hamiltonian over i.i.d. Gaussian random couplings J i 1 i 2 ···iq with zero mean and variance J 2 where J is a positive constant. The theory is solvable in the limit 1 βJ N , where an emergent conformal symmetry allows one to compute correlation functions of the theory [17,20].
For purposes of numerical computations, we can represent the N Majorana operators χ i for i = 1, ..., N in terms of Pauli strings, using the Jordan-Wigner transformation as an intermediate step. If we have a Hilbert space H N/2 j=1 C 2 with N/2 sites j = 1, ..., N/2, we can then define the N Majorana operators using the standard Clifford algebra representation as which are Hermitian and satisfy the desired anticommutation relation {χ i , χ j } = δ ij . Notice that each Majorana operator can be represented by a string of N/2 Pauli operators (including 2 × 2 identity matrices), and so can be expressed as a N/2 by N/2 complex Hermitian matrix. Thus, the Hilbert space of N Majoranas has dimension d = 2 N/2 , although it can be regarded as comprising N sites. Relevant for our discussion is that SYK has a particle-hole symmetry P which acts nontrivially on the spectrum [1,71,72]. Here, P is an anti-linear operator which can be written as complex conjugation times a product of Majorana fermions, specifically (5.9) P commutes with the Hamiltonian and squares to ±1. The N -dependence is [1] For N Majorana fermions, the statistics of the spectrum of H depends on N (mod 8), which can be understood by the action of P on the charge sectors. We understand that for time-reversal invariant systems, the ensemble E t may become k-invariant with respect to orthogonal or symplectic transformations. For fermionic systems with different charge-parity sectors, the non-trivial block diagonal structure of the Hamiltonian again necessitates a refinement of our definition of k-invariance. In what follows, we analyze how to correctly account for the particle-hole symmetry that acts nontrivially on the SYK spectrum, and determine the suitable version of symmetric k-invariance. SYK with N ≡ 2, 6 (mod 8) For N ≡ 2 or 6 (mod 8), P maps the even sector to the odd sector, i.e. even eigenstates are mapped to odd eigenstates. Therefore two sectors have the same spectra spec(H + ) = spec(H − ), and so the entire spectrum is doubly degenerate. We also note that P does not commute with either H + or H − and therefore regardless of whether P 2 = ±1, there is no symmetry within the sectors. This means each sector should have GUE statistics and become unitarily invariant.
The SYK Hamiltonian in this case has a block diagonal form H = diag(H + , H − ), where H + = H * − and so the diagonalizing unitaries U + and U − should not be taken as independent. We must take U + = U * − and then the expression for the symmetric 1-invariance in Eqn. (4.22) becomes where since the spectra of H + and H − are the same, we use H to denote a diagonalized single sector of the Hamiltonian. Integrating, we find where d = d/2 is the dimension of each sector, and the form factor for each sector is simply where we sum over only the eigenvalues in a sector, i.e. only the unique eigenvalues in the spectrum. Note that this is related to the naïve spectral form factor R 2 (t) computed from directly the Hamiltonian by a factor of four, R 2 (t) = 4 R 2 (t). We can write the Haar'ed frame potential for SYK in terms of the full form factor as (5.14) SYK with N ≡ 0 (mod 8) In the case of N ≡ 0 (mod 8), we have that the particle-hole operator maps the even sector to itself and the odd sector to itself. As P 2 = 1, the operator maps an eigenstate to itself and thus the entire spectrum of the Hamiltonian has no degeneracy. Moreover, the particle-hole operator is a symmetry of each sector and commutes with H + and H − . As P 2 = 1, we thus expect each sector to have GOE statistics and become orthogonally invariant.
To consider this case, we want to compute the Haar-invariant frame potential where the symmetry of the ensemble is a block diagonal orthogonal matrix of the form Et , the difference between the frame potential for the SYK model and its Haar-invariant ensemble at different number of Majoranas. Right: the difference between frame potentials using the invariant with respect to the appropriate symmetry class. Moreover, the blocks H + and H − are independent, so we need to integrate over Haar random O + and O − taken independently. The first two terms in Eqn. (4.22) for Haar random orthogonal matrices O + and O − , instead of Haar unitary matrices U + and U − , give our standard orthogonally-invariant expressions for each sector plus an interaction term: Now to compute the interaction term between the two sectors, we simply need to compute the coupled terms of the form and using the first orthogonal moment, we obtain In summary, for the SYK model with charge sectors with GOE statistics, we have the invariant frame potential with which to define symmetric 1-invariance In the case of SYK with N ≡ 4 (mod 8), again the particle-hole operator acts by mapping each charge-parity sector to itself. But as P 2 = −1, the operator cannot take an eigenstate to itself and thus must map each eigenstate within a sector to another eigenstate within the sector, meaning each charge-parity sector must be doubly degenerate. Furthermore, as P will commute with H + and H − , and squares to −1, each charge sector should have GSE statistics and achieve symplectic invariance at late times.
Et , the difference between the frame potential for the q = 6 SYK model and its Haar-invariant ensemble at different values of N . Right: the difference between frame potentials using invariance with respect to appropriate symmetry classes for the given N and q.
As the first symplectic moment is the same as the other Haar moments, we simply find the for the SYK model with charge sectors with GSE statistics, we have the invariant frame potential with which to define symmetric 1-invariance k-invariance for q ≡ 2 (mod 4) SYK Now we turn to considering k-invariance in SYK with q-body interactions with q ≡ 2 (mod 4). In this case, the particle-hole operator anticommutes with the Hamiltonian and thus the operator P maps each eigenstate to another eigenstate with the opposite signed energy [72]. The Hamiltonian can again be block-diagonalized into charge-parity sectors, where the statistics of the sectors depends on N .
SYK with q ≡ 2 (mod 4) and N ≡ 0, 4 (mod 8) In these cases the Hamiltonian is block diagonal, and each block is of Bogoliubov-de Gennes (BdG) type [72]. For the N ≡ 0 (mod 8) case, P 2 = 1 and is a bosonic operator, which means that the Hamiltonian in each sector is pure imaginary and skew symmetric. This is the extended ensemble BdG class D, meaning the Hamiltonian in each block lives in the tangent space of O(d). There is no special invariance class for these ensembles, and thus we simply achieve unitary 1-invariance in each sector independently.
Similarly, for the N ≡ 4 (mod 8) case, P 2 = −1 and each sector of the Hamiltonian belongs to BdG class C, meaning that each sector of the Hamiltonian lives in the tangent space of Sp(d). Again, there is no special invariance group and the invariant ensemble is simply given by two independent unitarily invariant sectors. The corresponding frame potential is: (5.21) SYK with q ≡ 2 (mod 4) and N ≡ 2, 6 (mod 8) In this case the particle-hole operator exchanges the two blocks H + and H − , relating the blocks to one another. However, the blocks are related with inverted spectra, H + = −H * − , where each block has GUE statistics. As the diagonalizing unitaries for each block are not independent, i.e. we have U + = U * − , the expression for the 1-invariance can be computed as where we find a contribution from the + and − sectors, but the interaction term must be computed with identical U + and U − . Computing this interaction term with the oppositely signed sectors, we find

Comments on k-invariance in JT gravity
It is difficult to analytically probe k-invariance due to the late timescales involved. In particular, there are few theories that enable precise, late time computations, for which the results are non-trivial. One notable exception is Jackiw-Teitelboim (JT) gravity [73][74][75]. We consider the version of JT gravity which includes the contributions of higher-genus spacetimes. Remarkably, this theory of quantum gravity in nearly-AdS 2 spacetimes admits a non-perturbative completion in terms of a matrix integral [76]. This means that the completion can be regarded as a theory of disorder-averaged Hamiltonians. There is a similar story for JT gravity in nearly-dS 2 [77] (see also [78]). We will presently focus on the case of nearly-AdS 2 gravity, for which Λ = −1.
Here, we briefly recount the relation between JT gravity and a disordered ensemble, and then explain the connection with k-invariance. The action for JT gravity on a surface M with boundaries ∂M can be written as Here, g µν is the metric, φ is a non-dynamical dilaton field, χ(M) is the Euler characteristic of M, R is the Ricci scalar, Λ is the cosmological constant, h is the induced metric on ∂M, and K is the induced scalar curvature on ∂M. The amplitude of a spacetime with boundary ∂M is given by the path integral where the sum over topologies entails summing over all surfaces with boundary ∂M. Since we are in two dimensions, ∂M is a union of circles, which are each parameterized by a (renormalized) length. If there are n such circles so that ∂M S 1 × · · · × S 1 n , then we can regard A(∂M) as a function of n lengths, which we denote by β 1 , ..., β n . Thus we write A(β 1 , ..., β n ). This amplitude can be written as an asymptotic series in the parameter e −S 0 , which is in fact a genus expansion of the spacetime(s) ending at the n boundaries.
Explicit details on all of these matters are given in [76] for the AdS 2 case (and in [77] for the dS 2 case; see also [78]). JT gravity on higher genus topologies can be expressed as a type of matrix integral. In particular, where we are integrating over d × d Hermitian matrices 6 (here, d can be regarded as the dimension of the Hilbert space) for a particular potential V , and taking a "double scaling limit" where d → ∞. Note that the potential V also contains an explicit dependence on d (i.e., it does not just depend on d via the matrix H), which enables the double scaling limit to be non-trivial. Also, we have used "∼" instead of "=" in Eqn. (5.27) to denote that the matrix integral agrees with the asymptotic series expansion in e −S 0 given by the JT path integral above, but contains additional non-perturbative corrections in ∼ e −e S 0 . It is intriguing that in JT gravity, the amplitudes for spacetimes are the same as spectral form factors in a certain random matrix ensemble. Since we are in the double scaling limit, the Hamiltonians comprising the ensemble are formally infinite-dimensional, although each has discrete level spacings. Nonetheless, amplitudes A still have an asymptotic expansion in e −S 0 , which essentially plays the role of 1/d as per our other analyses in this paper. For instance, we will have "large e S 0 factorization" like where the averages are taken with respect to the matrix integral. 6 In [76] and related works, L is often used instead of d. Now we consider coupling the gravitational theory to matter fields. At present, it is not clear how to introduce propagating matter fields into the matrix integral in Eqn. (5.27), but one can more modestly introduce probe matter fields in the path integral in Eqn. (5.26) and compute terms in the e −S 0 expansion. This was carried out in [79] for 1+1D conformal matter with conformal dimension ∆ (see [80] for earlier work). For times t S 0 , it was found that where the left-hand side is taken with respect to the gravitational path integral at finite temperature β, and [spectral quantity](β, t) is a purely spectral quantity depending only on t and β. This exemplifies spectral decoupling, analogous to 1-invariance.
In particular, the time-dependence is contained purely in the [spectral quantity](β, t) term, and the data of the operator O ∆ (packaged as a function f (∆) of the conformal dimension ∆) factors out. One can check that O ∆ (t)O ∆ (0) β does not factorize in this way at early times, and so the form of Eqn. (5.28) is indeed novel. The calculational techniques at hand are not sharp enough to determine the precise timescale of the onset of this (approximate) 1-invariance-like behavior -we only know that t 1-inv ≤ O(S 0 ).
There are some important differences with 1-invariance, however. The first is that our formulas for 1-invariance entail tr(O ∆ (0) O ∆ (0)) which is not well-defined for 1+1D conformal matter, and so our prediction for the form of a 2-point function approximately satisfying 1invariance (after a specified timescale) cannot literally hold. Relatedly, our derivation of the form of a 2-point function approximately satisfying 1-invariance is no longer valid, since the derivation involves manipulating traces of a product of operators which are not trace class. Nonetheless, if nearly-AdS 2 JT gravity coupled to conformal matter turns out to have a nonperturbative completion in terms of a (multi-)matrix integral, possibly with additional fields, then the 1+1D conformal matter theory may become resolved in some microscopic manner. In such a scenario, our k-invariance formulas may be more precisely realized.
Even in the absence of speculation, Eqn. (5.28) is already rather suggestive. The equation advances our point of view that the late time dynamics of correlation functions in highly chaotic many-body theories (such a gravity) are captured by purely spectral quantities. In our quantum gravity example here, it appears that the fluctuations of spacetime dominate the late time dynamics of matter correlators, since the [spectral quantity](β, t) in Eqn. (5.28) only knows about the spectrum of the pure gravitational theory. It would be interesting to understand spectral decoupling in variations of JT gravity with other symmetries, as in [81], and then find analogs with symmetric k-invariance.
Furthermore, it would be appealing if other theories of gravity, beyond JT in nearly-AdS 2 and nearly-dS 2 , can be expressed as disorder averaged theories. Then perhaps there is some appropriate analog of k-invariance in such theories.

Discussion
We have quantified spectral decoupling by a detailed analysis of k-invariance, and its refinement to systems with symmetry. Systems with quenched disorder which undergo approximate spectral decoupling adopt universal dynamics, captured by the joint eigenvalue statistics of the Hamiltonian ensemble. In such systems, the late time physics of OTOCs, Keldysh ordered correlation functions, and operator growth are captured by appropriate spectral forms factors. Thus spectral decoupling provides a bridge between recent diagnostics of quantum many-body chaos in terms of operators and correlation functions, and more traditional diagnostics in terms of random matrix eigenvalue statistics.
Our numerical and analytic evidence suggests that chaotic quantum many-body systems with non-local q-body interactions and quenched disorder become approximately k-invariant at late times, for k small relative to the number of sites N . The timescale of the onset of approximate k-invariance depends on the system, choice of norm, and value of k, but appears to be at most poly(N ). For random quantum circuits, the timescale of k-invariance can be quantified more precisely. Random circuits are exactly 1-invariant, and become 2-invariant just before the 2-design time. It would be desirable to have more precise calculations of the 1-invariance time in a chaotic system, either by analytical or numerical methods. Further study of k-invariance for higher k would also be interesting.
There are many natural questions and directions suggested by this work. Foremost is understanding the range and scope of Hamiltonian systems with quenched disorder which experience k-invariance, and types of spectral decoupling more broadly. This requires both more extensive numerical (and possibly analytic) understanding of spectral decoupling in systems with non-local interactions, as well as more detailed study of spectral decoupling in systems with geometrically local interactions. Quantifying the class of systems which undergo some form of spectral decoupling would provide a sharper boundary between stronger and weaker forms of quantum many-body chaos.
Although we have done an initial investigation of spectral decoupling in disordered manybody systems with geometrically local interactions and found positive results, there is more to be understood. The number of random variables in the ensemble may play a key role. Also, there may be a more refined notion of spectral decoupling and k-invariance in particular, that is better suited to geometrically local systems. This could enable more precise contact with disordered systems in nature. and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.

A Random unitaries and operator averages
In this appendix we provide a few definitions for the quantum information theoretic notions used in the paper, as well as derive some identities for averages over Pauli operators.

k-fold channels and unitary designs
The k-fold channel with respect to an ensemble of unitaries E U of an operator O was given in Eqn. (3.3) for a discrete ensemble. For a continuous ensemble, a subset of the unitary group equipped with a probability measure, the k-fold channel is Again, we say that an ensemble E U forms a unitary k-design if the k-fold channel of E U and the Haar ensemble are equal Φ (k) For p ≥ 1, the superoperator norm of a quantum channel is defined as where O p := tr|O| p 1/p is the Schatten p-norm of the operator O. We can now define the notion of distinguishability of quantum channels with the diamond norm, defined as where I d is the identity channel on an ancilla system of dimension d.
The diamond norm allows us to specify when an ensemble is close to being Haar random. We say that an ensemble of unitaries E U forms an -approximate k-design when the k-fold channels satisfy Φ (k) A review of different operator norms and the various definitions of approximate designs and their relation to one another is given in [30].
In the present paper, we focus on a weaker measure of approximate design called the frame potential, where the k-th frame potential for an ensemble E U defined by [31,32] F (k) which is lower bounded for any ensemble as F (k) The difference in frame potentials F (k) E Haar corresponds to the 2-norm distance to Haar randomness, specifically [32] Φ (k) where here we mean the operator 2-norm of the channels as operators acting on H ⊗2k , explicitly written as Φ (k) E U , the k-th moment operator of the ensemble: The frame potential is related to an approximate unitary k-design in Eqn. (A.4), where the difference in frame potentials bounds the difference in channels as [21] Φ (k) As we discussed, the notion of an approximate design is suited to capture the ergodicity of ensembles of time-dependent evolutions, such as random circuits or stochastic Hamiltonians, but for evolution by disordered Hamiltonians in Eqn. (3.1) we should instead focus on kinvariance. We defined k-invariance as the distance between the ensemble of time evolutions E t and the unitarily invariant ensemble E t . The frame potential for E t is lower bounded by that for E t as [4] F (k) Thus the difference in frame potentials quantifies a distance of the ensemble to unitary invariance. Similar to our discussion of approximate designs, we could also consider the distance between the k-fold channels over the two ensembles as a stronger notion of k-invariance but the difference in frame potentials bounds the difference in k-fold channels of the ensembles E t and E t as Φ Next, we establish that the difference in frame potentials is related to the 2-norm of the channels, by viewing the channel as an operator on H ⊗2k , i.e. the k-th moment operator of the ensemble in Eqn. (A.6). The difference in k-th moment operators of the two ensembles is Computing the 2-norm, and using the invariance of the E t ensemble, we find the first two moment of the orthogonal group The general expression for integration of monomials of Haar random symplectic matrices 7 of dimension d × d is given by [83,85] Sp(d) where ∆ J σ is the index contraction with respect to a pair partition, similar to that for the orthogonal group except with a symplectic J inserted in the contraction, and Wg Sp is the symplectic Weingarten function. Let us consider the canonical symplectic form J, defined by where 0 is the d/2 × d/2 matrix where every entry is zero, and 1 is the d/2 × d/2 identity matrix. For i, j = 1, ..., d, we can equivalently write J ij = δ i+d/2,j − δ i,j+d/2 . The first two moments of the symplectic ensemble are Here, S D := JS T J −1 , so that S D S = SS D = 1. 7 We note that the unitary symplectic group is often defined with even dimension 2d, as the intersection Sp(2d, C) ∩ U (2d). Here, for convenience and ease in comparing expressions, we denote the symplectic group as Sp(d) keeping in mind that d must be taken to be even.

k-point spectral form factors
For convenience, we explicitly define and write out the form factors that appear in many of the second moment calculations in the paper: Note that R 4,1 (t) is generically complex unless all eigenvalues of H come in pairs ±λ.
B More on k-invariance and correlators

k-invariant correlators
As we mentioned in Section 2, the central property of certain random matrix ensembles which allows us to compute correlation functions from spectral quantities is the invariance of the measure. For k-invariant Hamiltonians, the same property holds and we can compute correlation functions using the approximate invariance of the measure on our set of Hamiltonians. For instance, the invariant 2-point function A(t)B with A(t) = e iHt Ae −iHt and averaged over H ∈ E H . For traceless operators we have is the 2-point spectral form factor averaged over the ensemble. The above expression only requires the first moment and thus follows from 1-invariance of the ensemble.
We can also consider generic out-of-time-ordered correlation functions A(t)BC(t)D , where again we average over H ∈ E H . For traceless operators A, B, C, and D we find the full expression AD BC + AB CD ADCB , in terms of the spectral functions defined above. For OTOCs of the form A(t)BA(t)B , where both A and B are non-identity Pauli operators and assuming A = B, we find the above expression reduces to .

Approximate 2-invariance from OTOCs
We now review the calculation of the k = 2 invariant frame potential done in [4], and present a rederivation in terms of invariant correlation functions. As usual, consider an ensemble E t of unitary time evolutions to time t by an ensemble of disordered Hamiltonians. We can compute the second frame potential for the unitarily invariant ensemble E t as As the measure over E is defined to be unitarily invariant, we take H → U HU † and average over random unitaries using the fourth moment. We then find in terms of the spectral functions disorder-averaged over the ensemble E H . This expression was derived in Appendix C of [4]. We want to explore approximate 2-invariance from the perspective of the constituent 4-point functions. As described in Section 3.1, we understand the onset of approximate 1invariance as the decay of generic 2-point functions, and their closeness to the average 2-point function at late times. But as approximate 2-invariance is sufficient for several definitions of scrambling, we should attempt to precisely formulate its onset as the late time behavior of OTOCs.
Recall the for any ensemble of unitaries, we can write the frame potential as an average of OTOCs as F where we sum over all Pauli operators, and · Et is the ensemble averaged correlator. We want to evaluate this for the invariant ensemble E t . First we separately consider the terms in the sum with any of the Paulis being the identity. We find where the first term comes from the OTOC with four consituent identity operators, there is no contribution from OTOCs with three identity operators, the second term comes from OTOCs for which A = C = 1 or B = D = 1, the third term comes from the other 4 possible contributions from two of the operators in the OTOC equaling 1, and the fourth term comes from a single identity operator in the OTOC. The last term is summed over non-identity Paulis P . Here we have used that We explain how to compute these sums in the next part of the Appendix. Already, we see something interesting. Recall that at the dip time the GUE forms a k-design [4], for which the frame potential equals 2. This can be thought of as the contribution from the OTOCs when A, C = 1 and when B, D = 1. At late times, the GUE frame potential is equal to 10. We can already see where part of this arises, since at late times the double and single identity contributions give 2 and 4. However, there are still ∼ d 2 OTOCs that give a late time contribution in the final sum. But this already tells us that almost all of the (d 2 − 1) 4 terms in the remaining sum over non-identity Paulis, when summed over, contribute at subleading order in 1/d. We knew some of the contribution at late times had to come from 2-point functions and some from 4-point functions as the late time invariant frame potential is F The contributing correlators at the dip time are just AC and BD , i.e. just the timeindependent 2-point functions. At late times, the contributing OTOCs are these two 2-point functions, as well as the four 2-point functions of the form A(t)BC(t) . Moreover, each of the non-identity 4-point OTOCs gives a contribution at late times from the 4-point form factor as R 4 (t)/d 4 . At late times this goes as ∼ 1/d 2 . Squaring and accounting for the 1/d 2 out front, we get a 1/d 6 suppression. Over the d 8 correlation functions we sum over, d 6 of them contribute at leading order, thus giving the constant contribution to the 2-invariant frame potential. But we can also arrive at the answer explicitly. We derived the invariant 4-point function A(t)BC(t)D Et in Eqn. (B.2), expressed in terms of spectral form factors. Thus it remains to square the quantity in Eqn. (B.2) and explicitly sum over Paulis A, B, C, D ∈ P . This requires summing squared 4-point functions over non-identity Paulis, e.g. | ABCD | 2 as well as terms like ABCD DBCA . Making use of some explicit Pauli sums given in the next subsection, we can compute these and simplify the entire expression for the invariant frame potential in terms of the OTOCs. In doing so we recover the expression in for the full 2-invariant frame potential in Eqn. (B.7).

Some explicit Pauli sums
Consider a system with N qubits, and let d = 2 N . We would like to compute where σ 0 = 1. Then M j = 4 j−1 , and we have the recursion relation where we have used M N j = d 2(j−1) . Thus we have, for instance, and so on.

Relation to 2k-point correlation functions
To get a better understanding of symmetric k-invariance, we explain its relation to correlation functions. Consider the operator From which it follows that tr(S (k) † S (k) ) = F where W πcyc is the cyclic permutation operator on H ⊗2k . The operator W πcyc is unitary, and so W πcyc W † πcyc = W † πcyc W πcyc = 1. Using the above basis, we can rewrite S (k) as (0)) . Then we have 0 ≤ tr(S (k) † S (k) ) (B.21)  Therefore, if E t becomes approximately symmetric k-invariant, all 2k-OTOCs with respect to E t are close to all 2k-OTOCs with respect to E sym t .

C More on numerics
Here we give some additional details on how we numerically evaluate the frame potential for an ensemble of disordered Hamiltonians. As discussed in [4], for a finite ensemble E t , the k-th frame potential receives 'diagonal' contributions from the i = j terms. For uniform weights, each diagonal term contributes d 2k /|E H | 2 , and the sum of all such terms gives d 2k /|E H |. In the infinite ensemble size limit, these terms yield a vanishing contribution. But in performing numerics at finite |E H |, they can obscure interesting late time physics. Thus, in our numerics we subtract these diagonal terms and compute a modified frame potential with uniform weights p i = 1/(|E H |(|E H | − 1)). We numerically evaluate the invariant frame potentials F (k) Et using the defined invariance of the measure dH to relate the quantity to spectral functions, and simply numerically compute the spectral functions (such as R 2 (t)) for the ensemble of disordered Hamiltonians at some time t.
To characterize 1-invariance numerically, we generate an ensemble of Hamiltonians, compute the frame potential F (1) Et for the ensemble of time-evolutions at a time t, remove the 'diagonal' contributions, and then compute the invariant frame potential F (1) Et via the spectral form-factor R 2 (t). Then we calculate the distance to 1-invariance F

Constructing T -invariant bases
Here we describe a procedure for generating a T -invariant basis, which can be constructed solely from the antiunitary time-reversal operator T [15]. These are a class of bases in which time-reversal symmetric Hamiltonians commuting with T can be written as real matrices.
Consider a d×d Hermitian matrix which commutes with the antiunitary operator T , with T 2 = 1. The procedure is: (i) Generate d linearly independent vectors {|φ i }. It suffices to take a set of random vectors. (ii) Symmetrize the first vector by defining |ψ 1 = |φ 1 + T |φ 1 , such that T |ψ 1 = |ψ 1 , and then normalize the vector to unit norm ψ 1 |ψ 1 = 1. (iii) Use Gram-Schmidt orthogonalization and symmetrize with respect to T to construct a set of invariant orthonormal basis vectors. Take |φ 2 and define a vector orthogonal to |ψ 1 as | φ 2 = |φ 2 − ψ 1 |φ 2 |ψ 1 . Then T -symmetrize and define |ψ 2 = | φ 2 + T | φ 2 , and normalize the vector to unit norm. Repeat this procedure by taking the random vector |ψ j , orthogonalizing with respect to all i < j, symmetrizing so that T |ψ j = |ψ j , and rescaling to unit norm. The result is a T -invariant orthonormal basis {|ψ i }. With respect to this basis, the Hermitian matrix will be real, H = H * .