Operator hydrodynamics, OTOCs, and entanglement growth in systems without conservation laws

TL;DR

This paper studies operator spreading, OTOCs, and entanglement growth in 1D random unitary circuits without conservation laws, revealing a hydrodynamic description with a butterfly velocity smaller than the light cone velocity and diffusive broadening of operator fronts.

Contribution

It provides exact results on OTOC behavior and entanglement growth in non-conserving random circuits, introducing a hydrodynamic framework for understanding quantum information spreading.

Findings

Operator fronts broaden diffusively over time.

OTOCs do not exhibit a prolonged exponential growth regime.

Entanglement velocity can be smaller than butterfly velocity.

Abstract

Thermalization and scrambling are the subject of much recent study from the perspective of many-body quantum systems with locally bounded Hilbert spaces (`spin chains'), quantum field theory and holography. We tackle this problem in 1D spin-chains evolving under random local unitary circuits and prove a number of exact results on the behavior of out-of-time-ordered commutators (OTOCs), and entanglement growth in this setting. These results follow from the observation that the spreading of operators in random circuits is described by a `hydrodynamical' equation of motion, despite the fact that random unitary circuits do not have locally conserved quantities (e.g., no conserved energy). In this hydrodynamic picture quantum information travels in a front with a `butterfly velocity' $v_{\text{B}}$ that is smaller than the light cone velocity of the system, while the front itself broadens diffusively in time. The OTOC increases sharply after the arrival of the light cone, but we do \emph{not} observe a prolonged exponential regime of the form $\sim e^{\lambda_\text{L}(t-x/v)}$ for a fixed Lyapunov exponent $\lambda_\text{L}$. We find that the diffusive broadening of the front has important consequences for entanglement growth, leading to an entanglement velocity that can be significantly smaller than the butterfly velocity. We conjecture that the hydrodynamical description applies to more generic ergodic systems and support this by verifying numerically that the diffusive broadening of the operator wavefront also holds in a more traditional non-random Floquet spin-chain. We also compare our results to Clifford circuits, which have less rich hydrodynamics and consequently trivial OTOC behavior, but which can nevertheless exhibit linear entanglement growth and thermalization.

Equations

$$ZX$$

$$Z^q$$

$${\sigma }^{\mu }\left(\tau \right)\equiv U^{†}(\tau ){\sigma }^{\mu }U(\tau )=\sum _{\nu }{c}_{\nu }^{\mu }\left(\tau \right){\sigma }^{\nu }.$$

$${\rho }_R^{\mu }(s,\tau )\equiv \sum _{\nu }{\left|c_{\nu }^{\mu }\left(\tau \right)\right|}^2\delta (\text{RHS}(\nu )=s),$$

$$U(t)=\prod _{\tau =1}^{2t,\leftarrow }\prod _{x=1}^{L/2}W(2x-1+n_{\tau },\tau )$$

$$\overline{{\rho }_R^{\mu }}(s,\tau )\equiv \sum _{\nu }\overline{|c_{\nu }^{\mu }(\tau ){|}^2}\delta (\text{RHS}(\nu )=s).$$

$$\overline{{\rho }_R}(s,0)=\delta (\text{RHS}(\mu )=s),$$

$$\overline{{\rho }_R}(s,\tau +1)=(1-p)\left[\overline{{\rho }_R}(s,\tau )+\overline{{\rho }_R}(s+1,\tau )\right];$$

$$\overline{{\rho }_R}(s+1,\tau +1)=p\left[\overline{{\rho }_R}(s,\tau )+\overline{{\rho }_R}(s+1,\tau )\right],$$

$$\overline{{\rho }_R}(x,t)\equiv \overline{{\rho }_R}(s=2x-1,\tau =2t)+\overline{{\rho }_R}(s=2x,\tau =2t),$$

$$\begin{gathered} \overline{{\rho }_R}(x,t+1)=2p(1-p)\overline{{\rho }_R}(x,t)+ \\ +p^2\overline{{\rho }_R}(x-1,t)+(1-p{)}^2\overline{{\rho }_R}(x+1,t). \end{gathered}$$

$$\overline{{\rho }_R}(x,t)=\frac{q^{2(t+x)}}{(1+q^2{)}^{2t}}\left(\genfrac{}{}{0ex}{}{2t}{t+x}\right).$$

$$\overline{{\rho }_R}(x=v_{\text{B}}t+O(\sqrt{t}))=\frac{1}{\sqrt{\pi (1-v_{\text{B}}^2)t}}{e}^{-\frac{(x-v_{\text{B}}t{)}^2}{(1-v_{\text{B}}^2)t}},$$

$$R(x)$$

$$\overline{R}(x=v_{\text{B}}t+O(\sqrt{t}))\approx \frac{1}{2}\left[\text{erf}(\frac{x-t{v}_{\text{B}}}{\sqrt{t(1-v_{\text{B}}^2)}})+1\right],$$

$$\overline{R}(x)\approx (1-2\overline{{\rho }_R}(x))\Theta (x-v_{\text{B}}t)+\overline{{\rho }_R}(x),$$

$$\overline{{\rho }_R}(s=2x+n,\tau =2t)=\frac{q^{2(t+x-1+n)}}{(1+q^2{)}^{2t}}\left(\genfrac{}{}{0ex}{}{2t-1}{t+x-1}\right),$$

$$\mathcal{C}(s,\tau )$$

$$=1-\text{Re}⟨\psi \mid Z_0\left(\tau \right)Z_s{Z}_0^{-1}\left(\tau \right)Z_s^{-1}\mid \psi ⟩$$

$$1-\overline{\mathcal{C}}(s,\tau )$$

$$\overline{\mathcal{C}}(s,\tau )\approx 1-\overline{R}\left(s-1,\tau \right)+q^{-2}\overline{{\rho }_R}\left(s,\tau \right).$$

$$\overline{\mathcal{C}}(s,\tau )\approx e^{\frac{1}{2}\delta \log (\frac{\gamma s}{\delta })-\frac{1}{6\delta }}\times {\left(\frac{q^2}{q^2+1}\right)}^s\frac{(1+q^2)\sqrt{\delta }}{2s\sqrt{\pi }}.$$

$$\overline{\mathcal{C}}(s,\tau )\approx \frac{1}{2}\text{erfc}\left(\frac{s-v_{\text{B}}\tau }{\sqrt{2\tau (1-v_{\text{B}}^2)}}\right).$$

$$1-\overline{\mathcal{C}}(s,\tau )\approx \frac{(1+q^{-2})q^s}{\sqrt{8\pi \tau }}{\left(\frac{2q}{1+q^2}\right)}^{\tau }$$

$$\hat{\omega }=\frac{1}{q^L}\prod _{s=1}^L\left(\sum _{k=0}^{q-1}{Z}_s^k\right)=\frac{1}{q^L}\sum _{\nu \in \text{Z-strings}}{\sigma }^{\nu }.$$

$$e^{-S^{(2)}}=\frac{1}{q^{L_{\text{A}}}}\sum _{\nu ,{\nu }^{\prime }}\sum _{\mu \subset A}{c}_{\mu }^{\nu }{c_{\mu }^{{\nu }^{\prime }}}^{\ast }\approx \frac{1}{q^{L_{\text{A}}}}\sum _{\nu }\sum _{\mu \subset A}|c_{\mu }^{\nu }{|}^2$$

$$\overline{e^{-S^{(2)}(\tau )}}=\frac{1}{q^{L_{\text{A}}}}+\frac{q^2-1}{q^2}\sum _{x_0=1}^{L/2}\frac{\overline{R}(L_{\text{A}}/2-x_0,t=\tau /2)}{q^{L_{\text{A}}-2x_0}},$$

$$\overline{e^{-S^{(2)}(\tau )}}=q^{-L_{\text{A}}}+[1-q^{-L_{\text{A}}}]\frac{q^{\tau }}{{\left(1+q^2\right)}^{\tau }}\sum _{x=-\tau /2}^{L_{\text{A}}/2-1}\left(\genfrac{}{}{0ex}{}{\tau }{\frac{\tau }{2}+x}\right).$$

$$\overline{e^{-S^{(2)}(\tau )}}\approx {\left(\frac{2q}{1+q^2}\right)}^{\tau },$$

$$v_{\text{E}}\equiv \frac{1}{\log q}\frac{\text{d}S(\tau )}{\text{d}\tau }$$

$$\hat{U}=e^{-i\frac{T}{2}h{\sum }_s{X}_s}{e}^{-i\frac{T}{2}{\sum }_s[Z_s{Z}_{s+1}+g{Z}_s]}.$$

$${\partial }_t{\rho }_R(x,t)+{\partial }_{x}J(x,t)=0$$

Full text

I Introduction

The past decade has seen a great revival of interest in the foundations of quantum statistical mechanics. It has been driven by theoretical advances involving the long sought demonstration that many-body localization (MBL) exists Basko et al. (2006), ideas from quantum information theory and the study of integrable systems Essler and Fagotti (2016). It has been equally driven by experimental advances in the study of cold atomic gases which provide examples par excellence of closed macroscopic quantum systems for which the foundational questions of quantum statistical mechanics are especially acute. Bloch et al. (2008) Perhaps the broadest question has to do with identifying possible “ergodic universality classes” in quantum many-body systems and understanding their more detailed physics. Of such a potential classification much work has focused on fully MBL systems which are believed to exhibits a breakdown of statistical mechanics. Nandkishore and Huse (2015)

Most recently a great deal of attention has focused on the related question of how to quantify “scrambling” in many-body systemsSekino and Susskind (2008); Brown and Fawzi (2012); Lashkari et al. (2013); Shenker and Stanford (2014a); Hosur et al. (2016). In this work, we will use the word scrambling to denote those features of the spreading of quantum information which are quantified by the out-of-time-ordered commutator (OTOC), which has been studied in the SYK model and its descendantsSachdev and Ye (1993); Kitaev , as well as MBL systemsChen et al. (2016); Huang et al. (2016); Fan et al. (2017); Chen (2016); He and Lu (2017), in field theoretic settingsRoberts and Stanford (2015); Stanford (2016); Chowdhury and Swingle (2017); Aleiner et al. (2016); Dubail (2016), and numerically in interacting spin-chainsDóra and Moessner (2016); Bohrdt et al. (2016); Kukuljan et al. (2017); Che . We report some exact analytical results and supporting numerics for interacting non-integrable spin-chains that are interesting in the context of scrambling.

In the following we discuss OTOC behavior and entanglement growth (building on work in Refs. Nahum et al., 2016; Ho and Abanin, 2017; Mezei and Stanford, 2016) for the following three spin-chain models: I) a random circuit (see Eq. (4)) where the 2-site gates are randomly chosen; II) an ergodic Floquet system with nearest-neighbor interactions, defined in Eq. (28) and III) a periodic Clifford circuit defined in Eq. (33). Our approach relies on quantifying operator spreading, i.e., how the support of operators changes under Heisenberg picture evolution. We derive analytical formulas for operator spreading in model I), which we support with additional numerics, while for model II) we rely entirely on numerical calculations. Our numerical method is based on the matrix product operator (MPO) Verstraete et al. (2004) representation. Since all three types of time evolution we consider can be represented as a network of 2-site gates (see Fig. 1), the MPO can be time evolved straightforwardly by using the TEBD algorithm Vidal (2003). Our results for the three models are as follows:

I) In the random circuit model (Sec. III) we find that operator spreading can be described by a remarkably simple hydrodynamical picture, which gives rise to a biased diffusion equation. Using this, we find that the typical extent of an operator grows with butterfly velocity $v_{\text{B}}$ which is less than the light cone velocity $v_{\text{LC}}$, while the width of the front broadens diffusively in time (see Fig. 2). We use these results to derive exact formulas for the OTOC and entanglement growth. The OTOC travels with the same butterfly velocity $v_{\text{B}}$ and its behavior near the front is also sensitive to the diffusive broadening (Eq. (20)). At early times, before the arrival of the main front, the OTOC grows exponentially, with an exponent that increases with the initial separation of the two operators (Eq. (19)). At long times the OTOC saturates to $1$. This is summarized in Fig. 3. Our front propagation results also lead to an exact formula for the entanglement growth of an initial product state, from which we can extract an entanglement velocity $v_{\text{E}}$ (Eq. (27)). We find that the diffusive broadening of the operator front gives rise to the inequality $v_{\text{E}}<v_{\text{B}}$. This exact result is consistent with general non-rigorous argumentsHartman and Afkhami-Jeddi (2015); Mezei and Stanford (2016), the heuristic operator spreading model of Ref. Mezei and Stanford, 2016, numerous results in holographyHartman and Maldacena (2013); Liu and Suh (2014a, b), and the results derived for Clifford circuits in Ref. Nahum et al., 2016.

II) In Sec. IV we verify numerically that, for a family of ergodic Floquet circuits, there is a similar diffusively broadening front behavior as observed in the random circuit (see Fig. 7). This leads to the tentative conjecture that the diffusive front picture is valid for generic ergodic 1D spin-chains, along with the resulting consequences for OTOC and entanglement dynamics.

III) Finally in Sec. V we compare I) and II) to Clifford circuits. Within such circuits, strings of Pauli operators evolve to other particular strings, rather than superpositions of such; in particular, Clifford circuits do not exhibit a diffusively broadening operator front. We connect this fine tuned nature of local Clifford circuits to the fact that such circuits always exhibit trivial OTOC behavior. Nevertheless, they can still have linear entanglement growth and their local observables thermalize to infinite temperature. This demonstrates the broader point that the presence of ballistic entanglement growth and thermalization are not sufficient to predict scrambling behavior (i.e., the behavior of the OTOC).

The systems above do not have local conserved quantities (in particular, they do not conserve energy), so it is surprising that hydrodynamics arises in I) and II). We conjecture that such hydrodynamic behavior is universal, generically appearing in 1D ergodic systems with local unitary evolution and a bounded local Hilbert space. [In this connection, our work has obvious parallels with Ref. Nahum et al., 2016, although as we will discuss in Sec. III the hydrodynamics in our case has a different physical origin]. Moreover, our work makes a clear and precise connection between the spreading dynamics of operators, the scrambling behavior captured by the OTOC and other metrics of ergodicity such as entanglement entropy and the late time behavior of local correlation functions (see App. D).

II Quantifying operator spreading

Consider a one-dimensional chain of $L$ sites, for which the Hilbert space of a single site is $\mathcal{H}_{\text{site}}=\mathbb{C}^{q}$. There exist operators $X,Z$ on the single site Hilbert space obeying

[TABLE]

These generate a convenient complete basis for all operators on $\mathcal{H}_{\text{site}}$, namely $\{\sigma^{\mu}\equiv X^{\mu^{(1)}}Z^{\mu^{(2)}}:\mu\in\mathbb{Z}_{q}^{\otimes 2}\}$. Here $\mu$ is shorthand for the doublet $\mu^{(1)},\mu^{(2)}\in\left\{0,1,\ldots,q-1\right\}=\mathbb{Z}_{q}$. This basis is orthonormal, such that $\mathrm{tr}(\sigma^{\mu\dagger}\sigma^{\nu})/q=\delta_{\mu\nu}$. The operators $\sigma^{\mu}$ can be regarded as generalizations of Pauli matrices, where the usual Paulis correspond to the $q=2$ case. Generalizing this to the Hilbert space of a 1D chain, $\mathcal{H}_{\text{chain}}=\left(\mathbb{C}^{q}\right)^{\otimes L}$, a complete orthonormal basis of operators is given by the $q^{2L}$ Pauli strings, defined as $\sigma^{\mathbf{\mu}}\equiv\bigotimes_{r=1}^{L}\sigma_{r}^{\mu_{r}}$, where each string is indexed by a vector $\mathbf{\mu}\in\left(\mathbb{Z}_{q}^{\otimes 2}\right)^{\otimes L}$.

Our goal is to quantify how an initial Pauli string spreads over the space of all Pauli strings under local unitary time evolution. At time $\tau$ the Pauli string $\sigma^{\mathbf{\mu}}$ becomes

[TABLE]

This defines a set of ‘operator spread coefficients’ $c_{\mathbf{\nu}}^{\mathbf{\mu}}\left(\tau\right)\equiv\text{tr}\left(\sigma^{\mathbf{\nu}\dagger}U^{\dagger}(\tau)\sigma^{\mathbf{\mu}}U(\tau)\right)/q^{L}$. The full set of coefficients $\{c_{\mathbf{\nu}}^{\mathbf{\mu}}(\tau)\}$ encodes all information regarding the unitary time evolution. However, as we show below, for accessing most physically interesting quantities, such as entanglement entropiesHo and Abanin (2017); Mezei and Stanford (2016); Nahum et al. (2016), or out-of-time-order commutators, it is sufficient to consider a more coarse-grained description of the operator spreading. One particularly useful coarse-grained quantity is the total weight on all operators with right endpoint $s$ appearing in $\sigma^{\mu}(\tau)$, i.e.,

[TABLE]

where $\text{RHS}(\mathbf{\nu})$ denotes the rightmost site on which $\mathbf{\nu}$ is non-zero. 111A similar quantity was defined in Ref. Roberts et al., 2015 Note that the ‘density’ $\rho^{\mathbf{\mu}}_{R}$ is conserved, i.e., $\sum_{s}\rho^{\mathbf{\mu}}_{R}(s)=1$ at all times. Motivated by this, we refer to $\rho_{R}^{\mathbf{\mu}}$ as the operator density of the time evolved Pauli string $\sigma^{\mathbf{\mu}}$.

In this paper we consider systems where the time evolution can be represented as a circuit of 2-site unitary gates, arranged in the geometry shown in Fig. 1. Sites of the 1D chain are indexed by $s=1,\ldots,L$ while the layers of the circuit are indexed by the variable $\tau$. It is useful to introduce coarse-grained coordinates $x$ and $t$ that label pairs of sites and pairs of layers, respectively, as defined in Eq. (8) and illustrated in Fig. 1. Note that due to the geometry of the circuit all such models have a well defined light cone velocity $v_{\text{LC}}=\Delta s/\Delta\tau=\Delta x/\Delta t=1$, corresponding to the fact that with each successive time-step a local operator can spread at most one additional site in each direction.

In the following we investigate the three models I)-III) discussed in the introduction, all three of which can be represented by local circuits of the kind shown in Fig. 1. In Sec. III.1 we show that for model I) the average of $\rho^{\mathbf{\mu}}_{R}$ obeys a classical biased diffusion equation, and use this to derive exact formulas for the behavior of OTOCs and the dynamics of entanglement. In Sec. IV we investigate model II) and find that it shares many features with the random circuit model, such as a broadening of the propagating wavefront, which suggests that the random walk description obtained in Sec. III has applications in a wider class of ergodic systems. In Sec. V we contrast this with model III) where there is no diffusion and $\rho_{R}^{\mu}$ remains a delta-function at all times, which corresponds to a non-generic behavior of the OTOC.

III Random circuit model

