From operator statistics to wormholes

For a generic quantum many-body system, the quantum ergodic regime is defined as the limit in which the spectrum of the system resembles that of a random matrix theory (RMT) in the corresponding symmetry class. In this paper we analyse the time dependence of correlation functions of operators. We study them in the ergodic limit as well as their approach to the ergodic limit which is controlled by non-universal massive modes. An effective field theory (EFT) corresponding to the causal symmetry and its breaking describes the ergodic phase. We demonstrate that the resulting Goldstone-mode theory has a topological expansion, analogous to the one described in arXiv:2008.02271 with added operator sources, whose leading non-trivial topologies give rise to the universal ramp seen in correlation functions. The ergodic behaviour of operators in our EFT is seen to result from a combination of RMT-like spectral statistics and Haar averaging over wave-functions. Furthermore we analytically capture the plateau behaviour by taking into account the contribution of a second saddle point. Our main interest are quantum many-body systems with holographic duals and we explicitly establish the validity of the EFT description in the SYK-class of models, starting from their microscopic description. By studying the tower of massive modes above the Goldstone sector we get a detailed understanding of how the ergodic EFT phase is approached and derive the relevant Thouless time scales. We point out that the topological expansion can be reinterpreted in terms of contributions of bulk wormholes and baby-universes.


Introduction
A resolution of Hawking's information paradox, [1], is a cornerstone to our understanding of quantum gravity. While the statement of the paradox in terms of the seemingly featureless radiation that is emitted by a black hole is the most popularly known one, the underlying problem of loss of unitarity manifests itself in many other forms. For a gravitational theory on anti-de Sitter (AdS) backgrounds, where conventionally the radiation does not escape to infinity, a natural set of observables are the boundary correlation functions. In the semiclassical approximation of the bulk theory, such correlation functions decay indefinitely at late times, [2]. This is paradoxical when viewed in light of the AdS/CFT correspondence, according to which a gravitational theory is dual to a quantum theory living on its boundary in one lower dimension, since it can be argued [2,3] that the correlation functions cannot decay indefinitely in a quantum mechanical theory with a discrete spectrum. Recent works have discovered new geometrical configurations of the gravitational path integral, [5][6][7][8][9][10], which relieve the tension with unitarity in observables such as correlation function or entanglement entropies. These new configurations are associated with Euclidean geometries connecting several boundaries, also known as wormholes. Given that on each of these boundaries there resides a copy of the dual quantum mechanical theory, these gravitational configurations represent connected correlations in the spectrum of the dual theory. 1 Interestingly, these configurations also contribute to observables corresponding to single boundary geometries by introducing handles in these geometries. The type of correlations induced by these wormhole more typically arise in situations where one considers not a single theory and its observables, but instead an ensemble of theories and an associated averaging procedure, a canonical example being random matrix theory (RMT) [11][12][13][14][15]. This realisation has led to a shift in our understanding of holography in which one adopts the view that a semi-classical theory of gravity is dual to an ensemble of quantum theories, [16]. On the other hand, the authors of [4] have proposed an alternative paradigm in which the ensemble description in the quantum mechanical side of the holographic duality is an emergent phenomenon that arises in physical unitary theories, 2 more specifically as a consequence of quantum chaotic dynamics. This proposal is closely related to the quantum ergodic limit and late-time chaotic dynamics of quantum systems, as we now explain. 3

