Scrambling and quantum chaos indicators from long-time properties of operator distributions

Scrambling refers to the spreading of initially localized information to the rest of the degrees of freedom in a many-body system Hayden and Preskill (2007); Sekino and Susskind (2008); Lashkari et al. (2013); Hosur et al. (2016); Swingle et al. (2016). It plays an important role in describing diverse phenomena such as closed-system thermalization Lewis-Swan et al. (2019); Claeys and Lamacraft (2021), dynamical phase transitions Heyl et al. (2018); Dağ et al. (2019); Wang and Pérez-Bernal (2019), sampling hardness in random quantum circuits Boixo et al. (2018); Bouland et al. (2019), and information retrieval in black holes Sekino and Susskind (2008); Hayden and Preskill (2007); Basko et al. (2006). One of the most prominent quantifiers of scrambling is the out-of-time-ordered correlator (OTOC), which is a four-point correlation function of the form $\langle W^{\dagger}(t)V^{\dagger}(0)W(t)V(0)\rangle$, together with the closely related square commutator $\langle|[W(t),V(0)]|^{2}\rangle$ Maldacena et al. (2016); Cotler et al. (2017); Swingle (2018).

OTOCs and scrambling have become important actors in the dynamical characterization of chaos in many-body quantum systems Rozenbaum et al. (2017); García-Mata et al. (2018, 2022). Quantum chaos is typically defined in terms of kinematic features like statistical properties of energy spectra and their connection to random matrix theory Bohigas et al. (1984); Haake (1991); Wimberger (2014); Kos et al. (2018). However, a unifying dynamical description of chaos in general quantum systems remains an outstanding challenge. Many studies of scrambling in quantum systems have focused on understanding its early-time behavior, particularly through the decay of OTOCs, and have sought to define a quantum analogue of the Lyapunov exponent. Nevertheless, there are cases of quantum systems whose kinematic properties follow random matrix theory predictions (and are thus ‘quantum chaotic’), for which the OTOC does not decay exponentially Fortes et al. (2019); Roberts and Stanford (2015); Aleiner et al. (2016); Shenker and Stanford (2014, 2015). Conversely, some quantum systems with classically integrable counterparts showing unstable fixed points can lead to exponential decay of OTOCs Kidd et al. (2021); Xu et al. (2020); Pilatowsky-Cameo et al. (2020).

More generally, scrambling is a complex process which can be described by the way an initially simple operator evolves into a complicated superposition of configurations belonging to an exponentially large operator space Zhuang et al. (2019); Parker et al. (2019). As such, it is bound to require methods to characterize it that go beyond the short-time behavior of a single correlation function. Indeed, recent studies on operator growth have developed a more thorough characterization of scrambling by analyzing the dynamics of operators in operator space, most notably using the Krylov representation Parker et al. (2019); Noh (2021). Moreover, some studies have found important links between quantum chaos and the long-time behavior of the OTOC Fortes et al. (2019), rather than its initial decay. In this context it is useful to further probe the connection between the long-time properties of evolving operators and quantum chaos, and to explore other objects of interest beside the OTOC which allow us to construct diagnostics of scrambling in the long-time regime. An additional motivation is the inherent complexity of accessing OTOCs experimentally: even with the extraordinary control and isolation capabilities found in state-of-the-art quantum simulation experiments Landsman et al. (2019); Blok et al. (2021); Mi et al. (2021), accessing OTOCs require costly resources such as use of auxiliary systems or time-reversal operations Gärttner et al. (2017); Li et al. (2017); Swingle et al. (2016); Yao et al. (2016), among others, and thus require considerably more effort when compared to the usual quench-dynamics experiments (with notable exceptions, see Vermersch et al. (2019); Joshi et al. (2020); Blocher et al. (2022a)).

