Spectral Form Factor of a Quantum Spin Glass

It is widely expected that systems which fully thermalize are chaotic in the sense of exhibiting random-matrix statistics of their energy level spacings, whereas integrable systems exhibit Poissonian statistics. In this paper, we investigate a third class: spin glasses. These systems are partially chaotic but do not achieve full thermalization due to large free energy barriers. We examine the level spacing statistics of a canonical infinite-range quantum spin glass, the quantum $p$-spherical model, using an analytic path integral approach. We find statistics consistent with a direct sum of independent random matrices, and show that the number of such matrices is equal to the number of distinct metastable configurations -- the exponential of the spin glass"complexity"as obtained from the quantum Thouless-Anderson-Palmer equations. We also consider the statistical properties of the complexity itself and identify a set of contributions to the path integral which suggest a Poissonian distribution for the number of metastable configurations. Our results show that level spacing statistics can probe the ergodicity-breaking in quantum spin glasses and provide a way to generalize the notion of spin glass complexity beyond models with a semi-classical limit.


Introduction
An isolated quantum many-body system which reaches an effective thermal equilibrium state starting from an out-of-equilibrium initial state is often called "quantum chaotic." As commonly used, quantum chaos is a loose term referring to a family of phenomena that typically co-occur, including the ability of the system to serve as its own heat bath [1][2][3], hydrodynamic behavior of conserved quantities [4][5][6][7][8], and randommatrix-like energy eigenvalues [9][10][11][12]. Given this variety, it is crucial to understand the relationships between different manifestations of quantum chaos [13,14].
These relationships are complicated and interesting in large part because the systems in question have structure, such as locality and symmetry. For example, if the Hamiltonian has spatial locality, energy conservation implies the existence of slow hydrodynamic modes and an associated long time scale, the Thouless time, such that random-matrix behavior is only present for energy levels closer than the inverse Thouless time [15,16]. Similarly, if the Hamiltonian possesses a symmetry, then it can be organized into blocks labelled by irreducible representations of the symmetry. One finds random-matrix statistics within each individual block, but full ergodicity is broken because matrix elements between different blocks are forbidden [17][18][19][20][21].
It is natural to ask whether there are other ways in which ergodicity can be lost, and if so, what the resulting spectral statistics of the Hamiltonians are. In particular, we will better understand the relations between different measures of quantum chaos by understanding how they are lost and what replaces them.
Quantum spin glasses provide one well-established context to explore these questions, since they exhibit a rich phenomenology associated with the inability to fully thermalize [22][23][24][25][26][27][28] In this paper, we determine the spectral statistics of an analytically tractable spin glass model, the quantum p-spherical model. We find that up to times polynomial in the system size, the Hamiltonian can effectively be described as approximately block-diagonal. Each block behaves as a random matrix independent of the others, and the number of blocks depends on the energy per particle. At high energies, there is only one block and the system is ergodic. Below a critical energy density, the Hamiltonian breaks into exponentially many blocks -the average number of blocks jumps discontinuously from the high energy regime and then decreases as the energy density decreases further. We establish these results via a path integral computation of the spectral form factor (SFF), which measures correlations between pairs of energy levels [29][30][31][32][33].
In the remainder of the introduction, we give some physical context by reviewing the spectral form factor and mean-field spin glasses, and then summarize our results. In Sec. 2, we review the p-spherical model in detail. In Sec. 3, we calculate the SFF of this model in the high-temperature ergodic regime, and in Sec. 4, we do so in the non-ergodic regime. Finally, in Sec. 5, we investigate higher-moment analogues of the SFF. We then discuss implications of these results and directions for future work in Sec. 6.

Review of the spectral form factor
To study the spectral correlations of a Hamiltonian H, a standard tool is the spectral form factor (SFF) [34,35], defined as SFF(T ) ≡ Tre −iHT 2 .
In situations where the spectrum is unbounded, or when one wishes to concentrate on a portion of the spectrum, the trace in Eq. (1) is regulated by a filter function f (H): One common choice is f (H) = e −βH [29,36], and another is f (H) = e −c(H−E0) 2 . The latter allows one to study level statistics near a specified energy E 0 . For a single Hamiltonian, the SFF is an erratic function of time [35]. Thus one usually considers an ensemble of Hamiltonians and defines the SFF as the average of Eq. (2) over the ensemble. Throughout this paper, we use the notation E[ · ] to denote the ensemble average.
The SFF is closely related to the correlation function of the density of states. Formally, the (filtered) density of states is given by where n labels the eigenstate of H with eigenvalue E n , and its correlation function is The SFF is simply the Fourier transform of the correlation function with respect to ω, integrated over E (although the filter function allows one to concentrate on an arbitrary subset of the spectrum). It is conceptually useful to split the SFF into two contributions: The first term, the disconnected piece of the SFF, comes solely from the average density of states. It is the second term, the connected piece, that contains information on the correlation between energy levels. The assertion of "random matrix universality" [9,37] can be phrased as the statement that an ensemble of quantum chaotic Hamiltonians will generically have the same connected SFF as the canonical Gaussian ensembles of random matrix theory [11,38]. This conjectured universal behavior is illustrated in Fig. 2, which plots the disorder-averaged SFF of the Gaussian unitary ensemble (one of the aforementioned canonical ensembles). Note the three distinct regimes: • The "dip", occurring at short times, comes from the disconnected piece of the SFF (and thus its precise shape is non-universal). It reflects a loss of constructive interference -the different terms of Tre −iHT acquire different phase factors as T increases.
• The "ramp", occurring at intermediate times, is arguably the most interesting regime. In the canonical matrix ensembles, it is a consequence of the result [11] E ρ E + where b = 1, 2, 4 in the orthogonal, unitary, and sympletic ensembles respectively [11]. The right-hand side being negative is a reflection of the well-known level repulsion in quantum chaotic systems [39].
Taking the Fourier transform with respect to ω gives a term proportional to T for the connected SFF. Such a linear-in-T ramp is often taken as a defining signature of quantum chaos.
• The "plateau", occurring at late times, results from the discreteness of the spectrum. At times much larger than the inverse level spacing, one expects that all off-diagonal terms in the double-trace of the SFF sum to effectively zero, meaning that As the plateau regime is both challenging to access analytically and not particularly informative, we shall not consider it further in this work.
The bulk of our analysis in this paper is devoted to calculation of the ramp in a well-known quantum spin glass model, the p-spherical model (discussed below). The results can be understood via the elementary observation that when a Hamiltonian is block diagonal, then Tre −iHT = k Tre −iH k T . If the different blocks are independent, then the variance of Tre −iHT is the sum of the variance of each Tre −iH k T , i.e., the SFF is the sum of the SFF for each block. In particular, the coefficient of the universal linear-in-T ramp is multiplied by the number of independent blocks. Systems with only approximately block-diagonal Hamiltonians, for which there are small matrix elements between blocks, have this enhancement of the ramp up to the transition timescale between blocks. For a more detailed analysis, see Ref. [19].

Review of mean-field spin glasses
Broadly speaking, spin glasses are systems in which the magnetic moments σ i are frozen but disordered at low temperatures. However, this definition (much like that of "quantum chaos") encompasses a wide variety of phenomena which are in many ways quite distinct, as is made clear by the literature on the subject [22][23][24][25][26][27][28].
In the present paper, we focus on what are known as "one-step replica symmetry breaking" (1RSB) spin glass phases [27]. We are specifically interested in quantum spin glasses, but we first review the corresponding classical case, for which configurations are labelled by a list σ ≡ {σ 1 , · · · , σ N } and the Hamiltonian is simply a function of σ.
While the technical definition of 1RSB is somewhat involved, the qualitative physics is straightforward to understand and captured by the sketch in Fig. 3. The energy landscape, i.e., energy as a function of spin configuration, has many deep wells and steep barriers. In particular, the number of wells is e O(N ) and the heights of the energy barriers separating wells are O(N ), where N is the number of spins. As a result, below a certain energy density d , the system is extremely non-ergodic: it remains trapped within an exponentially small fraction of the thermodynamically relevant configuration space until exponentially long timescales. While the 1RSB phenomenon was originally studied in the context of stochastic classical dynamics [41][42][43][44], it has recently been shown to imply exponentially long tunneling timescales for isolated quantum dynamics as well [45][46][47][48][49].