Quantum chaos and wormholes
A key signature of quantum chaos is based on properties of the spectrum of the Hamiltonian of a quantum system. Berry and Tabor [21] and Bohigas-Giannoni-Schmit [22] conjectured a remarkably universal property of quantum chaotic systems: their spectral statistics coincides with that of an appropriate RMT, while that of integrable quantum systems obeys Poisson statistics. A distinguishing feature between these behaviours is the level repulsion between the neighbouring eigenstates. For systems that lie in the symmetry class of Gaussian unitary ensembles (GUE) translating this characteristic behaviour of level spacings into the time domain, in other words following the time evolution of correlation functions and similar probes, generates another universal signature: a linearly rising ramp, culminating in a late-time plateau, as further illustrated in Figure 1.
The onset of the linear growth, also referred to as ramp, is known as the ergodic/Thouless time, t er /t Th . This time scale marks the onset of the universal behaviour that can be attributed to the RMT-like spectral properties in a given chaotic quantum system. In a single realisation, the ramp-plateau are superimposed with fluctuations and the universal behaviour emerges as an average over these fluctuations. At the level of the spectrum of the theory, this universal behaviour manifests itself in the two-point correlation function of the density of states (DoS), The occurrence of this sine-kernel, sinc(ω) gives rise to the universal ramp-plateau behaviour. Here, ∆ is the mean level-spacing and is defined precisely below.
Operator expectation values in pure states or thermodynamic ensembles form a more physical set of observables in a theory as they can be measured directly in experiments. For example, twopoint functions of operators are the simplest observables measuring the linear response in a system. A perpetually decaying correlation function is characteristic of a coupling to an external system and hence non-unitarity. Therefore it is important to understand how the transition in the behaviour of a correlation function from an initial decay to a non-vanishing plateau at long times happens. Additionally, operators are more appropriate observables to understand how a system in a pure state might appear to be thermal,à la Eigenstate thermalisation hypothesis (ETH), [23][24][25][26][27][28]. The observable, the operator resolvent, that our work focuses on in this paper is related to many such 1 To see this the reader may associate each boundary with one copy of the partition function Z(β) of the boundary theory. These connected correlations between several partition functions are tantamount to connected correlations in the spectrum of the theory. 2 We often use the expression physical unitary theory to describe one in which an ensemble average over theories is not performed. 3 For some related discussion on emergence of ensemble averaging in theories also see [17][18][19][20]. . From the holographic perspective, the conventional gravitational contribution computes the decay in the non-ergodic regime. Region I corresponds to the fast decay determined by the smooth energy dependence of coarse-grained DoS. This corresponds to the contribution of the disk geometry to the correlation function. In region II, correlation function decay exponentially because of presence of non-ergodic and non-universal massive modes. The part of the curve in purple is the universal ramp (region III) and plateau (region IV) behaviour and is governed by the mean level statistics of the eigenstates. The timescales at which the linear ramp becomes conspicuous is known as t er or the ergodic time. Lastly, t H is the time at which the plateau begins and is of the order of Heisenberg time. The universal ramp-plateau corresponds to the contribution of the sum over disk geometries with handles around the standard and the Andreev-Altshuler saddle points. Also see Figure 5 for more accurate representation of operator correlation functions in the SYK model.
physical observables quite simply, as explained in subsection 2.1. One of the main results of this paper is to establish a similar sine-kernel behaviour for many-body operators, O, governed by the operator sine kernel It is interesting that the '1' in (1.1) corresponds to the disconnected part of the correlation function whereas in (1.2) it is part of the connected part of the correlation function. As we explain in detail in the bulk of this work, the traces on the right hand side are over the entire Hilbert space or parts of it. Fourier transforming into the time domain gives rise to the universal ramp-plateau behaviour of operator correlations illustrated in Figure 1. Another observation is that this result is consistent with the form of the ETH ansatz of [23,24], O αβ =Ō(E) δ αβ + e −S(E)/2 f (E, ω) R αβ , whereŌ is the expectation value of the operator in a microcanonical ensemble at the energy, E; S(E) is the corresponding entropy in the microcanonical ensemble, and R αβ is a random number of zero mean and unit variance. Our result here encapsulates universal contributions to ETH-like correlations in the ergodic phase of any chaotic quantum system, as well as ergodic and non-ergodic contributions for the Sachdev-Ye-Kitaev (SYK) model in particular. The emergence of the ergodic phase in a wide class of quantum many-body systems was recently studied in [4,29,30]. The relevant degrees of freedom in the ergodic phase correspond to certain pseudo-Goldstone modes that arise from the causal symmetry breaking in the study of correlation functions of the kind (1.1). This effective field theory (EFT) description applies in the late time (ergodic) limit, while at early times (nonergodic limit), observables decay at a rate determined by the smooth energy dependence of the coarse-grained DoS, ρ(E), of the system, and are insensitive to the finer details of the spectrum. In between these two regimes lies a non-universal intermediate regime which we analyse below for cases of interest in holographic duality.
In the bulk the initial coarse-grained decay corresponds to the prediction of semiclassical gravity. What then corresponds to the late-time universal ergodic behaviour? The naïve semiclassical gravity computation reproduces the initial decay of physical observables, but it does not reproduce the ramp-plateau behaviour; the observables instead decay to zero. Remarkably, including new wormhole configurations in the semi-classical gravity path integral reproduces the required linear ramp at late times. Our EFT description of the ergodic phase of a holographic many-body theory allows us to predict these contributions and to classify them systematically. We establish that the behaviour of the system in this regime is equivalent to that of a random-matrix theory where the matrix is interpreted as the Hamiltonian of the system. One may then interpret each diagram that arises in the EFT description in a topological expansion as constructing a bulk surface of given topology, along the lines of the old matrix models [31,32]. For the case of SYK and similar lowdimensional examples one can literally interpret this construction as the triangulation of the bulk geometry, but for higher-dimensional examples the bulk interpretation remains somewhat more mysterious.
At even later times, the plateau arises due to doubly non-perturbative corrections. From our discussion above, this means that such configurations in the semi-classical gravity contain information about the fine-grained spectrum of the quantum theory, albeit in an averaged sense. That is to say, they don't reproduce the noisy fluctuations that are observed in a physical unitary theory, but do reproduce their averaged behaviour.
We would like to highlight a fact that might be of particular holographic interest. Our work establishes that the wormhole geometries receive contributions not only from the spectral details of the theory as described in (1.1), but also from the Haar integration of the unitaries in an ensemble description. The causal-symmetry approach advocated here thus provides a unified description of the different field-theoretic descriptions of the gravitational wormhole configurations that have been put forward in [4,6,7] and [33]. Another important feature that has emerged and that we explicitly corroborate in this work is the fact that the gravitational configurations that contribute to the observables in the ergodic phase (at late times) depend upon the observables under study, in contrast to what was widely believed originally.

Going beyond the ergodic regime
In the pre-ergodic phase, the above-mentioned transition from an initial decay to the universal ergodic regime happens at intermediate timescales governed by the 'massive modes' of the theory. As discussed in detail in [4] and in the present work in the specific context of the SYK model, these massive modes correspond to Hilbert space non-homogeneous projections of the collective field. The ramp starts when the lightest of these modes decay at t er . Reflecting the non-universal nature of this regime, the ergodic time itself may vary for different observables, as is apparent by comparing the results of this work with those of [29].

