\section{Introduction}\label{s:intro}
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~\cite{srednicki1994,polkovnikov2011,rigol2008,gogolin2016,dalessio2016,deutsch2018}. Theoretical advances in recent years have highlighted the importance of quantum information scrambling~\cite{sekino2008,iyoda2018,hosur2016} 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~\cite{kitaev2015,maldacena2016,gu2020,chowdhury2022}, which received much interest in recent research due to its relevance in multiple branches of physics, ranging from the study of black holes~\cite{kitaev2015,davison2017,kitaev2018,sachdev2019,sarosi2018} to non-Fermi liquids~\cite{sachdev2015,song2017,cha2020,sachdev2010}. 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\geq 4$, the model manifests quantum chaotic behavior~\cite{kitaev2015,polchinski2016,garciagarcia2018}, diagnosed by the presence of a quantum Lyapunov exponent in out-of-time-order correlators (OTOCs)~\cite{maldacena2016bis,hashimoto2017,kobrin2021}. Specifically, this exponent is found to saturate its theoretical upper bound~\cite{maldacena2016bis}, implying that quantum scrambling occurs as rapidly as possible. While multiple properties can be derived exactly in the thermodynamic limit~\cite{kitaev2015,sachdev2015,chowdhury2022}, studying the model at finite size remains a challenging yet essential problem, especially in view of future experimental investigations that have recently been proposed~\cite{danshita2017,garciaalvarez2017,pikulin2017,chew2017,chen2018,luo2019,wei2021}.

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, leading to complete absence of operator growth on average. This corroborates the results of our previous work~\cite{bandyopadhyay2021}, 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 well-known results from the literature. While the lack of operator growth does a priori not apply to individual realizations of the model, it is 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, owing to the absence of operator growth, 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~\ref{s:model} presents the model and the superoperator framework adopted throughout the work. In Sec.~\ref{s:operator_size}, 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 Sec.~\ref{s:cumulant_expansion} we present the cumulant expansion technique that we implement to approximate the dynamics. Our theoretical findings are tested in Sec.~\ref{s:application}, 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 Sec.~\ref{s:conclusions}, 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. 
	
\section{Model and focus of our study}\label{s:model}

We focus on the family of SYK models of $N$ complex fermion modes on a fully connected lattice, defined by the Hamiltonian
\begin{subequations}\label{hamiltonian}
\begin{equation}\label{hamiltonian-a}
    \hat{H} = \mathcal{K}_q \sum_{\substack{i_1<\dots<i_{q/2}\\j_1<\dots<j_{q/2}}} J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}} \hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_{q/2}}\hat{c}_{j_1}\dots\hat{c}_{j_{q/2}},
    \end{equation}\begin{equation}
    \mathcal{K}_q = \sqrt{\frac{(q/2)!(q/2-1)!}{N^{q-1}}},
    \end{equation}
\end{subequations}
where $q$ is an even integer number, and $J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}$ are complex Gaussian independent random couplings with the following statistical properties ($\eave{...}$ denotes disorder averaging):
\begin{subequations}\label{disorder_properties}
\begin{align}
\label{disorder_properties-a}
& \eave{J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}}=0,\\
\label{disorder_properties-b}
& \eave{\left(J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}\right)^2}=0 \\
    & \quad \text{if } \{i_1,\dots,i_{q/2}\}\neq\{j_1,\dots,j_{q/2}\}, \nonumber\\
    \label{disorder_properties-c}
 & \eave{\left|J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}\right|^2} = J^2.
\end{align}
\end{subequations}
Notably, for $q>2$ the only integral of motion of any individual realization is the total charge $\hat{Q}=\sum_i \hat{n}_i$. 

Our analysis relies on the study of the dynamics of observable expectation values, namely
\begin{equation}
    \braket{\hat{W}(t)}  = \Tr\left(\hat{\rho}_0e^{i\hat{H}t} \hat{W} e^{-i\hat{H}t} \right),