In this work we purposely steer away from OTOCs and focus on studying scrambling in quantum systems by analyzing directly the dynamical properties of a suitably-defined probability distribution $\{P_{k}(t)\}$ over a coarse-grained operator basis. This distribution can be defined for arbitrary quantum systems, and here we focus on the cases of models of many spin-$1/2$ particles (including quantum circuits on qubits) and models of collective spins, which are effectively described by a single large spin $J$ Lerose and Pappalardi (2020). The dynamical transition from ‘simple’ to ‘complex’ operators is then encoded in different properties of the distribution such as its mean, variance, and (de)localization, which signify the growth and spreading of operators over the degrees of freedom of the system. We apply this framework to paradigmatic models of quantum chaos such as the ‘tilted’-field Ising model Karthik et al. (2007) and the quantum kicked top Haake et al. (1987), and show that analysis of the long-time averages of the distribution properties and their temporal fluctuations are good indicators of the onset of chaos in these systems. Furthermore, we discuss how some of these features can distinguish different properties of the nonergodic regimes, and show that both models can show very similar behavior in this picture even though their physical properties are quite different. We also apply these tools to the study of random quantum circuits Fisher et al. (2022); Harrow and Low (2009); Oliveira et al. (2007); Arute et al. (2019) with a tunable number of $T$-gates, and analyze how the properties of scrambling change as just a small fraction of non-Clifford gates are included in the dynamics. Finally the connection between averages of OTOCs and the moments of the distribution $\{P_{k}(t)\}$ is analyzed and we discuss the number of OTOC measurements needed to reconstruct the different measures we study.

Our work extends previous studies that have focused on the properties of the operator distribution in the study of scrambling Hosur et al. (2016); Qi et al. (2019); Roberts et al. (2018); Schuster and Yao (2022). Notably, these also include NMR experiments which routinely analyze the size of active clusters of spins, a quantity that is closely related to the mean operator size in the Heisenberg picture Álvarez and Suter (2010); Domínguez and Álvarez (2021). It also complements the approach of Ref. Fortes et al. (2019), which studied the connection between quantum chaos and the long-time properties of OTOCs by considering the properties of the operator distribution directly in a similar regime.

The rest of the work is organized as follows. In Sec. II we discuss scrambling for general quantum systems, and define the coarse-grained probability distribution as the object that characterizes it. In Sec. III the tilted-field Ising model is studied via numerical simulations, and we analyze the evolution of the probability distribution in its integrable and chaotic regimes. In Sec. IV we study the quantum kicked top, which is a collective spin model with a well-defined classical limit, and use it to relate the long-time properties of scrambling with the chaotic properties of the model. Then, in Sec. V we study a model of random Clifford circuits perturbed by a tunable number of $T$-gates, and study how the properties of the operator distribution change as the number of non-Clifford gates increases. Section VI discusses the connection between the operator distribution $\{P_{k}(t)\}$ and averages of OTOCs. Finally, we give an outlook and discuss potential future work in Sec. VII.

Consider a quantum system on a finite-dimensional Hilbert space of dimension $d$, with evolution from $t=0$ to some arbitrary time $t$ described by an unitary operator $\hat{U}(t)$. We will focus on the dynamics of a generic hermitian operator $\hat{O}$, which can be expressed as

where $\{\hat{\Lambda}_{j}\}$, $j=0,1,\ldots,d^{2}-1\equiv D$ is an operator basis which we take to be orthonormal, i.e. $\operatorname{tr}\left(\hat{\Lambda}_{i}^{\dagger}\hat{\Lambda}_{j}\right)=\delta_{ij}$. We also take $\hat{\Lambda}_{0}=\mathbb{I}/\sqrt{d}$ throughout, and consider $\operatorname{tr}\hat{O}=0$. The scalar coefficients $\{f[\hat{\Lambda}_{j},\hat{O}(t)]\}$ describe the dynamics of $\hat{O}(t)$ in the chosen basis, and at all times satisfy the normalization condition

Equation (2) allows us to treat the set of squared coefficient amplitudes as a probability distribution (after normalization) ¹¹1 The coefficients $\{f[\hat{\Lambda}_{j},\hat{O}(t)]\}$ can also be regarded as the complex amplitudes of an “operator wavefunction” Roberts et al. (2018), and thus Eqn. (2) represents the normalization of such wavefunction.. Additionally, in many situations of interest the operator basis admits a natural ordering related to some notion of complexity of the operator, which we schematically depict in Fig. 1. Qualitatively, we consider this ordering to induce a natural grouping of the basis, which has the form

where the complexity of operators is considered to grow as the index $k=1,2,\ldots,k_{\mathrm{max}}$ increases and $\sum_{k}\mathrm{dim}(C_{k})=D$. For instance, this structure could correspond to a multibody Pauli basis where $k$ corresponds to the weight or size of each operator, as we will discuss in Sect. II.1. Leveraging this structure we define the following coarse-grained probability distribution for the operator $\hat{O}(t)$