Ergodic regime in physical theories
Having argued for the interest and importance of spectral chaos in a quantum system as well as in a dual theory of gravity, we want to emphasise that demonstrating these behaviours in a given field theory is a non-trivial task. For example, quantum ergodicity has not been explicitly demonstrated for N = 4 SYM, which is a benchmark of holographic quantum field theories. In fact, the study of the spectral statistics of quantum systems has largely been a numerical endeavour. The SYK model, [34][35][36], is a toy-model of N Majorana fermions that facilitates explicit analytic computations. Hence, it provides an excellent laboratory to test and improve our understanding of many of these issues related to quantum chaos. Since it is dual to a two-dimensional theory of gravity, [37][38][39][40], it also sheds light on various aspects of AdS/CFT correspondence in general. In fact, it was the original observation by Kitaev, [35], namely that the OTOCs in this model grow exponentially with a Lyapunov index that saturates the bound proposed in [41] which has led to intensive study of this model in recent times. It might appear that the SYK model is already the same as a RMT since it is a theory with a Hamiltonian chosen from an ensemble, i.e. a theory of disorder averaging. However, for the SYK model with a (q/2)-body interaction and a Hilbert space of dimension, D = 2 N/2 , only N q ∼ N q matrix of elements out of a total of D 2 in the Hamiltonian are chosen randomly. This exponentially small amount of randomness indeed leads to the deviations from the RMT behaviour seen at early times. Quantum ergodicity in the SYK model was studied numerically in [42,43]; and subsequently, analytically in [29,44,45]. In these works, the late time behaviour of the spectral form factor is studied using two different techniques, namely the replica trick and supersymmetry, and the emergence of RMT-like statistics is derived analytically. The study of the spectral form factor as a fine-grained probe of the spectrum in mesoscopic physics is quite old, [46,47]; however, it was popularly introduced in the string theory community originally in [48] and more recently in [42]. In the present work, we generalise their results for the case of operator two-point functions. Before we proceed further, let us provide a brief structure of this paper along with a preview of our main results.
Main results and the structure of the paper We discuss the ergodic limit of the SYK model in section 2. In the process, we study the socalled operator resolvent which measures the behaviour of an operator two-point function in energy eigenstates as a function of the difference of their energies. This is closely related to finite time thermal two-point functions as discussed in subsection 2.1. In the rest of the section, we provide a path integral representation of this observable and derive the precise ergodic limit in the case of the SYK model. This is done by studying the model in the energy (or equivalently, time) regime when the universal σ-model physics becomes important. At late times, the average behaviour of the operator correlation functions has a ramp and the plateau behaviour that originates from the contribution of the Goldstone modes in the sigma model. Unlike the observables such as the spectral resolvent, in this case the relevant contribution comes not only from the spectral statistics but also from averaging over the SYK eigenstates, which in the ergodic limit becomes equivalent to the Haar averaging. Importantly, our result also clarifies the relationship between ETH and RMT predictions for the operator correlation functions. Not only do these observables depend on the behaviour of operators in energy eigenstates, they also contain the information about the distribution of the eigenstates themselves as a function of the energy.
In section 3, we study the slowest non-universal non-ergodic modes that mark the late time transition of observables from non-universal to universal ergodic behaviour, the timescales corresponding to the onset of the RMT behaviour. An important observation is that t er depends on the observables under study.
We finally conclude with some discussion in section 4. For a more in-depth discussion of the supersymmetric methods we point the reader to [4,[49][50][51]. In Appendix A, we discuss the dominance of the different saddle point solutions of the classical equation of motion. Various details of the derivation of the σ-model action and the computation of operator correlators in the ergodic limit that have been omitted in the main part of the paper are presented in Appendix B. In Appendix C we discuss the details of the non-ergodic modes that govern the physics at Thouless time.

Ergodic limit and σ-model
We begin this section with an introduction of the observables that we are interested in. The operator resolvents measure the probability of all quantum mechanical processes that exchange a given amount of energy in the presence of the operator insertion. In other words, they quantify the behaviour of the operator matrix elements in energy eigenbasis (state |α is an eigenstate with energy, E α ) corresponding to the states with a fixed energy difference across the spectrum. Below, we describe how the operator resolvents that we study in this work are related to the time dependent correlation functions in (micro-)canonical ensembles. In a general setup, we write a generating function as a path integral over some auxiliary variables to compute these observables.

Operator resolvent
The Fourier transform of the resolvent is related to the thermal two-point function of the operator O at infinite temperature, Note that the energy argument, ω, is dual to the real time parameter, t. Therefore, early/late times correspond to large/small values of ω, and are referred to as UV/IR, respectively. This is different from the other notion of UV/IR which refers to the low/high energy part of the spectrum of a theory. From the spectral resolvent one can obtain other more familiar quantities, such as The correlation function in the Canonical Ensemble: this is the Fourier transform of the spectral resolvent of an operator of the form O → e − β 2 H O at temperature, T = 1/β.
The correlation function in the Microcanonical Ensemble: this quantity can be obtained by an insertion of a projection operator, P W , that restricts the sum over the eigenstates to an appropriate microcanonical window, E α,β ∈ Ē − ∆E,Ē + ∆E =: W. This can be achieved by Note that the new insertions aren't Hermitian conjugates of each other. A possible choice of P W = e − β 2 H gives the regulated thermal two-point function, This regulated correlation function is related to the canonical correlation function by an analytic continuation, t → t + i β 2 , and is relevant for the observables in the holographic setting. Note that any operator can be written as a sum over a traceless and a trace component, The analysis of the term proportional to identity needs to be performed more carefully. Subsequently, for pedagogical reasons we show the details of only the traceless term. However, the full answer including the traceful and the traceless parts is presented in (2.38). The operator resolvent can be written in an integral representation that is more useful for our purpose, are the resolvents of the many-body Hamiltonian. 4 Henceforth, we concentrate on these two terms separately. The second term above doesn't contain the ergodic modes but needs to be kept in order to get the correct factors to obtain level-repulsion and is studied in Appendix B.4. The non-trivial physics of the σ-model arises from the first term. It can be rewritten as a path integral, is the generating function for the correlation functions. One can use either the replica trick [29] or supersymmetric techniques [44] to compute the generating function. A recent work by some of the authors, [4], discusses these techniques quite generally. We refer to Appendix B for the details of the computations presented in this section. Z[h] can be written as a graded Gaussian integral, Here,Ψ, Ψ are 4D dimensional auxiliary graded vectors with 2D Grassmann components and 2D c-number components acting in the tensor product R/A ⊗ b/f ⊗ V of three linear spaces corresponding to the retarded/advance, boson/fermion and the physical many-body Fock space of the SYK model, respectively. The retarded/advanced product space corresponds to the choice of the sign of the regulator in the kernel of the above integral, (2.6). This in turn ensures that the Gaussian integration over corresponding components of the graded vector, Ψ, gives a differently regulated Green's function, G ± (E) = (E ± i0 + − H) −1 . This is needed for generating both the factors of G ± in (2.4) using the integration of Ψ fields. Also,Ψ = Ψ † · σ RA 3 , which guarantees a convergence of the Gaussian path integral (2.6) in the advanced sector. The variable z is defined as, where the Pauli matrix σ RA 3 acts in the RA-space as indicated by the upper index. The kernel for the Gaussian integral contains the Hamiltonian of the theory under study, H, and appropriate sources, Here P b is the projector to the bosonic sector, while a block structure refers to the RA-space. It can be readily checked that (2.5) holds. The determinant factors arising from the Gaussian integral The dots represent the projectors on the bosonic sector as well as the insertion of the sources for the operators. It arises from computing the resolvent with + +-causality in (2.4) and therefore does not receive contributions from the ergodic modes, B,B. This diagram can be understood to arise from the integration over the eigenfunctions of the Hamiltonian. In this diagram we have refrained from explicitly demonstrating the disorder average. The figure (b) is a depiction of the diagram as a bulk geometry, which is a disk in this case. The bulk is understood to emerge out of disorder averaging.
cancel between the Grassmann and the c-number sector once the sources have been switched-off, h ± = 0, resulting in a normalised correlation function.
As described here, the above technique is quite general, although, it is a little formal and abstract. It provides little insight into how to obtain the spectral properties of a physical unitary system. An integral over the elements of the Hamiltonian in (2.6) provides an analytical handle to compute an average over such physical observables described above, [49][50][51][52]. The RMT was originally developed to model the strongly interacting Hamiltonians that arise in nuclear and atomic systems. The important question is to understand how it arises as an effective description in a physical unitary theory. The high frequency quantum noise that exists in physical unitary systems over and above the smooth averaged behaviour easily dephases under any form of averaging, not just that of the Hamiltonian. An averaging over even small energy scales leads to dephasing. This is reminiscent of the RG interpretation where coarse-graining over small timescales (equivalently, lengthscales) scales or large energy scales leads to an effective theory. In [4] it was argued that the emergence of this effective description is captured by a low-energy σ-model. Therefore, instead of searching for the randomness of the Hamiltonian in physical systems, one should try to understand in what physical regime the σ-model describes our strongly-interacting systems like N = 4 SYM theory. Demonstrating it for chaotic holographic field theories like SYM is still a daunting task, however there exist holographic many-body systems such as the SYK model which do provide a viable proving ground to demonstrate these ideas. As explained in the introduction, the SYK model, despite being defined as an ensemble of random Hamiltonians, is not a RMT. Its Hamiltonian has much smaller randomly chosen independent elements than the full RMT. We now specialise our calculations to the SYK model as a demonstration of our ideas. Some general lessons that we learn from the sigma model approach for the operator correlations functions in the ergodic limit of any theory can be easily summarised as follows. The leading contribution is due to planar diagrams with a single boundary as depicted in Figure 2. 5 We have drawn this diagram in the flavour space. Alternatively it can be drawn in the colour space where the blue lines that demonstrate the disorder average themselves are drawn in a double line notation corresponding to the colour-indices carried by the Hamiltonian. In this dual description the colour lines can be directly interpreted as triangulating a two-dimensional surface with a disk topology (see [4] for a more detailed discussion on the colour-flavour duality). This is reminiscent of the seminal work on the duality between a particular double-scaled matrix model and 2-dimensional gravity, [31,32]. 6 The difference between the present case and the models studied in the past is that the matrix models in the present case are only an effective description captured by the Gaussian matrix models in the late-time limit. Recent work of [57] might provide a more concrete connection in the case of the SYK model but we also believe that this observation holds for an arbitrary quantum theory which might not have an explicit matrix model description.
The subleading contribution to the operator correlation function is given by a torus with a boundary, depicted in Figure 3. This diagram contributes the sine kernel, (1.2), which gives rise to the ramp-plateau behaviour. It is important to note that in terms of the geometries the Euclidean wormholes directly contribute to the variance of the physical observables as discussed in [5,6]. More generally, the number of boundaries of the bulk geometry are related to the number of traces in the physical observables. However, wormholes still contribute as intermediate states (baby universes) in other observables. Within the sigma model this is manifested by the properties of the Wick contraction, (B.13). The first of these equations can be seen as a contribution of the wormhole, while the second one can be seen as an emission of baby universes.

SYK model
The SYK model, [34][35][36], is a (0+1)-dimensional theory of N Majorana fermions with the Hamiltonian, The coupling constants, J i1i2i3i4 , are chosen randomly from a Gaussian ensemble, Note that the SYK Hamiltonian is quite sparse, there are only n = N 4 independent matrix elements in it out of the total D 2 = 2 N −2 elements of the full Hamiltonian. This, in principle, is the reason behind the deviations from random matrix theory. For brevity, we denote the product of four Majorana fermions by where, a ≡ {i 1 , i 2 , i 3 , i 4 } is a condensed index labelling the q-fermions in the individual interaction terms of the Hamiltonian and satisfy the condition is the basis of matrices that span the space of sparse SYK Hamiltonians. We also denote a basis of arbitrary D × D matrices by a product of arbitrary number of Majorana fermions, The index µ k ≡ {i 1 , i 2 , . . . , i k } not only labels the individual fermions that enter the product but also the number of fermions (k) in the product itself. It therefore satisfies the constraints 1 ≤ i 1 < i 2 < . . . < i k ≤ N for k = 0, 1, . . . , N . The element X 0 = 1 corresponding to k = 0. In case of theories with an ensemble average, (2.6) is modified to include the disorder average, where, · dis denotes the disorder average. On substituting the SYK Hamiltonian in (2.12) and performing the disorder average, one generates a quartic interaction term in the graded super-vectorΨ, Ψ. This quartic term can be reduced to a quadratic term by introducing an auxiliary field, A a , using Hubbard-Stratonovic transformation. The remaining path integral is 7 We have introduced the energy scale 14) in the large N limit. This is a quadratic path integral in the vectors,Ψ, Ψ, which can be integrated out to obtain, We decompose the Hubbard-Stratanovic fields into average and difference fields, Since onlyĀ enters the logarithmic part of the action in (2.15), the path integral over difference fields A a is Gaussian and can be completed exactly. 8 The fieldĀ can be written using the basis (2.11),Ā = µā µ X µ . (2.17) Here,ā µ carries the 4 dimensional indices in the supersymmetry and adv./ret. space, while X µ carries the Hilbert space indices. 9 In terms of these fields, the generating function is where, K denotes the kernel of the quadratic integral over the superfields, Ψ andΨ, Note that the STr represents super-trace over the 4-dimensional space in the first term, but a super-trace over 4D-dimensional space in the second term. In future, we refrain from explicitly stating this fact as it will always be clear from the context. The permutation symbols s(µ), S(µ), count the sign when permuting X µ across themselves or X a . Specifically, s(µ) = ±1 is defined by (X µ ) † = s(µ)X µ and depends only on the number of Majoranas |µ| contained in µ. The sum S(µ) reads where sign factors s(µ, a) = ±1 stem from commutation relations X µ X a = s(µ, a)X a X µ . At this stage, we make an ansatz that will be justified retrospectively in section 3. We assumē i.e. the dominant saddle point solution(s) are homogeneous in the Hilbert space. We refer to this ansatz as homogeneous Hilbert space ansatz (HHA). With this ansatz, the path integral in (2.18) is reduced to an integral only over the homogenous modes, where 1 Fock stands for the identity operator in the physical Fock space. 8 The Jacobian of the integration over the difference fields is trivial because they are related to the original fields only through linear transformations. Additionally because the constraint on the sum of these fields in (2.16) is also linear, it doesn't give rise to any Jacobian. The one-loop determinant is trivial because these fields are graded fields and the integration is quadratic (this follows from Parisi-Sourlas-Efetov-Wegner theorem, [49,50]). For the details of these integrals we refer the reader to Appendix A of [29]. 9 We have simplified the notation so that µ collectively refers to the number of Majorana fermions in the basis, k, as well as all possible distinct choices of them for a fixed number.