Ergodic Non-ergodic, delocalized
Non-ergodic, localized Figure 3: (Left) Cartoon of the energy landscape in a 1RSB spin glass. The y-axis is energy per spin, E/N , where E is energy and N is the number of spins. Different points on the x-axis represent (very roughly, since the actual configuration space is N -dimensional) different spin configurations σ. The dashed line indicates the energy density d below which the system is non-ergodic. (Right) Sketch of the dynamical phase diagram for a quantum 1RSB spin glass. The x-axis represents parameters controlling the strength of quantum fluctuations, and the y-axis is energy density. Note that many other types of phase transitions are also present, in particular equilibrium transitions, but are not indicated here. See, e.g., Refs. [26,27,40] for more information.
TAP states (named after Thouless, Anderson, and Palmer [50]) provide a more quantitative description of such "deep wells". Arguably the most general definition (see Ref. [25] for others) is in terms of the Legendre transform of the free energy with respect to local fields: where H is the Hamiltonian of interest and the fields {h i } are chosen so that σ i = m i (where · indicates a thermal average). TAP states are simply the local minima of F ({m i }). Physically, each corresponds to a different "well" of the energy landscape, including thermal fluctuations around the lowest point (thus TAP states do generically depend on temperature). The partition function can be decomposed as a sum over TAP states: where α denotes a TAP state and δ σ∈α restricts the trace to only those states belonging to TAP state α. Note that in this discussion, σ can refer to any set of degrees of freedom: Ising spins, vector spins, continuous coordinates, etc. In all cases, Eqs. (10) and (11) can be interpreted accordingly. Quantum generalizations of spin glasses are usually obtained by adding non-commuting terms to the Hamiltonian. For example, with an Ising Hamiltonian, one often interprets σ i as the Pauli spin-z operator σ z i and includes an additional transverse field Γ i σ x i [51][52][53][54]. On the other hand, with systems having continuous degrees of freedom (including the one which we study in this paper), one can interpret σ i as a position coordinate and include the "kinetic energy" i π 2 i /2µ, where π i is the momentum operator conjugate to σ i [55,56]. Generically, the resulting system has a frozen spin glass phase at low energy and small quantum fluctuations (the latter being controlled by Γ and µ −1 respectively in the examples above), and has a paramagnetic phase at either high energy or large quantum fluctuations. A sketch of the typical phase diagram is shown in Fig. 3, with these two phases indicated by "non-ergodic" and "ergodic".
It has recently been noted that quantum 1RSB spin glasses can exhibit eigenstate phase transitions which are distinct from the above [57][58][59]. Qualitatively speaking, on the low energy/fluctuation side of the eigenstate phase boundary, each eigenstate of the Hamiltonian is localized on a single TAP state. This implies that under the system's internal dynamics alone (i.e., as given by the Schrodinger equation), the system cannot tunnel between TAP states on any timescale, even times exponential in the number of spins. On the other side of the phase boundary, each eigenstate is delocalized over many TAP states in accordance with random matrix behavior. As discussed in Ref. [48], while this implies that the system does tunnel between TAP states, the timescale for tunneling is necessarily exponential in system size, analogous to the activation times under open-system dynamics. Only when there exists a single TAP state can one identify the phase as genuinely thermalizing. As a result, one finds phase diagrams like that sketched in Fig. 3, with "non-ergodic"/"ergodic" indicating whether multiple TAP states exist and "localized"/"delocalized" referring to the eigenstate properties.

Summary of results
In this paper, we calculate the SFF for a particular ensemble of quantum spin glasses, the quantum pspherical model (PSM) [55,56,60]. We find that in the ergodic phase, the connected part of the SFF agrees with the expectation from random matrix theory (Eq. (62) below), while in the non-ergodic phase, it is enhanced by a factor which is precisely the number of TAP states (Eq. (109)). Given the discussion in Secs. 1.1 and 1.2, this makes precise and validates the idea that each metastable state (i.e., TAP state) corresponds to a block of the Hamiltonian that is quantum chaotic on its own but is nearly decoupled from all others, thus making the system as a whole non-ergodic [58]. This is the main result of the present work.
We also consider higher moments of the evolution operator and identify a set of saddle points (Eq. (134)) which, in addition to confirming the picture that different TAP states have independent level statistics, suggest that at least at low complexity, the number of TAP states at a given energy is Poisson-distributed and independent of other energies. Yet as we shall discuss, since our analysis does not consider the perturbative corrections around each saddle point, this does not constitute a complete calculation and serves more as motivation for future investigation.
2 Real-time dynamics of the quantum p-spherical model

