Unveiling operator growth in SYK quench dynamics

We study non-equilibrium dynamics induced by a sudden quench of strongly correlated Hamiltonians with all-to-all interactions. By relying on a Sachdev-Ye-Kitaev (SYK) based quench protocol, we show that the time evolution of simple spin-spin correlation functions is highly sensitive to the degree of locality of the corresponding operators, once an appropriate set of fundamental fields is identified. By tracking the time-evolution of specific spin-spin correlation functions and their decay, we argue that it is possible to distinguish between operator hopping and operator growth dynamics; the latter being a hallmark of quantum chaos in many-body quantum systems. Such observation, in turn, could constitute a promising tool to probe the emergence of chaotic behavior, rather accessible in state-of-the-art quench setups.

We study non-equilibrium dynamics induced by a sudden quench of strongly correlated Hamiltonians with all-to-all interactions. By relying on a Sachdev-Ye-Kitaev (SYK) based quench protocol, we show that the time evolution of simple spin-spin correlation functions is highly sensitive to the degree of locality of the corresponding operators, once an appropriate set of fundamental fields is identified. By tracking the time-evolution of specific spin-spin correlation functions and their decay, we argue that it is possible to distinguish between operator hopping and operator growth dynamics; the latter being a hallmark of quantum chaos in many-body quantum systems. Such observation, in turn, could constitute a promising tool to probe the emergence of chaotic behavior, rather accessible in state-of-the-art quench setups.

I. INTRODUCTION
The study of strongly correlated quantum systems dates back to the early stages of quantum mechanics, and it still represents one of the most intriguing and challenging subjects of research [1][2][3][4][5][6][7]. Recent technological advances both in condensed matter [8][9][10] and atomic physics [11][12][13][14][15][16] allowed the investigation of manybody physics at the nanoscale, or single atom level, both in equilibrium and non-equilibrium settings. In particular, state-of-the-art experiments with ultracold atoms and trapped ions [6,17,18] offer the possibility to engineer closed quantum systems with very high precision and to perform quantum quench protocols [6,19,20]; where non-equilibrium dynamics can be measured in real time after a sudden variation of some parameters of interacting Hamiltonians. Notable results have been already achieved in the context of relaxation dynamics and equilibration properties of quantum many-body systems [16,21]. Interestingly, optical lattice designs can be implemented to simulate quantum systems in various dimensions, under the presence of both local and q-body photon-mediated interactions [22][23][24][25][26] between different lattice sites.
The above-mentioned progresses have attracted attention [21] due to their potential in testing thermalization hypotheses and related conjectures in strongly correlated systems, which exhibit intriguing connections with quantum chaos and black-hole physics [27][28][29][30][31]. Popular observables in this context are out-of-time-ordered correlators (OTOCs) [32][33][34] which have been recently used to quantify the chaotic nature of a given quantum system [35]. In fact, due to the exponential growth of the OTOCs in time when dealing with chaotic dynamics, they have been linked to the Lyapunov exponents of classical chaos, thereby being proposed as a promising measure of quantum chaos [33,34]. They quantify scrambling, or fast spreading of an initial local perturbation across the system. The spread of information in chaotic systems is intimately related to the notion of operator growth [36][37][38][39][40], i.e. to the idea that, during time evolution, local operators develop into rather complicated operators with increased spatial support and non-locality, resembling the so-called butterfly effect.
Although a few pioneering cold atom experiments have reported the possibility to measure the OTOCs in specific systems [34], an implementation of the OTOCs remains very hard since it demands the ability of measuring different operators involving time-reversal of many-body dynamics. Several alternative diagnostics of quantum chaos [41][42][43][44] have been thus proposed to circumvent this difficulty. Among them, Ref. [41] proposed that quantum chaos can be inferred by inspecting the temporal evolution of the full probability distribution of randomly prepared initial states after a quantum quench.
Inspired by this approach, in this work we focus on a paradigmatic example of dynamical system equipped with all-to-all interactions, i.e. the Sachdev-Ye-Kitaev (SYK) model [45][46][47]. It describes a strongly correlated quantum system of Majorana or Dirac fermions with random, all-to-all, q-body interactions, with q being an integer larger than or equal to 2. In particular, when q > 2, both the energy-level statistics and time-evolution of the OTOCs indicate that non-local interactions and random disorders make the model highly chaotic [46][47][48][49].
Early-time quantum chaos in the SYK model was studied, from a slightly different perspective than the study of the OTOCs, also in [39]. The authors thereof studied the operator growth dynamics of the SYK model, with a particular emphasis on how simple operators, i.e. consisting of products between few Majorana fermions, evolve into more complicated operators that involve products between a stack of Majorana fermions.
In this paper, we aim to study how the operator growth dynamics, linked to the onset of quantum chaos, can be arXiv:2007.03551v2 [quant-ph] 24 Jul 2020 detected by quantum quench protocols [50,51]. Under the specific protocol where the quench Hamiltonian takes the form of SYK models, with different values of q, we investigate the dynamics of different spin correlation functions in 1D lattice spin systems. We argue that these correlators, which are quite accessible in state-of-the-art quench experiments by exploiting local imaging of quantum gas microscopes [52], contain useful information on the operator dynamics.
In particular, we demonstrate that the decay rate of these correlation functions can be traced back to the operator growth dynamics under the SYK q quench Hamiltonian, depending on the specific value of q.
The paper is organized as follows. In section II, we introduce the quench protocol and the models under study. In section III, we examine the quenched evolution of connected spin-spin correlation functions. Specifically, we observe that a seemingly innocent modification, i.e. of the orientation of a static external magnetic field coupled to the spins in an initial local Hamiltonian, strongly affects the quench dynamics. In section IV, we explain the above observation in terms of the dynamical evolution of the "fundamental" operators in the theory. We introduce the notion of operator hopping -to be contrasted to the operator growth -which controls the dynamics under the integrable SYK Hamiltonians at q = 2. In section V, instead, we consider the case of the chaotic SYK 4 models. We find that the dynamics is governed by operator growth instead of operator hopping, as expected for chaotic systems, demonstrating how temporal evolution of spin-spin correlators reflects the nature of operator dynamics. Section VI summarizes our main results and possible perspectives.