Causal symmetry breaking
The path integral in (2.21) can be studied in the classical saddle point approximation when D → ∞. The dominant solution of the saddle-point equations is A more detailed study of choosing the correct saddle point out of the 16 naïve possibilities can be found in any standard text on the subject like [49,50]. The procedure is also succinctly reviewed in Appendix A. Note that y → T · y · T −1 is the symmetry of the first term in (2.21) for all T ∈ U (1, 1|2). However, it is a symmetry of the full action only if T · z · T −1 = z. For ω = 0, ε → 0, the full T ∈ U (1, 1|2) remains the symmetry of the action which is weakly (and explicitly) broken by ω = 0. The saddle point solution, (2.22), breaks this symmetry spontaneously to U (1|1) × U (1|1). The coset space U (1, 1|2)/U (1|1) × U (1|1) is that of pseudo-Goldstones which have a finite action cost associated with them when ω = 0. However, this action cost is suppressed with respect to the saddle point action by a factor of ω/γ ∼ ω/(J √ N ).Thus the perturbative treatment is valid when where ∆ ∼ γ/D is the average many-body level spacing in the band center. This provides an explicit realisation of the universal causal symmetry breaking that was shown in [4] to give rise to the ergodic behaviour more generally. The sigma-model manifold generated by the T -transformations contains an alternate saddle point of the equations of motion corresponding to, [47], which can be arrived at from the dominant saddle point by a special transformation, where P b and P f denote the projectors to the boson/fermion sectors, respectively. We refer to this saddle point as the Andreev-Altshuler saddle point subsequently. The contribution of this saddle point is non-perturbative with respect to the contribution of the dominant saddle, (2.22). Ideally, one should integrate over the entire sigma-model manifold because these modes are 'soft'. However, because of the one-loop exactness of the theory, it is sufficient to sum the perturbative contributions around the individual saddle points, [49,50].