Most of the present paper is concerned with one dimensional local random unitary circuits, with the geometry shown in Fig. 1. Random circuits were also investigated in Ref. Nahum et al., 2016 with regards to the growth of entanglement from an initial product state, albeit with a different geometry where the location of the unitary gates is randomly chosen, instead of the regular arrangement used here. There it was argued that the evolution of entanglement obeys an equation belonging to the KPZ universality class, which determines certain universal exponents that appear in the average value and fluctuations of the entanglement entropy. Here we shift our focus from states to operators and derive exact results for their spreading, for arbitrary on-site Hilbert space dimension. In Sec. III.4 we relate our operator spreading results to the dynamics of bipartite entanglement, as captured by the second Rényi entropy, and find no sign of the universal fluctuations observed in Ref. Nahum et al., 2016. This suggests that the KPZ-like behavior is specific to circuits whose geometry (and not just the individual circuit elements) is random (compare Fig. 1 above to Fig. 6 of Ref. Nahum et al., 2016).

The random circuits we discuss are defined as follows. Consider a discrete time evolution, consisting of layers of two-site unitary gates acting on pairs of neighboring sites. Odd numbered layers act on all the odd bonds of the chain while even numbered layers act on even bonds. Each two-site gate is chosen independently from the Haar distribution over $q^{2}\times q^{2}$ unitary matrices. The time evolution after an even number of $2t$ layers is given by

[TABLE]

where $n_{\tau}=\frac{1+(-1)^{\tau}}{2}$ and $W(s,\tau)$ is a Haar random two site unitary acting on sites $s,\,s+1$. The product $\prod_{\tau=1}^{2t,\leftarrow}$ is defined to be time ordered. Such a circuit is graphically illustrated in Fig. 1.

The primary goal of this work is to quantify the spread of operators under random circuits Eq. (4), and to relate operator spread to entanglement growth. A related question is how correlation functions of local observables behave in this random circuit model in the thermodynamic limit. As we confirm in App. D, such correlations tend to their infinite temperature values at long times, similarly to the case of Floquet ergodic systemsLazarides et al. (2015, 2014); D’Alessio and Rigol (2014); Abanin et al. (2014); Ponte et al. (2015). This result holds for any random realization of the circuit.

Focusing on the problem of operator spreading in this random circuit model, we find that the average of the operator density $\rho_{R}$, defined in Eq. (3), performs a biased random walk, independent of the internal structure of the operator considered. Solving the random walk problem allows us to derive exact formulas for the averages of out-of-time-order commutators and entanglement growth, which we detail in Sec. III.2 and III.4, respectively. In Sec. III.3 we quantify numerically the fluctuations between different random realizations of the circuit.

III.1 Random walk dynamics of operator density

In the following we will quantify how operators spread under the time evolution generated by the random circuit defined above in Eq. (4). We focus on the average of the operator density, defined in Eq. (3), for which we derive an exact equation of motion. Upon solving this equation we find that the operator density moves in a front whose velocity $v_{\text{B}}$ is an increasing function of the on-site Hilbert space dimension and with a front width increasing diffusively in time.

We start by noting that under Haar averaging the operator spread coefficient $c_{\mathbf{\nu}}^{\mathbf{\mu}}(\tau)$ vanishes for any time $\tau\geq 1$, provided that $\mathbf{\mu}$ is non-trivial, since $c_{\mathbf{\nu}}^{\mathbf{\mu}}$ and $-c_{\mathbf{\nu}}^{\mathbf{\mu}}$ have equal probability. However, the average of its modulus squared, $\overline{|c_{\mathbf{\nu}}^{\mathbf{\mu}}(\tau)|^{2}}$, can be non-zero. (An explicit expression for this quantity is written down in App. C, using a mapping to a classical Ising model, but we will not require it for the subsequent discussion).

Following Eq. (3), we define the average operator density as

[TABLE]

As we show below, this quantity satisfies an equation of motion, Eq. (9), which does not depend explicitly on $\mathbf{\mu}$. Pre-empting this, we drop the explicit $\mathbf{\mu}$ dependence $\overline{\rho_{R}^{\mathbf{\mu}}}\rightarrow\overline{\rho_{R}}$ to declutter notation. In fact, $\mathbf{\mu}$ will enter considerations only as an initial condition on the operator density

[TABLE]

which is the same for all initial operators sharing the same right endpoint.

To understand how $\overline{\rho_{R}}$ evolves in time, consider the effect of applying a single two-site gate on sites $s$ and $s+1$. There are $q^{4}-1$ nontrivial operators acting on this two-site Hilbert space. Of these, $q^{2}-1$ contribute to $\overline{\rho_{R}}(s,\tau)$ (the ones that are trivial on site $s+1$), while the other $q^{2}(q^{2}-1)$ contribute to $\overline{\rho_{R}}(s+1,\tau)$. Under a two site Haar random unitary transformation all the possible transitions between any of these $q^{4}-1$ operators have, on average, the same probabilityBrown and Fawzi (2012). The upshot is that after the application of the unitary gate the density $\overline{\rho_{R}}$ evolves as

[TABLE]

with probabilities $p=\frac{q^{2}}{q^{2}+1}$ and $1-p=\frac{1}{q^{2}+1}$. To apply a similar argument for two subsequent layers of the circuit it is useful to redefine the density by grouping together the pairs of sites on which the first layer of the circuit acts. We abuse notation and denote this quantity as

[TABLE]

where we now only consider the value of the operator density at even time steps $\tau=2t$. Applying Eq. (7) for two layers we arrive at the equation

[TABLE]

Thus the right endpoints of Pauli strings perform a biased random walk on the lattice, where in each step they move to the right with probability $p^{2}$, to the left with probability $(1-p)^{2}$, and stay on the same site otherwise. A feature of the above equation is that the time evolution of $\overline{\rho_{R}}$ is independent of the internal structure of the operator and thus the solution $\overline{\rho_{R}}(x,t)$ will be the same for all initial Pauli strings, modulo a shift $x\to x-x_{0}$ where $x_{0}$ is defined by the right endpoint of the initial string.

The result of the random walk process outlined above is a front that propagates to the right from its initial position $x_{0}$ as $\langle x\rangle-x_{0}=v_{\text{B}}t$ with a butterfly velocity $v_{\text{B}}=p^{2}-(1-p)^{2}=\frac{q^{2}-1}{q^{2}+1}$. Thus the speed at which the operator weight travels is smaller than the light cone velocity: $v_{\text{B}}<v_{\text{LC}}=1$. This resonates somewhat with the result of Ref. Roberts and Swingle, 2016. The width of the front increases in time as $\langle x^{2}\rangle-\langle x\rangle^{2}=2Dt$ with diffusion constant $D=\sqrt{1-v_{\text{B}}^{2}}/4=\frac{q/2}{q^{2}+1}$. Note that in the limit $q\to\infty$ the ‘particle’ described by $\overline{\rho_{R}}(x,t)$ hops to the right with probability $1$ in each step, and consequently the front becomes infinitely sharp with velocity $v_{\text{B}}\to v_{\text{LC}}=1$.

The total weight of left endpoints, $\overline{\rho_{L}}(x,t)$, obeys a similar equation except that it propagates to the left with velocity $-v_{\text{B}}$, while diffusing at the same rate, as shown in Fig. 2. This means that at time $t$ the vast majority of quantum information initially stored in $\sigma^{\mathbf{\mu}}$ with left (right) endpoint $x_{l}$ ($x_{r}$) is carried by operators with support $[x_{l}-v_{\text{B}}t,x_{r}+v_{\text{B}}\tau]$, but where the precise position of either endpoint can be uncertain within a region of width $\Delta x\sim\sqrt{Dt}$.

We can find the full distribution of $\overline{\rho_{R}}(x,t)$ using a standard generating functional method. In the rest of this section we will use coordinates relative to the initial position of the front, i.e. $x-x_{0}\to x$. The solution to Eq. (9) than reads

[TABLE]

In the scaling limit $t,x\to\infty$ but keeping $x/t\approx v_{\text{B}}$ fixed this becomes (using Stirling’s approximation)

[TABLE]

so that the traveling front has the shape of a Gaussian, as one would expect from the solution of the continuum limit of the lattice diffusion equation Eq. (9).

As we shall see in the next section, it is also useful to compute the total weight of all Pauli strings contained entirely to the left of position $x$. This quantity, which we denote by $R(x),$ is given by

[TABLE]

Around the position of the front, where $x\approx v_{\text{B}}t$, we can integrate Eq. (11) to obtain

[TABLE]

where $\text{erf}(x)$ is the error function.

Later on we will also need an approximation for $R$ well away from the front. Using the fact that in the large $t$ limit $\overline{\rho_{R}}(x,t)$ increases sharply with $x$ for $x/t<v_{\text{B}}$, the sum in Eq. (12) is dominated by its largest term (i.e. $y=x$). Using the fact that $\overline{R}\left(x\geq v_{\text{LC}}t\right)=1$ we can similarly approximate $\overline{R}$ for $x/t>v_{\text{B}}$, to obtain

[TABLE]

where $\Theta$ is the Heaviside step function. This result is accurate up to multiplicative $O(1)$ constants when $|x-v_{\text{B}}t|/t=O(1)$ in the large $t$ limit. See the discussion in App. A for a precise statement and derivation of Eq. (14).

Using our results for the coarse-grained density $\overline{\rho_{R}}(x,t)$ we can also write a formula for the density in terms of the site coordinate $s$. Note that due to Eq. (7) the ratio $\overline{\rho_{R}}(s=2x+1)/\overline{\rho_{R}}(s=2x)=q^{2}$ is fixed at any time $\tau=2t$. Using this, the density on site $s$ after applying an even number of layers becomes

[TABLE]