The model
The classical p-spherical model (PSM) [61] is a disordered spin model with all-to-all p-body interactions. It is defined by the classical Hamiltonian where the couplings J i1...ip are independent Gaussian random variables with mean zero and variance Here and throughout, E indicates an average over couplings. The notation (i 1 · · · i p ) denotes sets of p indices such that 1 ≤ i 1 ≤ · · · ≤ i p ≤ N . The sum in Eq. (12) is over all such sets. Our treatment differs from the standard convention by including a parameter J for the overall strength of the disorder. To recover the standard expressions, simply set J 2 = p/2. We also include the combinatorial factor C i1···ip = 1≤i≤N n i !, where n i is the number of indices set equal to i. This term is almost always one, but its inclusion avoids 1/N corrections in the action. The σ i are real, continuous spin variables subject to the spherical constraint which ensures that the system has an extensive free energy. It is apparent that this is a mean-field model without any spatial structure. This allows for infinite free energy barriers around metastable states in the thermodynamic limit, making the model ideal for examining the impact of metastability on the spectral statistics of spin glasses.
In this work, we follow Refs. [55,56,60] in generalizing Eq. (12) to a quantum Hamiltonian H. We treat the σ i as commuting position operators, and define conjugate momentum operators π i which satisfy the commutation relations The quantum PSM simply includes a kinetic energy term in the Hamiltonian: The mass µ is an additional parameter controlling the strength of quantum fluctuations. To incorporate the spherical constraint, we take the Hilbert space to be the subspace in which i σ 2 i has eigenvalue N . The quantum PSM may be interpreted as a soft-spin version of the Ising p-spin model in an external transverse field -itself the subject of much study [53,[62][63][64][65] -where µ −1 is analogous to the transverse field. Alternatively, if we think of σ ≡ {σ 1 , · · · , σ N } as a position vector in N -dimensional space, the quantum PSM has a natural interpretation as a particle of mass µ moving on a hypersphere of radius √ N . This particle experiences the Gaussian random potential whose correlation function is Note that there is a very important difference between p = 2 and p > 2: the former is a Gaussian model, essentially (but for the spherical constraint) a system of linearly coupled harmonic oscillators. It therefore has qualitatively different behavior than the p > 2 models, which are genuinely interacting and serve as reasonable toy models for rugged energy landscapes. In this work, we exclusively consider p > 2.

Schwinger-Keldysh path integral
Just as other all-to-all models have a saddle-point/mean-field description at large N , so too does the PSM. We start with the disorder-averaged (i.e., "annealed") path integral on the Schwinger-Keldysh contour at inverse temperature β, illustrated in the left column of Fig. 4. While it is in general incorrect (often grossly) to disorder-average the path integral itself, it is known that the annealed approximation is accurate in the PSM as long as β is less than a critical value β s [26,66]. We shall assume that this is true throughout. The annealed path integral is where For brevity, we use C to denote the entire contour. Thus C dt indicates a contour integral within the complex-t plane. The Lagrange multiplier z(t) is included to enforce the spherical constraint. It can be interpreted as a time-dependent harmonic potential whose value is chosen such that i σ i (t) 2 = N at all times. Thus the measure Dσ N is simply the product measure over each σ i independently. From here, the same manipulations used to get Schwinger-Dyson equations for the SYK model will give us equations of motion for the PSM. One can immediately perform the Gaussian integrals over the couplings to obtain Spectral form factor Schwinger-Keldysh

Real-time correlation functions
Integration contour

Equations of motion
Re t < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 + U Q M R / C Z j C Q g f 7 y y J L E 2 9 I 2 t s U = " > A A A C P X i c b Z D P S x t B F M d n o 1 W b 1 h r 1 2 M t g K C Z g w 2 4 p 6 E U Q P O j R Q n 5 I k 3 R 5 O 3 m J g z O 7 y 8 x b M S z 7 j 3 n x f + i t N y 8 e F P H a q 5 O Y g x o f D H z 4 f t 9 j 3 v t G q Z K W f P + f V 1 p Y / L C 0 v P K x / O n z 6 p e 1 y v p G 2 y a Z E d g S i U r M a Q Q W l Y y x R Z I U n q Y G Q U c K O 9 H 5 4 c T v X K C x M o m b N E 6 x r 2 E U y 6 E U Q E 4 K K 8 2 j M M + y o k Y 7 n L b r f J / / / p N / D w r e I 7 w k o / O m K T g 6 q R c h A T 9 2 h g Y 6 E 6 D y p m M r R x p C W a P 6 C 9 6 u h 5 W q 3 / C n x e c h m E G V z e o k r P z t D R K R a Y x J K L C 2 G / g p 9 X M w J I X C o t z L L K Y g z m G E X Y c x a L T 9 f H p 9 w b 8 5 Z c C H i X E v J j 5 V X 0 7 k o K 0 d 6 8 h 1 T n a 3 b 7 2 J + J 7 X z W i 4 1 8 9 l n G a E s X j + a J g p T g m f R M k H 0 q A g N X Y A w k i 3 K x d n Y E C Q C 7 z s Q g j e n j w P 7 R + N w G 8 E v 3 5 W D / x Z H C v s K 9 t i N R a w X X b A j t k J a z H B r t g N u 2 P 3 3 r V 3 6 z 1 4 j 8 + t J W 8 2 s 8 l e l f f / C V 7 i r O M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 + U Q M R / C Z j C Q g f 7 y y J L E 2 9 I 2 t s U = " > A A A C P X i c b Z D P S x t B F M d n o 1 W b 1 h r 1 2 M t g K C Z g w 2 4 p 6 E U Q P O j R Q n 5 I k 3 R 5 O 3 m J g z O 7 y 8 x b M S z 7 j 3 n x f + i t N y 8 e F P H a q 5 O Y g x o f D H z 4 f t 9 j 3 v t G q Z K W f P + f V 1 p Y / L C 0 v P K x / O n z 6 p e 1 y v p G 2 y a Z E d g S i U r M a Q Q W l Y y x R Z I U n q Y G Q U c K O 9 H 5 4 c T v X K C x M o m b N E 6 x r 2 E U y 6 E U Q E 4 K K 8 2 j M M + y o k Y 7 n L b r f J / / / p N / D w r e I 7 w k o / O m K T g 6 q R c h A T 9 2 h g Y 6 E 6 D y p m M r R x p C W a P 6 C 9 6 u h 5 W q 3 / C n x e c h m E G V z e o k r P z t D R K R a Y x J K L C 2 G / g p 9 X M w J I X C o t z L L K Y g z m G E X Y c x a L T 9 f H p 9 w b 8 5 Z c C H i X E v J j 5 V X 0 7 k o K 0 d 6 8 h 1 T n a 3 b 7 2 J + J 7 X z W i 4 1 8 9 l n G a E s X j + a J g p T g m f R M k H 0 q A g N X Y A w k i 3 K x d n Y E C Q C 7 z s Q g j e n j w P 7 R + N w G 8 E v 3 5 W D / x Z H C v s K 9 t i N R a w X X b A j t k J a z H B r t g N u 2 P 3 3 r V 3 6 z 1 4 j 8 + t J W 8 2 s 8 l e l f f / C V 7 i r O M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 + U Q M R / C Z j C Q g f 7 y y J L E 2 9 I 2 t s U = " > A A A C P X i c b Z D P S x t B F M d n o 1 W b 1 h r 1 2 M t g K C Z g w 2 4 p 6 E U Q P O j R Q n 5 I k 3 R 5 O 3 m J g z O 7 y 8 x b M S z 7 j 3 n x f + i t N y 8 e F P H a q 5 O Y g x o f D H z 4 f t 9 j 3 v t G q Z K W f P + f V 1 p Y / L C 0 v P K x / O n z 6 p e 1 y v p G 2 y a Z E d g S i U r M a Q Q W l Y y x R Z I U n q Y G Q U c K O 9 H 5 4 c T v X K C x M o m b N E 6 x r 2 E U y 6 E U Q E 4 K K 8 2 j M M + y o k Y 7 n L b r f J / / / p N / D w r e I 7 w k o / O m K T g 6 q R c h A T 9 2 h g Y 6 E 6 D y p m M r R x p C W a P 6 C 9 6 u h 5 W q 3 / C n x e c h m E G V z e o k r P z t D R K R a Y x J K L C 2 G / g p 9 X M w J I X C o t z L L K Y g z m G E X Y c x a L T 9 f H p 9 w b 8 5 Z c C H i X E v J j 5 V X 0 7 k o K 0 d 6 8 h 1 T n a 3 b 7 2 J + J 7 X z W i 4 1 8 9 l n G a E s X j + a J g p T g m f R M k H 0 q A g N X Y A w k i 3 K x d n Y E C Q C 7 z s Q g j e n j w P 7 R + N w G 8 E v 3 5 W D / x Z H C v s K 9 t i N R a w X X b A j t k J a z H B r t g N u 2 P 3 3 r V 3 6 z 1 4 j 8 + t J W 8 2 s 8 l e l f f / C V 7 i r O M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 + U Q M R / C Z j C Q g f 7 y y J L E 2 9 I 2 t s U = " > A A A C P X i c b Z D P S x t B F M d n o 1 W b 1 h r 1 2 M t g K C Z g w 2 4 p 6 E U Q P O j R Q n 5 I k 3 R 5 O 3 m J g z O 7 y 8 x b M S z 7 j 3 n x f + i t N y 8 e F P H a q 5 O Y g x o f D H z 4 f t 9 j 3 v t G q Z K W f P + f V 1 p Y / L C 0 v P K x / O n z 6 p e 1 y v p G 2 y a Z E d g S i U r M a Q Q W l Y y x R Z I U n q Y G Q U c K O 9 H 5 4 c T v X K C x M o m b N E 6 x r 2 E U y 6 E U Q E 4 K K 8 2 j M M + y o k Y 7 n L b r f J / / / p N / D w r e I 7 w k o / O m K T g 6 q R c h A T 9 2 h g Y 6 E 6 D y p m M r R x p C W a P 6 C 9 6 u h 5 W q 3 / C n x e c h m E G V z e o k r P z t D R K R a Y x J K L C 2 G / g p 9 X M w J I X C o t z L L K Y g z m G E X Y c x a L T 9 f H p 9 w b 8 5 Z c C H i X E v J j 5 V X 0 7 k o K 0 d 6 8 h 1 T n a 3 b 7 2 J + J 7 X z W i 4 1 8 9 l n G a E s X j + a J g p T g m f R M k H 0 q A g N X Y A w k i 3 K x d n Y E C Q C 7 z s Q g j e n j w P 7 R + N w G 8 E v 3 5 W D / x Z H C v s K 9 t i N R a w X X b A j t k J a z H B r t g N u 2 P 3 3 r V 3 6 z 1 4 j 8 + t J W 8 2 s 8 l e l f f / C V 7 i r O M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " y 3 n p l V w b B r r e s T Q d v q i 9 e B u 5 m j g = "     where

g Q Y 9 I K R C H v Z E u c k 0 x a A v M 8 = " > A A A C G n i c b V D L S g M x F M 3 U V 6 2 v q k s 3 w S J W q G V S B N 0 I B U H F V Q X 7 g L 7 I p G k b m n k 0 u S O U o d / h x l 9 x 4 0 I R d + L G v z F t Z 6 H V A 4 G T c 8 4 l u c c J p N B g 2 1 9 W Y m F x a X k l u Z p a W 9 / Y 3 E p v 7 1 S 0 H y r G y 8 y X v q o 5 V H M p P F 4 G A Z L X A s W p 6 0 h e d Q Y X E 7 9 6 z 5 U W v n c H o 4 A 3 X d r z R F c w C k Z q p 8 l l F n I Y D o / w O b 5 p F f B V f G 1 F w T E Z 5 3 B j O A x p J 5 Y n I d J O Z + y 8 P Q X + S 0 h M M i h G q Z 3 + a H R 8 F r r c A y a p 1 n V i B 9 C M q A L B J B + n G q H m A W U D 2 u N 1 Q z 3 q c t 2 M p q u N 8 Y F R O r j r K 3 M 8 w F P 1 5 0 R E X a 1 H r m O S L o W + n v c m 4 n 9 e P Y T u W T M S X h A C 9 9 j s o W 4 o M f h 4 0 h P u C M U Z y J E h l C l h / o p Z n y r K w L S Z M i W Q + Z X / k k o h T + w 8 u T 3 J F O 2 4 j i T a Q / s o i w g 6 R U V 0 j U q o j B h 6 Q E / o B b 1 a j 9 a z 9 W a 9 z 6 I J K 5 7 Z R b 9 g f X 4 D T C 2 b Y Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 y a k 3 g Q Y 9 I K R C H v Z E u c k 0 x a A v M 8 = " > A A A C G n i c b V D L S g M x F M 3 U V 6 2 v q k s 3 w S J W q G V S B N 0 I B U H F V Q X 7 g L 7 I p G k b m n k 0 u S O U o d / h x l 9 x 4 0 I R d + L G v z F t Z 6 H V A 4 G T c 8 4 l u c c J p N B g 2 1 9 W Y m F x a X k l u Z p a W 9 / Y 3 E p v 7 1 S 0 H y r G y 8 y X v q o 5 V H M p P F 4 G A Z L X A s W p 6 0 h e d Q Y X E 7 9 6 z 5 U W v n c H o 4 A 3 X d r z R F c w C k Z q p 8 l l F n I Y D o / w O b 5 p F f B V f G 1 F w T E Z 5 3 B j O A x p J 5 Y n I d J O Z + y 8 P Q X + S 0 h M M i h G q Z 3 + a H R 8 F r r c A y a p 1 n V i B 9 C M q A L B J B + n G q H m A W U D 2 u N 1 Q z 3 q c t 2 M p q u N 8 Y F R O r j r K 3 M 8 w F P 1 5 0 R E X a 1 H r m O S L o W + n v c m 4 n 9 e P Y T u W T M S X h A C 9 9 j s o W 4 o M f h 4 0 h P u C M U Z y J E h l C l h / o p Z n y r K w L S Z M i W Q + Z X / k k o h T + w 8 u T 3 J F O 2 4 j i T a Q / s o i w g 6 R U V 0 j U q o j B h 6 Q E
Next introduce a "fat unity", The integral over the self-energy F(t, t ) runs along the imaginary axis, making the second line simply the identity dpe ipx = 2πδ(x) (we absorb factors of 2π into the measure DF). However, when we ultimately evaluate the path integral by saddle point, we shall find that the saddle point value of F(t, t ) is real. Inserting Eq. (23) into the path integral gives where We can now perform the integral over σ i , resulting in where At large N , the remaining path integral can be evaluated within the saddle point approximation. The locations of the saddle points are determined by setting to zero the functional derivatives of Eq. (27): Keep in mind that the time arguments in Eq. (28) are complex and range over the entire Schwinger-Keldysh contour. In particular, although it is hidden in this compact notation, the infinitesimals dt acquire different phases depending on the branch of the contour: dt is a positive real infinitesimal on the upper ("forward") real-time branch, a negative real infinitesimal on the lower ("backward") real-time branch, and a negative imaginary infinitesimal on the thermal branch. G(t, t ) is the order parameter of this theory. As is clear from the manner by which it was introduced (top line of Eq. (23)), expectation values of G(t, t ) within the path integral are equivalent to expectation values of N −1 i σ i (t)σ i (t ). The latter are simply time-ordered correlation functions. We shall focus on the real-time correlation functions, for which it is more transparent to explicitly indicate the branches by α ∈ {u, l} and have t be simply a real variable. Formally, we have that where T denotes time ordering and T denotes time anti-ordering. Note that we can omit the sum over i because the different spins (upon disorder-averaging) have equivalent behavior. A number of formal properties of G αα (t, t ) are evident from Eq. (29). For one thing, G αα (t, t ) clearly depends only on the time difference t − t , and we shall often write G αα (t) with t = 0. Since the four components differ only in time ordering, we see that for any function f (x), We can further express all four components in terms of a single complex-valued function (equivalently two real-valued functions). For example, write G lu (t) in terms of its real and imaginary parts as G R (t) + iG I (t).
Since G lu (t) * = G lu (−t), G R (t) is even and G I (t) is odd. One can easily confirm that One of the most important features of G αα (t, t ) is the limiting behavior at large |t − t |, as a function of the inverse temperature β. Numerical solution of Eq. (28) demonstrates that there is a critical value β d (which is less than β s ): • For β < β d , lim |t−t |→∞ G αα (t, t ) = 0. We call this the "ergodic" phase (E σ i (t) = 0 by symmetry regardless of temperature, and so in this phase • For β d < β < β s , lim |t−t |→∞ G αα (t, t ) = q EA > 0. We call this the "non-ergodic" phase. The quantity q EA is referred to as the "Edwards-Anderson" order parameter.
• For β s < β, our initial annealed approximation is no longer valid. The replica trick is required to obtain accurate results [23,24], but (at least for finite-time dynamical properties) the behavior is qualitatively similar to that of the non-ergodic phase.

TAP equations on the Schwinger-Keldysh contour
The dynamical calculation described above only hints at the complexity of the non-ergodic phase. A more complete picture emerges from a generalization in the spirit of the TAP equations. Our treatment follows that of Ref. [67], which derived TAP equations on the thermal circle for the quantum PSM. While the extension to real-time dynamics is straightforward, we are not aware of any explicit calculation in the literature. Thus we present a detailed derivation of the following equations in App. A. As discussed in Sec. 1.2, the TAP free energy (or Gibbs potential) is the Legendre transform of the free energy with respect to local fields. It is therefore a function of the magnetization m i of each spin. For the free energy of quantum systems, the magnetization should also have an imaginary time index m i (τ ). The imaginary-time correlation function G(τ, τ ) becomes an additional order parameter.
We define the TAP action on the Schwinger-Keldysh contour analogously. It is a function of the magnetizations m i (t) and the correlation function G(t, t ), with t again being complex-valued and ranging over the entire contour. Specifically, where C denotes the Schwinger-Keldysh contour and The fields h i (t) and Λ(t, t ) are not independent parameters. They are instead chosen so that σ i (t) = m i (t) and N −1 i σ i (t)σ i (t ) = G(t, t ), just as z(t) is again chosen to enforce N −1 i σ i (t) 2 = 1, where the expectation value is with respect to the action in Eq. (32).
Due to the Legendre-transform structure of S TAP , we have that The TAP equations are those for m i (t) and G(t, t ) which one gets by setting the right-hand sides of Eq. (34) to zero. The solutions are therefore the values of magnetization and correlation function which the system can consistently possess "on its own," without any external fields. In this sense, each solution corresponds to a distinct metastable state. There is no reason why there cannot be many self-consistent solutions, and indeed, spin glass models such as the PSM do have many at sufficiently low temperature. We calculate the TAP equations in App. A. They are simplified by the fact that we can take m i (t) = m and z(t) = z. We also define q EA ≡ N −1 i m 2 i . The equations come out to be (together with G(t, t) = 1) Note that Eq. (36) is N equations, one for each spin i, and that it holds equally for any value of t due to time translation invariance. Defining F(t, t ) ≡ J 2 G(t, t ) p−1 , Eq. (35) is quite similar to Eq. (28). The only difference is that Eq. (35) uses ∆G(t, t ) ≡ G(t, t ) − q EA and ∆F(t, t ) ≡ F(t, t ) − J 2 q p−1 EA , which decay to zero at large |t − t |, rather than G(t, t ) and F(t, t ) themselves.
Despite the more involved derivation, G(t, t ) remains a contour-ordered expectation value. Thus, returning to the notation in which α ∈ {u, l} labels branches and t is real, G αα (t − t ) possesses the same formal properties as discussed in the previous subsection (Eqs. (30) and (31)). Of particular importance will be the Fourier transform of ∆G αα (t) at zero frequency, denoted ∆ G αα (0), as well as its (matrix) inverse, ∆ G −1 αα (0). Also define L ≡ ∞ −∞ dt∆G R (t) and Λ ≡ ∞ 0 dt∆G I (t). Then from Eq. (31), we see that The multiplicity of solutions to the TAP equations comes from Eq. (36). By use of Eqs. (35), (37), and (38), it can be written (associating u with 0 and l with 1) Eq. (39) is identical to that which appears and has been well-studied for the classical PSM [26,42,61,68]. Thus we simply quote the following results. In addition to the inverse temperature β, solutions to Eq. (39) are parametrized by the quantity which can be interpreted as a "normalized" potential energy density: each magnetization has a value which is (very roughly) comparable to q 1/2 EA , and thus the natural scale for the interaction energy is Jq EA . The value of q EA for a given E < 0 is given by the largest solution to − 2Jq where Λ depends on q EA through Eq. (35).
The exponent Σ(E) is referred to as the "complexity" in the spin glass literature. The connection between this TAP approach and the conventional Schwinger-Keldysh path integral lies in the fact that: i) the inverse temperature β d at which TAP states with non-zero magnetization appear is identical to that at which the autocorrelation function acquires a non-zero late-time limit; ii) the overlap determined by Eq. (41) is identical to the late-time value of the autocorrelation function. This strongly suggests the following picture: • For β < β d (the "ergodic" phase), there exists a single equilibrium state with zero magnetization, and the correlation function decays to zero on a finite timescale.
• For β d < β (the "non-ergodic" phase), there exist exponentially many metastable states having non-zero magnetization. The number of states is given by the exponential of the complexity Σ(E). Dynamically, in the N → ∞ limit, a system prepared in one metastable state will remain in that state for all time. At finite N , it is only on a timescale exponential in N that the system can transition between states.
Much more can be said about these phases (in particular how the replica-symmetry-breaking transition at β s appears within the TAP approach), and a large body of literature is devoted to this topic. We refer in particular to Ref. [26] as an excellent starting point.
In the present work, we determine the spectral statistics of the PSM in both the ergodic and non-ergodic phase. Those of the former can be computed very much along the lines of Ref. [29], which we do in Sec. 3. Those of the latter, however, require novel calculations which we present in Sec. 4. Unsurprisingly, the properties of TAP states shall play an essential role.