Integral on the σ-model manifold
Following the above discussion, when the integration over the y-variable is reduced to the σ-model manifold, the generating function is given by the following path integral, Q is a point on the coset space U (1, 1|2)/U (1|1)×U (1|1). We refer the reader to Appendix B for the details of the computation on how to derive this equation from (2.21). Here, ρ(E) is the mean-field density of states given by, This is a good approximation of the density of states only near the centre of the band of the SYK spectrum. For modifications in the present analysis to give the correct DoS away from the band-centre see [58]. The action S src (h) contains all terms up to quadratic order in the sources, which are obtained by expanding the log-term in (2.21) w.r.t. h ± ., see details in Appendix B. On differentiating the generating function (2.26) twice, one gets the building block of the correlation function (2.5), (2.28) The path integral above features the semi-classical exactness and thus can be evaluated by performing the Gaussian integration around the two saddle points, Λ and Λ AA . In the first case the T matrix within the perturbation theory can be parametrised as where we introduced the matrices Here x, y ∈ C, while µ,μ, ν andν are independent Grassmann variables. In analogy with low-energy QCD, one may think of W as the pions as described for example in [4]. For the contribution around the Andreev-Altshuler saddle point Λ AA one parametrise T = T 0 exp(−W ), where W is defined as before.

Perturbations around the dominant saddle point
It can be verified that it is sufficient to expand the action in (2.28) to quadratic order in W and the sources in the pre-exponent to quartic order. The action then becomes It can be regarded as the "quark mass" for the pion field W . An astute reader might point out that it is not consistent to expand to different orders. However, it can be checked that the contribution arising from the quartic interaction term in the exponent evaluates to zero. This is related to the non-renormalisation theorem for the GUE symmetry class, [49], the only one we treat explicitly in this paper. The expression for sources expanded to higher order in W can be found in Appendix B.
Performing the integrals on the pion fields and taking the real part to compute the resolvent, we obtain We have used the dimensionless frequency, This concludes our evaluation of the contribution from the standard saddle point to the operator resolvent. To this we need to add the contribution from the Andreev-Altshuler saddle point, which we shall do now.

Perturbations around the Andreev-Altshuler saddle point
The Gaussian action around the Andreev-Altshuler saddle point (2.24) evaluates to Note that the subleading saddle point Λ AA corresponds to the Weyl symmetry, which is equivalent to exchanging the two eigenvalues in the fermionic sector of the standard saddle point Λ.We refer the reader to [4] for a detailed discussion on this symmetry in terms of exchanging the north and the south poles of the S 2 factor of the coset manifold. The evaluation of the resolvent gives Adding together the contributions of both the saddle points, the total resolvent is To get the final answer we need to subtract from the above the contribution of R + + computed in (B.25), However, recall that so far we have been working with traceless operators. Generalising to the case of traceful operators requires a more careful analysis and we present only the final result here, This is the operator resolvent in the ergodic limit, (1.2), advertised in the introduction. Figure 4 compares the numerical computation of the operator resolvent in the SYK model with this analytical porediction. Before we consider the contributions of the non-ergodic modes in the next section, let us look at the behaviour of the correlators as a function of time.

