Absence of operator growth for average equal-time observables in charge-conserved sectors of the Sachdev-Ye-Kitaev model

Quantum scrambling plays an important role in understanding thermalization in closed quantum systems. By this effect, quantum information spreads throughout the system and becomes hidden in the form of non-local correlations. Alternatively, it can be described in terms of the increase in complexity and spatial support of operators in the Heisenberg picture, a phenomenon known as operator growth. In this work, we study the disordered fully-connected Sachdev-Ye-Kitaev (SYK) model, and we demonstrate that scrambling is absent for disorder-averaged expectation values of observables. In detail, we adopt a formalism typical of open quantum systems to show that, on average and within charge-conserved sectors, operators evolve in a relatively simple way which is governed by their operator size. This feature only affects single-time correlation functions, and in particular it does not hold for out-of-time-order correlators, which are well-known to show scrambling behavior. Making use of these findings, we develop a cumulant expansion approach to approximate the evolution of equal-time observables. We employ this scheme to obtain analytic results that apply to arbitrary system size, and we benchmark its effectiveness by exact numerics. Our findings shed light on the structure of the dynamics of observables in the SYK model, and provide an approximate numerical description that overcomes the limitation to small systems of standard methods.


Introduction
The investigation of out-of-equilibrium quantum systems is an extremely rich and active topic in modern statistical and condensed matter physics, which aims at clarifying fundamental concepts like thermalization and the emergence of the ensemble description [1][2][3][4][5][6]. Theoretical advances in recent years have highlighted the importance of quantum information scrambling [7][8][9] as a general feature associated to the relaxation of observables: Locally encoded information spreads throughout the system during the dynamics, and, while formally it is not lost, it cannot be recovered through any local measurement at late times. A paradigmatic example of a system exhibiting scrambling behavior is the Sachdev-Ye-Kitaev (SYK) model [10][11][12][13], which received much interest in recent research due to its relevance in multiple branches of physics, ranging from the study of black holes [10,[14][15][16][17] to non-Fermi liquids [18][19][20][21]. There are multiple versions of the SYK model, parameterized by an even number q and referred to as SYK q . Each one consists of fermionic modes coupled via all-to-all disordered q/2-body interaction amplitudes. For q ≥ 4, the model manifests quantum chaotic behavior [10,22,23], diagnosed by the presence of a quantum Lyapunov exponent in out-of-time-order correlators (OTOCs) [24][25][26]. Specifically, this exponent is found to saturate its theoretical upper bound [24], implying that quantum scrambling occurs as rapidly as possible. While multiple properties can be derived exactly in the thermodynamic limit [10,13,18], studying the model at finite size remains a challenging yet essential problem, especially in view of future experimental investigations that have recently been proposed [27][28][29][30][31][32][33].
In this work, we investigate the disorder-averaged dynamics of observables in the SYK model, highlighting the manifestation of a counter-intuitive symmetry. In contrast to the scrambling nature of the system, we prove rigorously that disorder-averaged operators do not increase in complexity throughout the time-evolution. This corroborates the results of our previous work [34], which demonstrated numerically a non-trivial universality with respect to the choice of the initial state in the dynamics of some observables. The argument only applies to single-time observables, and, in particular, it does not hold for OTOCs: All correlators that can be used to diagnose scrambling and chaotic behavior are unaffected by our discussion, which resolves the apparent contradiction of our analysis with the wellknown results from the literature. While our findings do not apply a priori to individual realizations of the model, they are expected to manifest in self-averaging single-time physical quantities for sufficiently large system sizes. In addition, working in a framework typical of open quantum systems, we leverage these results to show how to approximate the average dynamics by means of a cumulant expansion of the effective dynamical map. Specifically, the functions approximating the dynamics of observables can be computed analytically for arbitrary system sizes. Finally, we benchmark the performance of our approximation scheme through comparison with results found using exact diagonalization.
The paper is organized as follows. Section 2 presents the model and the superoperator framework adopted throughout the work. In Section 3, we introduce the concept of operator size, and we prove its conservation throughout the dynamics in ensemble average. Moreover, we discuss how this result does not conflict with well-known scrambling properties of the model. Then, in Section 4 we present the cumulant expansion technique that we implement to approximate the dynamics. Our theoretical findings are tested in Section 5, in which we investigate operator growth, and we probe the performance of the cumulant expansion using exact diagonalization simulations for the SYK 4 model. We summarize our findings in Section 6, and we discuss their possible applications, as well as future lines of research. Appendices expand the discussion of the main work, both by providing detailed descriptions of some proofs, and by showing explicit analytic results.

Model and focus of our study
We focus on the family of SYK models of N complex fermion modes on a fully connected lattice, defined by the Hamiltonian where q is an even integer number, and J i 1 ,...,i q/2 ;j 1 ,...,j q/2 are complex Gaussian independent random couplings with the following statistical properties (E [...] denotes disorder averaging): Notably, for q > 2 the only integral of motion of any individual realization is the total chargeQ = in i . Our analysis relies on the study of the dynamics of observable expectation values, namely where typically (but not necessarily)ρ 0 = |ψ 0 ψ 0 | is a pure-state density matrix: This corresponds to preparing the system in a disorder-independent initial state |ψ 0 , performing a quantum quench to the SYK q Hamiltonian, and tracking the time-evolution of physical quantities. Specifically, since the Hamiltonian involves randomness, we focus on the disorder-averaged evolution E Ŵ (t) . A well-known feature of disordered systems is that disorder-averaging introduces mixing, making the average dynamics non-unitary [35]. As a consequence, it becomes natural to adopt a superoperator-based framework, typical of dissipative open quantum systems [36,37]. We introduce the Liouvillian L• = −i Ĥ , • , as well as the superoperators α • = −i ĥ α , • , where α = {i 1 , . . . , i q/2 ; j 1 , . . . , j q/2 } is a multiindex, andĥ α =ĉ † i 1 . . .ĉ † i q/2ĉ j 1 . . .ĉ j q/2 appears in the Hamiltonian of Eq. (2.1a). With this notation, the ensemble-averaged dynamics reads In what follows, we provide a detailed analysis of the time evolution of such disorderaveraged expectation values.

Operator size symmetry in the disorder-averaged ensemble
In this section, we introduce the notion of operator size, and we show that it plays the role of a conserved quantum number for the average dynamics of operators. In this sense, when the total charge is fixed, no operator growth takes place for expectation values of observables. As we discuss, this result does not hold for more sophisticated correlators, e.g., OTOCs, which are able to diagnose quantum scrambling.

Operator size symmetry
An analysis of the spectral properties of the disorder-averaged time-evolution superoperator allows us to determine some exact features of the dynamics. Specifically, we will prove further below that the average dynamics distinguishes different operator sizes. Loosely speaking, we say that an operator consisting of the product of m creation and n annihilation operators has operator size (m, n) [38,39]. This definition is appropriate if all fermionic operators refer to different lattice sites, but we must analyze more carefully the cases involving one or more number operatorsn i . In any charge-conserved sector of the Hilbert space, the operatorQ = in i is proportional to the identity, and thus its size is (0, 0). 1 As a consequence, the size of the number operatorsn i is not simply (1, 1), because each of them has a finite overlap with the identity operator. Despite this complication, it is possible to define the concept of operator size for the most general combination of creation and annihilation operators (see the discussion of Appendix A). Any linear combination of operators sharing a single common size also maintains that same size. This allows us to define an orthonormal basis of operators with well-defined sizes, over which any operator can be uniquely decomposed. We now show that the dynamics preserves the sizes of operators on average, in the sense that it does not introduce components with operator sizes different from those present at t = 0. We briefly present the main reasoning here, whilst a more detailed description of the following procedure can be found in Appendix B.
Taking the disorder average of a simple annihilation operatorĉ i (t) in the Heisenberg picture gives where the superoperators α i are defined at the end of Section 2, and the sum over repeated multi-indices is implicit. After using Wick's theorem, the previous equation results in a sum of nested commutators involving the operatorsĥ α 1 , . . . ,ĥ αn , as well as their Hermitian conjugatesĥ † α 1 , . . . ,ĥ † αn , appearing in a certain order. In particular, each multi-index α k is shared precisely byĥ α k andĥ † α k , as follows from Eqs. (2.2b) and (2.2c). If we expand each nested commutator in sequences of fermionic operators, for eachĉ † k appearing in any string there will always be exactly oneĉ k , with the same index, somewhere within that same string. For each chain of fermionic operators, we may use the anticommutation relations to moveĉ i towards one of the edges, for example the left one. Each new sub-string originated by this process involves only the operatorĉ i and multiple sums over all sites of paired creation and annihilation operators. Following this argument, any initial string is equal to the original operatorĉ i multiplied by a complex combination of fermionic operators that, however, cannot depend on any lattice site in particular, as all indices involved are summed over: As a consequence, the latter will be a function of the number of sites N and the charge operatorQ. In conclusion, we have proved that where the dependence on a function of t 2 arises from the absence of odd terms in Eq. (3.1). The above result shows that, if we limit its action to a subset of the Hilbert space with fixed Q,ĉ i is an eigenoperator of the generator of the ensemble-averaged dynamics. In other words,ĉ i (t) manifests no operator growth (nor hopping) on average when the total charge is fixed. Our above argument is readily extended to all fermionic strings of typeŜ m,n = c † i 1 . . .ĉ † imĉ j 1 . . .ĉ jn in which all indices are distinct, i.e., E e −Lt Ŝ m,n =Ŝ m,nfm,n (t 2 ; N,Q).
These results apply also to strings containing repeated indices, i.e., number operators, and thus they hold for any generic operator with fixed size (m, n): The proof is presented in Appendix B. Notice that the functionsf m,n are independent of the specific lattice sites appearing inŜ m,n , consistently with the S N symmetry of the SYK model under ensemble average. In summary, we know how all operators with a well-defined size evolve on average: The operator size is conserved, and it determines the evolution in time. We point out that the derivation of this conclusion is model-independent, up to the requirement that the couplings in the Hamiltonian bear no dependence on lattice indices. For this reason, not only does this proof apply to any SYK q model, but also to linear combinations of SYK Hamiltonians with different values of q, as well as to other fully-connected homogeneous models. Our results apply also to the Majorana version of the SYK model [11], for which no charge operatorQ can be defined; in that case, the functions f m,n will depend only on t 2 and N . 2 In addition, we expect the result to extend also to fermionic and spin SYK-like models with non-Gaussian disorder, as well as bosonic versions of the system with Gaussian randomness. 3

Implications for correlation functions and OTOCs
The previous results might misleadingly be interpreted as total absence of operator growth in the SYK model, which would be an incorrect conclusion. We stress that operator size symmetry holds only under disorder average, whereas individual realizations of the system do clearly present operator growth (see Section 5.1 for a detailed discussion). Indeed, operator growth is related to scrambling of quantum information, which is one key characteristic of the SYK model.
While it is important to distinguish physical features computed in single realizations and on average, there are cases in which the two coincide, namely for self-averaging observables. At suitably large N , such quantities approach their ensemble-averaged values for any single instance of the model, and thus operator size symmetry may manifest also for individual disorder realizations. Typical self-averaging quantities are expectation values of operators with extensive support [40,41], which are precisely in the form of Eq. (2.4). As a consequence, if, for instance,Ŵ =Ŝ m,m has a fixed operator size (m, m), andρ(0) has fixed number of particles Q, we have for large N . This implies that the functional form of the time-evolution is independent of the choice of the initial state, and the latter sets only the amplitude of the curve, as was observed in our previous work [34]. As anticipated previously, operator size conservation on average does not imply absence of scrambling in single physical realizations. An important aspect of the previous discussion is the limitation of our framework to the study of quantities that are linear in the timeevolution superoperator. In particular, OTOCs of the form which diagnose quantum chaotic behavior, are beyond the scope of our investigation, because they involve correlations between two superoperators e −Lt , as well as the thermal density matrixρ β = e −βĤ /Z. We conclude that quantum information scrambling must be encoded in these correlations, which arise because operators themselves are not selfaveraging, namelyŴ (t) = E Ŵ (t) even if their quantum-mechanical expectation values become equal at large N . The same argument applies to the so-called Krylov complexity [42][43][44], which is used to quantify operator growth. This is consistent with the belief that simple equal-time correlators are unable to capture the chaotic properties of the model, e.g., the quantum Lyapunov exponent. Such features are manifested only by more complicated quantities, such as, for instance, OTOCs.
To conclude this section, we point out that multiple studies in the literature focus on operator growth from the perspective of a distribution for the operator size [39,[45][46][47]. In such studies, operators are expanded over a basis of operators with fixed sizes, and a size distribution is obtained from the squared moduli of the expansion coefficients. Eventually, this framework allows one to compute the average size (see also [44,48]) and its dynamics, which can be directly related to OTOCs. Nevertheless, this bears no contradiction with our findings. The distribution has quadratic dependence on the time-dependent expansion coefficients, and thus its disorder-average probes correlations between two time-evolution superoperators; as discussed previously, our results do not apply to such a scenario.

Cumulant expansion method
In the previous section, we discussed the structure of the average time-evolution of operators. However, there is no straightforward way to explicitly determine the functionŝ f m,n (t 2 ; N,Q): The model is characterized by a single energy scale, and thus we cannot set up conventional perturbation theory starting from a known non-interacting solution. Instead, we now develop a cumulant expansion scheme to approximate the disorder-averaged dynamics of equal-time observables.
The numerical investigation presented in Ref. [34] shows that the quench dynamics of an operator can be well approximated by a Gaussian, even though the agreement is not perfect. Motivated by this result, we look for an exponential representation of the disorder-averaged time-evolution superoperator, in such a way as to recover the Gaussian shape as the lowest order approximation, and to implement additional corrections. We formally write where C(t) will be referred to as a cumulant generating superoperator, in agreement with usual nomenclature in statistics. Accordingly, the superoperators C k will be referred to as cumulants. Since we are considering Gaussian disorder, if L was just a scalar function instead of a superoperator, then all cumulants with k > 2 would vanish. In the present case, however, the infinite series on the right-hand side of the previous equation does not terminate at finite order, because the operators that multiply different disordered couplings do not commute.
A convenient way to determine the cumulants is by expanding the right-hand side of Eq. (4.1), and comparing it to Eq. (3.1) by matching equal powers of t. We immediately see that all superoperators C k with odd k vanish. The identification of C k requires the previous determination of all C l with l < k, and thus the procedure is iterative. Here, we present the first two non-vanishing terms: Leading order truncation of the cumulant expansion yields precisely the time-evolution superoperator used in Ref. [34], which was derived alternatively by writing an effective master equation for the ensemble-averaged density matrix. Higher-order cumulants provide corrections to the results of this previous work, allowing to achieve a better approximation of the dynamics. For the SYK model, all cumulant superoperators commute. This property is proved by observing that each superoperator of type α 1 . . . α 2k , in which the indices are contracted in pairs and summed over, preserves operators with well-defined size (see Appendix B for the details), and thus all cumulants share the same eigenoperators. In particular, each individual cumulant fulfills an eigenvalue equation analogous to Eq. (3.3). As a consequence, we can differentiate Eq. (4.1) to obtain an effective master equation for the disorder-averaged density matrix in the Schrödinger picture, 4 yielding The sum within the round brackets is the effective Liouvillian that generates the disorderaveraged dynamics.
For practical purposes, we can only evaluate the cumulant expansion up to some given finite order. Even if all C k were known, finding an analytic expression for C(t) remains as hard as computing the exact disorder-averaged time-evolution superoperator; as a consequence, the series must be truncated. We point out, however, that even if only a finite number of terms is known, the resulting approximate description of the dynamics can be quite accurate over a wide time window. Indeed, using a truncated cumulant expansion scheme to characterize the average dynamics already proved successful in systems with non-Markovian noise [49]. Evaluating the series up to any given order is only as expensive as computing a short-time expansion up to that same order. Nevertheless, while the latter is limited to early times only, the former is potentially able to reasonably reproduce the evolution at arbitrary times. The reason lies in the exponential form of the time-evolution superoperator: Even when C(t) is truncated, Eq. (4.1) still involves an infinite series of powers of t, and thus it can represent non-polynomial time-dependence. In addition, we can argue that the impact of higher-order cumulants is relevant only on longer timescales, as they are suppressed by a factor (2k)!. 5 Appendix C presents a numerical check of this statement, which indeed supports its validity. Suppose that the disorder-averaged density matrix quickly approaches a steady state, so that In this case, neglecting high-order cumulants that become relevant only at times greater than τ is a good approximation, because, as seen from Eq. (4.3), their action on the steady state density matrix is practically zero (notice that all cumulants share the same steady states, because they commute). In conclusion, for a system that thermalizes quickly, a finite number of cumulants is sufficient to obtain a valid approximation of the exact dynamics. As shown in Ref. [34], the quench dynamics of observables in the SYK model fall within this situation, as they manifest super-exponential relaxation to stationary values.
The above cumulant expansion allows us to approximate the functionsf m,n (t 2 ; N,Q). , which is smaller than (2k −1)!! (it being the total number of contractions of order 2k). By irreducible, we refer to those contractions that cannot be written as products of lower-order ones. For example, referring to Eq. (4.2b), the superoperators appearing with the plus sign are irreducible, whereas the other one is reducible. In conclusion, not only 2M k < 2(2k − 1)!! < (2k)!, but also different superoperators sum up in a destructive way. As argued before, a finite number of functionsλ (2k) m,n can be sufficient to obtain a very good approximation of the exact result. The advantage of this approach is that eigenoperators of cumulants are known, and thus we can evaluate their eigenvalues analytically. Some exact expressions for the SYK 4 model are presented in Appendix C. From a practical point of view, this greatly simplifies the numerical application of our formalism: We can approximate Eq. (4.5) for arbitrarily large system sizes without resorting to exact diagonalization of cumulant superoperators.

Numerical results
This Section presents numerical results that complement the previous theoretical findings. For a given operator, we first investigate the dynamics of those components which have a well-defined size, and then we compare the exact time-evolution to the approximate one produced by the cumulant expansion. The exact dynamics is obtained by exact diagonalization of the SYK 4 Hamiltonian within a sector of fixed total charge Q. We evolve an initial disorder-independent pure state |ψ 0 , and then we evaluate the expectation value Ŵ (t) of an observableŴ that commutes with the chargeQ. In our numerics, we choose |ψ 0 to be the Néel state |1010 . . . , where the ordering of lattice sites is arbitrary. Regarding the choice of the operator, we study the staggered magnetizationR = N i=1 (−1) i+1n i (the indices are assigned with the same convention used for the Néel state). Disorder-averaged quantities are computed by iterating this procedure multiple times for independent realizations of the disordered Hamiltonian given in Eq. (2.1a), and finally taking a statistical average over the sample.

Operator size throughout the dynamics
We now proceed to study the dynamics of individual components with well-defined size of a given operator. As mentioned previously, it is possible to define an orthonormal (with respect to the Hilbert-Schmidt product) basis of operators such that each basis element has a well-defined operator size. Using the notation introduced in Appendix A, we denote these basis operators asT (k) m k ,n k , and each of them has size (m k , n k ). We can extract their amplitudes in the expansion ofŴ (t) by taking the Hilbert-Schmidt product where Tr Q denotes the trace over states with fixed charge Q, as defined in Appendix A. For W =R, we use exact diagonalization to study the dynamics of some coefficients a (k) m k ,m k (t), representative of others. We focus on m = 1, 2, and for each operator size we design two basis operators, one diagonal and the other off-diagonal in the basis of Fock states. Following the prescription of Eq. (A.4) and fixing Q = N/2, we introducê where the prefactors N (k) m k ,m k are such that these operators have unitary Hilbert-Schmidt norm. Figure 1 presents the disorder-averaged dynamics of the coefficients a (k) m k ,m k (t), as well as the evolution for some individual disorder realizations. As expected from the exact theoretical discussion in Section 3.1, the operatorR has no overlap with the size (2, 2). The coefficient a (off-diag) 1,1 (t) is also vanishing because the staggered magnetization only involves number operatorsn i . In contrast, a (diag) 1,1 (t) manifests a quick decay to zero, mirroring the super-exponential time-dependence of the expectation value R (t) itself. This is substantiated by the inset in Figure 1a, which shows that the time-dependence of E a diag 1,1 (t) coincides with that of E R (t) apart from fluctuations that can be attributed to the statistical error of the averages. If we consider the dynamics of a single realization of the SYK model, in generalR(t) acquires components of size different from (1, 1), and it does not remain diagonal in the basis of Fock states. Still, Figure 1a shows that the typical time-evolution is close to the disorder-averaged one. This effect is a consequence of self-averaging: In fact, we observe larger deviations from the average curve as we decrease the system size.

Application of the cumulant expansion
We now test the performance of the cumulant expansion by comparing it to the dynamics obtained through numerically exact simulations. As discussed in Section 3, each constituent of an operatorŴ with well-defined size will evolve differently. The disorder-averaged dynamics thus reads E where w m,m = Tr Q ρ 0 P m,mŴ , and P m,m is a projection superoperator over the subspace with operator size (m, m). Notice that no components of size (m, n) with m = n can appear becauseŴ preserves the total charge. We explicitly see that the only role of the initial state is to set the amplitudes w m,m . In our previous work [34], we provide numerical evidence that, for some observables, the time-dependence of E Ŵ (t) manifests universality with respect to the choice of the initial state. This is readily explained by Eq. operator size evolve according to a single dynamical function, and the initial state affects only the amplitude. In contrast, if two or more coefficients (excluding w 0,0 ) are non-zero, then changingρ 0 modifies the relative weights of different dynamical functions, and no universality is found. To summarize, the shape of the dynamics of an operator with welldefined size is independent of both the initial state and the precise definition of the operator itself.
We study the dynamics of the staggered magnetizationR and of its squareR 2 . The former has operator size (1, 1), while the latter contains both sizes (0, 0) and (2,2). Figure 2 presents the disorder-averaged dynamics of the operators for N = 8, 12, 16, comparing it to the results of the cumulant expansion method truncated at different orders. Focusing on the operatorR 2 , for all considered system sizes, the quality of the approximation is found to improve with the number of cumulants considered. In addition, the plots suggest that the curves converge rapidly to the exact result as more cumulants are included; even though this claim cannot be rigorously confirmed with only three cumulant superoperators, it agrees with the discussion of Section 4, where we argued that high-order corrections provide negligible contributions as compared to low-order ones. For the operatorR, instead, some curves manifest unphysical divergencies at late times. This happens because a cumulant eigenvalue λ . The same is true in the case of R when adding C 4 for N = 8, 12, whereas long-time divergences appear when adding C 6 for all considered system sizes, and already when adding C 4 for N = 16. These pathologies are due to an inappropriate truncation of the cumulant at some system sizes. We also point out that generally (in the time windows preceding divergencies, if present) the variation induced by introducing the fourth cumulant is larger than that due to the sixth cumulant, which suggests quick convergence of the cumulant series.
to divergent long-time behavior (see Eq. (4.5)). Instead, one should evaluate higher-order cumulants, until a negative eigenvalue is found, and the series can be truncated safely.
We point out that also λ (4) 2,2 and λ (6) 2,2 , which characterize the dynamics ofR 2 , acquire positive values for larger system sizes than those considered here. Any divergence must be compensated by the presence of a higher-order cumulant with a negative eigenvalue that restores the relaxation to the steady state. As previously exposed in Ref. [34], the steady state of the ensemble-averaged dynamics corresponds to the infinite temperature state. Some deviations are observed for small system sizes, but these can be attributed to finite-size effects. This result is reproduced by the cumulant expansion approach, as the identity is an eigenoperator of all cumulants, and is unique for all cases which we studied.

Discussion and conclusions
In this work, we have studied the out-of-equilibrium properties of the SYK model in terms of the spectral features of the effective time-evolution superoperator for the disorder-averaged dynamics. Owing to the absence of any spatial ordering, it is possible to identify the exact eigenoperators of E e −Lt , and to prove that its (time-dependent) spectrum manifests high degeneracies, corresponding to different operator sizes. As a result, for all system sizes, the dynamics is characterized by a non-trivial symmetry, which hides the scrambling properties of the model when probing the evolution of standard expectation values of observables in the ensemble average. This conclusion is limited to correlation functions that involve a single time-evolution superoperator, and thus it does not apply to multiple-time functions nor, specifically, OTOCs. Our findings prove rigorously that the dynamics of observables with well-defined operator size manifest universal features with respect to the initial state, as highlighted by our previous numerical study [34].
The proof of operator size symmetry implies the existence of dynamical functions, obtained from the diagonalization of E e −Lt , that completely characterize the dynamics. Nevertheless, the demonstration does not indicate a simple way to compute them explicitly, and thus we developed a cumulant expansion scheme to achieve an approximation thereof. By construction, cumulant superoperators inherit the spectral structure of the exact timeevolution generator, and thus their eigenvalues can be determined analytically, as functions of the system size N and total charge Q. A direct comparison between this approach and exact diagonalization simulations reveals that the method is successful in reproducing the evolution to good accuracy in all time regimes. In addition, numerical results corroborate theoretical arguments that the cumulant expansion converges quickly with the number of terms included. Although this approximation scheme appears to be very effective for the system considered here and, potentially, variants thereof, we expect its usefulness to be limited for generic models: In the absence of operator size symmetry, it is impossible to determine the eigenoperators of the cumulants analytically, which implies one would need to resort to computationally expensive explicit numerical diagonalization. This difficulty is overcome for the SYK model, thus enabling the cumulant expansion method to effectively compute the dynamical functions.
We believe our investigation enriches the current understanding of the SYK model by showcasing the presence of an unexpected symmetry, and by illustrating that scrambling can be completely absent in some physical quantities even if the model is quantum chaotic. While the main focus of our work is the complex SYK model, our findings directly extend to its Majorana counterpart, as argued in Section 3.1. We expect our results to be relevant for future experimental implementations of the SYK model, as regular equal-time expectation values of observables are typical measurable quantities, further highlighting the need for measurement schemes that access more intricate observables [50]. In addition, the cumulant expansion approach that we developed manages to address arbitrary system sizes, thus enabling to access a regime that is out of reach of practical numerical simulations. We stress that a large enough system is fundamental to exploit the self-averaging property, which allows one to compare individual measurements to disorder-averaged quantities.
A natural topic of future investigation is whether it is possible to generalize our open system framework to study OTOCs, and thus to potentially reveal generic features also for such quantities. Moreover, as mentioned previously, the present work bears implications for other systems: We believe that similar features to the ones discussed here are also manifested by other fully-connected disordered models. In particular, the proof of operator size symmetry does not rely on the details of the SYK model, and can be extended to fermionic and spin systems with other disorder distributions, as well as to bosonic models with Gaussian disorder.

A Definition of operator size
As mentioned in Section 3.1, strings of creation and annihilation operators acting on different lattice sites allow for an unambiguous definition of operator size. In contrast, number operators do not have a well-defined size, because they contain a contribution proportional to the total chargeQ, which is proportional to the identity when acting on any chargeconserved subspace. This argument suggests that the definition of operator size, in the most general case, should take into account the specific value of the charge Q we are fixing. Therefore, throughout the following discussion, we limit the action of operators on kets with a fixed total charge Q.
The linear space of operators has the structure of a Hilbert space after introducing the Hilbert-Schmidt inner product Â B = Tr Â †B . (A.1) Any operatorÂ is represented as Â using a ket notation. Focusing on a specific charge sector, the overlap between two operators can be quantified through the charge-constrained Hilbert-Schmidt product, namely where Tr Q denotes the trace over the Hilbert subspace with fixed charge Q. Notice that, in general, the charge-constrained trace does not have the cyclic property of the standard trace.
Having introduced the necessary tools, we now proceed to the definition of operators with fixed size. In order to understand the general case, it is instructive to first focus on those operators for which we are already able to provide a rigorous definition, namely strings of typeŜ m,n =ĉ † i 1 . . .ĉ † imĉ j 1 . . .ĉ jn , where all indices are distinct. As discussed in Section 3.1, such an operator has size (m, n). LetŜ m ,n be another operator of this type, with (m, n) = (m , n ). It is easily checked that these operators are orthogonal with respect to the charge-constrained Hilbert-Schmidt inner product Ŝ m,n Ŝ m ,n 3) can be used to generalize the notion of size to operators involving one or moren i . Specifically, we may extract a component with well-defined size from an operator of typeŜ m,nni 1 . . .n i k by orthogonalizing it with respect to all operators of lower sizes. We now show how to perform this procedure in an iterative way. Suppose that all operators with well-defined sizes (m , n ) with m ≤ m and n ≤ n are known. It is always possible to choose suitable linear combinations of them to obtain an orthonormal set. Let us denote the basis elements asT (k) m k ,n k , where k is an index that counts them, and (m k , n k ) labels the size. We now pick a generic operatorŜ m,n with welldefined size (m, n). Since the operatorn i is a combination of terms with sizes (1, 1) and (0, 0), the productŜ m,nni contains a component of size (m + 1, n + 1), which is given bŷ It is easily checked that, with this definition,Ŝ m+1,n+1 is orthogonal to all basis elements, and thus it does not belong to the manifold with operator sizes below or equal to (m, n). We conclude that this operator has size (m + 1, n + 1). Proceeding in this way, we can formally determine all operators of size (m + 1, n + 1) that are required, together with the trivial ones that do not involve anyn i , to build an orthonormal set that spans all operators in the size sector (m + 1, n + 1). Once these are known, we can iterate this procedure to generate operators with even larger sizes. Finally, replacing Q →Q provides a general definition ofŜ m,n , independent of the Q-sector. We point out that the actual computation of operators with large sizes can be quite expensive. Still, for practical purposes, one is typically interested in studying only operators with small sizes, where the implementation of the previous method is viable.

B Proof of operator size symmetry
This Appendix complements the proof of operator size symmetry presented in Section 3 by providing a more detailed derivation. Let us definê which is (up to a constant) of the form of the summands entering Eq. (3.1). We consider the first non-trivial term, namelŷ where we made use of the disorder properties given in Eq. (2.2). Each of the four terms in the previous equation can be expanded in strings of fermionic operators. For example, For each of them, we can bringĉ i to the left by using the anticommutation relations. In particular: • if there is an annihilation operatorĉ k on the immediate left ofĉ i , we directly swap them, which yields a minus sign; • if, instead, there is a creation operatorĉ † k , we haveĉ † kĉ i = δ i,k −ĉ iĉ † k . The Kronecker delta restores the presence ofĉ i by constrainingĉ k , present somewhere else on the fermionic string, to have k = i. Both terms originated by the anticommutation relation still involve a singleĉ i and sums of pairedĉ † k andĉ k operators. Eventually, the operatorĉ i is brought to the left of all fermionic strings. What remains on the right are homogeneous sums of operators, in the sense that they do not depend on the lattice index i, nor any other lattice index in particular: These can always be written in terms of the lattice size N and the chargeQ. Finally, the discussion can be generalized to each operatorÂ n , and thus Eq. (3.2) follows.
The result we just proved forĉ i is immediately generalized to any string of fermionic operators with distinct lattice indices, as defined above Eq. (3.3). In fact, we can repeat the same procedure for eachĉ ( †) k independently, as all creation and annihilation operators anticommute; this leads to Eq. (3.3). In contrast, the validity of this result is not obvious for operators of the form given in Eq. (A.4). We can, however, generalize the proof using induction. Consider a generic operatorŜ m,n of size (m, n) (possibly involving also number operators), and assume that Eq. (3.3) holds for all sizes below or equal to (m, n). Under this inductive hypothesis, we want to prove that the eigenvalue equation is also valid for all operators of size (m + 1, n + 1). Specifically, we need to show that the operatorŜ m+1,n+1 , defined as in Eq. (A.4), also satisfies the eigenvalue equation. For this purpose, let i and j be lattice indices that do not appear in the definition ofŜ m,n , and consider the action of the disorder-averaged time-evolution superoperator onŜ m,nĉ † iĉ j . We put no constraint on the values of i and j themselves, they may be equal or different. Throughout the procedure of moving the operators to the left, due to the appearance of Kronecker deltas when using the anticommutation relations, it is not guaranteed thatĉ † i andĉ j remain in this order. In addition, some operators of sizes below or equal to (m, n) may appear ifŜ m,n itself has been obtained through the orthogonalization procedure described in the previous section. 6 In general, we find whereÂ,B,Â k , andB k are functions of t 2 , N , andQ, but not of i nor j, and the operatorsT (k) m k ,n k form an orthonormal basis for all sizes (m , n ) with m ≤ m and n ≤ n, as introduced in Appendix A. For i = j, the operatorŜ m,nĉ † iĉ j is guaranteed to have size (m + 1, n + 1). In particular, sinceŜ m,n satisfies Eq. (3.3) by assumption, thenŜ m,nĉ † iĉ j must fulfill the same equation (with (m, n) → (m + 1, n + 1)) because we are simply adding unpaired creation and annihilation operators. Requiring that the eigenvalue equation is recovered for i = j yields the conditionsÂ −B =f m+1,n+1 andÂ k =B k . It follows that for i = j we obtain Notice thatŜ m,n belongs to the space spanned by the basis elementsT (k) m k ,n k , and thus we can absorbŜ m,nB in the last sum of Eq. (B.5) by redefining the coefficientsÂ k . For brevity, we do not change their notation, and we have We will now relateÂ k to the dynamical functionsf m ,n . Adopting the bra-ket notation for operators, introduced in Appendix A, the action of the disorder-averaged time-evolution superoperator on the bra T (k) m k ,n k is defined as We consider the following matrix element, which can be written in two ways by using where we usedf † m,n =f m,n , 7 as well as the cyclicity of the trace. SinceÂ k ,f m+1,n+1 , and f m,n all depend onQ, the trace is conveniently decomposed as a sum of traces on subspaces with fixed charge Q, leading to We now argue that each element of the sum over Q must vanish individually. This is a consequence of charge conservation. Suppose we apply the disorder-averaged time-evolution superoperator toŜ m,nniĝ , whereĝ =ĝ(Q) is an arbitrary function of the total charge operator. Charge conservation implies thatĝ can be extracted from each commutator of Eq. (3.1), so that the analogue to Eq. (B.6) reads E e −Lt Ŝ m,nniĝ =Ŝ m,nnifm+1,n+1ĝ + kT (k) m k ,n kÂ kĝ . (B.11) Repeating the same calculations done previously results in a modified version of Eq. (B.10), in which additional weights g(Q) appear in the sum. Since the function is arbitrary, each element of the sum must be zero, and we finally obtain

C Analytic eigenvalues of cumulant superoperators
Cumulant superoperators can be analytically diagonalized by applying them to fixed-size operatorsŜ m,n , and performing the procedure described in Appendix B manually. The calculation can be carried out in an algorithmic way. Here, we present some exact results valid for the SYK 4 model, and for the operators of smallest sizes. For the first non-vanishing cumulant C 2 , we have:λ  For higher-order cumulants, the calculation can be performed analogously but becomes increasingly cumbersome. All these expressions are rational functions of the number of  and thus some eigenvalues are not independent. We numerically verified the exactness of some of the previous formulae (as well as some others that we do not provide here explicitly) by building matrix representations of the cumulant superoperators, and diagonalizing them in charge-conserved sectors of the Hilbert space. Specifically, we verifiedλ 1,1 ,λ 2,2 ,λ