which naturally satisfies $\sum_{k}P_{k}(t)=1$. We are interested in situations where the initial operator $\hat{O}$ is included within one of the groups of low complexity (say $k=1$), and our goal is to characterize the scrambling process in the system via the evolution of the distribution $P_{k}(t)$ over time, as depicted in Fig. 1. To characterize the distribution we will consider the first two cumulants of the distribution,

as well as its inverse participation ratio (IPR),

Here $\eta_{\mathrm{IPR}}\simeq 1$ indicates a very localized distribution (“low participation”), while $\eta_{\mathrm{IPR}}\simeq 1/k_{\mathrm{max}}$ shows a delocalized distribution (“high participation”). The IPR has been extensively used to measure delocalization of eigenstates in quantum chaos Sieberer et al. (2019), and a similar object has been used in the context of resource theories of Magic Leone et al. (2022). Here we will use it to assess the delocalization of the operator distribution in the scrambling process. We point out that higher-order cumulants of the full distribution in the Krylov basis have been studied recently in Bhattacharjee et al. (2022).

In Appendix A we show that for Haar-random evolution the operator $\hat{O}(t)$ will be uniformly spread in any operator basis. The resulting distribution thus has the form

and it is straightforward to compute the indicators introduced once the sets $\{C_{k}\}$ are defined. This limiting case will be helpful in order to study the onset of chaos as dictated by the randomization of the evolution. In the following subsections we investigate two cases of interest: systems of $N$ spin-$\frac{1}{2}$ particles or qubits, and collective spin systems described by a single large collective spin $J=N/2$.

We begin by studying the properties of the operator distribution $\{P_{k}(t)\}$ in the different regimes of the Ising model, a standard paradigm in the study of many-body quantum systems Edwards et al. (2010); Bernien et al. (2017); Zhang et al. (2017). This model describes a set of $N$ spin-$\frac{1}{2}$ particles interacting in $1$D via nearest-neighbor interactions and in the presence of an external magnetic field with a transverse and a longitudinal component. The Hamiltonian can be written as

where we take $0\leq\theta\leq\frac{\pi}{2}$. Here $\sigma_{n}^{\alpha}$ are the usual Pauli operators on site $n$ with $\alpha=x,y,z$. For $\theta=0$ Eq. (18) is the transverse-field Ising model (TIM), whose equilibrium and nonequilibrium properties have been studied extensively Sachdev (2011). This model is integrable since it can be mapped to a noninteracting system of fermions via the Jordan-Wigner transformation Jordan and Wigner (1993); Sachdev (2011). The case of pure longitudinal field $\theta=\frac{\pi}{2}$ is diagonal in the computational basis and thus also trivially integrable. For other values of $\theta$ (and generic choices of $B/J$), this “tilted-field” Ising model is quantum chaotic as revealed by several of the usual metrics, which have been studied in previous works Prosen (2002); Mejía-Monasterio et al. (2005); Karthik et al. (2007). For completeness, we revisit some of those results here. First, the eigenenergies of $H_{\mathrm{Ising}}$ in this chaotic regime display level repulsion as predicted by random matrix theory, and this can be quantified by the average adjacent spacing ratio $\overline{r}$, the details of which we present in Appendix B. We plot a normalized version of this quantity, $\overline{r}_{\mathrm{norm}}$, as a function of $\theta$ in Fig. 2 (a). Values close to $1$ indicate agreement with the Gaussian Orthogonal Ensemble (GOE) predictions and thus quantum chaotic behavior, while deviations towards $0$ indicate uncorrelated level statistics typical of integrable systems. As an additional metric, we also consider the average entanglement entropy of the excited states of $H_{\mathrm{Ising}}$ for an equal bipartition of the chain Karthik et al. (2007). This is defined as

where $\rho^{(i)}_{\frac{N}{2}}$ is the state resulting from tracing out half of the particles from $|\phi_{i}\rangle\langle\phi_{i}|$, $\Ket{\phi_{i}}$ is an eigenstate of $H_{\mathrm{Ising}}$, and the sum is carried out over the bulk the spectrum (i.e. avoiding the ground and low energy states of $H_{\mathrm{Ising}}$ and $-H_{\mathrm{Ising}}$). Regimes of maximum average bipartite entanglement are then associated with quantum chaos, as can be verified from the results shown in Fig. 2 (b).