\end{equation}
where typically (but not necessarily) $\hat{\rho}_0 = \ket{\psi_0}\bra{\psi_0}$ is a pure-state density matrix:
This corresponds to preparing the system in a disorder-independent initial state $\ket{\psi_0}$, performing a quantum quench to the SYK$_q$ Hamiltonian, and tracking the time-evolution of physical quantities, as carried out, e.g., in~\cite{bandyopadhyay2021}.
Specifically, since the Hamiltonian involves randomness, we focus on the disorder-averaged evolution $\eave{\braket{\hat{W}(t)}}$. A well-known feature of disordered systems is that disorder-averaging introduces mixing, making the average dynamics non-unitary~\cite{kropf2016}. As a consequence, it becomes natural to adopt a superoperator-based framework, typical of dissipative open quantum systems~\cite{minganti2018,manzano2020}. We introduce the Liouvillian $\mathcal{L} \bullet = -i\comm{\hat{H}}{\bullet}$, as well as the superoperators $\ell_\alpha \bullet = -i \comm{\hat{h}_\alpha}{\bullet}$, where $\alpha = \{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}\}$ is a multi-index, and $\hat{h}_\alpha = \hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_{q/2}}\hat{c}_{j_1}\dots\hat{c}_{j_{q/2}}$ appears in the Hamiltonian of Eq.~\eqref{hamiltonian-a}. With this notation, the ensemble-averaged dynamics reads
\begin{equation}\label{observable_dyn}
    \eave{\braket{\hat{W}(t)}} = \Tr\left(\hat{\rho}_0 \eave{e^{-\mathcal{L} t}}\hat{W} \right).
\end{equation}
In what follows, we provide a detailed analysis of the time evolution of such disorder-averaged expectation values. 

\section{Operator size symmetry in the disorder-averaged ensemble}\label{s:operator_size}

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. This leads us to the conclusion that 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. 

\subsection{Operator size symmetry}\label{ss:size_sym}

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)$~\cite{carrega2021,roberts2018}. 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 operators $\hat{n}_i$. In any charge-conserved sector of the Hilbert space, the operator $\hat{Q}=\sum_i \hat{n}_i$ is proportional to the identity, and thus its size is $(0,0)$. As a consequence, the size of the number operators $\hat{n}_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 App.~\ref{a:op_size}). 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 App.~\ref{a:proof_op_size}. 

Taking the disorder average of a simple annihilation operator $\hat{c}_i(t)$ in the Heisenberg picture gives
\begin{equation}\label{expand_exp}
    \eave{e^{-\mathcal{L} t}}\hat{c}_i = \sum_{n=0}^\infty \frac{(\mathcal{K}_q t)^{2n}}{(2n)!}\eave{J_{\alpha_1}\dots J_{\alpha_{2n}}}\ell_{\alpha_1} \dots\ell_{\alpha_{2n}}\hat{c}_i,
\end{equation}
where the superoperators $\ell_{\alpha_i}$ are defined at the end of Sec.~\ref{s:model}, 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 $\hat{h}_{\alpha_1},\dots, \hat{h}_{\alpha_n}$, as well as their Hermitian conjugates $\hat{h}^\dag_{\alpha_1},\dots, \hat{h}^\dag_{\alpha_n}$, appearing in a certain order. In particular, each multi-index $\alpha_k$ is shared precisely by $\hat{h}_{\alpha_k}$ and $\hat{h}^\dag_{\alpha_k}$, as follows from Eqs.~\eqref{disorder_properties-b} and~\eqref{disorder_properties-c}. If we expand each nested commutator in sequences of fermionic operators, for each $\hat{c}^\dag_k$ appearing in any string there will always be exactly one $\hat{c}_k$, with the same index, somewhere within that same string.
For each chain of fermionic operators, we may use the anticommutation relations to move $\hat{c}_i$ towards one of the edges, for example the left one. Each new sub-string originated by this process involves only the operator $\hat{c}_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 $\hat{c}_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 operator $\hat{Q}$. In conclusion, we have proved that
\begin{equation}\label{absence_growth}
    \eave{e^{-\mathcal{L} t}}\hat{c}_i = \hat{c}_i \hat{f}(t^2;N,\hat{Q}),
