Many-body quantum chaos: Analytic connection to random matrix theory

A key goal of quantum chaos is to establish a relationship between widely observed universal spectral fluctuations of clean quantum systems and random matrix theory (RMT). For single particle systems with fully chaotic classical counterparts, the problem has been partly solved by Berry (1985) within the so-called diagonal approximation of semiclassical periodic-orbit sums. Derivation of the full RMT spectral form factor $K(t)$ from semiclassics has been completed only much later in a tour de force by Mueller et al (2004). In recent years, the questions of long-time dynamics at high energies, for which the full many-body energy spectrum becomes relevant, are coming at the forefront even for simple many-body quantum systems, such as locally interacting spin chains. Such systems display two universal types of behaviour which are termed as `many-body localized phase' and `ergodic phase'. In the ergodic phase, the spectral fluctuations are excellently described by RMT, even for very simple interactions and in the absence of any external source of disorder. Here we provide the first theoretical explanation for these observations. We compute $K(t)$ explicitly in the leading two orders in $t$ and show its agreement with RMT for non-integrable, time-reversal invariant many-body systems without classical counterparts, a generic example of which are Ising spin 1/2 models in a periodically kicking transverse field.


I. INTRODUCTION
RMT was introduced into physics in the 1950s by Wigner [1] for providing a statistical description of nuclear resonance/excitation spectra. It should be intuitively clear that a system consisting of a few tens of nucleons coupled via short and long-range interactions is complicated enough that a successful description of experimental spectral fluctuations in terms of an ensemble of random Hamiltonians with independent stochastic matrix elements is not that surprising. An example of a robust phenomenological measure of fluctuations is the statistical variance of the number of energy levels in an interval of fixed length ∆E which, in RMT and experimental nuclear spectra [2], grows as ∼ log |∆E/ρ| (known as spectral stiffness), rather than ∼ ∆E/ρ as in the Poissonian random spectrum (ρ is the average density of states). The atomic spectra observed already by 1960 exhibited the so-called 'level repulsion' which can be quantitatively explained [3] with Wigner's RMT. However, in the early 1980s a much more surprising fact has been revealed, namely that RMT also works extremely well for capturing spectral fluctuations of simple singleparticle systems whose corresponding classical dynamics are completely chaotic, such as dispersive (Sinai) billiards or hydrogen/Rydberg atoms in external magnetic or microwave fields. These observations [4][5][6], termed as the quantum chaos conjecture which has been concisely * Corresponding author. email: tomaz.prosen@fmf.uni-lj.si stated in [7], have driven the field of quantum chaos for decades. The first, partial explanation for the success of RMT in simple chaotic systems came from Berry's semiclassical (small effective ) calculation [8] of the spectral form factor K(t) in terms of a double sum over classical unstable periodic orbits, which we shall explain below. K(t) is defined as a Fourier transformation of the two-point correlation function of the spectral density ρ(E) = j δ(E − E j ), {E j } being the energy spectrum (1) whereρ = ρ(E) E and . . . E represents local energy average over an energy shell (say of width ∆E) containing many levelsρ ∆E 1, in case of autonomous (timeindependent) systems. In case of periodically driven, i.e. Floquet systems that we shall discuss later in this paper, the average over the full range [0, 2π) of quasienergies -eigenvalues of unitary Floquet (one-period) propagator U is normally considered as physical properties are not expected to depend on the particular value of quasi-energy. A fruitful intuition stems from observation that K(t) characterizes all pair-correlation properties including the level repulsion and spectral stiffness, since in RMT {E j } can be considered as a fictitious one dimensional (Dyson's) gas with a logarithmic pairwise interactions [9].
For integrable systems, possessing a complete set of conserved quantities, the energy spectrum {E j } is conjectured [10] to represent a Poisson random uncorrelated sequence, so the spectral form factor (1) can be exactly computed as K(t) ≡ t H = 2π ρ = const for all t > 0, and thus provides a clear discriminator between integrable and chaotic systems, since for the latter K(t) ∝ t, in agreement with explicit predictions of RMT, as we shall explain in Subsect. I A below.
A clear heuristic derivation of RMT spectral form factor K(t) for classically strongly chaotic (hyperbolic) systems from semiclassical periodic orbit theory, starting from Berry's diagonal approximation [8], upgraded to second order in t by Sieber and Richter [11,12], and finally completed to all orders in a tour de force by Müller et al. [13,14], has been arguably the main accomplishment of the field of quantum chaos of single or few particle systems. Nevertheless, a rigorous proof of the quantum chaos conjecture has so far only been possible for a much more abstract class of single-particle systems, specifically for mixing quantum graphs [15,16]. These semiclassical periodic-orbit approaches have a natural generalisation to a quantum many-body problem for bosons when the number of quanta per mode is large [17][18][19][20].
However, RMT has also been found to excellently describe spectral fluctuations in the simplest, say low-dimensional and locally interacting, non-integrable many-body systems where local degrees of freedom have no classical limit at all, such as spins 1/2, qubits, fermions etc. [21][22][23][24] and where no semiclassical or mean-field approach can be applied. Due to such phenomenological success, RMT statistics of level spacings is nowadays used essentially as a definition of the so-called quantum chaotic or ergodic phase (see e.g. [25][26][27][28][29][30]). Moreover, the ergodic phase has been intensively theoretically investigated in recent years and its most concise characterization is provided by the so-called eigenstate thermalization hypothesis (see e.g., [25,31]). Nevertheless, there has so far been no proposition of the underlying dynamical (or microscopic) mechanism (such as unstable periodic orbit pairings in semiclassical chaotic systems discussed above). Recent studies of out-of-time-ordered correlations in many-body systems, some of which establish exponential growth (in particular in 0+1 dimensional systems such as the Sachdev-Ye-Kitaev model), have no clear connection to Lyapunov instability as it is understood in classical dynamical systems theory, and is the only mathematically meaningful definition of chaos. This has to do with the lack of the concept of classical orbits and the corresponding unstable (non-linear) equations of motion which result in sensitive dependence on initial conditions (e.g. the butterfly effect). In short, the concept of orbits and Lyapunov chaos does not make sense at ' ∼ 1'. One thus urgently needs alternative concepts which would enable one to explain the surprising success of RMT in simple many-qubit systems.
Providing one such concept is the main objective of this article. We identify a coherent structure, in a class of generic many-body quantum systems with the lowest, two-dimensional local Hilbert space (qubits, or spins 1/2), which is responsible for building up level (spectral) correlations. Expanding K(t), which is written as the product of two traces of the quantum mechanical propagator (see Subsect.I A), in the computational spin basis and writing it in a discrete-path-integral like fashion, we find that the leading contribution comes from constructive interference which corresponds to a partition function of a classical one-dimensional Ising model. Furthermore, sub-leading contributions can be interpreted as a family of partition functions of so-called twisted Ising models which are classified using a novel diagrammatic technique. In terms of this expansion the leading contributions to K(t) are shown to exactly correspond to RMT results K(t) 2t for times longer than a certain crossover time t * , while the non-universal results at t < t * are shown to reproduce numerical data extremely well. The time-scale t * , which scales logarithmically with the system size, can be interpreted as a quantum many-body analogue of the Ehrenfest (or Thouless) time. Finally, we identify the non-semiclassical analog of the Sieber-Richter pairing mechanism [11] and exactly reproduce the sub-leading RMT term −2t 2 /t H as well.
A. Spectral form factor in Floquet systems and periodic orbit theory In order to address the setup with a minimal amount of inessential technical complications we decide to study periodically driven (Floquet) many-body systems in the absence of any conserved charges or unitary symmetries. Even in single-particle context these are the minimal models of quantum chaos [9] and correspond to Dyson's circular ensembles of random unitary matrices [32].
In this subsection we define the main object of our study, namely the spectral form factor for Floquet systems, and for comparison with the main derivation in Sect. II of our paper, outline the key steps of historical semiclassical derivation of the RMT form factor in terms of periodic orbit theory. For a unitary one-period Floquet propagator U we write the eigenphases ϕ n and eigenvectors |n as U |n = e −iϕn |n , n = 1, . . . , N , where N denotes the dimension of the Hilbert space. The spectral density (1-point function) is now defined as and is normalized to a unit mean level density Locally averaged density is expected to be ϕ-independent which makes Floquet systems particularly appealing for studying spectral fluctuations. These are encapsulated in the connected (2-point) spectral correlation function which is, again, expected to be homogeneous (ϕindependent). An equivalent, and very convenient quantity is the spectral form factor K(t), t ∈ Z, defined as an appropriately scaled Fourier transform Finally, one writes where ... represents an appropriate additional averaging, either over local windows of time t (moving time average) or over an ensemble of similar systems, which is needed since the spectral form factor (5) is not a selfaveraging quantity [33]. For circular random matrix ensembles [32] which are expected to model Floquet systems in RMT (orthogonal/unitary ensemble (OE/UE) for systems with/without time-reversal or more general anti-unitary symmetry) the spectral form factor up to Heisenberg time, t < N , reads Note that exactly the same expressions hold as well for Gaussian ensembles of RMT which model timeindependent systems. For Floquet systems with a well defined classical limit, where the motion is hyperbolic (chaotic) everywhere in the phase-space, one can write tr U t in terms of a Feynman path integral and evaluate it by the method of stationary phase in terms of a finite sum over all periodic orbits p of length t, with classical actions S p and amplitudes A p which are proportional to inverse square root of stability exponents Note that we chose to work at fixed t rather than at fixed energy E as is customary in semiclassical analysis of timeindependent systems. Here one assumes that the effective is small, i.e. S p for all p, which is justified for large Hilbert space dimensions N 1. A semiclassical representation of the spectral form factor can then be written as Berry identified the leading RMT contribution K(t) t from the diagonal terms of paired orbits, arguing that the non-diagonal terms of unequal orbit pairs average out in the leading order due to random phases. The diagonal contribution then results from Hannay-Ozorido de Almeida sum rule [34] p |A p | 2 = t, which is just a restatement of classical ergodicity. For systems with time-reversal invariance (TRI), the leading order of RMT result (7), K(t) 2t, then simply follows by pairing each orbit with itself p = p and its time-reversed partner p =p, noting that Sp = S p and Ap = A p . However, Berry's result only holds on time-scales much shorter than the Heisenberg time t t H which translates to spectral correlations on quasi-energy ranges much larger than the mean level spacing. That result can thus be considered as the leading order of a power series expansion (in t) of the RMT expression for K(t). Further progress came only 16 years later when Sieber and Richter [11,12] correctly identified the next-to-leading RMT term of H ) for TRI systems via the self-encountering periodic orbit doublets. Specifically, they decomposed the periodic orbit sum into two parts, the first containing a majority of orbits which never come close to themselves before the full period and in the second part, they considered orbits which experience a close self-encounter. They argued that the orbits from the second group form doublets with very similar actions S p and amplitudes A p which thus coherently interfere in the double sum (10) and result, after careful bookkeeping, exactly in the second order term −2t 2 /t H of RMT. It took another few years of efforts until this endeavour has finally been completed in [13,14] (see also [35] for the analysis of unitary-to-orthogonal ensemble crossover) by correctly identifying all the terms in the power-series expansion of K(t) from sums over chaotic periodic orbits with an arbitrary number of selfencounters.

II. PARTITION FUNCTION EXPANSION OF THE SPECTRAL FORM FACTOR
Here, however, we consider an interacting many-body system of quantum excitations without any meaningful classical limit, so the semiclassical periodic orbit theory is not applicable. We consider a system of spins 1/2 (qubits) described by Pauli spin operators σ where the time-evolution is given by the following two-step unitary Floquet propagator of a periodically pulse-driven Hamiltonian (time is measured in units of pulse period) where v is a 2×2 matrix with elements v 00 = v 11 = cos h, v 01 = v 10 = i sin h. In the basis of N = 2 joint eigenstates of σ x |s = (−1) sx |s , labelled by classical spin configurations s = (s 1 , . . . , s ), s x ∈ {0, 1}, W acts as a pure phase factor while the matrix elements of V factorize The propagator (12) defines a generic family of Ising models periodically kicked with a uniform transverse field (generalizing the kicked Ising chain [36,37], where RMT spectral fluctuations have been verified to a high precision [38]). See also Ref. [39] for a related discussion of transfer-matrix evaluation of the many-body propagator. In more abstract terms, one can also view H 0 as a generic integrable or many-body localized system with l-bits σ x and V as a global perturbation. The model is time-reversal invariant as the matrices of H 0,1 are real, i.e., V, W are symmetric.
We note immediately that the method that shall be developed below can be used as well to study a continuoustime version of the transverse field Ising model, where the kicked model (11) represents its trotterization (via the Trotter formula) by substituting J k x... → (∆t)J k x... , h → (∆t)h and carefully performing double scaling → ∞ and ∆t → 0 where the thermodynamic limit should be considered first. Further, more general forms of off-diagonal perturbations V can be considered by allowing an arbitrary spatial dependence of the magnetic field h → h x . Nonetheless, the particular system that we choose to study in the present paper represents a minimal generic model of many-body quantum chaos at ∼ 1.
We start by considering an expression (6) for the spectral form factor of Floquet systems K(t) = (tr U t )(tr U −t ) , defined for positive integer time t. Inserting multiple identities s τ |s τ s τ | = 1 in tr U t and s τ |s τ s τ | = 1 in tr U −t , we obtain Note that taking the trace implies periodic boundary conditions in time t + 1 ≡ 1. Assuming pseudo-randomness of the phases θ s one has  where s 1 , s 2 , . . . s t represents a lexicographically ordered string of words s 1 , s 2 , . . . , s t . In the ideal case, where all 2 phases θ s can be assumed to be independent random and uniform in [0, 2π) (which is equivalent to the assumption that all the coupling constants J k x,x ... are i.i.d.), the fluctuation term in (17) exactly vanishes. We shall refer to such an ideal situation in which fluctuations in (17) are set to zero as a random phase model (RPM). Below we show how to compute K(t) for the RPM and demonstrate that the result describes both the universal (RMT) and non-universal (short time) regimes of large families of clean kicked Ising models excellently.
For times much shorter than the Heisenberg time t H = 2 one may assume that all configurations s τ in the string s 1 , s 2 , . . . , s t are different. Then, eq. (17) implies that there exists a permutation π ∈ S t : τ → π(τ ), such that s τ = s π(τ ) , therefore (1 − δ sτ ,sτ+1 ) a half-number of domain-walls in a periodic spin sequence s (which is always an integer), Z π can be rewritten as Note that for the identity permutation, Z Id is a partition function of a classical one-dimensional Ising model (on a ring of circumference t) which can be calculated via a 2 × 2 transfer matrix T ss = |v ss | 2 , namely Z Id = tr T t = 1 + (cos 2h) t . Z π equals Z Id for any other permutation which does not change any neighbours in the string s, i.e. conserves the wall counting function w(s). These are exactly the t cyclic permutations and t anti-cyclic permutations -compositions of cyclic permutations with inversion τ → t + 1 − τ . For all other permutations π which contain at least one pair of neighbour-changes, one can show that the twisted partition functions Z π =Id are strictly smaller, and can be systematically computed using a diagrammatic technique (see Appendix A). Upon approaching the thermodynamic limit → ∞ at fixed t one thus finds an exact asymptotic result t * = − ln ln cos 2h .
This result can be interpreted as an analogue of Berry's diagonal approximation and yields the first order of RMT [Eq. (7)], while exact non-universal behaviour is predicted for times t t * = O(h −2 ln ). The time scale t * which separates the universal from non-universal behaviour can be interpreted as a kind of quantum manybody Ehrenfest time. Note that in the limit h → 0, the spectral form factor grows to the saturation value N at t ∼ 1, and one obtains the expected Poissonian behaviour which is typical for completely integrable systems [10]. At the two points, h = 0 and h = π 2 , where the crossover time t * [Eq. (22)] formally diverges (and close to them, for finite ) our expansions in t brakes down (see appendix A and Fig. 4). The case h = π/2 actually corresponds to a generic realization of a Floquet time crystal [29,40] which is non-ergodic and where the discrete translational invariance in time is spontaneously broken. This corresponds to a staggered behaviour of spectral form factor K h= π 2 (t) = N , t even, 0, t odd. (24) We note that in the range h ∈ π 4 , 3π 4 the RPM spectral form factor (21) in the non-RMT regime, t < t * , still displays characteristic period-2 oscillations.
Carefully subtracting double counted terms where exactly one configuration (word) in the string ('orbit') s 1 , s 2 , . . . , s t appears twice one obtains exactly the RMT result Eq.(7) up to the second order (see Appendix B). This can be considered as a manyqubit analogy of the Sieber-Richter pairing mechanism [11,12]. We conjecture that it should be possible to obtain RMT result to all orders by implementing the multiple-counting technique generalizing the diagrammatics sketched in Appendixes A,B.

III. KICKED TRANSVERSE FIELD ISING MODEL: THEORY EXPLAINS NUMERICS
We compare analytic results for the RPM to exact numerical computations of the spectral form factor in the following family of kicked Ising models x,x ... ≡ 0, (26) with normalization constants defined as and interaction effectively being short range N 1,2 = O( 0 ) for α > 1. Power-law decaying 1-spin terms and 2spin interactions are motivated by a requirement for the spectrum of H 0 to be non-degenerate and free from any other discrete symmetry, which should be a generic situation. Our results are not sensitive to the exact choice of α, as long as we are sufficiently far away from, either strictly local interactions α = ∞, or very long-ranged interactions α ≈ 0, where the model becomes mean-field like and describable by a single semi-classical degree of freedom.
In order to avoid the need of ensemble averaging we define a time-integrated spectral form factor as shown in Fig. 1, which is indeed a self averaging quantity as demonstrated in the inset. We observe very good agreement with the RPM. As K int (t) propagates deviations at short times to longer times, we also show the non-self averaging K(t) at short times in Fig. 2 and again observe very good agreement with RPM upon averaging over an ensemble of values of parameter J or taking a moving time-average over a short window of time for fixed J (and a, b, α). Even the fluctuations of the two averages around the theoretical prediction (RPM) for similar statistical sample size (n = 10) look quantitatively comparable. The data reported in Figs. 1,2 were obtained for locality exponent α = 3/2. In Fig. 3 we investigate the role of α in more detail. When the model becomes increasingly short ranged, i.e. increasing α, the fluctuations part in Eq. 17 can no longer be neglected and the model develops deviations from RPM. The crossover time t * at which K(t) begins to follow RMT becomes larger and starts to depend on the parameter J. But even for strictly local, nearest-neighbor interactions (α → ∞) at fixed J this time seems to scale polynomially, perhaps like ∝ 2 , and is still much smaller than the Heisenberg time t H = 2 .
In order to quantify a transition between RPM and non-RPM physics we define the following order parame-ter: Since the RPM prediction for long times becomes equivalent to RMT and is expected to match K(t) well, provided the model is non-integrable and ergodic, the order parameter ψ becomes independent of t max as long as t * t max t H . As shown in Fig. 3, we indeed find a phase-transition-like behavior of ψ(α) with two, short-range and long-range critical points, α * long ∼ 0.5 and α * short ∼ 2.5 : For α * long < α < α * short , agreement with RPM is excellent and ψ(α) is small and of the order of expected statistical fluctuation due to averaging over J, while outside the range ψ(α) quickly grows.

IV. CONCLUSION
Our work discloses the first theoretical mechanism which connects RMT to simple many-qubit systems in (effectively) low dimensions. There are many immediate further questions and generalizations which are to be studied: i) The assumption of the pure δ-correlator of phases (17) is on the same level of rigour as the randomphase approximation in standard semiclassics, but one may hope to find a more rigorous justification here. ii) The interesting case of local Ising interactions (α = ∞) also obeys RMT physics [38] but needs to be studied separately as r.h.s. of eq. (17) then acquires extra systematic contributions. iii) One may generalize our technique to study universal behaviour of dynamical correlation functions (i.e. spin structure factors, etc) in the quantum chaotic regime, iv) Furthermore, one may expand our methods by introducing quenched disorder, say in the transverse field, and attempt to approach the manybody localization transition [25] from the ergodic side, for instance, by tuning the locality exponent α.
Our results have a direct relevance for understanding the vast body of numerical experiments, simulations, and in the near future possibly also experimental spectra of highly excited simple many-body systems, which correspond to ever longer accessible observation times of perfectly coherent out-of-equilibrium quantum systems (see e.g. [30]). The ideas of many-body quantum chaos and random matrix theory are also vividly debated in the context of high energy physics and holography [41][42][43], where our results and methods could also be applied.
After our work has been submitted for publication, we have learned of a series of related works [44,45], where the RMT spectral form factor has been computed for local Haar-random unitary nearest-neighbor quantum circuit propagator, which corresponds to UE universality class of RMT, in the limit of large local Hilbert space dimension q. It is remarkable that in Ref. [45], where the authors consider a related variant of RPM, but insisting on the strict locality of the (nearest-neighbour) interaction at the expense of having to consider a large q limit, a b  (18), we have to sum over all permutations π ∈ S t . As we already noted, the (anti-)cyclic permutations together with the identity permutation, which form a subgroup of S t , yield identical leading contributions which become expo-nentially (in ) dominant in the thermodynamic limit.
Here we identify and explicitly calculate the contributions of the permutations which yield the leading (first order) corrections. From eq. (20) it follows that each contribution to Z π depends only on the number of domain walls in periodic strings s = (s 1 , s 2 , . . . , s t ) and π(s) = (s π1 , s π2 , . . . , s πt ). Therefore Z π depends solely on a diagram obtained by plotting a directed graph of sequentially arranged nodes (s 1 , s 2 , . . . , s t ) with the links s π1 → s π2 , s π2 → s π3 , . . . , s πt → s π1 . For example, the (anti-)cyclic permutations are then represented as circular loops, meaning that they preserve sequential order. The next to leading order Z X comes from the so called X-diagrams (shown in the diagrammatic expression below and illustrated numerically in fig. 4), where all connections apart from two (order changes) are kept intact (sequential). The diagrammatic expression for the case where the first sequential stretch has length τ , and the second length t − τ − 2, reads: Circular black arcs represent summations over stretches of sequential spins [τ , τ ], namely over all s τ , τ < τ < τ . These are given by the powers of the transfer matrix T τ −τ +1 s τ ,s τ , specifically T p ss = 1 2 1 + (−1) s−s λ p where λ = cos 2h plays the role of the coupling constant. Red dots correspond to remaining spins which one still needs to sum over while putting matrix element v ss for each broken sequential link s → s (dotted) and v * ss for each crossed link s → s .
Of similar importance are the XX-diagrams with two sequential crossings:  which are straightforwardly evaluated (31) One can show that contributions of all other diagrams, starting with two crossings separated by sequential stretches, triple crossings, etc., are of the form Z other (τ ) = 2 −n (1 + O(λ k )) with n ≥ 2 and k ≥ 1, so they contribute to Z only beyond the second order in t/2 and will be ignored here.
Each X and XX diagram (for fixed τ ) has multiplicity t 2 (or t 2 /2 for τ = t/2), since any Z π is invariant under 2t cyclic and anti-cyclic permutations and we can start drawing the diagram at t/2 inequivalent points.
We now see that the smallest gap of Z Id − Z π > 0 comes from Z X (1) for λ > 0 and Z X (2) for λ < 0 (see fig. 4). Since these terms enter to power , the sub-leading contributions to K(t) are exponentially suppressed by a factor of the order of ∼ t((1 + λ)/2) for λ > 0 and t((1 + λ 2 )/2) for λ < 0. When we approach the Heisenberg time t ∼ t H = 2 these contributions become important, as we will see in the next section.
Appendix B: Second order term in t/t H .-RMT predicts that the next term in the expansion is −2t 2 /2 , eq. (7). We now show that the contributions of the X and XX diagrams and the possible repetitions of the spin configurations almost cancel, yielding exactly the RMT result.
On the right, the corresponding diagrams are shown for clarity, where the blue diamond sites connected with a dashed line depict the repeating spin. The first two sums come from enumerating all X and XX diagrams (as explained in Appendix C). The sum from the repeating spin configurations comes next, and is written in two parts. The combinatorial factor of the diagrams is 4t when τ = 1, 2, t − 1, t − 2 and 8t otherwise, which we take into account by writing two sums with different starting and final value of τ . In the last line of eq. (34) we note a remarkable cancellation of all terms apart from the RMT result for times t > t = O( log λ) where t H /t is still exponentially large in . This could be viewed as a quantum many-body analogy of the Sieber-Richter selfencountering-orbits mechanism.