The semiclassical ramp in the ergodic phase
To reiterate, we are evaluating where H is the PSM Hamiltonian (Eq. (16)) and f is a filter function as discussed in Sec. 1.1. Here we consider the ergodic phase, for which the results are analogous to those of SYK [29]. We then consider the non-ergodic phase in Sec. 4.

Effective action
The calculation begins by retracing the steps described in Sec. 2.2, only on a modified contour. We still have upper and lower branches indicated by α ∈ {u, l} (with u = 0 and l = 1), but now each is separately periodic. Furthermore, we no longer have a thermal branch. See the right column of Fig. 4, as compared to the left column. While some care is required to account for the filter functions (as discussed in Appendix C), we ultimately arrive at an expression analogous to Eq. (27): where the "energy density" α [G] is defined as (0 + denotes a positive infinitesimal) See App. C for details. The saddle point of S eff is given by the equations (compare to Eq. (28)) Denoting averages with respect to the path integral of Eq. (44) by · , the expectation value of G is related to the original degrees of freedom as follows (we omit the filter functions here for brevity): where T denotes time ordering and T denotes time anti-ordering. One immediately sees from Eq. (48) that: i) all components of G αα (t, t ) are time-translation invariant and have period T ; ii) G uu (t, t ) and G ll (t, t ) are even functions of t − t ; iii) G ul (t, t ) and G lu (t, t ) are in fact independent of both time arguments; Solutions to Eq. (47) do not necessarily share all these properties, since some of the symmetries may be spontaneously broken. However, one simple solution that obeys all of the above is to take G ul (t, t ) = G lu (t, t ) = 0. The resulting action is precisely what one would get from averaging each factor of Tre −iHT separately, i.e., this solution gives the disconnected contribution to the SFF: where · · · denotes the contribution to the path integral from non-zero G ul and/or G lu . Eq. (49) holds equally well in the non-ergodic phase, and thus the remainder of this paper will be concerned with determining those additional contributions.