We now study the scrambling process in the different regimes of the Ising model by analyzing the probability distribution $P_{k}(t)$ defined in Eq. (11). In Fig. 2 (c) and (d) we display exact numerical results of the time-dependent distribution corresponding to a chain of $N=6$ particles where the initial operator sits either (c) in middle of the chain $\hat{O}(0)=\hat{\sigma}_{y}^{(N/2)}/\sqrt{2}$, or (d) at the edge $\hat{O}(0)=\hat{\sigma}_{y}^{(1)}/\sqrt{2}$. In both cases, the initial distribution is initially concentrated in $P_{1}(0)=1$ and then evolves in time displaying different features depending on the system parameters. All cases with $0\leq\theta\lesssim\frac{\pi}{4}$ show a rapid decrease of the initial component $P_{1}$ as the operator spreads into a superposition of larger-size configurations. Crucially, however, the chaotic case $\theta=\pi/6$ shows a fast equilibration to the values corresponding to the random distribution shown as dashed lines, cf. Eq. (8). As $\theta\rightarrow 0$ and the model becomes integrable, the distribution shows further oscillations and fails to equilibrate completely in the timescale shown. The deviations from ergodicity are enhanced when the initial operator sits at the edge of the chain, a situation in which the initial configuration has less options to equilibrate to since the site has a single neighbor instead of two due to the open boundary conditions.

The other integrable regime, occurring at $\theta=\pi/2$ but already noticeable for $\theta=\pi/3$, corresponds to a very different type of evolution. In the diagonal case the external field commutes with the interaction, and thus a single site operator like $\hat{\sigma}_{y}^{(l)}$ evolves to only two- and three- site operators, independent of the length of the chain. The distribution then shows very little spreading as most of the elements are never populated. The situation remains roughly the same even in the presence of a small transverse field, as can be seen in the cases corresponding to $\theta=0.45\pi$ shown in Fig. 2.

While the operator probability distribution $\{P_{k}(t)\}$ already reduces the description of observable evolution from the exponentially large basis to a set of only $N$ numbers, it is still helpful to analyze measures which describe particular aspects of the distribution at each time. We thus turn to study the quantities introduced in Sec. II, namely the mean $\mu(t)$, variance $\sigma(t)$, and IPR $\eta_{\mathrm{IPR}}(t)$ of the distribution. Figure 3 (a) shows the evolution of these three quantities for different values of $\theta$ for the initial operator located in the middle of the $N=6$ particle chain (for completeness, the case where the initial operator sits at the edge is shown in Appendix C). The evolution of the mean $\mu(t)$ shows a increase from $\mu(0)=1$ towards $\sim N$ with different features depending on the value of $\theta$. The most chaotic case, $\theta=\pi/6$, grows until reaching the random value $3N/4$ (see Table 1), while most other cases show oscillations. Note that the TIM case ($\theta=0$) grows beyond the random value, meaning that during certain periods, the integrable model leads to mean operator sizes which are larger than those found for Haar-random evolution. Finally, as the other integrable limit is approached for $\theta\sim\pi/2$, the distribution stops shifting beyond $N=3$, as discussed before, and shows indefinite large-amplitude oscillations.

For the evolution of the variance of the distribution $\sigma^{2}(t)$, we observe that cases far from the diagonal case (i.e. $\theta\ll\pi/2$) show an initial increase from $\sigma^{2}(0)=0$ which shoots up significantly above the Haar-random prediction $\sim 3N/16$. Then, most cases equilibrate to that value, albeit in different timescales. Interestingly, as $\theta$ reaches $\pi/3$ the variances become consistently larger than the Haar prediction. This indicates a nontrivial behavior in which not necessarily the most chaotic evolution leads to the largest width of the operator distribution. A similar trend is observed in the IPR of the distribution, $\eta_{\mathrm{IPR}}(t)$, shown also in Fig. 3 (a). We point out that $\mu(t)$ and $\sigma(t)$ as a function of time for the chaotic case $\theta=\pi/6$ show a remarkable similarity to the ones observed for the SYK model by Roberts et al. in Ref. Roberts et al. (2018).

On the other hand, we observe that the studied properties of the distribution have a harder time distinguishing the integrable model at $\theta\rightarrow 0$ from the ergodic case. For instance, the time-averaged values of the mean, variance, and IPR stay very close to the random predictions as $\theta\rightarrow 0$, indicating that the integrable transverse Ising model leads to significant, random-like operator spreading at long times. This is true for most choices of the initial operator; however, we find that choosing $\hat{O}(0)$ at the edge leads to somewhat different features, particularly in the $\theta\simeq 0$ regime. For this regime and choice of initial operator, we see that the time-averaged mean drops and the time-averaged variance rises, signaling a clear deviation from ergodicity. Interestingly, while the behavior of the mean is similar to the other integrable limit, the case of the variance is opposite – the width of the distribution increases above the chaotic case in the transverse-field regime.