where $n=0,1$. We can use Eq. (15) to derive the total operator weight left of site $s$, i.e. $R(s)=\sum_{r\leq s}\rho(r)$, which, as we will see in the next section, is closely related to the OTOC between sites [math] and $s$.

III.2 Behavior of out-of-time-order commutators

We relate our results for the time evolution of operator weights to another oft used measure of information spreading in many-body systems, the so-called out-of-time-order commutator (OTOC)Shenker and Stanford (2014b, a, 2015); Roberts et al. (2015); Roberts and Stanford (2015); Maldacena et al. (2016); Hosur et al. (2016). For concreteness consider the following OTOC between two Pauli operators separated by distance $s$ (in this section we work in a shifted co-ordinate system where one of the Pauli operators resides at site [math])

[TABLE]

where $s$ and $\tau$ are the original time/lattice co-ordinates (as opposed to the coarse grained co-ordinates $t,x$ below Eq. (8)). We will show how the OTOC $\mathcal{C}(s,\tau)$ behaves in the scaling limit $\tau\rightarrow\infty$, with $\kappa\equiv s/\tau$ held fixed. We we detail below, for $s$ outside of the light cone ($1<\kappa$) it is zero. As $s$ enters the light-cone ($\kappa<1$ and close to $1$) it increases exponentially. When $s$ is near the operator front ($\kappa=v_{\text{B}}<1$) the OTOC becomes $O(1)$. After the front has passed ($\kappa<v_{\text{B}}$) the front has passed, the OTOC exponentially saturates to the value $1$ with an exponent that is independent of $s$. See Fig. 3 for a summary.

Let $\overline{\mathcal{C}}$ denote the average of the OTOC over all unitary circuits with geometry shown in Fig. 1. Note that due to the averaging this quantity is independent of the choice of Pauli operator, i.e., it would be the same if we replaced either or both operators in Eq. (III.2) with another local Pauli different from $Z$. We will be concerned with the second term in Eq. (III.2) which equals

[TABLE]

where $e^{\theta_{\mu,Z_{s}}}$ is a $q^{\text{th}}$ root of unity arising from commuting $\sigma^{\mu}$ past $Z_{s}$, and $c_{0}^{\mu}(\tau)$ are the operator spreading coefficients of $Z_{0}(\tau)$. Notice that the Haar average forces $\mu=\nu$ which removes all dependence on the particular initial state $\psi$222To see that $\overline{c^{\mu}c^{\nu*}(\tau)}\propto\delta_{\mu\nu}$: If $\mu\neq\nu$ there exists a single site Pauli $\sigma^{\alpha}_{r}$ such that the group commutators $[\sigma^{\alpha}_{r}:\sigma^{\mu,\nu}]=e^{i\theta_{\mu,\nu}}$ are unequal. Using the invariance of the two-site Haar measure under multiplication by single site operators, this implies $\overline{c^{\mu}c^{\nu*}(\tau)}=e^{i(\theta_{\mu}-\theta_{\nu})}\overline{c^{\mu}c^{\nu*}(\tau)}$ from whence the result follows.. In particular, the average OTOC value in any state, pure or mixed, will be identical to the average OTOC value at infinite temperature, i.e., $\text{tr}(\frac{1}{2}|[Z_{0}(\tau),Z_{s}]|^{2})/2^{L}$. At this point, we can use Eq. (17) and Eq. (53) to write an exact closed from expression for the OTOC. However, instead of doing that, we will write a more manageable asymptotic expression for Eq. (17) using simpler results from Sec. III.1.

To perform the sum over Pauli strings in Eq. (17), we first need to prove the following statement: $\overline{|c_{0}^{\mathbf{\mu}}(\tau)|^{2}}$ depends only on the position of the two endpoints of the string $\mathbf{\mu}$. The proof goes as follows. First, it is easy to verify that under Haar averaging $|c_{\mathbf{\nu}}^{\mathbf{\mu}}(\tau)|^{2}=|c_{\mathbf{\mu}}^{\mathbf{\nu}}(\tau)|^{2}$, which means that the average probability of the one-site operator $Z_{0}$ evolving into a specific string $\mathbf{\mu}$ is the same as the probability of string $\mathbf{\mu}$ evolving into $Z_{0}$. In the random walk picture this latter process corresponds to both left and right endpoints ending up on site [math] at time $\tau$ during their respective random walks. As we noted previously, these random walks are independent of the internal structure of the initial string. Thus $\overline{|c_{0}^{\mathbf{\mu}}(\tau)|^{2}}$ depends only on where the two endpoints of $\mathbf{\mu}$ are located. We confirm this argument more concretely with an explicit expression for such operator spread coefficients in App. C.

The above statement has important consequences for the OTOC. If site $s$ lives in the support of $\mathbf{\mu}$ then the contribution to Eq. (17) coming from the strings with the same support as $\mathbf{\mu}$ have an equal weight for each possible value $\theta_{\mu,Z_{s}}\in\frac{2\pi}{q}\{1,\ldots,q\}$, so that the cosine term averages to zero. The remaining part is the total weight due to Pauli strings which are supported on intervals that do not contain site $s$ (along with some corrections for Pauli strings which border on site $s$). Deferring the full justification App. B, the upshot is that provided $\kappa>0$ in the $\tau\rightarrow\infty$ limit, the OTOC behaves as

[TABLE]

up to exponentially small corrections in $\tau$. Hence the OTOC physics is directly related to the operator density, and changes appreciably at the operator front $s=v_{\text{B}}t$, as we show in Fig. 3.

Let us now summarize the behaviour of the OTOC as a function of space and time, as parameterized by the ratio $\kappa\equiv s/\tau$ and taken in the limit $\tau\rightarrow\infty$. We distinguish four regimes of OTOC behaviour which we illustrate in Fig. 3.

OTOC trivial at early times ($1<\kappa$): In this regime the events $(\tau,0),(0,s)$ are causally disconnected, so the commutator in Eq. (III.2) (and hence the OTOC) is exactly zero. 2. 2.

Early OTOC growth ($v_{\text{B}}<\kappa<1$): This regime describes the behaviour after site $s$ has entered the light cone, but before it encounters the main operator front. Here we approximate the OTOC using Eq. (14), so that $\overline{\mathcal{C}}(s,\tau)\approx c_{1}\overline{\rho_{R}}(s-1,\tau)+\overline{\rho_{R}}(s,\tau)$, where $c_{1}>1$ is bounded in the $s,\tau\rightarrow\infty$ limit. Fortunately a simple closed form expression already exists for $\overline{\rho_{R}}$, namely Eq. (15). We obtain a more convenient expression for the initial OTOC growth by expanding Eq. (15) near the light cone in the $\delta^{2}/s\rightarrow 0$ limit, where $\delta\equiv\tau-s$ is the distance between $s$ and the light cone

[TABLE]

up to multiplicative O(1) constants, where $\gamma=e\left(1-v_{\text{B}}^{2}\right)/2$. This formula demonstrates that the OTOC will increase with an exponent $\lambda\sim\log s$ for $0<\delta\ll s$. Due to its dependence on $s$, and its limited range of validity, it is unclear whether this should be viewed as a Lyapunov exponent as in Ref. Maldacena et al., 2016. Note that the exponential increase occurs in a regime where the overall scale of the OTOC is still exponentially small in $s$. In the regime where the OTOC increases to an $O(1)$ value (i.e., when the operator front hits, see next point) its behavior is not exponential. Furthermore, we note that $\gamma\sim 1/q^{2}$ for large $q$, such that the regime in which the exponential behavior can be observed becomes smaller in the large $q$ limit. 3. 3.

Near the front ($|\kappa-v_{\text{B}}|=O(1/\sqrt{\tau})$): As mentioned, the above approximation breaks down when the main front, which we recall travels at speed $v_{\text{B}}$ and has width $\sim\sqrt{\tau}$, arrives at site $s$. In this intermediate regime, we estimate the OTOC by combining Eq. (18) and Eq. (13)

[TABLE]

This formula describes the behaviour of the OTOC in the regime when it increases from a value exponentially small in $s$ to an $O(1)$ number. 4. 4.

Late times ($0<\kappa<v_{\text{B}}$): After the main front has passed the OTOC relaxes exponentially to $1$. Expanding Eq. (14) for fixed $s-v_{\text{LC}}\tau$ and large $\tau$ we find that the OTOC in this late time regime is

[TABLE]

Thus the OTOC decays to its equilibrium value with an exponent $\log(\frac{1+q^{2}}{2q})$.

III.3 Fluctuations from circuit to circuit

The results discussed above concern quantities averaged over many different random circuits with the same geometry but different choices of two-site gates. The question remains regarding whether these average quantities are also ‘typical’, i.e. how large are the fluctuations between different realizations of the random circuit. In this section we investigate this problem numerically. Our numerical method relies on representing the operator $Z_{0}(t)$ as a matrix product operator (MPO), which allows us to apply the two-site unitary gates efficiently. Two layers of the random circuit can be applied by just a single step of the TEBD algorithm, which allows us to go up to bond dimension $\chi=20000$. Both the infinite temperature OTOC and the total operator weight contained in an arbitrary subregion can be extracted straightforwardly from the MPO representation (both calculations are similar to computing the overlap of two matrix product states, but in the computation of $R(s)$ only the legs corresponding to sites $\leq s$ are contracted).