\end{equation}
where the dependence on a function of $t^2$ arises from the absence of odd terms in Eq.~\eqref{expand_exp}. The above result shows that, if we limit its action to a subset of the Hilbert space with fixed $Q$, $\hat{c}_i$ is an eigen-operator of the generator of the ensemble-averaged dynamics. In other words, $\hat{c}_i(t)$ manifests no operator growth (nor hopping) on average. 

Our above argument is readily extended to all fermionic strings of type $\hat{S}_{m,n}=\hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_m}\hat{c}_{j_1}\dots\hat{c}_{j_n}$ in which all indices are distinct, i.e.,
\begin{equation}\label{absence_growth_general}
    \eave{e^{-\mathcal{L} t}}\hat{S}_{m,n} = \hat{S}_{m,n} \hat{f}_{m,n}(t^2;N,\hat{Q}).
\end{equation}
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 App.~\ref{a:proof_op_size}. Notice that the functions $\hat{f}_{m,n}$ are independent of the specific lattice sites appearing in $\hat{S}_{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: There is no growth, and the time-dependence of the evolution depends only on the operator size. 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. In addition, we expect the result to extend also to fermionic and spin SYK-like models with non-Gaussian disorder, as well as a bosonic versions of the system with Gaussian randomness~\footnote{Gaussian disorder implies that the superoperators $\ell_\alpha$ are contracted in pairs. For fermionic systems, this is not fundamental, because the square of a creation or annihilation operator vanishes: Even if multiple operators with the same lattice index appear within a fermionic string, eventually the expression simplifies leaving at most one of them. A similar argument applies to spin models. In contrast, more-than-pairwise contractions generate new distinct operators for bosons, and the proof of operator size symmetry does not hold.}.


\subsection{Implications for correlation functions and OTOCs}\label{ss:single_realizations}

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 manifest present operator growth (see Sec.~\ref{ss:single_realizations} 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~\cite{castellani2005,torresherrera2020}, which are precisely in the form of Eq.~\eqref{observable_dyn}. As a consequence, if, for instance, $\hat{W}=\hat{S}_{m,m}$ has a fixed operator size $(m,m)$, and $\hat{\rho}(0)$ has fixed number of particles $Q$, we have
\begin{equation}
    \braket{\hat{W}(t)} \approx \eave{\braket{\hat{W}(t)}} = f_{m,m}(t^2;N,Q)\Tr\left(\hat{W} \hat{\rho}_0\right),
\end{equation}
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~\cite{bandyopadhyay2021}.

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 time-evolution superoperator. In particular, OTOCs of the form
\begin{equation}
    \eave{\Tr\left(\comm{\hat{W}(t)}{\hat{V}(0)}^2\right)},
\end{equation}
which diagnose quantum chaotic behavior, are beyond the scope of our investigation, because they involve correlations between two superoperators $e^{-\mathcal{L} t}$. We conclude that quantum information scrambling must be encoded in these correlations, which arise because operators themselves are not self-averaging, namely $\hat{W}(t)\neq\eave{\hat{W}(t)}$ even if their quantum-mechanical expectation values become equal at large $N$. The same argument applies to the so-called Krylov complexity~\cite{parker2019,rabinovici2022,jian2021}, 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~\cite{roberts2018,lensky2020,qi2019,lucas2020}. 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~\cite{lin2019,jian2021}) 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.


\section{Cumulant Expansion Method}\label{s:cumulant_expansion}

In the previous section, we discussed how operators with well-defined size do not grow on ensemble average. Nevertheless, there is no straightforward way to explicitly determine the functions $\hat{f}_{m,n}(t^2;N,\hat{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.~\cite{bandyopadhyay2021} 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
\begin{equation}\label{cumulant_series}
    \eave{e^{-\mathcal{L} t}} = e^{\mathcal{C}(t)} = \exp\left(\sum_{k=1}^\infty \frac{(Jt)^k}{k!}\mathcal{C}_k\right),
\end{equation}
where $\mathcal{C}(t)$ will be referred to as a cumulant generating superoperator, in agreement with usual nomenclature in statistics. Accordingly, the superoperators $\mathcal{C}_k$ will be referred to as cumulants. Since we are considering Gaussian disorder, if $\mathcal{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.~\eqref{cumulant_series}, and comparing it to Eq.~\eqref{expand_exp} by matching equal powers of $t$. We immediately see that all superoperators $\mathcal{C}_k$ with odd $k$ vanish. The identification of $\mathcal{C}_k$ requires the previous determination of all $\mathcal{C}_l$ with $l<k$, and thus the procedure is iterative. Here, we present the first two non-vanishing terms:
\begin{subequations}
\begin{equation}\label{c2}
    \mathcal{C}_2 = \mathcal{K}_q^2 \sum_\alpha \ell_\alpha^2,
\end{equation}
\begin{equation}\label{c4}
     \mathcal{C}_4 = \mathcal{K}_q^4 \sum_{\alpha,\beta}\left(\ell_\alpha\ell_\beta^2\ell_\alpha + \ell_\alpha\ell_\beta\ell_\alpha\ell_\beta-2\ell_\alpha^2\ell_\beta^2\right).
\end{equation}
\end{subequations}
Leading order truncation of the cumulant expansion yields precisely the time-evolution superoperator used in Ref.~\cite{bandyopadhyay2021}, 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 $\ell_{\alpha_1}\dots\ell_{\alpha_{2k}}$, in which the indices are contracted in pairs and summed over, preserves operators with well-defined size, exactly as the full $\eave{e^{-\mathcal{L} t}}$ in Eq.~\eqref{absence_growth_general} (see App.~\ref{a:proof_op_size} for the details). As a consequence, we can differentiate Eq.~\eqref{cumulant_series} to obtain an effective master equation for the disorder-averaged density matrix in the Schrödinger picture~\footnote{Notice that $e^{\mathcal{C}(t)} = \eave{e^{-\mathcal{L} t}} = \eave{e^{\mathcal{L} t}}$ because all odd powers in the series expansion vanish.}, yielding
\begin{equation}\label{me}
    \dv{}{t}\eave{\hat{\rho}(t)} = \dv{}{t}e^{\mathcal{C}(t)}\hat{\rho}_0= \left(J\sum_{k=1}^\infty \frac{(Jt)^{2k-1}}{(2k-1)!}\mathcal{C}_{2k}\right)\eave{\hat{\rho}(t)}.
\end{equation}
The sum within the round brackets is the effective Liouvillian that generates the disorder-averaged dynamics. Equation~\eqref{me} is a generalization of the master equation proposed in Ref.~\cite{bandyopadhyay2021}, and it is formally exact.

For practical purposes, we can only evaluate the cumulant expansion up to some given finite order. Even if all $\mathcal{C}_k$ were known, finding an analytic expression for $\mathcal{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~\cite{groszkowski2022}. 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 $\mathcal{C}(t)$ is truncated, Eq.~\eqref{cumulant_series} 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)!$~\footnote{The timescale at which a cumulant becomes relevant is long if its magnitude in Eq.~\eqref{expand_exp} is small. Any $\mathcal{C}_{2k}$ with $k>1$ should vanish if all superoperators $\ell_\alpha$ commuted, and thus each cumulant must contain a certain number $2M_{2k}$ of superoperators of type $\ell_{\alpha_1}\dots\ell_{\alpha_{2k}}$, half of them with the positive sign, half with the negative one (for example, see Eq.~\eqref{c4}): We do not expect $\mathcal{C}_{2k}$ to grow as fast as the factor $(2k)!$, appearing at the denominator. On top of that, $M_{2k}$ is the number of irreducible Wick contractions arising from the $2k$th term of Eq.~\eqref{expand_exp}, 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.~\eqref{c4}, 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.}. Appendix~\ref{a:exact_evals} 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 $\eave{\hat{\rho}(t)}\approx\eave{\hat{\rho}(\infty)}$ for $t \gtrsim\tau$. In this case, neglecting high-order cumulants that become relevant only at times greater than $\tau$ is a good approximation, because, as seen from Eq.~\eqref{me}, 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.~\cite{bandyopadhyay2021}, 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 functions $\hat{f}_{m,n}(t^2;N,\hat{Q})$. First, as observed previously, an equation analogous to Eq.~\eqref{absence_growth_general} can be written for each individual cumulant, so that
\begin{equation}
    \mathcal{C}_{2k} \hat{S}_{m,n} = \hat{S}_{m,n}\hat{\lambda}^{(2k)}_{m,n}(N,\hat{Q}),
\end{equation}
where $\hat{S}_{m,n}$ is a generic operator in the form of Eq.~\eqref{size_operators}. The dynamical functions $\hat{f}_{m,n}(t^2;N,\hat{Q})$ are then expanded as
\begin{equation}\label{expand_fmn}
    \hat{f}_{m,n}(t^2;N,\hat{Q}) = \exp\left(\sum_{k=1}^\infty \frac{(J t)^{2k}}{(2k)!}\hat{\lambda}^{(2k)}_{m,n}(N,\hat{Q})\right).
\end{equation}
As argued before, a finite number of functions $\hat{\lambda}^{(2k)}_{m,n}$ can be sufficient to obtain a very good approximation of the exact result. The advantage of this approach is that eigen-operators of cumulants are known, and thus we can evaluate their eigenvalues analytically. Some exact expressions for the SYK$_4$ model are presented in App.~\ref{a:exact_evals}. From a practical point of view, this greatly simplifies the numerical application of our formalism: We can approximate Eq.~\eqref{expand_fmn} for arbitrarily large system sizes without resorting to exact diagonalization of cumulant superoperators.
 

\section{Numerical results}\label{s:application}

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 $\ket{\psi_0}$, and then we evaluate the expectation value $\braket{\hat{W}(t)}$ of an observable $\hat{W}$ that commutes with the charge $\hat{Q}$. In our numerics, we choose $\ket{\psi_0}$ to be the Néel state $\ket{1010\dots}$, where the ordering of lattice sites is arbitrary. Regarding the choice of the operator, we study the staggered magnetization $\hat{R} = \sum_{i=1}^N(-1)^{i+1}\hat{n}_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.~\eqref{hamiltonian-a}, and finally taking a statistical average over the sample.

\subsection{Operator size throughout the dynamics}\label{ss:operator_growth}

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 App.~\ref{a:op_size}, we denote these basis operators as $\hat{T}^{(k)}_{m_k,n_k}$, and each of them has size $(m_k,n_k)$. We can extract their amplitudes in the expansion of $\hat{W}(t)$ by taking the Hilbert--Schmidt product
\begin{equation}\label{c_mn}
    c^{(k)}_{m_k,n_k}(t) = \Tr_Q\left(\left(\hat{T}^{(k)}_{m_k,n_k}\right)^\dag \hat{W}(t)\right),
\end{equation}
where $\Tr_Q$ denotes the trace over states with fixed charge $Q$, as defined in App.~\ref{a:op_size}. For $\hat{W}=\hat{R}$, we use exact diagonalization to study the dynamics of some coefficients $c^{(k)}_{m_k,m_k}(t)$, representative of others. We focus on $m=1,2$, and for each operator size we design two operators, one diagonal and the other off-diagonal in the basis of Fock states. Following the prescription of Eq.~\eqref{size_operators}, we introduce
\begin{subequations}\label{test_operators}
\begin{align}
\label{test_operators-a}
& \hat{T}_{1,1}^{(\text{diag})} = \mathcal{N}^{(\text{diag})}_{1,1}\left(\hat{n}_1-\frac{1}{2}\right),
\\
\label{test_operators-b}
& \hat{T}_{2,2}^{(\text{diag})} = \mathcal{N}^{(\text{diag})}_{2,2}\left(\hat{n}_1\hat{n}_2-\frac{\hat{n}_1+\hat{n}_2}{2} + \frac{N}{4(N-1)}\right),
\\
\label{test_operators-c}
& \hat{T}_{1,1}^{(\text{off-diag})} = \mathcal{N}^{(\text{off-diag})}_{1,1}\left(\hat{c}^\dag_1\hat{c}_2+\mathrm{h.c.}\right),
\\
\label{test_operators-c}
& \hat{T}_{2,2}^{(\text{off-diag})} = \mathcal{N}^{(\text{off-diag})}_{2,2}\left(\hat{c}^\dag_1\hat{c}^\dag_2\hat{c}_3\hat{c}_4+\mathrm{h.c.}\right),
\end{align}
\end{subequations}
where the prefactors $\mathcal{N}^{(k)}_{m_k,m_k}$ are such that these operators have unitary Hilbert--Schmidt norm. Figure~\ref{f:operator_growth} presents the disorder-averaged dynamics of the coefficients $c^{(k)}_{m_k,m_k}(t)$, as well as the evolution for some individual disorder realizations. As expected from the exact theoretical discussion in Sec.~\ref{ss:size_sym}, the operator $\hat{R}$ has no overlap with the size $(2,2)$. The coefficient $c^{(\text{off-diag})}_{1,1}(t)$ is also vanishing because the staggered magnetization only involves number operators $\hat{n}_i$. In contrast, $c^{(\text{diag})}_{1,1}(t)$ manifests a quick decay to zero, mirroring the super-exponential time-dependence of the expectation value $\braket{\hat{R}(t)}$ itself. This is substantiated by the inset in Fig.~\ref{f:operator_growth}a, which shows that the time-dependence of $\eave{c_{1,1}^\text{diag}(t)}$ coincides with that of $\eave{\braket{\hat{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 general $\hat{R}(t)$ acquires components of size different from $(1,1)$, and it does not remain diagonal in the basis of Fock states. Still, Fig.~\ref{f:operator_growth}a 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.
\begin{figure}[ht]
    \includegraphics[width =   
    \linewidth]{operator_growth.pdf}
    \caption{Dynamics of the coefficients $c^{(k)}_{m_k,m_k}(t)$ computed in single disorder realizations (red curves) and in ensemble average (blue curves) for the operator $\hat{R}$. The inset of panel (a) portrays the difference $\Delta(t) = \left|\eave{ c_{1,1}^\text{diag}(t)/c_{1,1}^\text{diag}(0)} - \eave{\braket{\hat{R}(t)}/\braket{\hat{R}(0)}}\right|$. Data is for $N=10$ at half filling $Q=N/2$. The ensemble averages are estimated by performing statistical averages over $10000$ independent disorder realizations. The terms of size $(2,2)$ vanish on average, in agreement with the predicted absence of operator growth. In contrast, there is a finite overlap with the diagonal operator of size $(1,1)$, and the time-evolution matches the full dynamics of $\braket{\hat{R}(t)}$, as can be appreciated by the inset of panel (a). Curves for individual realizations show that components with sizes different from $(1,1)$ are in general developed throughout the evolution.}
    \label{f:operator_growth}
\end{figure}

\subsection{Application of the cumulant expansion}\label{ss:cumulant_exp}

We now test the performance of the cumulant expansion by comparing it to the dynamics obtained through numerically exact simulations. As discussed in Sec.~\ref{s:operator_size}, each constituent of an operator $\hat{W}$ with well-defined size will evolve differently. The disorder-averaged dynamics thus reads
\begin{equation}\label{decompose_dynamics}
    \eave{\braket{\hat{W}(t)}} = \sum_m c_{m,m} f_{m,m}(t^2;N,Q),
\end{equation}
where 
\begin{equation}
    c_{m,m} = \Tr_Q\left(\hat{\rho}_0\mathcal{P}_{m,m}\hat{W}\right),
\end{equation}
and $\mathcal{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\neq n$ can appear because $\hat{W}$ preserves the total charge. We explicitly see that the only role of the initial state is to set the amplitudes $c_{m,m}$. In our previous work~\cite{bandyopadhyay2021}, we provide numerical evidence that, for some observables, the time-dependence of $\eave{\braket{\hat{W}(t)}}$ manifests universality with respect to the choice of the initial state. This is readily explained by Eq.~\eqref{decompose_dynamics}: Operators that overlap with a unique 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 $c_{0,0}$) are non-zero, then changing $\hat{\rho}_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 well-defined size is independent of both the initial state and the precise definition of the operator itself.

\begin{figure*}[ht]
    \includegraphics[width =   
    \linewidth]{ed_vs_ce_dynamics.pdf}
    \caption{ Comparison of disorder-averaged dynamics of $\hat{R}$ and $\hat{R}^2$ obtained using exact diagonalization (solid blue curve), and the cumulant expansion truncated at second (dashed orange curve), fourth (dashed red curve), and sixth (dashed purple curve) order. Data is for $(\mathrm{a})$ and $(\mathrm{d})$ $N=8$ with 100000 samples, $(\mathrm{b})$ and $(\mathrm{e})$ $N=12$ with 3000 samples, $(\mathrm{c})$ and $(\mathrm{f})$ $N=16$ with 500 samples, at half filling $Q=N/2$.
    For $\hat{R}^2$, in all cases the addition of an extra cumulant manifestly improves the approximation of the exact dynamics (see insets). The same is true in the case of $\hat{R}$ when adding $\mathcal{C}_4$ for $N=8,12$, whereas long-time divergences appear when adding $\mathcal{C}_6$ for all considered system sizes, and already when adding $\mathcal{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.
    }
    \label{f:numerics}
\end{figure*}

We study the dynamics of the staggered magnetization $\hat{R}$ and of its square $\hat{R}^2$. The former has operator size $(1,1)$, while the latter contains both sizes $(0,0)$ and $(2,2)$. Figure~\ref{f:numerics} 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 operator $\hat{R}^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 Sec.~\ref{s:cumulant_expansion}, where we argued that high-order corrections provide negligible contributions as compared to low-order ones. These conclusions hold also for the operator $\hat{R}$, with the clear exception of some curves that manifest unphysical divergencies at late times. This happens because a cumulant eigenvalue $\lambda_{m,n}^{(2k)}$ can have a positive sign for some values of $N$, which makes the truncation of the cumulant series to that specific order insufficient. We point out that also $\lambda_{2,2}^{(4)}$ and $\lambda_{2,2}^{(6)}$, which characterize the dynamics of $\hat{R}^2$, acquire negative 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.~\cite{bandyopadhyay2021}, 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.

\section{Discussion and conclusions}\label{s: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 $\eave{e^{-\mathcal{L} t}}$, 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~\cite{bandyopadhyay2021}.



The proof of operator size symmetry implies the existence of dynamical functions, obtained from the diagonalization of $\eave{e^{-\mathcal{L} t}}$, 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 time-evolution 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. 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~\cite{torlai2020}. 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.

\emph{Acknowledgements.---}We thank Andrea Legramandi for reading the manuscript and giving useful comments. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 804305), the Bundesministerium für Wirtschaft und Energie through the project ``EnerQuant” (project ID 03EI1025C), Provincia Autonoma di Trento, and by Q@TN, the joint lab between University of Trento, FBK-Fondazione Bruno Kessler, INFN-National Institute for Nuclear Physics and CNR-National Research Council. 
S.B. acknowledges CINECA for the use of HPC resources under ISCRA-C project ISSYK-2 (HP10CP8XXF). 