Late time behaviour
As discussed in the introduction, the important feature of the observables in energy space that leads to the universal ramp-plateau behaviour at late times is the appearance of the sine-kernel. However, to see this behaviour explicitly one needs to perform additional integrals of the resolvent, (2.38), over the energies E, ω, The time at which the ramp transitions into the plateau depends on the value of energy, E, of the microcanonical window because ∆(E) is a function of E. Therefore, the above integral is a convolution of various ramps and plateaus defined for each microcanonical energy window. The last ramp ends at the Heisenberg time t H = 2D/γ corresponding to the band center. When working  with the infinite-temperature canonical ensemble, the integration with the DoS (2.27) leads to the answer where we have used scaled time,t = t/t H . While this is not the same linear ramp-plateau that one obtains in a microcanonical ensemble and is more familiar from the literature, the long-time behaviour of the observables still has a monotonic growth followed by a plateau, Figure 5. The initial monotonic growth comes from the contributions of the correlations between the energy levels that are separated by a few mean level spacings, ω ∆. In this regime, sin 2 (s) ∼ 1/2, can be averaged over a few level spacings, and the sine-kernel behaves like 1 − 1/(2s 2 ). Note that while the ramp in Figure 5 only starts at ergodic time, t er , as we will discuss in detail in the next section. At even longer times, the regime s 1 is explored giving rise to the plateau. Having explained the ergodic limit in the SYK model and the behaviour of the operators in this limit, we next move on to understand the non-universal deviations of the observables from the RMT-like behaviour at earlier times in the following section.