To quantify the fluctuations we look at an ensemble of 100 random circuit realizations (which is enough to reliable reproduce the exact average quantities, see Fig. 4) with on-site Hilbert space dimension $q=2$ and compute a) the OTOC $\mathcal{C}(s,\tau)=1-\mathrm{Tr}[(Z_{0}(\tau)Z_{s})^{2}]/2^{L}$ and b) the total operator weight $R(s,\tau)$ of $Z_{0}(\tau)$ contained within the region left of site $s$. Both $R(s,\tau)$ and $\mathcal{C}(s,\tau)$ are functions of the distance $s$ and the number of layers $\tau$. We find that for both quantities, the circuit-to-circuit fluctuations are largest at the traveling wavefront and become smaller deep behind it. This is shown in Fig. 4. This also shows that there is a well defined front for the information propagation in each individual circuit.

We also find that the the fluctuations decrease in time. Fig. 4 (c) shows the standard deviation of the weight $R(s)$ for different times. We find that the maximum of this standard deviation over all values of $s$ decreases in time, approximately as $\propto\tau^{-\beta}$ with an exponent $0.4<\beta<0.5$.

III.4 Relationship to entanglement spreading

Another question closely related to operator hydrodynamics is the problem of entanglement growth. An initial product state develops spatial entanglement during time evolution. In systems without quenched disorder the entanglement is expected to grow linearly in time, with a growth rate characterized by the ‘entanglement velocity’ $v_{\text{E}}$. In this section we use our results for operator spreading to compute the typical value of the second Rényi entropy between two sides of a spatial entanglement cut and extract the entanglement velocity from it. We find that this velocity is smaller than the butterfly velocity for any finite $q$ and approaches the light cone velocity logarithmically slowly, so that $v_{\text{LC}}-v_{\text{E}}\propto 1/\log{q}$ for large $q$. At long times the Rényi entropy saturates to its maximal value with the saturation becoming increasingly sharp as $q$ is increased.

Consider an initial ferromagnetic product state of the 1D chain where the state on site $s$ is an eigenstate of the local Pauli operator $Z_{s}$ with eigenvalue $+1$ (Note that for the average behavior of the random circuit the choice of initial product state is unimportant). The density matrix $\hat{\omega}$ corresponding to this state is then a sum over all possible $Z$-strings, i.e. Pauli strings that only contain powers of the operator $Z$ on each site:

[TABLE]

The density matrix at time $t$ is obtained by replacing each Pauli string $\sigma^{\nu}$ in Eq. (22) with its time evolved counterpart $\sigma^{\mathbf{\nu}}(t)$.

Let us now divide the system into two regions, $A$ and $B$, the first of which corresponds to sites $1,\ldots,L_{\text{A}}$. Generalizing the formula of Refs. Ho and Abanin, 2017; Mezei and Stanford, 2016, the second Rényi entropy $S^{(2)}=-\log\text{tr}(\hat{\omega}_{\text{A}}^{2})$ of the reduced density matrix $\hat{\omega}_{\text{A}}=\text{tr}_{\text{B}}(\hat{\omega})$ is related to the operator spreading coefficients by

[TABLE]

where the strings $\mathbf{\nu}$ and $\mathbf{\nu^{\prime}}$ are both $Z$-strings and $\mu$ has support entirely in subsystem $A$. In the last equality of Eq. (23) we assumed that the off-diagonal contributions are negligible, which becomes exactly true in the random circuit model once we average over different realizations.

Let us assume that $L_{\text{A}}$ is even. Reverting back to the coarse-grained position $x$ (see Eq. (8)) we recognize Eq. (23) as the total operator weight in region $A$, $R(x=L_{\text{A}}/2,t)$, as defined in Eq. (12), summed over all initial $Z$-strings. As we noted previously, this quantity is on average the same for all initial strings with the same endpoints $x_{0}$. The number of different $Z$-strings with right endpoint $x_{0}$ is $q^{2(x_{0}-1)}(q^{2}-1)$. After averaging over random circuits, and assuming an even number of layers, Eq. (23) thus becomes

[TABLE]

where we have used that $\overline{\rho_{R}}(x)$ (and consequently $\overline{R}(x)$) only depends on the position $x$ relative to the initial endpoint $x_{0}$. The first term in Eq. (24) is the contribution coming from the identity operator, which is responsible for the saturation of the entanglement at long times.

Using the exact solution Eq. (10) one can perform the sum over initial positions to find

[TABLE]

The sum over binomial coefficients can be expressed in terms of a hypergeometric function.

Eq. (25) describes an entanglement that initially increases linearly with time and saturates to the maximum value $L_{\text{A}}\log{q}$ at long times. For $\tau\ll L_{\text{A}}$ we find

[TABLE]

from which we can identify the entanglement velocity

[TABLE]

Note that the entanglement velocity approaches $1$ logarithmically slowly for large $q$, i.e., $v_{E}\sim 1-\log(2)/\log(q)$. This is a separate velocity scale, distinct from, and smaller than the front speed $v_{\text{E}}<v_{\text{B}}$. This difference comes from the diffusive broadening of the operator wavefront. First, it is straightforward to verify that if the wavefront is sharp, i.e., $\overline{R}(x,t)=\Theta(x-v_{\text{B}}t)$ then Eq. (24) gives $v_{\text{E}}=v_{\text{B}}$. Second, we have checked that Eq. (24) gives $v_{\text{E}}=v_{\text{B}}$ even if the wavefront has a width which is finite but independent of time333We checked this for front profiles arising from $\rho_{R}(x,t)\propto e^{-(x-vt)^{2}/\sigma^{2}}$ and $\rho_{R}(x,t)\propto e^{-|x-vt|/\sigma}$.. Hence, we attribute the difference between $v_{\text{B}}$ and $v_{\text{E}}$ to the fact that the operator front broadens in time.

In the right panel of Fig. 5 we compare the exact formula Eq. (25) to the second Rényi entropy as computed numerically (using a matrix product state representation), averaging over $100$ realizations of the circuit and find extremely good agreement. Moreover, the numerical calculation allows us to compare the typical and average values of the Rényi entropy, defined as $S^{(2)}_{\text{typ}}=-\log{\overline{e^{-S^{(2)}}}}$ and $S^{(2)}_{\text{avg}}=\overline{S^{(2)}}$, respectively. We find no significant difference between the two values, showing that there are no strong circuit-to-circuit fluctuations in the entropy and both are captured well by our exact formula. We also found numerically that replacing the Rényi entropy with the von Neumann entropy leads to a slightly larger entanglement velocity.

The entanglement saturates when the contribution of the identity becomes significant (i.e. when all other operators have essentially left the subsystem). Note that the saturation softens, compared to the prediction of the simple operator spreading model of Ref. Ho and Abanin, 2017, which is another consequence of the diffusive broadening of the front. This intermediate saturation regime becomes smaller with increasing $q$, as shown in the left panel of Fig. 5.

IV Comparison with the kicked Ising model

A natural question that emerges in relation to the results stated above, is to what extent are they representative of other, more generic thermalizing quantum many-body systems. To address this question we investigate a system with a periodically driven nearest neighbor Hamiltonian. Our model has the same geometry as the random circuit, shown in Fig. 1 and it similarly does not conserve energy. However, unlike random circuits, it is periodic in time and its two-site (and one-site) gates take a specific form, rather than randomly chosen. Despite these two significant differences, we find that several details of the operator spreading described in Sec. III, such as the diffusive broadening of the wavefront, remain approximately valid.

For concreteness, we consider a model with on-site Hilbert space dimension $q=2$ that consists of switching back and forth between two Hamiltonians, such that one period of the time evolution (with period time $T$) is given by

[TABLE]

This system  which we refer to as the ‘kicked Ising model’  is known to be ergodic, provided that both $g$ and $h$ are sufficiently large. Since at any given time the terms in the Hamiltonian all commute with each other, the time evolution can be represented as a circuit of two-site unitaries (with the one-site rotations included in the two-site gates) with the same geometry as in Fig. 1. As such, it is in fact contained among the ensemble of random circuits considered before. The question is to what extent do the properties of this specific circuit coincide with the average quantities discussed above.

At first, operator spreading in the Floquet system seems very different from the case of the random circuit. There is no inherent randomness and the evolution is completely determined by the internal structure of the initial operator $\sigma^{\mu}$, while for the random circuit the average behavior was independent of the internal structure. The random unitary and Floquet systems do have something in common, however: both involve evolution under a local unitary circuit. This puts strong constraints on the evolution of $\rho_{R}$: not only does it obey global conservation law $\int\text{d}x\rho_{R}(x,t)=const.$ (we revert to a continuum notation for ease of presentation), it should also obey a local conservation law

[TABLE]

for some local current density $J(x,t)$. This conservation law puts severe restrictions on the equation of motion of $\rho_{R}$. For example, one can imagine that in a coarse grained picture, on long enough time scales, the constitutive equation $J\sim v\rho_{R}+D\partial_{x}\rho_{R}+\ldots$ becomes valid; note that the discretized version of this constitutive relation is exactly the random walk equation we derived for the random circuit averaged $\rho_{R}$, see Eq. (9). Therefore it seems plausible that in a sufficiently coarse grained picture, the dynamics might be well approximated by a biased diffusion similar to the one described in Sec. III for the Floquet circuit, with hopping probabilities depending on the microscopic couplings. Here we present numerical evidence in support of this conjecture. Our results can be summarized in three points:

• •

The butterfly velocity $v_{\text{B}}$ depends strongly on the coupling $g$ and can be tuned to be much smaller than the light cone velocity $v_{\text{LC}}$

• •

When tuning the couplings to decrease $v_{\text{B}}$ from its maximal value $v_{\text{B}}\approx v_{\text{LC}}$, the front also becomes wider, as expected for a random walk when increasing the probability of hopping to the left at the expense of the probability of hopping right.

• •

The operator wavefront gets wider during time evolution, with the width increasing in time as $\sim t^{\alpha}$, with an exponent $0.5\lesssim\alpha\lesssim 0.6$