Connected solutions
Following Ref. [29], we construct approximate solutions to Eq. (47) which become accurate at large T . We first present the solutions and justify them afterwards. Take G αα (t, t ) to be the Schwinger-Keldysh correlation function at inverse temperature β aux , exactly as given in Sec. 2.2 (Eq. (29) in particular). Again define F αα (t, t ) ≡ J 2 G αα (t, t ) p−1 . A solution to the SFF saddle point equations (up to terms which vanish at large T ) is Here ∆ can be any real number between 0 and T . Thus Eqs. (50) and (51) constitute a two-parameter family of solutions, the parameters being β aux and ∆. Every such solution contributes to the SFF.
As for the Lagrange multipliers z α (t), they are independent of t due to time translation invariance. We further have that z u = z l ≡ z: both equal the value of the chemical potential needed to satisfy the equilibrium spherical constraint, i.e., N −1 i TrZ −1 SK e −βauxH σ 2 i = 1 (time translation invariance then implies that G αα (t, t) = 1 for all times and both branches).
To justify Eqs. (50) and (51), it is essential that G αα (t − t ) decay exponentially to zero as |t − t | → ∞. Thus these solutions only apply in the ergodic phase. With this in mind, the following comments together establish their validity: • The sum over n ensures that G αα (t − t ) has period T , even though G αα (t − t ) does not.
• Since G αα (t − t ) decays exponentially, G αα (0) ∼ 1 up to terms which are exponentially small in T .
• The equation F αα (t, t ) = J 2 G αα (t, t ) p−1 is satisfied up to exponentially small terms because, when raising Eq. (50) to the p − 1'th power, all cross terms are exponentially small (as is the sum over them).
In other words, • G αα (t, t ) obeys Eq. (28), written explicitly in terms of components as where v denotes the thermal branch of the contour. For t, t 1 (which still allows t − t to take any value), G αv (t + iτ ) is exponentially small for all τ and the last term on the left-hand side can be neglected. We can also take the lower limit of the t integral to −∞. Thus when checking whether Eq. (50) satisfies Eq. (47), we have that again making use of the fact that G αα (t − t ) is exponentially small when |t − t | is large. The equation is indeed satisfied.
• Finally, the off-diagonal components G ul (t, t ) and G lu (t, t ) contain the parameter ∆ because they break the separate time translation symmetries in t and t (see property iii above). Thus if any choice of ∆ solves Eq. (47), so do all choices of ∆ ∈ [0, T ).
As noted above, we have thus identified a two-parameter family of solutions to the SFF saddle point equations. In what follows it will be more convenient to parametrize the solutions by the equilibrium energy density (β aux ) corresponding to inverse temperature β aux . We can express (β) in terms of G (and thus G) by inserting a factor of H into the Schwinger-Keldysh contour. Since H clearly commutes with the evolution operator e −βH e iHt e −iHt , it can be inserted at any point, in particular at a late time for which (again because G αα (t − t ) decays exponentially) the thermal branch can be neglected. By following the same steps as in Appendix C, we find that (β) is given precisely by Eq. (46), evaluated on either branch:

Contribution of connected solutions
Having demonstrated that Eqs. (50) and (51) where ω ∈ 2πZ/T and G αα (ω) ≡ T 0 dte iωt G αα (t). Note that the Lagrange multiplier terms have dropped out since z u = z l . Furthermore, since dtG αα (t) p ∼ dtG αα (t) p , the general relation in Eq. (30) implies that the first term of Eq. (56) in fact vanishes.
For the second term, note that by Eq. (50), The exponential decay of G αα (t) implies that G αα (ω) (and thus G αα (ω)) is an infinitely differentiable function of ω. Strictly speaking, since the path integral is regularized by a timestep ∆t → 0, G αα (ω) is furthermore periodic with period 2π/∆t. The same is true of log Det G αα (ω). Thus the Euler-Maclaurin formula [69] gives up to terms which vanish faster than any polynomial in T −1 . Thus S eff is proportional to T , and we only need to evaluate the proportionality constant. Rather than calculate the integral directly, we follow Ref. [29] and evaluate the derivative dS eff /dT starting from Eq. (45). It is convenient to rescale time as t → T t, so that T becomes simply another parameter: Note that, since S eff is evaluated at a solution of the saddle point equations, we only need to differentiate the explicit factors of T : (60) Returning to unscaled time and using Eq. (47), we have that again using Eqs. (30) and (55). Thus the proportionality constant is in fact zero, i.e., S eff = 0.