II. MODEL AND QUENCH PROTOCOL
The main focus of this work is the study of nonequilibrium properties of strongly correlated quantum many-body systems. To this end, we consider a quench protocol, where the interactions between single entities of a quantum system are suddenly switched on at time t = 0 asĤ where θ(t) is a step function and the dynamics for t > 0 is entirely governed byĤ 1 in a scale invariant fashion. Such quantum quench protocols can be engineered, and have been realized, for example in cold-atoms setups [6]. The system is initially prepared (for t < 0) in the ground state of a generic free and localĤ 0 , i.e. a quantum system on a spin lattice of length L, and then it evolves under the action of a strongly interacting Hamiltonian with all-to-all couplings. For the latter, we consider the SYK q Hamiltonian that represents a paradigmatic example of all-to-all interacting systems [27,[45][46][47]. These models recently gained a lot of attention due to intriguing connections with black-hole physics and since they can show quantum chaotic behavior [27]. This family of Hamiltonians, classified in terms of an integer parameter q, can be written in terms of 2L Majorana fermions, interacting via q-body and all-to-all coupling terms as [45][46][47] where the Majorana fermion operators,γ i , satisfy the following Clifford algebra relations: and where the coupling constants J i1...iq are extracted from a Gaussian distribution, with null mean value and variance Throughout the paper, we set = 1 and the overline will denote the Gaussian average over all the SYK coupling constants J i1···iq . Unless otherwise specified, we also set J 2 = 1. While at q = 2 the SYK model has been known to show integrable behavior [53], the situation gets completely different for higher values of q > 2: although sharing all-to-all correlations, in this case, the model turns out to be chaotic based on the study of the OTOCs [46,47] as well as the energy-level statistics [48,49]. It is thus interesting to take both the q = 2 and q > 2 models as the quench Hamiltonian, investigating how they differently affect the relaxation dynamics, looking for features related to the presence of quantum chaos.
For sake of simplicity, we consider the free and local H 0 as an ensemble of non-interacting spin-1 2 variables on a lattice of length L immersed on a transverse magnetic field oriented along the a = x, y, z direction whereσ a i denotes the a-th Pauli matrix at site i and ω is the magnetic field strength, hereafter assumed to be the same for all directions (and equal to 1) without loss of generality. The spin Hamiltonian defined in Eq. (5) can be mapped into a system of 2L Majorana fermions via the following Jordan-Wigner (JW) map [54], which defines the duality between L spins and 2L Majorana fermions. Obviously, also the SYK q Hamiltonians can be mapped to the spin chain variables via the JW map (6), but this would result in expressions that are highly non-local and not easy to put in a compact form. The lattice spin Hamiltonians (5) can be realized in several settings, such as cold atoms, trapped ions or solid state devices. For example, one can engineer a system of neutral atoms with hyperfine interactions and coupling them to Rydberg states [55][56][57][58]. A transverse field can be then introduced by applying a resonant microwave or Raman coupling between two hyperfine states, and several type of correlations, ranging from local to all-to-all interactions, have been proposed and recently realized, exploiting photon-mediated correlations in optical cavities as well [22,23]. In passing, we mention that some recent proposals have suggested the possibility to realize SYKlike Hamiltonians in all these settings [27,30,59,60].
We are interested in tracing the time-evolution of the ground states, |0 a , after the action of the quench protocol of Eq. (1). In order to make fair comparisons of timeevolution dynamics across different spin models with various sizes and interaction types, the quench Hamiltonian can be properly normalized asĤ →Ĥ to have a unit bandwidth. This can be done [50] by renormalizing all coupling constants by the energy bandwidth of the interaction Hamiltonian, ∆Ĥ(q) 1 , defined as the energy gap between the maximum/minimum eigenvalues, i.e.
In the following, in line with the recent results of [41], we show that the SYK quench Hamiltonian can produce markedly distinct effects on the dynamics starting from different initial states and depending on the operator dynamics under investigation. To this end we will focus on two possible choices, a = x or a = z, forĤ a 0 (the case a = y being completely equivalent to the case a = x).
At first sight, the distinction between the two models may look very minor, since they can be related by a simple rotation of the external magnetic field. However, they have different dynamical features: for example, the global symmetries of the constant Hamiltonians (5) are broken by the quench term in different ways. Indeed, H x 0 preserves the total spin along the x axes and the quench term completely breaks this symmetry (the exactly same pattern happens for the a = y Hamiltonian). On the other hand,Ĥ z 0 preserves the total spin along the z axes and this symmetry is solely partially broken by the quench term down to the parity symmetry.

III. DYNAMICAL SPIN-SPIN CORRELATION FUNCTIONS
Time-resolved detection of propagating correlations in an interacting quantum many-body system after a quantum quench have been recently reported [12,16,61,62]. In addition, the ability offered by quantum gas microscopy [52,63,64] of single-site imaging in an optical lattice [65] allows the spatial resolution and sensitivity to reveal the real-time evolution of a many-body system at the single-particle level. Motivated by these progresses, we numerically investigate, by means of discrete timeevolution simulations, the dynamics of the evolved states after the quantum quench by studying the time evolution of the connected part of the two-point function or dynamical susceptibility, where the expectation value · · · is taken over the evolved ket, defined as For the sake of simplicity, we start discussing the case of quench Hamiltonian of SYK 2 type, while larger values q > 2 are considered in later sections. We recall that, when q = 2, the SYK Hamiltonian is non chaotic, thus providing an interesting example of disordered, all-to-all, integrable dynamics. The time evolution of χ a (t) for both a = x, z is reported in Fig. 1. In both cases, this quantity shows an initial rise, up to a maximum value χ a max , followed by a decrease toward 0. Interestingly, the peak is drastically larger in the a = x model, with the difference getting parametrically enlarged by increasing L. More in details, we have observed that χ x max grows quadratically with L, while χ z max grows linearly with L (See the inset in Fig 1). Such difference is not obvious to be explained, given that the two pre-quench models as well as their correlators just differ by a global rotation.
To better clarify these behaviors, it is instructive to inspect the single spin-spin correlation functions appearing in the sums of Eq. (8), with i, j indicating two lattice sites. Across all i, j combinations, the timeevolution of χ a ij (t) exhibits the same qualitative pattern as of χ a (t), which is an initial rise followed by a decay. This is because both σ a iσ a j and σ a i σ a j decay monotonically from one to zero but with different rates, such that χ a ij (t) and χ a (t) develop a peaked structure. To understand the differences in the peaks between the a = x and a = z models, we study the height of the peaks P a ij ≡ max t (χ a ij (t)) as a function of i and j. Naively, one would expect P a ij to be homogeneous and independent of i and j, because the initial Hamiltonian (5) treats all the spins equally and does not involve spin-spin couplings, and the SYK term (2) is all-to-all. Hence, one would expect that the models do not have any notion of neighboring sites, therefore P a ij would be independent of the lattice sites after averaging over the Gaussian random couplings.
Numerical results on P a ij are displayed in Fig. 2. In agreement with the motivations just explained above, P z ij is independent of i, j. On the other hand, P x ij exhibits a very intriguing behavior; it is highly dependent on i and j. In addition, the following patterns clearly emerge; the peaks are more pronounced for i and j being closer to each other, and they are further magnified when i and j approach to the center of the lattice.
Notice that these differences are not very easy to understand solely from the symmetry breaking perspective. For example, it is not obvious why P x ij is site-dependent. Moreover, in App. A we provide further evidences that the symmetry breaking pattern alone cannot explain the dynamical differences observed betweenĤ x 0 andĤ z 0 , considering hard-core bosonic versions of the SYK models. Although the bosonic variants show the same symmetry breaking patterns as in their fermionic counterparts, they exhibit clearly distinct dynamical properties.

IV. LOCALITY IN THE OPERATOR SPACE AND QUANTUM DYNAMICS
Here we demonstrate how the difference between P x ij and P z ij presented above can be explained by inspecting the temporal decay properties of the correlation functions σ a iσ a j and σ a i . Our approach is similar to the one discussed in [39] and based on the notion of operator size.
It is fairly simple to convert the Pauli matricesσ a i as well as the productσ a iσ a j of any two Pauli pairs into a product of Majorana fermionsγ i [66], by making use of the JW maps defined in (6). For this reason, we find more convenient to treat the Majorana fieldsγ i , rather than the Pauli operatorsσ a i , as the fundamental operators [39] to get a more transparent tracking of the time-evolution of the spin correlators.
Given the set of fundamental operators, one can introduce the notion of size of operators: an operatorÔ is said to have size k if it can be written as a product of k fundamental operators. More generally, an operator is k-local if, once rewritten as a sum of operators with definite size, the maximum size of its constituents is k. We will denote the locality of an operatorÔ in the superscript, i.e.Ô (k) . It has been emphasized in [39] that a proper identification of the set of fundamental operators, and the associated notion of operator size, allow a clear description of the dynamics of a strongly correlated quantum system. Let us start considering the a = z case, whose description is simpler. By making use of the JW map (6), the operatorsσ z i andσ z iσ z j can be re-written in terms of the Majorana variables as follows: thus showing that they are, respectively, of size 2 and 4 in the space of the Majorana operators. It is noticeable that the sizes of the operatorsσ z i andσ z iσ z j are independent of the lattice sites i and j.
The situation is more involved in the a = x case. Here, the spin operators, σ x i , show varying sizes which depend on the lattice site. More precisely, we find thatσ x i is an operator of size (2i − 1), i.e.
The application of the Clifford algebra (3) shows similarly that the product operatorσ x iσ x j is of size 2|j − i|. We summarize the size of the above spin operators as follows: The post-quench evolution of a generic spin operator O (k) of size k follows the Heisenberg equation of motion, which involves the commutator between the SYK 2 Hamiltonian and the operator itself. An advantage of the Majorana representation comes from the fact that the above commutator can be computed in a very simple way that manifestly preserves the operator size, by using Hence, the operator dynamics, in the space of Majorana operators, under the SYK 2 Hamiltonian is a kind of operator hopping: an operator O (k) , under time evolution, moves along the space of the operators of the same size k. The operator size does not grow in time. We underline that this characteristic is very distinct from the operator dynamics under SYK q Hamiltonians with q > 2, studied in [39] which will be discussed in the later section.
We are now at the position to explain the observed differences in the time evolution of the spin-spin correlation functions χ a ij (t), in terms of σ a iσ a j and σ a i , for the a = x and a = z cases. First of all, after averaging over the SYK coupling constants, the above correlators turn out to be return amplitudes, i.e. they compute the amplitude that at t > 0 the evolved operator takes exactly the same form of the initial operator. This is a consequence of the Gaussian averaging, as we show explicitly in App. B. [67] Secondly, when t sufficiently increases, the averaged correlators should decay to 0, since the initial operator is of measure 0 in the entire space of size k operators.
Given these facts, we argue that the vanishing of the spin-spin correlator happens faster when the size k of the corresponding operator is larger, if k ≤ L, and when the size k is smaller, if k > L. Indeed, as k becomes bigger, the dimension of the space of all size k operators increases monotonically until it reaches k = L, i.e. exactly half of the maximum possible size for an operator, and then decreases up to k = 2L. In short,Ô (k) (t) has a larger space to explore as |k − L| decreases. This causes the corresponding return amplitude Ô(k) (t) to vanish faster, since the chance to go back to the original opera-torÔ (k) will get reduced in a shorter amount of time. To numerically confirm the above argument, we have studied the time evolution of various spin operators O k (t) in the a = x model. We have checked (not shown) that the evolution of the return amplitude turns out to be solely characterized by the operator size k. For instance, the decay patterns of σ x iσ x i+ (t) are essentially identical across all 1 ≤ i ≤ L − , while they strongly depend on the value of . Furthermore, Fig. 3 exhibits the pattern that the operators with higher k ≤ L decay faster to zero, in agreement with our prediction. We also confirmed (not shown) that a greater value of k > L leads to a slower decay rate, since the dimension (16) of the operator space shrinks as the size k > L gets bigger.
The relation between the size of the operator and the decay rate suffices to explain how the connected part of the spin-spin correlation functions (10) evolve in time.
In the a = z model, the first term, which is a correlator of size 4, vanishes slightly slower than the second term, which is the product of two size 2 operators. The difference between the two terms is small and independent of the choice of i and j, in agreement with Fig. 2.
On the contrary, in the a = x model, χ x ij (t) is the difference between a correlator of size 2(j − i) and the product of two correlators of respective sizes (2j − 1) and (2i − 1). The peak height P x ij is therefore maximized when i and j are adjacent, i.e. j = i ± 1, and i ∼ j ∼ L/2. Furthermore, since the product of correlators diminishes much more rapidly than the single correlator term, P x ij is generically larger than P z ij , thereby contributing to the bigger bump of χ x (t), as visualized in Fig 1. It is remarkable that, although the initial Hamiltonians, (5) and (2), do not explicitly involve any couplings between neighboring spins, the notion of locality on the spin lattice emerges from the Majorana representation of the spin operators, due to the operator hopping dynamics under the SYK 2 quench and by averaging over the Gaussian couplings.
An important point to stress is that the fast decay of the averaged correlator is a consequence of the property that the interactions are all-to-all, which remains true for any SYK q models irrespectively of the value of q. In the current case of the SYK 2 model, thanks to all-toall interactions, the operator O k (t) can sweep across the full space of size k operators, leading to the fast decay of the return amplitudes. More generally, the condition for return amplitudes to decay is that the space that O k can explore under the unitary time evolution is large enough, which can be satisfied whenever the operator O k sits on a connected (hyper)-graph [68,69].

V. FROM OPERATOR HOPPING TO OPERATOR GROWTH
In this section, we inspect how the evolution of O k (t) after a quantum quench protocol can discriminate between two different types of operator dynamics: operator hopping, discussed in section IV, versus operator growth, discussed, for example, in [39]. To this end, we turn now to the case in which the quench Hamiltonian is the SYK q Hamiltonian with q > 2 and we focus on the case q = 4.
Contrary to the q = 2 algebra, (15), which preserves the operator size, the commutation relation at q = 4, realizes an example of the operator growth dynamics; the evolution of a fundamental operatorγ i is not confined in the space of size 1 operators, but its degree of locality continues to increase with time until it saturates L. Notice that the notion of the operator growth is a hallmark of the early-time quantum chaos. The underlying idea is that the operator growth dynamics can develop simple (even fundamental) operators into more complex and extended ones, thereby realizing a quantum analogue of the well-known "butterfly effect" [36,70]. We recall that the decay rate of O k (t) is governed by the dimension of the operator space reached out from an initial k-local observable O k through the time evolution. This implies that the decay rate of O k (t) should be less sensitive to the value of k under the operator growth, compared to the hopping dynamics. As explained before, the operator hopping does not change the operator size, therefore a set of possible trajectories is also confined in the space with the fixed dimension (16). Under operator growth, however, the degree of locality of time-evolved operator quickly saturates to L regardless the initial value of k. Since the dimension (16) of the operator space with a fixed size k is maximized at k = L, the Gaussian averaged correlator O k (t) should vanish in a shorter time span, as well as being less sensitive to k.
To validate the above reasoning, we have computed the averaged correlators under the SYK 4 quench and plotted them into Fig. 4. By contrasting Fig. 4 with Fig. 3, which displays the averaged correlators under the SYK 2 quench, we recognize that the decay curves start to overlap, becoming indistinguishable from each other, at a smaller values of k in the case of SYK 4 quench. This is consistent with the argument just presented.
To better distinguish the decay pattern of the correlators, O k (t) , under the SYK 4 vs SYK 2 quench dynamics, one can directly look at the decay rate, focusing on its maximum value, max t (D k (t)), the highest speed of correlation decay over time. As we compare the maximum decay rate across different SYK q models and various sizes L of the spin lattice, it is also convenient to normalize the degree k of locality as k rel ≡ k/L, and the maximum decay rate, max t (D k (t)), as In Fig. 5, we display the normalized maximum value (19) of the decay rate as a function of k rel for different SYK q models. We have found that these plots are insensitive to the lattice size L, only depending on the value of q. We observe that R(k rel ) saturates at a smaller value of k rel under the SYK 4 quench than in the SYK 2 case. This is in agreement with the expected difference between operator growth and hopping: since the hopping dynamics preserves the size of operators, the associated decay rate must be more sensitive as explained above.
For comparison, we also depict the values of R(k rel ) for the SYK 6 model, for which we observe that a saturating k rel is even smaller than the one for the SYK 4 dynamics. The smaller saturating k rel for the SYK 6 model is easy to explain, since the sextic SYK Hamiltonian is faster in increasing the operator size than the quartic SYK Hamiltonian. Also, it is worth to notice that the difference between the SYK 4 and SYK 6 models is much less pronounced than that between the SYK 4 and SYK 2 models.
The above numerical results strongly suggest that the maximum value R(k rel ) of the decay rate of the averaged spin correlators can effectively distinguish the dynamics of operator growth versus hopping. However, from its definition, R(k rel ) is not a convenient quantity for direct experimental access. Indeed, to distinguish between hopping and growth, one still should in principle measure the correlators O k (t) for all the possible values of k, from k = 1 to k = L. Considering large spin-chains, this would imply the necessity of measuring a very large number of averaged correlators, thus making the requirements in terms of number of measurements very demanding.
Given the considerations, we ask whether there is a simple correlation function, which alone can discriminate between operator hopping and growth. As it can be inferred from Fig. 5, the averaged correlators at large degree of locality are not very useful since both hopping and growth are in the saturation regime at large k and the resulting dynamics will be indistinguishable. Instead, the difference between two kinds of dynamics must be evident at small degree of locality, such as k = 1 or k = 2. It is therefore interesting to study the correlator χ x 1L , involving the difference between an operator of size 2(L − 1), whose dynamics is identical to the operator of initial size 2, and the product of two operators of effective size 1. Hence, we expect to observe the maximum difference between hopping and growth using this probe. Fig. 6 contrasts the time evolution of the correlator χ x 1L (t) in the q = 2 and the q = 4 models. The ensemble averages are taken over 100 sets of random couplings. It is immediate to notice that the behaviors are rather different. In SYK 4 model, the correlation function exhibits a clear maximum, higher than the noise amplitude at late times by a few order of magnitudes; i.e. the peak is by far larger than the late time fluctuations, still present after the average over the Gaussian couplings, although highly suppressed. On the other hand, in SYK 2 model, the height of the initial peak is of the same order of the fluctuation amplitude at late times; the SYK 2 fluctuations are more pronounced than the SYK4 noises, since in the case of operator hopping the evolution trajectory of the given operator is confined to the relatively smaller space of operators of size 2(L − 1).
We have extensively checked that these marked differences are robust by increasing the number of ensemble realizations over which the Gaussian average is performed. Moreover, they become parametrically more evident by increasing lattice size L. These observations suggest that the evolution of χ x 1L (t) in time can be considered as a useful diagnostics of operator growth versus hopping, being able to identify the nature of the operator size dynamics.

VI. CONCLUSIONS
In this paper, we have investigated whether quantum chaos, and specifically operator growth, can be revealed by performing quantum quench protocols on systems defined over spin lattices.
By using the celebrated SYK q model as the quench Hamiltonian with all-to-all interactions, we have established that the time evolution of the spin-spin correlation functions can be used as a probe of operator growth.
Mapping the spin variables to the Majorana fields, which here constitute the fundamental operators, the as-sociated size of different spin-spin correlation functions have been identified.
We have demonstrated how the decay rate of the averaged spin correlators O k (t) is controlled by their initial sizes. Moreover, the relative decay rate, R(k rel ), can distinguish operator growth from operator hopping; the former being a hallmark of early-time quantum chaos while the latter shows a rather trivial dynamics in operator space. Finally, we have discussed that the difference in operator dynamics strongly affects the particular averaged spin correlator χ x 1L (t): the amplitude of the latetime fluctuations is comparable to the height of an initial peak under operator hopping, but significantly suppressed under operator growth. Such marked distinction can be useful to detect quantum chaos in strongly correlated systems, especially given the fact that the evolution of the spin correlators is experimentally traceable by using state-of-the-art imaging techniques and quantum gas microscopy [52]. We believe that a precise analysis of the connections between the peaked structure that emerges in the dynamics of spin-spin correlation functions χ x 1L and the OTOCs is worth of future investigations. In particular it would be extremely interesting to understand whether a quantitative evaluation of the Lyapunov exponents can be extracted from χ x 1L . Also, it would be intriguing to analyze the evolution of averaged spin correlators in other chaotic systems, e.g. random circuits [71,72], from the perspective of operator dynamics. ACKNOWLEDGEMENT We thank F. Haehl, T. Nosaka, V. Rosenhaus for related discussions. JK acknowledges the support from the NSF grant PHY-1911298. DR is supported by a KIAS Individual Grant PG059602 at Korea Institute for Advanced Study. Most numerical computations were done thanks to the computing resources provided by the KIAS Center for Advanced Computation (Abacus System) and the Institute for Advanced Study. Some of the numerical results have been obtained by making use of the Wolfram Mathematica package QuantumManyBody, freely available on GitHub. The rendering of the plots has been realized using the the Wolfram Mathematica package MaTeX, freely available on the Wolfram Library Archive.

Appendix A: Bosonic models
In this appendix, we investigate the real hard-core boson variants of the SYK models. We study how their dynamics differ from their fermionic counterparts. We recall that the complex hard-core boson SYK model was introduced and studied in [73,74]. Despite the same symmetry breaking pattern as in the fermionic models, the bosonic dynamics turns out to be completely different. This difference can be understood again in terms of the size of the operators involved.
To begin with, we define the real hard-core boson operators,χ i , with i = 1, . . . , 2L, as the operators satisfying the following algebra: The hard-core bosonic operators can be mapped to operators defined on a spin lattice of length L via the following JW-like mapχ Similarly to the fermionic case, one can define the following SYK-like Hamiltonian for the operators {χ i }, where the coupling constants J i<j<k<l are sampled from the same Gaussian distribution as in (4). The factor s(i 1 , . . . , i q ) reads and it must be introduced to ensure thatĤ q 1, B is Hermitian under the bosonic algebra (A1).
Here, for sake of brevity, we focus on the q = 4 model. We contrast the temporal evolution of averaged correlators under the quench protocol (1) with the Hamiltonian either being the bosonic SYK 4 (A3) or its fermionic counterpart (2). Notice that the symmetry breaking pattern in the bosonic models is the same as in the fermionic systems; in the a = x model the spin symmetry along the x axes is fully broken by the quench term, while in the a = z model the spin symmetry along the z axes is partially broken to the chiral symmetry.
Let us start by considering the connected part of the two-point function, defined in (8). As we have shown in the main text, the time evolution of χ a (t) as well as individual summands χ a ij (t) is controlled by the initial size of the operators in σ a i σ a j and σ a i . A rapid inspection, by making use of the JW-like map (A2), shows that in the a = x case, the bosonic σ x i σ x j operators have size 2, while the bosonic σ x i operators have size 1. For the a = z case, however, the operator size is exactly the same as in the fermionic case: the bosonic σ z i σ z j operators have size 4, while the bosonic σ z i operators have size 2. It is important to note that, in the bosonic case, all the operator sizes are independent of i and j.
The correlators χ a (t), for bosonic/fermionic SYK 4 models are compared in Fig. 7. We clearly observe that, for a = z, the evolution curves for both bosonic and fermionic models are barely distinguishable; this shows a perfect agreement with the expectation based on the size of the operators involved, since when a = z the size is independent of the bosonic/fermionic nature of fundamental variables. On the other hand, the situation is completely different for a = x case; the characteristics for the fermionic and bosonic models are rather opposite, where the fermionic/bosonic models exhibit a huge/tiny value of the peak, respectively. Such distinction is again in perfect agreement with the expectation based on the operator size, because the fermionic a = x model is the only case for which we have very large differences in the operator sizes between σ x i σ x j and σ x i . It is also interesting to observe that the fermionic a = x model is the only model for which the height of the peak, χ x max (L), scales quadratically as a function of L (Fig 1). We also display the height P a ij of the peak for all bosonic/fermionic SYK 4 variants in Figs 8 and 9.
Appendix B: The average over the Gaussian coupling Here we present an analytic argument, based on a short time expansion, to establish a relation between the degree of operator locality and the decay rate of correlation functions. As a product, we show that the numerical results presented in Fig. 3 and Fig. 4 indicate that the correlators O k (t) are computing a return amplitude, as the effect of the average over the Gaussian couplings.
Thanks to the Majorana representation, (11) and (12), of the spin variables, it is easy to consider the timeevolution of spin operators. One needs to resolve the commutator [Ĥ 1 , O k (t)] that appears in the Heisenberg between the SYK q coupling (2) and a Majorana operator. For notational convenience, we assume the quench HamiltonianĤ 1 to be the standard SYK q Hamiltonian H (q) 1 without normalizing the bandwidth. Instead, the bandwidth normalization has been directly taken into account as a proper rescaling of the time variable t in (B2).
Since the post-quench Hamiltonian is static, the solution O k (t) of the Heisenberg equation can be written where O k ≡ O k (0). Let us replace all the iterated commutators in (B2) by using the commutation relation (B1). Averaging over all Gaussian coupling constants, J i1···iq , leads to huge simplification. All odd-order moments of Gaussian couplings drop out, implying O k (t) is an even polynomial in t. The remaining even-order moments can be handled with Wick's theorem, enforcing the sum of surviving operators on the right-hand side to be proportional to O k (0). This shows that the averaged correlator O k (t) is indeed a return amplitude, i.e. it computes the amplitude that O k (t) returns to itself, O k (0). Also, it manifests that the early-time scale, which determines the regime of validity of the expansion (B2), is inversely proportional to the standard deviation J of the SYK coupling constants.
We specialize to a simple case with q = 2 and k = 1 as an illustrative example. The commutator action [Ĥ (2) J ab J cd J ef J gh δ bn δ da δ f c δ he ·γ g , (B4) and so on. The right-hand side expressions are further simplified by taking the coupling constant average,: J ab J cd J ef J gh = J ab J cd · J ef J gh + J ab J ef · J cd J gh (B6) Just like (B6), the higher-order moments can be obtained by Wick contractions. Inserting them back, the summations appearing in iterative commutators, e.g., (B3)-(B4), can be easily carried out. [Ĥ 1 , [Ĥ The above expressions can be simply generalized to the followings that hold for spin operators with k-locality k ≥ 1: [Ĥ 1 , [Ĥ [Ĥ 1 , [Ĥ and they constitute the solution (B2) of Heisenberg equation. In Fig. 10 we compare some representative correlation functions, O k (t) , against their analytical quartic predictions based on (B3) and (B4). From the figure we clearly see that the short time agreement is extremely good, thus confirming that the numerical average over the Gaussian couplings, performed to get the numerical plots of O k (t) , is sufficiently precise to convince that the numerical correlators O k (t) are effectively computing the return amplitudes of the operators O k . The extension to the q > 2 case is also straightforward in principle. The main impediment for an actual q > 2 calculation is the rapid growth of the number of operators that appear in iterative applications of the commutator action (B1). Since the right-hand side of the commutator (B1) contains C(2L − 1, q − 1) distinct terms, the number of operators appearing in the n-times commutator action on O k grows very quickly, for q > 2, as observed in the following L = 5 example: × 2 5 (3k + 1)L 5 − 2 4 (15k 2 + 38k)L 4 + 2 3 (36k 3 + 122k 2 + 207k − 99)L 3 − 2 2 (48k 4 + 168k 3 + 543k 2 + 352k − 510)L 2 + 2 3 (9k 5 + 21k 4 + 168k 3 + 88k 2 + 156k − 304)L − 12k 6 − 336k 4 − 624k 2 + 864 , which, compared with the analogous formulas for the SYK 2 model, (B8) and (B9), suggest how involved can be the dynamics of the model for q > 2. This is, of course, a manifestation of the fact that for q > 2 the dynamics turns from operator hopping to operator growth, and this change is reflected in the cumbersome formulas, (B10) and (B11). Anyway, we can make use of (B10) and (B11) to compare the short time decay of the numerical correlators in the SYK 4 model, against their analytical predictions based on (B10) and (B11). The results are reported in Fig. 11. Also in this case, an excellent agreement is found at short times. In particular, it is possible to note that already at the 4th perturbative order the averaged correlators at large values of k do overlap with each other.