More generally we observe that the integrability of the model at $\theta=0$ consistently leads to increased long-time temporal fluctuations in all quantities, independent of the choice of the initial operator. Enhanced oscillations in the mean, variance, and IPR are seen in both integrable regimes as compared to the chaotic case. We thus find that these temporal fluctuations can be used to distinguish chaotic and nonchaotic regimes of the model. These results are aligned with the findings of Ref. Fortes et al. (2019), who showed that the long-time behavior (particularly the properties of the frequency spectra) of OTOCs serves as a good indicator to distinguish quantum chaos for integrability in various models, including the Ising model considered here. In Sec. VI we will further discuss the connection between the quantities studied here and OTOCs.

Finally, we point out that the results shown here for a fixed system size of $N=6$ are representative of other cases, which we show in Appendix C. In particular we show in Fig. 9 that cases with $N=5,6,7$ behave very similarly and essentially coincide (when properly normalized) in the chaotic regime. We also observe that the magnitude of temporal fluctuations decay with increasing system size, a behavior typical of ergodic systems.

We now turn our attention to the study of operator spreading and scrambling in collective spin systems, which we introduced in Sec. II.2. We will consider the dynamics of a quantum kicked top (QKT), a paradigmatic model of quantum chaos first introduced by Haake, Kuś, and Scharf Haake et al. (1987), which has been the subject of many theoretical Schack et al. (1994); Sieberer et al. (2019); Ghose et al. (2008) and experimental Chaudhury et al. (2009); Neill et al. (2016); Lysne et al. (2020) studies. The QKT time evolution operator of interest for this study can be written as

where $\mu=x,y,z$ and the total angular momentum is $J=2N$. The model in Eq. (22) is constructed in such a way as to avoid parity and time-reversal symmetry, which are present in the original QKT of Ref. Haake et al. (1987); see also Kus et al. (1988) for the study of similar models. Each of the unitaries $\hat{U}_{\mu}$ can be regarded as generated by a “twisting and turning” collective Hamiltonian Sieberer et al. (2019), composed of a rotation term $\hat{J}_{\mu}$ and a twisting or interaction term $\hat{J}^{2}_{\mu}=\sum_{ij}\hat{\sigma}_{i}^{(\mu)}\hat{\sigma}_{i}^{(\mu)}/4$. Since the symmetric subspace dimension scales linearly with the number of particles $d=N+1$, large values of $N$ can be accessed numerically in these types of systems.

The quantum chaotic properties of the QKT can be fully understood by studying the associated classical kicked top, which can be recovered as the mean-field limit of the map generated by Eq. (22) Haake et al. (1987). The resulting classical area-preserving map acts on a spherical phase space whose coordinates are $\mathbf{R}\equiv(X,Y,Z)=\lim_{J\to\infty}\langle\hat{\mathbf{J}}\rangle/J$. In Fig. 4 (a) we show the Poincaré sections corresponding to this map, where we have chosen the system parameters to be $\alpha_{x}=1.7$, $\alpha_{y}=1$, $\alpha_{z}=0.8$ and $\gamma_{x}=0.85\gamma$, $\gamma_{y}=0.9\gamma$ and $\gamma_{z}=\gamma$. For $\gamma=0$ the system is trivially integrable as $\hat{U}_{\mathrm{QKT}}$ generates only rotations, and as $\gamma$ is increased the classical phase space becomes mixed, with islands of regular motion separated by areas of chaos. For $\gamma\gtrsim 2$, most of the phase space becomes chaotic. This transition to chaos can be clearly observed from the normalized average adjacent spacing ratio $\overline{r}_{\mathrm{norm}}$, introduced in the previous section (see also Appendix B), which we show in Fig. 4 (b).