Evaluation of the SFF
To finally compute the SFF, we simply need to sum over all connected solutions, i.e., integrate over aux and ∆. However, there are additional discrete symmetries which give further solutions: i) we can time-reverse the off-diagonal components, i.e., take G ul (t) = G ul (−t) and G lu (t) = G lu (−t); ii) if p is even, we can take G ul (t) = −G ul (t) and G lu (t) = −G lu (t). These must be summed over as well, giving an additional factor of 2(1 + δ p even ), where δ p even is the indicator function on p being even (1 if true, 0 if false). Thus our final expression is The measure 1/2π can be derived using hydrodynamic methods [19,29], but its precise value is not essential for our purposes. The key feature is simply that the linear-in-T ramp has emerged. However, keep in mind that Eq. (62) is only valid if the filter function is such that all contributing values of aux lie in the ergodic phase. In the following section we modify this analysis to hold in the non-ergodic phase as well. We shall see that it is necessary to incorporate the structure of multiple TAP states.

The semiclassical ramp in the non-ergodic phase
As we have stressed repeatedly, the results of Sec. 3 rely heavily on having an equilibrium correlation function which decays to zero at late times. Thus a new approach is needed to calculate the SFF in the non-ergodic phase, where G αα (t − t ) → q EA = 0 as |t − t | → ∞. More specifically, we can no longer neglect the integral over the thermal branch in Eq. (28), and G αα (t−t ) no longer solves the SFF equations of motion (Eq. (47)).
However, in the TAP equations of motion, Eq. (35), we can neglect the thermal branch since G(t) − q EA does decay to zero exponentially quickly. This suggests that a viable strategy is to construct solutions for the SFF using the TAP correlation function. Since TAP states are parametrized by the quantity E in Eq. (40), it will be necessary to first modify the SFF path integral so as to involve E. We associate the magnetizations and overlap from the TAP approach with time-averaged functions of the spin configuration, namely The choice to use only the upper contour in defining m i [σ] and q[σ] will become convenient in Sec. 5, but for now one could equally well use any other combination of branches, say the average of σ i (t) over the lower branch or over both branches symmetrically. With these definitions, we introduce E via Eq. (40).

Effective action
To begin, insert an additional fat unity into the path integral: With this addition, the full path integral is Proceeding as usual -averaging over disorder, introducing G αα (t, t ) and F αα (t, t ) as before, integrating out spins -we arrive at where we are denoting q[G] ≡ T −2 T 0 dtdt G uu (t, t ). The argument of the filter function is modified as well; it is now Note that E aux enters linearly into the action. Thus if we were to integrate over E aux at this point, we would obtain a δ-function forcing λ = 0. The action would then reduce to the ergodic-phase expression, Eq. (45). While reassuring, this would not have accomplished anything, so we instead treat E aux as a fixed parameter for now. We obtain saddle point equations by differentiating Eq. (67) only with respect to λ, z, G, and F .
The saddle point equations, assuming time translation invariance from the outset, are as well as the usual requirement G αα (0) = 1. Here G αα (ω) is the Fourier transform of G αα (t).

Connected solutions
With E aux fixed, let G αα (t) be the solution to the TAP equation of motion (Eq. (35)) corresponding to inverse temperature β aux . Denote the Edwards-Anderson order parameter at E aux and β aux by q EA . Also recall the various auxiliary quantities we defined in Sec. 2.3: the self-energy F αα (t) ≡ J 2 G αα (t) p−1 , the deviations ∆G αα (t) ≡ G αα (t) − q EA and ∆F αα (t) ≡ F αα (t) − J 2 q p−1 EA , and the quantity Λ ≡ ∞ 0 dt∆G I (t). We have that G αα (t) and F αα (t) obey Eq. (35), which by taking t and t to be far from the thermal branch can be written We also have, as a result of Eq. (39), the relationship Finally, recall the expression for the complexity Σ(E), the logarithm of the number of solutions to the TAP magnetization equations at E: where E 2 th = 4(p − 1)/p 2 . Using Eq. (73), we can express Σ(E aux ) in terms of Λ and q EA : Our solution to the SFF saddle point equations, Eqs. (69) through (71), is best written in the frequency domain (tildes denote Fourier transforms): We again take z α to be the equilibrium value corresponding to β aux . The precise form of the correction terms g αα (ω), f αα (ω), and δ is largely unimportant -the essential feature is simply that they are O(1) and the corrections are thus O(T −1 ). Note that in the time domain, this solution amounts to The sums are convergent because ∆G αα (t) and ∆F αα (t) decay rapidly to zero as |t| → ∞.
Although we have omitted it for notational simplicity, we can add a term δ α =α ∆ to the time arguments of ∆G αα (t + nT ) and ∆F αα (t + nT ) for any ∆ ∈ [0, T ), exactly as in Sec. 3. Due to the separate time translation symmetry on each branch of the SFF contour, all such solutions are equally valid and contribute the same action. Thus we shall demonstrate the validity of Eqs. (76) through (78) and evaluate the action only for ∆ = 0, but then integrate over all ∆ ∈ [0, T ) in the final expression for the SFF.
Let us first confirm that our solution satisfies the saddle point equation for λ, Eq. (71). Referring to Eq. (37), we have that ∆ G αα (0) = L + (−1) α 2iΛδ αα . Thus and Eq. (71) becomes Solving for λ indeed gives Eq. (78). The O(T −1 ) terms determine δ as a function of the other quantities. Now turn to Eq. (70). In the frequency domain, the right-hand side evaluates to 2 (83) Along the lines of Eq. (52), we have that which agrees at O(1) due to Eq. (72). At zero frequency, we instead have The O(T ) terms come out to be Yet from the TAP magnetization equations, Eq. (36), it follows that Thus Eq. (87) evaluates to It remains only to evaluate the action, Eq. (67), at the above solution. The action can be written as Interestingly, we can determine S eff up to a single additive constant simply by noting that dS eff /dE aux = −iλ (recall that S eff is stationary with respect to variations in all quantities other than E aux ). With λ given by Eq. (78) and q p/2−1 EA Λ given by Eq. (41), we can carry out the integral to obtain that for some unknown constant C. Comparing to Eq. (42), this is highly suggestive that S eff = −Σ(E aux ). Of course, we do need to determine the remaining constant, and so we now turn to a more elaborate calculation. Rather than substitute Eqs. (76) and (77) into Eq. (90), we instead use the simpler functions and show that the error incurred in doing so vanishes at large T . Let us demonstrate that the error is negligible first. At any non-zero frequency, we have that The partial derivatives of S eff at nonzero ω are which vanish when evaluated at G αα (ω) = ∆ G αα (ω) and F αα (ω) = ∆ F αα (ω). Thus the O(T −1 ) difference between G αα (ω) and G αα (ω), as with F αα (ω) and F αα (ω), translates only to an O(T −2 ) difference in the action. Even after summing over all ω = 0, the total error 3 is only O(T −1 ). Neglecting non-zero frequencies, F αα (ω) is identical to F αα (ω) and G αα (ω) differs only by g αα (0)δ ω0 /T . In the time domain, the latter corresponds to Yet When evaluated at G αα (t) and F αα (t), the O(T ) contribution vanishes (see Eq. (84)). Thus ∂S eff /∂G αα (t) is O(1), and an O(T −2 ) change to G αα (t) leads only to an O(T −1 ) change in the action even after integrating over t.
Since all errors are O(T −1 ), we can safely evaluate S eff at Eqs. (92) and (93) rather than the full solution (we still use Eq. (78) for λ). The first line of Eq. (90) can be computed straightforwardly. It comes out to be p E aux − 2Jq where we used Eq. (73) to obtain the right-hand side. Next consider the bottom line. Since F αα (ω) = ∆ F αα (ω) for ω = 0, while F αα (0) = F αα (0), we can write the determinant term as The top line vanishes by exactly the same reasoning as in Sec. 3.3: it is proportional to T by the Euler-Maclaurin formula, and then must be zero since the derivative with respect to T vanishes. Given Eq. (37), the bottom line is simply For the middle line we take an indirect approach. We have that iz(−1) α δ αα + (−1) α+α F αα (0) is the matrix inverse to G αα (0) (using the full solution for the latter, Eq. (76)). Written out, Rather than this (u, l) basis, express Eq. (102) in the (u + l, u − l) basis (called "classical"/"quantum" in the Keldysh literature), denoted (+, −): We can read off that Det G(0) −1 = F ++ (0)/ G ++ (0). Note that we only need G ++ (0) and F ++ (0) to O(T ) in order to calculate the determinant to O(1). Thus the middle line of Eq. (100) evaluates to (log J 2 q p−2 EA )/2, and the total contribution of the determinant term is Lastly consider the middle line of Eq. (90). Since G αα (t) = G αα (t) (up to exponentially small corrections), αα (−1) α+α G αα (t) p = 0 by virtue of Eq. (30). We are left with The first term is again proportional to αα (−1) α+α G αα (t) p = 0. The second term would appear to be more problematic, since f αα (0) (for which we have not given an explicit expression) contributes at O(1) due to G αα (0) being O(T ). However, we only need the component f −− (0)/T = F −− (0), and from Eq. (103) we see that Eq. (105) evaluates to again using Eq. (73). We finally have the large-T limit of the action, given by the sum of Eqs. (99), (104), and (107): Comparing to the complexity Σ(E) given in Eq. (75), we see that S eff is precisely −Σ(E aux ).