Non-ergodic modes
While the SYK model is defined as a disordered Hamiltonian, that is via an average over coupling constants, it is still far from being a random matrix model where every independent entry of the Hamiltonian is a random variable, as emphasised previously in the introduction. Nevertheless after the ergodic time (c.f. Figure 5) every ergodic chaotic quantum system is thought to be well described by random matrix theory.
This includes, as a special case, the quantum chaotic SYK models, i.e. (2.8). Up to this point, we have concerned ourselves with understanding how the full RMT physics originates in this system and what are its phenomenological effects on the operator correlation functions. However, in a generic quantum chaotic system one expects to see deviations from the universal RMT behaviour at early enough times. This is no different in the SYK class of models, with the added benefit that we can analytically study the non-ergodic regime as well. While, as we saw, the RMT physics comes from Hilbert space homogeneous modes, the non-universal modes correspond to the projection of the fieldĀ in the non-homogeneous directions,ā µ =0 . Previously, we had made an ansatz, (2.20), to study the saddle point equations arising from (2.18). Fluctuations around the saddle point solution are captured byĀ = y 1 + These fluctuations contribute to both the terms of the action (2.18). The first quadratic term contributes following additional term to the quadratic action of (2.21), We expand the second term in (2.18) perturbatively in the fieldsā µ k . We do this by writing the matrix K that appears in (2.18) as a sum of various components, representing the homogeneous term, frequency ω-dependent terms, sources and non-homogeneous terms, respectively. The last three terms are treated perturbatively in the expansion of the STr[ln K] term. The tracelessness and orthogonality of the X µ k matrices imply that only terms that contain even numbers of each of these matrices are non-zero. Thus the action for these modes is even and there are no tadpole diagrams that contribute perturbatively (see Appendix C for more details, especially the discussion preceding (C.6)). Therefore, these modes can be treated perturbatively around the saddle point, (2.22), as well as the entire σ-model manifold. This provides a selfconsistent justification of the ansatz. Following the detailed analysis of these massive modes presented in Appendix C, one finds that for generic observables the two point function is given by, where, C RMT (t) is given by (2.40). The function, f (x), decays exponentially for x 1. Therefore, individual terms in (3.3) decay at a rate determined by the mass, k (0). The sum in the above expression runs over all the basis elements of the D × D Hermitian matrices, (2.11). A given typical light operator has non-zero projection only over a few of these basis elements. The ergodic time that marks the transition between the universal late time physics and the non-universal early time physics corresponds to the minimum of the C(t). Quite importantly, it depends on the operator of interest. For a k-local operator, O = X ν with k = |ν|, as long as |N/2 − k| 1. In the Figure 5, we have plotted C(t) for the operator O = ψ i ψ j (i = j). These results complements the results of known literature, [29,[59][60][61], where the ergodic time for the spectral form factor was studied. Numerical studies of the ergodic time for operator correlation functions is beyond the scope of the present work and we leave it for future investigations, [62].

Summary & Discussion
Over the years, random-matrix physics has been thought to be an important ingredient in understanding the quantum physics of black holes. In recent years this connection has been reinforced, with the surprising realisation that semi-classical gravity itself appears to capture moments of the probability distribution governing this random-matrix like behaviour [4,5,33,[63][64][65]. In the present work we extended our understanding of the emergence of random matrix behaviour beyond the study of spectral indicators, in particular in operator correlation functions. Using the causal symmetry breaking approach of [4], we established that at late times two point functions of typical non-extensive operators demonstrate universality in their time dependence, succinctly captured in the operator sine kernel, (2.37). While the EFT description of the ergodic behaviourpost Thouless time -should apply more widely to ergodic theories (e.g. also to higher-dimensional holographic field theories), the approach to and exact value of the Thouless time are not universal. As an example of a holographic quantum system over which we have microscopic control, the SYK model provides a valuable case study that allows us to go further. Indeed we systematically derived the deviations from the universal ergodic behaviour at earlier times. These deviations depend on the theory under consideration, as well as on the operators under study. The techniques discussed in this work hold as long as the operator sources can be treated perturbatively as the analysis of the saddle point and the fluctuations around it is independent of the sources. As we have shown the deviations are controlled by a set of massive modes, inhomogeneous in the Hilbert space, and whose mass is directly related to the amount of inhomogeneity. The technology developed in our work to study the dominance of homogeneous solutions as well as the non-homogeneous mode that are non-universal is generalisable to other strongly interacting Hamiltonians. However, the exact implementation of this technology to more complicated theories is an interesting future challenge. Lastly, it is important to note that the result (2.38) is consistent with the prediction of ETH for the two point function: the off-diagonal contribution is suppressed by an e −S factor with respect to the diagonal contribution. However, we emphasise once more that this result encapsulates the information of ETH as well as spectral statistics in it. Comparing our results directly with other proposals of how quantum chaotic physics makes a contribution to gravity is therefore subtle, and we now undertake a somewhat detailed comparison to other related results that have appeared in the literature, before moving on to a brief outlook of interesting open issues and future directions.

Comparison with known results
The problem of understanding and characterising the onset of random matrix behaviour in manybody systems was addressed numerically in [43]. In the Jackiw-Teitelboim theory, [66,67], or dilaton-gravity in two dimensions, the emergence of RMT has been studied in [5,6]. In these studies it was argued that t er ∼ log(S) for the SYK model, where S is the entropy. This argument is supported by explicitly constructing states that contribute to a linear growth of the spectral form factor. However, in [29], as well as in the present work, it was demonstrated that the existence of non-universal modes masks this universal linear growth, making it not apparent in observables. Consequently, the linear ramp becomes detectable only after the last/lightest of these non-universal modes have decayed completely. The time scales corresponding to this transition depend on the operators under consideration. For example, t er ∼ √ S log(S) for the spectral form factor and t er ∼ S 3/2 for the kind of k-local operators considered in this work. In the double-scaled SYK model, this also agrees with the analytic studies of [61], where the more conventional collective field, or the "G, Σ" approach, is used to arrive at the same time scales. Comparisons with numerics in [59,60] lend support to this result. Our results are also consistent with [68,69]. Numerical studies for the operator correlation functions will be performed in [62]. Having addressed some subtle and perhaps technical issues that have arisen in our work and the literature, we would now like to discuss some of the more qualitative issues implied by the operator statistics derived in this paper.
In a recent work, [33], it has been proposed that the universal late time behaviour of observables in semiclassical gravity can be understood to arise from the Haar-averaging over states in a microcanonical energy window. As it turns out this proposal is quite similar to the physics dictated by our EFT, in particular in the case in which one projects the observables in a microcanonical window as described in subsection 2.1. Our results imply the picture that wormhole-like contributions in quantum gravity arise from quantum chaos via Haar averaging over eigenstates, but crucially the detailed analysis we presented demonstrates that observables receive contributions from both the spectral statistics as well as the Haar-averaging over the eigenstates.
Finally, let us comment that recently the reference [70] addressed a similar question of finding RMT statistics within physical theories by looking for the signatures of random-matrix behaviour at small frequencies (equivalently, at late times) in operator correlation functions. Our approach can also be used to analytically address the questions studied in this work.

Outlook and generalisations
The question of generalising the picture we have presented in this work to higher dimensions for theories of physical interest, such as N = 4 SYM remains a non-trivial one. We hope to address this question in future work. While the optimistic vision would be to use our techniques to understand thermalisation in N = 4 SYM to understand the questions about apparent thermalisation in high energy states as well as the corresponding holographic statements in the theory of gravity, a more prudent approach might be to apply them to the next simplest example on the ladder of complexity: the AdS 3 /CFT 2 correspondence. Combined with the results of [4], especially those pertaining to the minimal string theory discussed in that paper, we have a powerful tool promising to unravel the holographic duality with a finer scrutiny at late times.
In relation to the SYK model itself, a relevant question that is studied in [52] concerns the relation between the more widely familiar collective field treatment of the SYK model (originally developed by [35,71,72]) to the σ-model approach discussed in the present work. A unified approach will be quite useful to develop a coherent picture. The recent work, [57], proposes a string theoretic origin of the SYK model. It would be interesting to compare the matrix model described in their work with the emergent RMT discussed in our present work.
It would also be very interesting as well as important to study simple models akin to the SYK model which lie in other random matrix universality classes. Some of these can be achieved within the SYK model itself with different number of Majorana fermions, [42]. We hope to address these directions of research in the future. How a chaotic quantum system approaches the behaviour of a full dimH sized random matrix is an extremely important question in general and in particular in the context of holography, where such an understanding promises to microscopically resolve the unitary late time behaviour of black holes. It appears thus of crucial importance to obtain a more complete understanding of causal symmetry breaking in the bulk gravity, the minimal string example worked out in [4] being a low-dimensional prototype.

A Saddle point analysis
The discussion presented in this appendix reproduces the well known results of standard references, [49][50][51], for the reference of the reader. The path integral after taking the homogenous Hilbert-space ansatz in (2.21) is given by, where we have dropped the terms that we treat perturbatively in our analysis (the sources and the ω = 0 terms). The saddle point equation in the large D limit is, The solution of this equation is given by, where,Λ is a diagonal matrix with entries ±1. However, only a few of these solutions can be reached by contour deformations starting from the initial path integral. Let us denote by λ a,r b,f the eigenvalues of the y-matrix in the boson/fermion and advanced/retarded sectors, respectively. When integrating along λ a,r b in the bosonic sector of the y-path integral, one encounters poles at Note, that only one saddle point, either can be reached by a contour deformation in the retarded/advanced sectors, respectively (see Figure 6). worth noting here that STr(P f y 2 ) = −tr(P f y 2 ), thus the original integration contour runs along the imaginary axis in the complex plane of variables λ r/a f . It can be shown, see monograph [50], that the leading saddle point corresponds to, The other saddle point solution corresponds to, [47], T 0 ∈ U (1, 1|2)/U (1|1) ⊗ U (1|1), which transforms the standard saddle point, Λ to the sub-leading saddle Λ AA is given by, B Perturbation analysis around the saddle points

B.1 The σ-model with source terms
We begin our analysis by expanding the trace-log piece of the action, (A.1), containing the terms The terms in the above expansion are (1) at most linear in ω, and (2) at most quadratic in sources.
In addition, we use The saddle solution of y is here we have defined, τ 3 = σ RA 3 ⊗ 1 bf , as the Pauli matrix in the adv./ret. space. The perturbations of the saddle point by T ∈ U (1, 1|2)/U (1|1) × U (1|1) can be parametrised by, T matrix is perturbatively written as W is the pion field,

B.2 Gaussian integration around the standard saddle point
In this case the perturbations are parametrised by By defining the diagonal matrix where the bosonic eigenvalue e 1 and the fermionic one, e 3 , form the retarded sector, one can use it to write the individual terms as The first three terms are the kinetic terms for our perturbation theory and the last term as well as the higher-order terms that are denoted by the ellipsis can be treated perturbatively. However, we find that these higher order terms are not relevant due to non-renormalisation theorems that are applicable for the GUE. It is useful to write the propagator for the B,B fields as (using e 3 = e 1 , e 4 = e 2 ) where, X, Y are arbitrary supermatrices. Subsequently, we expand the source terms to quartic order in W . 10 However, to avoid clutter we explicitly write only the terms that are non-zero using the above contraction rules. and, So the path integral that we are interested in computing is given by (B.16) where, S source is the shorthand for all the source, h, dependent terms. Note that only the term linear in h contributes to the term that involves a single derivative with respect to h and the term quadratic in h to the second term. However, since we are working with traceless operators (B.14) is zero. The vertex in (B.16) can be diagrammatically represented in the double-line notation as shown in Figure 8. The coloured fat-line propagators demonstrate matrix fields B,B respectively in the double line notation. These B,B fields have one index each in advanced and the retarded sector as is apparent from its definition in (B.8), which is demonstrated by different colours of the double lines. The fields B,B have a concealed trace over the Fock-space, generated by the integration over the Hamiltonian (see (2.13)), which is represented by the blue dotted lines in this diagram. The insertion of projectors in the vertex is demonstrated by the black dots. The projectors in the definition of the source, (2.7), define how they change the causalities in this diagram. Evaluating the path integral one gets,

B.3 Gaussian integration around Andreev-Altshuler saddle point
In general, our analysis remains the same as above. However, there are some important qualitative differences between the two path integrals. Around the subleading saddle point, that can be obtained by a Weyl transformation of the standard saddle point as discussed in (A.8), the kinetic term for the Grassmann variables, µ,μ, ν,ν vanishes. Thus, for the Grassmann integrals to be nonvanishing, all four factors of the Grassmann variables should arise from the pre-exponent terms. This means that the first non-trivial contribution to the path integral comes from terms that are quartic in the pion fields W in the pre-exponential factors. Explicitly, where, once again the higher order terms don't contribute. Of all the terms arising in the expansion of the source term only the following terms contribute, Note that diagrammatically this vertex corresponds to the same Feynman diagram as given in Figure 8, however, this time the fields are expanded around the Andreev-Altshuler saddle point. All the above terms contribute equally to the path integral giving a total contribution, Adding this to the above answer, (B.17), and taking the real part we obtain: To this answer, we still need to add the contribution of the correlation function involving identical causalities.

B.4 Correlation function with two advanced Green's function
The correlation function that appears in the second term on the RHS of (2.4) can be defined as a path integral analogous to the first term, Note that in this case, the supervectorsΨ, Ψ comprise two copies of retarded space as opposed to a copy of advanced and retarded space each. It is not hard to see, following the discussion of Appendix A, that the saddle point for the path integral that arises in this case is The expansion in small ω to linear order and in sources to quadratic order is given by,

C Perturbation analysis of the massive modes
Let's go back to equation (2.18), where, where the quadratic term can we rewritten more explicitly as, Here, we have explicitly split the 'notational' µ-sum in (C.1) as a sum over cardinality of |µ| = k, as well as, a sum over N k choices of fermions for each k. We have also separated out the k = 0 (HHA) term from the remaining terms. We refer to the non-homogeneous modes as the massive modes and treat them perturbatively. The following analysis provides a self-consistent justification for this approach as well as the HHA ansatz that we had taken in section 2. We write K as a sum over various pieces, where the rest of the definitions are the same as in (B.1). To treat K M perturbatively, we need to assume K M K 0 . The linear and quadratic terms in the expansion of the logarithm are, where the linear term is zero because the matrices X µ are traceless. In fact, this implies that all odd powers in the logarithmic expansion are zero. The potential for the non-homogeneous fields is even with a saddle point solution at zero. This is a stable saddle point which justifies our perturbative treatment as well as homogeneous Hilbert space ansatz. Using the identity the quadratic action is It is useful to break the fieldsā µ k into diagonal and off-diagonal components, because they don't interact.
In what follows we proceed with the analysis of fluctuations around the standard saddle point. In this case the diagonal part ofā µ k commutes with Λ while the off-diagonal part anti-commutes. Thus we have an action in terms of these fields where, y * 0 = − E 2γ 1 − iΛ 1 − E 2 4γ 2 and y 0 y * 0 = 1. We simplify the action further to get, The quadratic action for the massive mode additionally receives contributions from the ω = 0 corrections, (C.10) The kinetic action for the massive modes is sum of (C.10) and (C.9), STr S(µ) −1 − 1 w µ k · w µ k − y 0 ·ê π 2 γρ(E) · w µ k · w µ k +STr S(µ) −1 1 − y 2 0 · v µ k · v µ k − y 0 ·ê · y 2 0 π 2 γρ(E) · v µ k · v µ k . (C.11) Note that the quadratic integral for massive modes, (C.11), are proportional to STr [W +− W −+ ], STr V ++ 2 and STr V −− 2 . It is easy (although cumbersome) to check that the terms that are quadratic in source, h, don't have non-zero massive propagators and therefore don't contribute to the order that we are interested in. The terms linear in h are, (C.12b) These terms represent the interaction of the sources with the massive modes and contribute to the pre-exponential terms ∂ h+ S src ∂ h− S src .
Tr O † X µ k 1 X µ k 2 Tr [OX µ k 3 X µ k 4 ] + · · · , (C. 13) where, "· · · " represents terms that evaluate to zero. An explicit computation shows that the terms proportional to Tr O † X µ k 1 X µ k 2 Tr [OX µ k 3 X µ k 4 ] evaluate to zero and the contribution of the massive modes to the correlation function is given by where, is the spectrum of the massive modes in the SYK model at energy E. Let us estimate this propagator for various massive modes. For all the above cases (µ = 0), assuming |N/2 − k| 1, 11 16) This shows that all modes are strongly overdamped, with their masses being much larger than the many-body level spacing. Performing a Fourier transform over the frequency, ω, to obtain the correlation function as a function of time, (C.17) One can integrate over the spectrum to obtain the contribution of the massive modes to the correlation function,C The last relation is a property of the 'Fourier' transform in the Hilbert space which is based on the orthogonality relation Tr[X µ X † ν ] = Dδ µν . Since, from (C.20), we know that the lightest modes correspond to k = 1, N − 1.
Contribution of the massive mode along with the ergodic answer, (2.40), gives the full contribution to the correlation function, For the two-body hopping operators, O = γ i γ j , typical plots of the full correlation function C(t) are shown in Figure 5.

C.1 Ergodic time
Let's further estimate the ergodic time, t er , at which the minimum of C(t) is reached. Assuming that O = X ν with k = |ν| and |N/2 − k| 1 along witht 1 the correlation function expressed via scaled timet = t/t H becomes, The minimum of this function is reached at timet er satisfying Omitting further all constants of order unity and using the asymptotics of Bessel function at large arguments, the above equation is reduced to

C.2 Evaluation of energy integral
Taking into account the energy dependence of ρ(E) and k (E), the energy integral which follows from (C.18) becomes where the dimensionless function f (x) reads