Being able to access large system sizes also means that, even when coarse grained, it can be hard to visualize the evolution of each component of the probability distribution $\{P_{k}(t)\}$ defined in Eq. (17) in a manner similiar to what was done for a small instance of the Ising model in Fig. 2. In Fig. 5 (a) we present density plots of the distribution at short times for four representative values of $\gamma$ and an initial choice of operator $\hat{O}(0)=\hat{J}_{z}$. Each horizontal slice corresponds to a snapshot of the distribution at a given time, while vertical slices show the evolution of individual components. The plots illustrate how the distribution, which is concentrated at rank $k=1$ at $t=0$, spreads onto higher ranks faster as $\gamma$ is increased. In order to capture the long-time properties of this evolution, we display in Fig. 5 (b) the mean $\mu(t)$, variance $\sigma(t)$, and IPR $\eta_{\mathrm{IPR}}(t)$ of the operator distribution as a function of time. As expected, we observe how the values predicted for Haar-random evolution (dashed lines) are readily attained as $\gamma$ increases and the QKT becomes chaotic. We also observe an interesting similarity between the behavior of both the mean and the variance when compared to the Ising case, cf. Fig. 6 (a). In the ergodic cases the mean $\mu(t)$ increases steadily and saturates at the prediction from Table 1, but the variance temporarily ‘shoots up’ before it equilibrates. We emphasize that this same behavior has been observed for the SYK model in a previous work Roberts et al. (2018). The fact that the QKT, which is essentially a quantized classical system, also displays similar features is indicative of the presence of unifying features in the evolution of the operator distribution, even when considering quantum chaotic models of very different natures. Interestingly, we observe in both the QKT and the Ising model that the variance of the distribution can be systematically larger in the non-chaotic regime with respect to the chaotic case.

Finally, we study the long-time averages Eq. (20) and time-fluctuations Eq. (21) of the different measures of the operator probability distribution as a function of the $\gamma$ parameter in the QKT, analogous to the results presented for the Ising model in Fig. 3 (b). Results for the QKT are shown in Fig. 6 for two choices of operators: $\hat{O}(0)=\hat{J}_{z}$ (dark lines) and $\hat{O}(0)=\hat{J}_{y}$ (light lines). Similar to the results obtained for the Ising model, in the chaotic regime of the QKT the time-averaged mean, variance, and IPR closely match the predictions from Haar-random evolution, with vanishing temporal fluctuations. In the opposite (integrable) regime, the trivial dynamics at $\gamma\sim 0$ shows mostly localized distributions and little spreading (and, correspondingly, small fluctuations), akin to the regime of $\theta\sim\pi/2$ of the Ising model.

In the transition to chaos, for $0\lesssim\gamma\lesssim 2$, the properties of the distribution show some unifying features, but is overall operator-dependent. We observe that the initial operator $\hat{O}(0)=\hat{J}_{z}$ shows significantly more spreading than $\hat{J}_{y}$. Analyzing the classical phase spaces in Fig. 4 (a), one readily observes that stable islands tend to be localized on the equator of the spherical phase space, with unstable areas around the poles. Around these unstable fixed points is where chaos first appears already at $\gamma=1.0$ Reichl (1992); Chinni et al. (2021). We attribute the enhanced growth of $\hat{J}_{z}$ to these instabilities, in a phenomenon closely related to the previously studied saddle-point scrambling Kidd et al. (2021). Aside from this we find for both operators that the transition to chaos is characterized by an enhanced variance of the distribution. Interestingly, the shape of $\overline{\sigma(t)^{2}}$ for the QKT in the regime $0\leq\gamma\leq 2$ is very similar to that of the Ising model in the equivalent regime $\pi/4\lesssim\theta\leq\pi/2$. This similarity once again hints at a unified behavior of the operator growth in the transition from nonergodic to ergodic behavior.

Finally, we find that temporal fluctuations (bottom row of Fig. 6) in the QKT are a good proxy for quantum chaos in the model, similar to what we observed earlier for the Ising model. However, in this case the correspondence does not extend to very small values of $\gamma$, since the trivial dynamics of the QKT leads to a roughly constant operator distribution. As in the Ising model we find for the chaotic case of the QKT that the temporal fluctuations decrease with system size while the time-averaged measures are roughly independent of the values of $N$ (see Appendix C).

As a final case study, we analyze how the operator distribution $\{P_{k}(t)\}$ behaves for a quantum circuit akin to the ones studied in quantum computing Nielsen and Chuang (2010); Brylinski and Brylinski (2002). We consider circuits formed by gates in the universal set $\{H,T,S,CX\}$, where

It is known that combinations of $\{H,S,CX\}$ create a Clifford circuit Gottesman (1998), the action of which takes a Pauli operator to another Pauli operator. However, the $T$ gate is not a Clifford gate and thus transforms a Pauli operator to a sum of Pauli operators. From the Gottesman-Knill theorem, the evolution of a Clifford circuit can be simulated efficiently on a classical computer Gottesman (1998), whereas the presence of $T$ gates make the simulation of the non-Clifford circuit inefficient Gottesman (1998); Aaronson and Gottesman (2004); Bravyi and Gosset (2016). Universal quantum computing requires non-Clifford operations such as the $T$ gates, and the role of these operations in information scrambling has been studied in recent works Leone et al. (2021).