Evaluation of the SFF
We have shown that, at a given E aux and for each value of inverse temperature β aux , there is a solution to the SFF saddle point equations with S eff = −Σ(E aux ). The full (connected) SFF is obtained by integrating over all E aux and β aux , as well as the symmetry-broken order parameter ∆ (which contributes an overall factor of T ) and an additional factor 2(1 + δ p even ) from the discrete symmetries. As in Sec. 3, it is more convenient to integrate over the energy density (E aux , β aux ). We show in App. B that (E, β) comes out to be precisely the argument of the filter function, Eq. (68), when evaluated at the saddle point solution.
Our final result is that 4 SFF(T, f ) = ETrf (H)e −iHT 2 + 2 1 + δ p even T pN 2π dE aux e N Σ(Eaux) where the inner integral runs only over the range [ − (E aux ), + (E aux )] in which solutions to the TAP equations exist. Furthermore, one can easily generalize Eq. (109) by making the filter function E-dependent, i.e., f (E, aux ). The resulting quantity is the SFF for the projection of the system into certain TAP states. Compare Eq. (109) for the non-ergodic phase to Eq. (62) for the ergodic phase, and recall the discussion of block-diagonal Hamiltonians in Sec. 1.1. Our result demonstrates that each metastable (i.e., TAP) state can be thought of its own quantum chaotic subspace, one which is independent of any others. This is the central result of our paper. While the qualitative idea has been proposed in previous work [58], the present analysis both makes it precise and proves it.

Higher moments of the evolution operator
In this final section we consider higher moments of Tre −iHT , i.e., the quantities The saddle points of these higher moments exhibit an interesting structure that will shed further light on the distribution of TAP states, although care must be taken in interpreting the results. We first present the calculation and discuss afterwards.  which we evaluated in Sec. 4. In other words, this contribution to the n'th moment is simply SFF(T, f ) n . However, by the permutation symmetry described above, we actually have n! such contributions: where the ellipses denote additional solutions.

Connected solutions
In general, for arbitrary values of E aux,a , we have been unable to find any further saddle points. However, when some replicas have equal values of E aux , we can construct additional solutions. Pick any set of inverse temperatures β a (not necessarily equal), and suppose that the replicas {1, · · · , n} partition into groups A ≡ {a 1 , · · · , a |A| }, such that E aux,a equals a common value E aux,A for all a ∈ A. We again take G aα,aα (t−t ) to be the solution from Sec. 4. For a and a in different groups, we still set G aα,a α = 0. For a and a in the same group A, however, we now set where q EA,a is the Edwards-Anderson order parameter corresponding to E aux,A and β a . This corresponds to the replicas lying within the same TAP state (see Fig. 5). We can write this compactly as G aα,a α (t) = q EA,a q EA,a 1/2 + δ aa ∆G a,αα (t) + O(T −1 ) . (122) It remains only to check that Eq. (114) can be satisfied. It is automatically solved at non-zero frequencies, since then G(ω) and F (ω) reduce to δ aa ∆ G(ω) and δ aa ∆ F(ω) respectively. At zero frequency we confirm that the equation is solved to O(T ) (the O(1) terms only determine subleading corrections). Following the same steps as in Sec. 4.2, the left-hand side of Eq. (114) simplifies to JT q p−1 EA,a q EA,a a ∈A λ a + 2i(p − 1)Jq as desired.
Note that in this solution, only the sum a λ a is determined -all orthogonal components of the vector λ are free to take any values. This does not imply that there are multiple such solutions, however. Returning to the effective action in Eq. (112), the fact that the saddle point equations determine only G, F , and a λ a means that, if we first integrate over them, the resulting λ-dependent action is of the form for some function S of the single quantity a λ a (as well as all E aux ) and for any choice of orthonormal basis vectors u ab orthogonal to the all-1 vector. When we integrate over a λ a u ab , we thus get a δ-function forcing a u a b E aux,a = 0. Together, the δ-functions force all E aux,a to equal a common value E aux,A . Not only is this consistent with our original assumption, it shows that our construction cannot work for any other values of E aux,a .

Contribution of connected solutions
To evaluate the action, note first of all that since the numbers β a define a continuous family of solutions, and since the action is by definition stationary at these solutions, all choices of β a must give the same value of the action. We thus take all β a to equal a common value β for simplicity. The action evaluated at this solution still decomposes into a sum over groups, but now the contribution of a single group A is Note that now E aux , q EA , and Λ are all independent of the replica a (within a given group A). We are also free to set all λ a = λ, meaning that our saddle point solution simplifies to (in frequency space) F aα,a α (ω) = T J 2 q p−1 EA + ipJq p/2−1 EA and we again must determine certain components of F (0) and Det G(0). As before, it is expedient to use the (+, −) basis with respect to contour indices. We also switch to the Fourier basis with respect to factor indices: from Eq. (126), We see that Det G 0 (0) = G 0+,0+ (0)/ F 0+,0+ (0) ∼ 1/J 2 q p−2 EA , and F 0−,0− (0) (which is in fact the only element of F (0) needed in Eq. (129)) is given by which is again precisely −Σ(E aux ).

Evaluation of the Higher SFF
In the above calculation, note that we get a single contribution of complexity for the entire group A. However, there is still a factor (2T ) |A| (1 + δ p even ) 2|A|−1 due to the separate time translation, time reversal, and reflection symmetries of each replica 5 . Finally, the sum over all connected solutions amounts to a sum over the possible ways of partitioning n elements, in addition to the n! ways of pairing upper and lower contours. Using P ≡ {A 1 , · · · , A |P | } to denote a partition, we have that In particular, suppose the filter function is chosen so as to have a small width ∆E 1/N around a certain value E (as in Sec. 4.4, the above calculation can easily be modified to allow for E-dependent filter functions). Then the n'th moment simplifies to Eq. (134) has a nice interpretation as the n'th moment of a sum over a Poisson-distributed number of Gaussians. To be precise, suppose we have an infinite sequence of i.i.d. complex Gaussians, where M is itself a Poisson-distributed random variable with mean µ. The n'th moment of SS * , averaging over both Gaussians and M , can be written where p µ (m) denotes the Poisson distribution of mean µ, and Wick's theorem is used for the latter equality. It is known that the n'th moment of a Poisson distribution is P µ |P | , where the sum is again over all partitions of n elements. Thus If we associate σ 2 with the SFF of a single TAP state at E, and associate µ with the number of TAP states, then Eqs. (134) and (136) are identical. It is quite tempting to interpret this as saying that the number of TAP states at E is Poisson-distributed with mean given by Eq. (138), and that each TAP state has a Gaussian-distributed value of Tre −iHT with variance (i.e., SFF) given by Eq. (137). We are not aware of any results in the literature which would contradict such a claim. However, keep in mind that SFF (n) has perturbative corrections around each saddle point which are suppressed by powers of N , whereas every connected partition in Eq. (134) is suppressed exponentially relative to the fully disconnected one, whose contribution is given by Eq. (118) 6 . Thus we cannot claim to have rigorously computed the n'th moment to any level of accuracy beyond the disconnected piece. Nonetheless, the structure of saddle points which we have identified is highly suggestive and warrants further investigation.