We find numerically a linearly propagating wavefront for the time evolved operator $Z_{0}(t)$, which shows up in both the OTOC and the weight $R(s)$, with the OTOC $\mathcal{C}(s,t)$ saturating to 1 behind the front. While for the random circuit the speed of the front was set by the on-site Hilbert space dimension $q$, for the kicked Ising model we find that this speed can be tuned continuously by changing the value of the transverse field $g$ 444we found that the butterfly velocity does not depend significantly on the value of $h$, provided that it is not too small, as shown in Fig. 6. Note that changing $g$ does not affect the light cone velocity, which is $\Delta s/\Delta T=1$ due to the geometry of the circuit that represents the Floquet time evolution. For $g=0$, an initial operator $Z_{0}$ remains localized on the same site. As we make $g$ larger, the butterfly velocity gradually increases and it reaches $v_{\text{B}}\approx v_{\text{LC}}$ for $g\approx 0.9$ with period time $T=1.6$ 555we note that this set of parameters, for which $v_{\text{B}}$ is near maximal, is close to the parameter choice found in Ref. Zhang et al., 2016 to exhibit very fast thermalization. Looking at Fig. 6 we notice that decreasing $v_{\text{B}}$ from its maximum corresponds to an increased front width at any given time. This is consistent with a coarse-grained random walk description, wherein increasing the probability of hopping to the right results in both a larger butterfly velocity and a suppression of the diffusion constant.

The most important evidence in support of a hydrodynamic description comes from examining the front width as a function of time. Similarly to the random circuit model, we find the wavefront of the operator spreading broadens as we go to longer times. To quantify the width by looking at the standard deviation of $\rho_{R}(s)=R(s)-R(s-1)$, i.e.

[TABLE]

As shown in Fig. 7, at long times the width grows algebraically in time as $\sigma(t)\propto t^{\alpha}$ with an exponent $0.5\lesssim\alpha\lesssim 0.6$. This is roughly consistent with the random walk description of operator spreading put forward in Sec. III. This diffusive broadening is expected to result in the strict inequality $v_{\text{E}}<v_{\text{B}}$ for the entanglement velocity, according to the arguments put forward in Sec. III.4. We confirmed numerically that this indeed holds in this model for various values of $g$.

Finally, one might wonder whether the above story continues to hold when we consider a system with energy conservation i.e., a time independent Hamiltonian. We have been informedChe that there is evidence of a diffusively growing front in a family of energy conserving ergodic spin-chains.

V Fractal Clifford circuits

In this section we compare the results obtained for the random circuit model of Sec. III to another set of circuits which do not in general exhibit energy conservation, namely Clifford circuits. We show that despite the fact that Clifford circuits can be ‘ergodic’ in certain senses – they can exhibit linear entanglement growth and correlations heating up to infinite temperature – both their spectrum and their OTOC behavior is anomalous and non-chaotic.

General Clifford circuits have a particularly simple structure to their operator spreading coefficients: Under time evolution by $t$ steps, a simple Pauli string becomes another Pauli string

[TABLE]

where $M_{t}$ is a linear endomorphism acting on the set of strings $\left(\mathbb{Z}_{q}^{\otimes 2}\right)^{\otimes L}$. Thus, a Pauli operator evolves to a single Pauli operator, rather than the superposition of Pauli operators allowed by Eq. (2). $M_{t}$ has to obey a number of constraints. In particular, time evolution should preserve the commutation relations amongst the Pauli stringsSchlingemann et al. (2008); Gütschow (2010); Gütschow et al. (2010). (Incidentally, using these constraints, it is possible to classify all translation invariant Clifford circuits into three types called fractal, glider, and periodicGütschow (2010); Gütschow et al. (2010).)

In line with the stringency of the constraint Eq. (31), it is unsurprising that Clifford circuits have pathological properties distinguishing them from more generic ergodic systems. Calculating the OTOC Eq. (III.2) for initial Pauli strings $\sigma^{\mu},\sigma^{\nu}$ for a Clifford circuit gives

[TABLE]

where $\theta$ is the phase obtained commuting $\sigma^{M_{t}(\mu)}$ through $\sigma^{\nu}$, and the final result is independent of the state $\psi$. For generic clifford circuits, this result shows persistent oscillations at late times, and does not settle to a specific limit. This is in contrast to the OTOC decay seen in Sec. III.2 and the OTOC decay expected in more generic thermalizing systems (see Sec. IV).

Additionally, Clifford circuits tend to have pathological spectral properties not associated with ergodic systems. For instance, one can prove that translation invariant Clifford circuits have exact recurrences $U^{t_{\text{rec}}}_{\text{f}}\propto 1$ on time scales linear in system size $t_{\text{rec}}=O(L)$ (see App. F). This directly implies that the eigenvalues of $U^{t_{\text{rec}}}_{\text{f}}$ are $t_{\text{rec}}$-th roots of unity, which in particular implies that the average level degeneracy is $O(2^{L}/L)$. This spectral structure does not exhibit the level repulsion we expect in systems with ETH and no apparent global symmetries.

Although operators obey the stringent condition Eq. (31), clifford circuits can still exhibit many ergodic properties usually associated with ‘ergodicity’. Indeed, it can be proved rigorously that the above mentioned ‘fractal’ clifford circuits exhibit: (a) Linear entanglement growth starting from certain so-called stabilizer initial states666Theorem IV.8 in Ref. Gütschow et al., 2010. and relatedly (b) starting from an initial product state, at long times all local observables tend towards their infinite temperature expectation values777Prop III.2 in Ref. Gütschow et al., 2010.. These fractal clifford circuits are periodic in time, in addition to being spatially translation invariant. An explicit example of such a circuit has $q=2$ and takes form $U(t)=U_{\text{f}}^{t}$ where

[TABLE]

Note that the resulting circuit has the geometry illustrated in Fig. 1 (with the one-site rotations merged into the 2-site gates), and the circuit elements repeat every $2$ layers. Operators evolve under this circuit in a particularly simple way:

[TABLE]

Note that for certain strings there is a possibility of cancelation, e.g. $Y_{s}Z_{s+1}\to Y_{s-1}X_{s}1\!\!1_{s+1}$.

Fig. 8 shows that for this particular circuit even more generic initial states show near linear entanglement growth, with a rate that is mostly independent of the initial state at long times. We can explain this using the operator spreading picture of entanglement growth, discussed in Sec. III.4. Looking at (34) we notice that a string with a Pauli operator $X$ or $Y$ at its right endpoint will keep moving to the right at a fixed speed of 1 site / 1 period forever, leading to the fixed rate linear entanglement growth seen for stabilizer states. The deviations from this behavior for a random product initial state come from strings that, up to time $t$, failed to start growing to the right. However, such operators need to have a very specific structure and their number is exponentially small in $t$ (going roughly as $4^{-t}$). Consequently, after the first few periods, most product states settle to the same entanglement velocity exhibited by stabilizer states. This is in contrast with the behavior of the operator entanglement of any Pauli string $\sigma^{\mu}(t)$ which remains 0 at all times due to Eq. (31).

In summary, Clifford circuits can be ‘ergodic’ in the sense that they can exhibit linear entanglement growth and thermalization of local observables, but do not have the spectral or OTOC behavior expected of generic ergodic systems. This suggests that linear entanglement growth is a rather coarse measure of quantum information spreading, sensitive only to the fact that operators tend to grow in extent over time. The OTOC, on the other hand, is sensitive to both the fact that operators grow in extent, and the fact that they become complicated superpositions of many Pauli strings. Note also that the entanglement entropy tests the behavior of a large ensemble of different Pauli strings, while the OTOC characterizes the evolution of a single initial Pauli operator, thus giving more detailed information on the dynamics.

VI Conclusions

We considered the spreading of quantum information in one-dimensional systems with local unitary time evolution but no conservation laws. Our key results are as follows. For random circuit systems we derived an exact hydrodynamical description of operator spreading. According to this description, operators grow into superpositions of Pauli strings which tend to be supported at a front, propagating with velocity $v_{\text{B}}$ where $v_{\text{B}}$ is a velocity scale distinct from the light cone velocity (see Fig. 2). An important consequence of the hydrodynamic equation is that the front itself undergoes a diffusive broadening in time. We proved that the velocity $v_{\text{B}}$ also determines the characteristic scale of change of the out-ot-time-order commutator (OTOC) between two local operators, while at very early times, long before the arrival of the front, we find a regime of exponential growth for the OTOC, with an exponent that depends on the initial separation of the two operators – it remains to be determined whether this early time exponential growth represents a quantum analogue of the Lyapunov behavior present in classical chaotic systems. The exact description of operator spreading in this model also allows us to give an exact formula for evolution of (typical) entanglement entropy across a cut, starting from an initial product state. We find that the entanglement grows with a third distinct velocity scale, the entanglement velocity $v_{\text{E}}$, which is strongly affected by the diffusive behavior, making $v_{\text{E}}<v_{\text{B}}$.

Comparing our exact results for the random circuit to an ergodic Floquet spin chain, we verified the presence of a diffusively broadening operator front in the kicked Ising model. This leads us to a tentative conjecture that this behavior is present for generic ergodic interacting quantum systems in 1D, evolving under local unitary evolution, allowing for a universal coarse grained hydrodynamic description in these systems. We contrasted this with the fine tuned behavior observed in Clifford circuits, for which operators can spread ballistically but do not become entangled superpositions of many operators. The ballistic spread of operators leads to the linear entanglement growth and the thermalization of local observables seen in certain Clifford circuits, but the fine tuned nature of the operator growth shows up in the pathological behavior of OTOCs. This demonstrates that linear entanglement growth and thermalization are not good predictors of the scrambling behavior captured by the OTOC.