Here we focus on how the probability distribution $P_{k}(t)$ changes as one considers Clifford versus non-Clifford circuits. We study a random circuit model where each layer has a random arrangement of Clifford gates as in Fig. 7 (a), and where, with probability $p_{T}$, a $T$ gate is applied in the layer. The numerical results are obtained after randomly sampling 40 instances of the circuit for an initial operator $\sigma_{z}^{(1)}$ for $N=6$ qubits. In Fig. 7 (c) each probability $P_{k}(t)$ is plotted as a function of the circuit depth for different values of $p_{T}$. From our analysis it is clear that if the probability $p_{T}$ is not too small, the operator distribution becomes the one predicted from the Haar-random evolution after some depth of the circuit, and the required depth gets lower as we increase the $p_{T}$. This behavior is in line with the fact that including the non-Clifford gates makes the gate set universal, and thus random circuit instances are able to explore uniformly the space of unitaries. On the other hand it is interesting to analyze the behavior at small $p_{T}$, where the system behaves (roughly) as a Clifford circuit. As can be seen from the evolution of the probabilities $P_{k}(t)$, this does not preclude the spreading of operators into higher weights. However, since Pauli operators are approximately mapped into Pauli operators, the distribution is very localized at all times. This is the signature of quasi-scrambling (as opposed to genuine scrambling Zhuang et al. (2019)), also termed operator spreading (as opposed to operator entanglement) Mi et al. (2021).

Finally, in Fig. 7 (d), we analyze the mean, variance, and IPR of the operator distribution for this random circuit model, similarly to the analysis of tilted field Ising model and QKT in the previous sections. For sufficiently large $T$ gate probabilities $p_{T}$, yielding a large number of $T$ gates, the behavior of these are identical in nature to the one found for other models in the chaotic regime, see Fig. 5 (b) and tilted field Ising model in Fig. 3 (a): a sharp increase of the mean and variance and quick equilibration to the random predictions, including the shoot-up of the variance at intermediate times before also equilibrating. For small $p_{T}$, the quasi-scrambling behavior is clearly seen in these measures: while the mean of the distribution increases in a manner very similar to the other cases (albeit with enhanced temporal fluctuations), the variance and IPR show very different behaviors as the distributions localize when $p_{T}\rightarrow 0$. This analysis shows how the breaking down of classical tractability at $p_{T}\neq 0$ can be witnessed, dynamically, by the early-time growth of the operator distribution variance (or, conversely, the decay of its IPR).

In Sec. II we argued that the properties of scrambling, understood as delocalization of quantum information along the degrees of freedom of a system, are encoded in the coarse-grained operator distribution $\{P_{k}(t)\}$ defined over a particular partitioning of the operator basis, cf. Eqs. (3) and (4). Studies of scrambling in the literature, however, are often focused on the analysis of out-of-time-ordered correlators (OTOCs) of the form

where here we consider the OTOC to be evaluated for a thermal state at infinite temperature. In this section we discuss the mathematical connection between the operator distribution and the OTOC by summarizing some previous results in the literature Roberts et al. (2018); Qi et al. (2019) and showing novel relations.

However, one might argue that for many situations of interest, merely obtaining lower order moments like $\mu_{1}$ and $\mu_{2}$ would be sufficient. In an experimental setup, one is then confronted with the fact that OTOCs are intrinsically hard to access and typically require auxiliary systems Landsman et al. (2019), time-reversal operations Gärttner et al. (2017); Li et al. (2017); Swingle et al. (2016), or statistical correlations through randomized measurements Vermersch et al. (2019); Joshi et al. (2020), and clearly measuring $\sim N^{2}$ hard objects is undesirable. Notice that one might also invoke arguments of self-averaging, i.e., measuring a single choice of $\hat{R}$ in Eq. (27) might give a satisfactory indicator of the behavior of the averaged OTOC (and thus of the mean operator size) if the system is sufficiently ergodic. However, as we have shown in this work, the properties of the operator distribution are able to distinguish interesting features of the behavior of the system even in nonergodic regimes, where self-averaging might not work. For instance, both the variance ‘shoot-up’ at short times and the enhanced time-averaged spreading of the distribution observed in Figs. 2 and 5 are only noticeable beyond the chaotic regime.