Concluding Remarks
The focus of this paper was the derivation of Eq. (109), which gives the connected SFF of the quantum PSM in the non-ergodic phase. Our result demonstrates that each metastable state (i.e., TAP state) can be considered as an independent chaotic phase with an independent RMT-like Hamiltonian, at least as far as level statistics are concerned.
It is interesting to compare our result for the PSM with Ref. [29], which performs a similar calculation for the SYK model (the latter is structurally quite similar but with fermionic degrees of freedom). In the ergodic phase, we find an essentially identical result using extremely similar methods. Below the dynamical transition, however, the PSM displays an enhanced ramp quite different from that of the SYK model, which has no analogous phase.
Since we only calculate the SFF up to times polynomial in system size, our results are consistent with but do not test the distinction between localized and delocalized phases shown in Fig. 3, which is only relevant beyond the exponentially long timescale corresponding to tunneling between TAP states. We leave it for future work to incorporate such instanton effects into the path integral, expecting that they will reduce the SFF to the random matrix result precisely in the non-ergodic delocalized phase (and even then only beyond the exponential tunneling timescale). At the same time, consideration of exponential scales can allow one to identify a plethora of additional dynamical phases (see Ref. [47] for an example), and so the structure of instantons in these spin glass models may be quite rich.
In addition to the SFF, we have also considered higher moments of the evolution operator. We identified an important family of saddle points and evaluated its contribution to these higher SFFs (Eq. (134)). The results suggest that: i) the number of TAP states at a given energy is Poisson-distributed, and ii) the numbers of TAP states at different energies are independent. However, since we have not evaluated the perturbative corrections around each saddle point, which at finite complexity would generically dominate over any subleading saddle points, we cannot claim to have an accurate calculation. It is another direction for future work to study the distribution of TAP states more systematically.
Our results can be further understood by comparing to Refs. [18] and [19] on one hand and Refs. [70] and [71] on the other. The first set of papers argues that for a system which separates into weakly coupled sectors, the SFF enhancement is the sum of return probabilities over all configurations. If the time evolution can be considered as an effective Markov process with transfer rates between sectors given by some matrix M , then the SFF enhancement factor is Tre M T . The second set of papers argues that for a classical spin glass undergoing Markovian stochastic dynamics with generator M , the number of TAP states can be calculated -and perhaps even defined -as Tre M T . In this sense, the present paper can be considered as a "missing link" that extends the results of Refs. [70] and [71] to quantum systems.
The fact that SFF enhancement is related to return probabilities suggests that the spectral statistics of spin glasses may contain information on aging dynamics as well. Another open question is whether the equilibrium replica-symmetry-breaking transition has any consequences for spectral statistics. These, as well as those already mentioned, are all promising directions for future work.
Finally, let us briefly comment on the case p = 2, which -being a Gaussian model but for the global spherical constraint (and still integrable in any case [72]) -exhibits very different behavior than the p > 2 models considered here. The spectral statistics of the analogous SYK model with two-point interactions among fermions have been studied in the mean-field limit [31,73], and found to have ramps growing exponentially in time. One explanation for these ramps is the spontaneous breaking of a hidden SU (2) k symmetry in the saddle point equations, since for p = 2 the matrix G has a separate SU (2) conjugation symmetry at each frequency. Inspections show that the p = 2 spherical model has the same symmetry at tree level, although it is broken by higher loop effects. A full analysis of this system would be yet another excellent topic for further research.

A Derivation of Schwinger-Keldysh TAP equations
Here we derive the TAP equations on the Schwinger-Keldysh contour, given by Eqs. (35) and (36) of the main text. Our derivation is a straightforward generalization of that in Ref. [67] for the thermal circle, but since we are not aware of it appearing in the literature, we aim to make this section as self-contained as possible.
We begin with the expression for the TAP action, given by Eq. (32) and reproduced here: where C denotes the Schwinger-Keldysh contour and Note that we have included an additional parameter η governing the strength of interactions. We calculate S TAP by expanding in η, ultimately setting η = 1. One can show that all terms beyond second order vanish in the thermodynamic limit [67,74], so we need only calculate S TAP and its first two derivatives at η = 0. Then We proceed term by term.

Zeroth order
At zeroth order in η, we have that The first term evaluates to We further have that Note that the left-hand sides are constrained to be G(t, t ) − Q(t, t ) and m i (t) (summing over i for the former), where Q(t, t ) ≡ N −1 i m i (t)m i (t ). Thus we can use Eqs. (144) and (145) to express the action in terms of them. The first term of Eq. (140) becomes and the second line can be written t ). (147) Combining the two, we have that (up to an unimportant constant)

First order
Strictly speaking, since h i (t) and Λ(t, t ) are chosen so that σ i (t) = m i (t) and i σ i (t)σ i (t ) = N G(t, t ) at every value of η, they themselves are functions of η (as is z(t)). However, the Legendre-transform structure of the TAP action ensures that the total derivative with respect to η equals the partial derivative holding fields fixed (for example, the contribution to ∂S TAP /∂η from ∂h i (t)/∂η comes with a factor σ i (t) − m i (t) = 0). Thus we have simply that At η = 0, the spins are non-interacting and the expectation value factors:

Second order
The second derivative requires a bit more work. From Eq. (149), we have that where for brevity, we use I to denote (i 1 · · · i p ) and σ I (t) to denote σ i1 (t) · · · σ ip (t). To evaluate the partial derivatives on the right-hand side, which we only need at η = 0, note that where we set i σ i (t) 2 = N ρ(t) (even though ultimately ρ(t) = 1). Thus, using Eq. (150), The second derivative simplifies to Just as one can show that higher-order η derivatives vanish in the thermodynamic limit, one can also show that it is safe to replace J I J I by its average in Eq. (154) [67,74]. Again using that expectation values factor at η = 0, we therefore have that Full TAP action Inserting Eqs. (148), (150), and (155) into Eq. (141), we find that As discussed in the main text, the TAP equations come from setting to zero the derivatives with respect to m i (t) and G(t, t ). However, since G(t, t) is not free to vary (it must equal 1 by the spherical constraint), we include a Lagrange multiplier that we again denote by z(t). The TAP equations are thus Time translation invariance allows us to further simplify by setting m i (t) = m i and z(t) = z. Note that then Q(t, t ) = N −1 i m 2 i ≡ q EA . After some rearrangement, the TAP equations can be written These are precisely Eqs. (35) and (36) from the main text.

B Energy of a TAP state
Here we derive an expression for the total energy density in terms of a solution G(t, t ) to the TAP equations, Eqs. (159) and (160). Note that, despite the complicated derivation, G(t, t ) is simply the contour-ordered expectation value of N −1 i σ i (t)σ i (t ). Thus, as we discuss below in App. C, the kinetic energy per spin is given by µ 2N ∆t 2 i σ i (t + 2∆t) − σ i (t + ∆t) σ i (t + ∆t) − σ i (t) = − µ 2 ∂ 2 t G(t + , t), where t + denotes t + ∆t. To evaluate the potential energy, we replace J i1···ip by (1 + ζ(t))J i1···ip in the original path integral. On the one hand, ∂ ∂ζ(t) S TAP On the other hand, it is easy to see how the explicit calculation of S TAP is modified by ζ(t): (163) Thus and the total energy per spin is given by we need to be careful in how we treat the filter functions, since f (H) does depend on the specific disorder