We note that within existing field theoretic calculations Stanford (2016); Aleiner et al. (2016); Chowdhury and Swingle (2017) the broadening of the operator wavefront is not present, in contrast with our results. This discrepancy could be due to the unbounded local Hilbert spaces implicit in the field theoretic treatment; in any case, it should be resolved. Another possible direction for future work is finding the exact range of validity of the hydrodynamic operator spreading picture proposed here. It is known for example that the propagation of the OTOC becomes slower than ballistic in strongly disordered (but still thermalizing) systemsLuitz and Bar Lev (2017) and it would it is an interesting question whether this effect can be incorporated into some form of hydrodynamic description.

Related work:

Shortly before completing this manuscript, we became aware of related work888A. Nahum, S. Vijay, and J. Haah, Operator spreading in Random Unitary circuits., which should appear in the same arxiv posting. We have also been informed of forthcoming related work999V. Khemani et al, forthcoming..

Acknowledgements.

This work was supported by the Princeton Center for Theoretical Science (CvK), Research Unit FOR 1807 through grants no. PO 1370/2-1 (TR and FP), and DOE Grant No DE-SC/0016244 (SLS). We acknowledge useful conversations with Mark Mezei, Douglas Stanford and David Huse. SLS would like to thank Lauren McGough and Roderich Moessner for early discussions about OTOCs.

Appendix A Validity of Eq. (14)

It suffices to work in coarse grained co-ordinates $t,x$. Recall that

[TABLE]

[TABLE]

Our task is to justify the formula Eq. (14), which can be more carefully phrased as follows: In the limit, $t\rightarrow\infty$, with $|x-v_{\text{B}}t|/t=\varepsilon$ held to be a fixed nonzero number, the integrated operator density obeys

[TABLE]

where $c_{0,1}$ are positive numbers, bounded in the $t\rightarrow\infty$ limit, and we defined $\kappa\equiv x/t$ and work in units where $v_{\text{LC}}=1$. Define

[TABLE]

The quotient $Q_{x}$ is always positive. It is easy to verify that for $-1<\kappa<v_{\text{B}}$ the quotient is greater than 1, and an increasing function of $x$. On the other hand for $v_{\text{B}}<\kappa<1$, this quotient is less than $1$, and a decreasing function of $x$.

When $-1<\kappa<v_{\text{B}}$ we use these facts to bound

[TABLE]

On the other hand, it is immediate that $\overline{R}(x)\geq\overline{\rho_{R}}(x,t)$. Noting that in the large $t$ limit $Q_{x}=\frac{(1+v_{\text{B}})(1-\kappa)}{(1-v_{\text{B}})(1+\kappa)}$ we find

[TABLE]

Hence $c_{0}\leq 1+\frac{\left(1+\kappa\right)\left(1-v_{\text{B}}\right)}{2\varepsilon}$is an $O(1)$ constant.

Similarly we consider $v_{\text{B}}<\kappa<1$. Using $1=\sum_{y}\rho_{R}(y,t)$ it follows that $\overline{R}(x)=1-\sum_{y>x}\rho_{R}(y,t)$. In the present case $Q_{x}<1$ and is straightforward to derive a similar bound

[TABLE]

Hence $c_{1}\leq 1+\frac{\left(1-\kappa\right)\left(1+v_{\text{B}}\right)}{2\varepsilon}$ is an $O(1)$ number in the large $t$ limit. Note that near the edge of the light cone $\kappa=\mp 1$ respectively, the results Eq. (37) become increasingly exact as each of the bounds Eq. (A) and Eq. (43) become tighter. On the other hand, as we approach the front $\varepsilon\rightarrow 0$ the bounds become looser and Eq. (37) is less reliable – in this regime, the near front expansion Eq. (12) becomes more useful.

Appendix B Application to OTOC

Following Sec. III.2, the quantity $f(s,\tau)\equiv 1-\overline{\mathcal{C}}(s,\tau)$ can be written as

[TABLE]

where $e^{i\theta_{\mu,Z_{s}}}$ is a $q^{\text{th}}$ root of unity acquired by commuting $\sigma^{\mu}$ past $Z_{s}$. We can reparameterize the sum by summing over the left right endpoints of the Pauli string $\mu$

[TABLE]

where $h(l,r)$is simply the average $\overline{\left|c^{\mu}(\tau)\right|^{2}}$ for a Pauli string $\mu$ with left/right endpoint – recall that this value does not depend on the internal structure of $\mu$(see below eqnref or above ), only on the endpoints of $\mu$. Here $\mathfrak{F}(l,r)\equiv\{\mu:\text{supp}(\mu)=[l,r]\}$. Those intervals $\left[l,r\right]$ such that $s\in(l,r)$ do not on net contribute to this sum, because there results a sum over $q^{\text{th}}$ roots of unity disappears

[TABLE]

There are also contributions to (45) which arises when $s$ is on the left and/or right edge of an interval i.e., $s=l$ or $s=r$.

First, perform the sum over $\mu$ in (45) for the case $r=s>l$

[TABLE]

Next, perform the sum over $\mu$ in (45) for the case $l=s<r$

[TABLE]

Combining these sum identities Eq. (46), Eq. (47) back into Eq. (45) yields (after some rearrangement)

[TABLE]

The first line is

[TABLE]

It is readily verified that the second line disappears exponentially quickly in $\tau$ as $(2q/(q^{2}+1))^{2\tau}$, for any $s$ using Eq. (57). Hence, in the $\tau\rightarrow\infty$ fixed $s/\tau$ limit, we can approximate

[TABLE]

where we used the fact $\overline{L}(s+1)=\overline{R}(-s-2)$ and $\overline{\rho_{L}}(s+1)=\overline{\rho_{R}}(-s-2)$. For $0<s<\tau/v_{\text{LC}}$, and for any $q$, we can use Eq. (12) and Eq. (15) to argue that the second line is suppressed by a factor of $\sim q^{-s}$ relative to first line. Therefore, provided $\kappa>0$ in the $\tau\rightarrow\infty$ limit, the OTOC behaves as

[TABLE]

up to exponentially small corrections in $\tau$.

Appendix C Exact operator spread coefficients

In what follows we give an exact expression for the averaged operator spread coefficients, and a sketch of the derivation. We leave a more detailed derivation to future work. The averaged operator spread coefficients are defined as

[TABLE]

The Haar-averaging can be performed explicitly using the identity (62). After averaging, each two-site gate can be represented by a classical, Ising-like variable taking only two possible values. Due to the geometry of the circuit, these Ising variables form a triangular lattice. Eq. (52) becomes a classical partition function, i.e. a sum over all possible spin configurations. The partition function consists of two edge parts, which depend on the configurations on the first (last) layer and the Pauli strings $\mathbf{\nu}$ ($\mathbf{\mu}$), and a bulk part which is independent of the Pauli strings in question. Due to the Haar-averaging, the only information that remains in the partition function about the strings $\mathbf{\mu}$ and $\mathbf{\nu}$ is which sites they act on non-trivially. The bulk transfer matrix enforces a light cone structure on the spin variables. A light cone with velocity $v_{\text{LC}}=1$ emanates from the two endpoints of the string $\mathbf{\nu}$ such that all spins outside of the light cone have to point up (otherwise the configuration has zero weight in the partition function for $\overline{|c_{\mathbf{\mu}}^{\mathbf{\nu}}(t)|^{2}}$).

In the case of an initial Pauli operator acting on a single site the partition function for the operator spreading can be evaluated exactly. In this case the fact that $\mathbf{\mu}$ acts on one site only yields a boundary condition for the partition function wherein in the first row there is a single spin pointing down while all others point up. The partition function then becomes a sum over all possible ways this initially one-site domain can spread within the light cone, as shown in Fig. 9. Furthermore, the bulk interaction terms are only non-trivial at the boundary between the two domains and consequently give the same contribution for all domain configurations with the same depth. Thus the calculation simplifies to counting the possible domain configurations which can be done by considering it as a two-particle random walk for the two endpoints of the domain.

The calculation outlined above yields all average squared coefficients $\overline{|c_{\mathbf{\mu}}^{\mathbf{\nu}}(t)|^{2}}$ where $\mu_{j}=0$ if $j\neq 0$. In the following we simplify the notation by dropping the first index and denoting these as $\overline{|c^{\mathbf{\nu}}(t)|^{2}}$. The exact formula for these squared coefficients is given in Eq. (53). A surprising property of this formula only depends on the positions of the right and left endpoints $l,r$ of the Pauli string $\mathbf{\nu}$, and not on more detailed information concerning the internal structure of the string i.e., $\overline{|c^{\mathbf{\nu}}(t)|^{2}}=h(l,r)$. This is a consequence of the Haar averaging which in each step washes out all memory of the internal structure.

Form a circuit with even number $D$ layers. Number the spaces between the two site unitaries in the last layer $-D/2,1-D/2,\ldots D/2$. The support of an operator string $\mathbf{\nu}$ can be represented by $x,y\in\{-D/2,1-D/2,\ldots+D/2\}$ with $x<y$. The average square coefficient is obtained by plugging $x,y$ into the formula.

[TABLE]

Note that this expression depends only on $x,y$ and $q$ so we also denote it as $\overline{|c_{\mathbf{\nu}}^{\mu_{1}}|^{2}}=h(x,y,D)$ where we drop the $q$ dependence for simplicity. These expressions are of-course complicated. We note as an aside that this formula has a slightly neater expression in terms of hypergeometric functions.

C.1 Useful limits

Let us calculate the weight on an operator with endpoints $x,y$, in the large $D$ limit (in circuit co-ordinates) starting from Eq. (53). Re-express the weight as