The present discussion thus shows that while OTOCs allow access to some aspects of the operator distribution, the quantitative connection between the two is not straightforward in practice and might become hard to probe even at the level of the first moments. It is thus desirable to think about other potential methods to probe the distribution more directly, a subject which has attracted considerable attention recently Qi et al. (2019); Schuster and Yao (2022) and which will be the focus of an upcoming work by the authors Blocher et al. (2022b).

In this work we studied scrambling in quantum systems by analyzing the spreading of initially simple operators on a coarse-grained basis, a process which we describe via the operator distribution $\{P_{k}(t)\}$. We considered systems of spin-$1/2$ particles (qubits) in the basis of Pauli operators ordered by size, and kicked collective spin systems in the basis of spherical tensor operators ordered by rank.

We presented a numerical analysis of two paradigmatic models of quantum chaos in both the many-body and few-body setting: the ‘tilted-field’ Ising model and the quantum kicked top. For both cases we computed the evolution of the operator distribution and studied its properties via computing standard distribution measures such as the mean, variance, and localization. Focusing on long-time properties, we showed that in the chaotic regimes both models evolve to spread-out distributions whose properties match the predictions for Haar-random evolution. In particular, the mean operator size (rank) in these cases is proportional to the system size $N$ while temporal fluctuations are suppressed with increasing system size, and thus the dynamics essentially equilibrates to the random distribution. In the different nonergodic regimes of these models the behavior becomes nongeneric as expected. However, several interesting unifying features are observed, like the enhancement of temporal fluctuations and the increase of the distribution variance to values above the random prediction. In the trivially integrable regimes, the distributions remain localized and show long-lived oscillations. In all the studied models we have seen that the long-time properties of the operator distribution allows one to reconstruct the integrability-to-chaos transition in the studied models. Chaotic regimes are characterized by i) the distribution mean and variance matching with Haar-random predictions and ii) the suppression of temporal oscillations. We have found that deviations from one of these two conditions indicate some degree of nonergodicity.

We also applied the operator distribution framework to a random circuit model and showed how the different properties of the operator distribution change as a Clifford circuit is turned into a universal circuit containing non-Clifford gates. Finally, we studied the connection of the different properties of the distribution to averages of out-of-time-ordered correlators (OTOCs). The ideas and results presented in this work build on previous works and set a path where scrambling could be studied directly from the operator distribution, which can be more physically transparent than OTOCs. An important challenge is to devise experimental protocols that allow one to probe these distributions without resorting to the experimentally challenging OTOCs. Another important aspect that we leave for future work is the study of the short-time behavior of the operator distribution. In several cases, it has been observed that single-site OTOCs show exponential decay, and the corresponding exponent is associated with a quantum Lyapunov exponent. The connection discussed in Sec. VI entails that, if all single-site OTOCs behave in this way (or rather, if the average single-site OTOC does so), then we expect an exponential increase of the mean operator size. Interestingly, this does not say anything about the behavior of other properties of the operator distribution, and it is in principle possible to devise models in which the timescale associated with, e.g. the distribution variance is different than that given by the Lyapunov exponent.

Along these lines, an important path forward is to study how to apply the picture presented in Fig. 1 to systems beyond spin-$1/2$ particles. The analysis of collective spin models studied in this work represents an advance in this direction, since these collective models are completely equivalent to single particles of a fixed total spin $J$. There is thus a notion of scrambling in a single particle for any $J>1/2$ , similar to the manner in which scrambling can be defined for a single bosonic mode Zhuang et al. (2019) (note that the maximum rank for $J=1/2$ is $N=1$ and so the operator distribution is trivial). When considering, for instance, chains of spin-$1$ particles (i.e. circuits on qutrits Blok et al. (2021)) or systems of interacting bosons in lattices Bohrdt et al. (2017), one should think about defining a coarse-graining of the operator basis that takes into account operator spreading within one subsystem, as well as among different subsystems.

Finally, the analysis presented in this work highlights that interesting features about the nonergodic regimes of many-body systems may be studied from the operator distribution. While this work has studied a system which is integrable by mapping to noninteracting particles, there are other systems which are integrable by other mechanisms, like the Heisenberg model Santos et al. (2012). More generally, some many-body systems show dynamical phase transitions defined from their out-of-equilibrium properties Heyl (2018). An exciting path forward is to elucidate whether some of these dynamical transitions can be viewed as a transition in the operator distribution, which in turn is initial-state independent.