[TABLE]

Defining:

[TABLE]

where

[TABLE]

we find that

[TABLE]

Putting this altogether, we have

[TABLE]

For fixed $x,y$ taking the large $D$ limit, the $v_{D}$ functions can be approximated as (using the same reasoning as in Eq. (14) ):

[TABLE]

Plugging in this approximation we find.

[TABLE]

So, fixing $x,y$, the weight on an operator with endpoints $x,y$ decays exponentially quickly $~{}\left(\frac{2q}{q^{2}+1}\right)^{2D}$ with large $D$.

Appendix D Long time correlations

In Sec. I we anticipated that random unitary circuits should ‘heat to infinite temperature’, much like ergodic Floquet circuits. Here we back up this claim by examining the long time behaviour of various correlations functions, demonstrating that they relax to their expected infinite temperature values. For simplicity we consider 1 and 2 point functions for specific 1-site Pauli operators of form $\sigma_{0}^{\alpha}$ – although most of the results below follow for more general operators as well. Fix any initial state $\omega$ and times $t_{2},t_{1}$. We find that

[TABLE]

This follows from two observations. First, the operator in the expectation value can be written $U_{01}U_{12}\sigma_{0}^{\alpha}U_{21}\sigma_{x}^{\beta}U_{10}$ where $U_{ij}=U^{-1}_{ji}$ is shorthand for the unitary evolutions between times $t_{i}$ and $t_{j}$. Now $U_{12}$ is statistically independent from $U_{01}$, so we can average over these disjoint circuits independently. Provided $t_{1}\neq t_{2}$ and $\alpha\neq 0$ it is straightforward to see that the Haar average $\overline{U_{12}\sigma_{0}^{\alpha}U_{21}}=0$. This can be re-expressed succinctly as $\overline{c^{\alpha}_{\mu}(t)}=0$ for all $t=t_{2}-t_{1}\neq 0,\alpha\neq 0$. Note that this result is independent of the initial state $\omega$ and the value $\beta$ – in particular we can recover the behavior of 1 point functions by setting $\beta=0$.

While the above correlation functions disappear on average, we may also quantify how their variance behaves at long times. Indeed, the variance decays exponentially in time, at least for random initial product states $\psi$. The variance is

[TABLE]

We argue that this variance vanishes as we increase $t\rightarrow\infty$ while fixing $t_{1}$. We can show this rigorously for an infinite temperature state in App. D.1.

For a random product state we have a less rigorous argument which proceeds as follows. First, given any $\epsilon>0$, it will be true that for sufficiently long times $t$, all of the weight of $\alpha$ is invested in strings $\mu$ with left/right endpoints $l<-(v_{\text{B}}-\epsilon)t,r>(v_{\text{B}}-\epsilon)t$ respectively, up to exponentially small corrections in $t$. These statements follow from Eq. (11). Second, as $\alpha$ is a 1-site operator, $\overline{|c^{\alpha}_{\mu}(t)|^{2}}$ only depends on the endpoints of $\mu$ rather than the detailed internal structure (See the discussion under Eq. (17)). Hence, up to exponentially small corrections in time, $\sigma_{0}^{\alpha}(t)$ is made up of an equal amplitude superposition of all operators $\mu$ with left/right endpoints near $\mp v_{\text{B}}t$ respectively. The vast majority of such strings contain an extensive $O(v_{\text{B}}t)$ number of Pauli operators. The expectation values of such strings on a random product states is exponentially decaying in the number of Pauli operators. As a typical $\mu$ string contains $O(v_{\text{B}}t)$ Pauli operators, we find $|\langle{\sigma^{\mu}\sigma_{z}^{\beta}}\rangle_{\omega}|^{2}\sim e^{-\zeta v_{\text{B}}t}$ for some constant $\zeta$. Hence for $t_{1}$=0, $\overline{\langle{\sigma_{0}^{\alpha}(t_{2})\sigma_{z}^{\beta}(0)}\rangle^{2}_{\psi}}\sim e^{-\zeta v_{\text{B}}t}$. For $t_{1}$ nonzero, we expect even more marked decay $\overline{\langle{\sigma_{0}^{\alpha}(t_{2})\sigma_{z}^{\beta}(t_{1})}\rangle^{2}_{\psi}}\sim e^{-\zeta v_{\text{B}}(t+t_{1})}$ because the support of $\sigma^{\mu}\sigma_{z}^{\beta}$ is further increased under time evolution.

D.1 Infinite temperature results

Consider the variance of the infinite temperature expectation value function

[TABLE]

This is a Haar averaged single site weight. In the large $t$, $z/t\rightarrow 0$ limit, we have an expression for this quantity (in the coarse grained lattice basis). It is approximately equal to Eq. (57) using $x=z/2,y=x+1,D=2t$. Asymptotically then the variance in the infinite temperature average decays exponentially as

[TABLE]

Appendix E Haar Identities

Behold the following Haar moments for $d\times d$ random unitary matrices.

[TABLE]

The following higher moment is useful for deriving the average square operator spread coefficients.

[TABLE]

This can be more elegantly expressed as a sum over elements of permutation group $S_{2}$

[TABLE]

Appendix F Recurrence times in Translation invariant Clifford quantum cellular

automata (CQCA)

The aim of this section is to show that translation invariant Clifford circuits have linear in system size recurrence times (for a certain family of system sizes). We utilize the technology and formalism of Ref. Gütschow, 2010.

Fact F.1**.**

(Ref. Gütschow, 2010 Theorem II.5) Clifford quantum circuits with translation symmetry (unit cell size 1) are in correspondence with the set of $2\times 2$ matrices $t$ with elements in $\mathbb{Z}/2\mathbb{Z}\left[u,u^{-1}\right]$ (polynomials in $u,u^{-1}$ over ring $\mathbb{Z}/2\mathbb{Z}$) obeying

[TABLE]

where by coprime, we mean the polynomials over ring $\mathbb{Z}/2\mathbb{Z}$ do not possess any common non-trivial factors. We say a clifford quantum circuit is centered if $a=0$.

Fact F.2**.**

From Ref. Gütschow, 2010 Proposition II.13). A centered CSCA $t$ is periodic with period c+2 if $\text{tr }t=c$ for$c\in\mathbb{Z}/2\mathbb{Z}$ (so the period is either $2$ or $3$).

Proof.

This is actually quite straightforward. From Fact F.1, we have $t^{2}=\text{tr}\left(t\right)t+1$. If $\text{tr}(t)=0$ then $t^{2}=1$. Else $t^{3}=t(t+1)=t+1+t=1$. ∎

Theorem 1**.**

Translation invariant centered CQCAs $U$ have linear in system size recurrence times (at most $t_{\text{rec}}=12L$), at least for system sizes $L=2^{n}$.

Proof.

Let $t$ be the $2\times 2$ matrix corresponding to clifford unitary $U$. By the Cayley Hamilton (CH) theorem (which holds for matrices over arbitrary rings), and the fact that $\det t=1$ we have

[TABLE]

where the matrices are written over the ring of Laurent polynomials with $\mathbb{Z}/2\mathbb{Z}$ coefficients i.e., $\mathbb{Z}/2\mathbb{Z}\left[u,u^{-1}\right]$. Recall from Fact F.1 that for a centered circuit, the elements of $t$ are symmetric Laurent polynomials. Note that the CH theorem can be iterated

[TABLE]

where for convenience we denote $p\equiv\text{tr}\left(t\right)$, which is of course a symmetric Laurent polynomial. Proceeding inductively, and squaring successive equations, one can show

[TABLE]

Now consider the trace of $t^{2^{k}}$:

[TABLE]

Now we can write

[TABLE]

where $c_{r}=c_{-r}\in\mathbb{Z}/2\mathbb{Z}$. Now iteratively square this expression. The cross terms disappear because the ring is $\mathbb{Z}/2\mathbb{Z}$

[TABLE]

For a system of size $L$, and periodic boundary conditions, the constraint $u^{L}=u^{0}=1$ is imposed on our polynomial ring. Hence, setting $L=2^{n}$ and $k=n$ we get

[TABLE]

But $c_{0}$ is just a constant. Using Fact F.2, and the fact $\text{tr}t^{2^{n}}$ a constant, we have $U^{L=2^{n}}$ is a periodic circuit. By this, Gütschow et al. mean that $t^{2^{n}(c_{0}+2)}=\text{1}$ is the identity matrix. This implies that $U^{2^{n}(c_{0}+2)}$ does not permute Pauli matrices – it only multiplies them by phases (which can only be $\pm 1$ in order to preserve generating relations for Pauli matrices). This in turn implies that $U^{2^{n+1}(c_{0}+2)}\propto 1$, the identity matrix on the many-body Hilbert space. Hence, as $c_{0}+2=2,3$ both divide $6$, we certainly have $U^{2^{n+1}\times 6}\propto 1$. In other words $U^{12L}\propto 1$ for $L=2^{n}$, and any clifford circuit obeying the conditions of this theorem. The upshot is that all such clifford circuits obeying the conditions of the theorem have a linear in system size recurrence time $t_{\text{rec}}\leq 12L$ , for $L=2^{n}$.∎

Corollary 1.1**.**

Translation invariant CQCAs have, on average, and exponential in system size level degeneracy.

Proof.

Straightforward. There are $2^{L}$ states in the Hilbert space. The linear in system size $\kappa L$ recurrence time (e.g., $\kappa=12$) means the eigenvalues are WLOG $\kappa L$-th roots of unity. Hence, the unitary has at most $\kappa L$ eigenvalues. Hence, the average level degeneracy is $2^{L}/\kappa L$ – exponentially large in system size. ∎