Weak Quantum Chaos

TL;DR

This paper introduces the concept of weak quantum chaos, characterized by polynomial growth of OTOC density in non-integrable quantum systems, and demonstrates its properties through numerical and analytical methods.

Contribution

It proposes the OTOC density as a more effective measure of quantum chaos in bounded local operators and establishes its polynomial growth in non-integrable systems.

Findings

OTOC density grows at most polynomially in time in non-integrable systems

In integrable cases, OTOC density saturates to a plateau

Numerical and analytical results support the weak quantum chaos concept

Abstract

Out-of-time-ordered correlation functions (OTOC's) are presently being extensively debated as quantifiers of dynamical chaos in interacting quantum many-body systems. We argue that in quantum spin and fermionic systems, where all local operators are bounded, an OTOC of local observables is bounded as well and thus its exponential growth is merely transient. As a better measure of quantum chaos in such systems, we propose, and study, the density of the OTOC of extensive sums of local observables, which can exhibit indefinite growth in the thermodynamic limit. We demonstrate this for the kicked quantum Ising model by using large-scale numerical results and an analytic solution in the integrable regime. In a generic case, we observe the growth of the OTOC density to be linear in time. We prove that this density in general, locally interacting, non-integrable quantum spin and fermionic dynamical systems exhibits growth that is at most polynomial in time---a phenomenon, which we term weak quantum chaos. In the special case of the model being integrable and the observables under consideration quadratic, the OTOC density saturates to a plateau.

Equations

$$C\left(x,t\right)=-⟨[w_x(t),v_0(0){]}^2{⟩}_{\beta },$$

$$C\left(x,t\right)\propto e^{{\lambda }_L\left(t-|x|/v_B\right)}.$$

$$C\left(x,t\right)\le 4{\Vert v\Vert }^2{\Vert w\Vert }^2.$$

$$C\left(x,t\right)\le 4{\Vert v\Vert }^2{\Vert w\Vert }^2{e}^{-\mu \max \{0,|x|-v_{LR}t\}}.$$

$$c^{(N)}(t):=-\frac{1}{N}\left({⟨[W(t),V(0){]}^2⟩}_{\beta }-{⟨[W(t),V(0)]⟩}_{\beta }^2\right).$$

$$c^{(N)}(t)\le A{t}^{3d},$$

$$H(t)=H_{\text{Ising}}+H_{\text{kick}}\sum _{n\in \mathbb{Z}}\delta \left(t-n\right),$$

$$U$$

$$w_{2j}=\left(\prod _{k<j}{\sigma }_k^z\right){\sigma }_j^x,$$

$$U=e^{-J{\sum }_j{w}_{2j-1}{w}_{2j}}{e}^{-h{\sum }_j{w}_{2j}{w}_{2j+1}}.$$

$$M_z=-\mathrm{i}\sum _{j\in \mathbb{Z}}{w}_{2j}{w}_{2j+1},$$

$$U$$

$$\cdot \prod _k\left(\cos \left(h\right)-w_{2k}{w}_{2k+1}\sin \left(h\right)\right)=U_{\text{Ising}}{U}_{\text{kick}}.$$

$$w(\theta )=\sum _j{w}_{2j}{e}^{\mathrm{i}\theta j},$$

$$\mathcal{U}(J,h,\theta )={\mathcal{U}}_{\text{kick}}(J,h,\theta ){\mathcal{U}}_{\text{Ising}}(J,h,\theta )$$

$$=\left(\begin{array}{cc}\cos (2h)& -\sin (2h)\\ \sin (2h)& \cos (2h)\end{array}\right)\left(\begin{array}{cc}\cos (2J)& e^{\mathrm{i}\theta }\sin (2J)\\ -e^{-\mathrm{i}\theta }\sin (2J)& \cos (2J)\end{array}\right).$$

$$\mathcal{U}(J,h,\theta )=V^{†}(J,h,\theta )\left(\begin{array}{cc}{e}^{\mathrm{i}\kappa (J,h,\theta )}& \\ & e^{-\mathrm{i}\kappa (J,h,\theta )}\end{array}\right)V(J,h,\theta ),$$

$$\kappa (J,h,\theta )=$$

$$+\cos (\theta )\sin (2J)\sin (2h)],$$

$$K_{ab}^{kj}(t):=⟨w_{2k+a-1}{w}_{2j+b-1}(t)⟩,$$

$$K^{kj}(t):=K^{j-k}(t)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\mathrm{d}\theta e^{-\mathrm{i}\theta \left(j-k\right)}{\mathcal{U}}^t(\theta ).$$

$$c_z(t)$$

$$\Vert \left[a(t),b\right]\Vert \le \xi \min \left\{\left|\text{supp}(a)\right|,\left|\text{supp}(b)\right|\right\}\Vert a\Vert \Vert b\Vert e^{-\mu \max \left\{0,\mathrm{d}\left(\text{supp}(a),\text{supp}(b)\right)-v_{LR}t\right\}}.$$

$${a|}_{\Gamma }:=\frac{{\text{tr}}_{{\Gamma }^C}a}{\text{tr}\text{mathbbm}{1}_{{\Gamma }^C}}\otimes \text{mathbbm}{1}_{{\Gamma }^C},$$

$$\Vert a(t)-{a(t)|}_{\Gamma }\Vert \le \xi \left|\text{supp}(a)\right|\Vert a\Vert e^{-\mu \max \left\{0,\mathrm{d}\left(\text{supp}(a),{\Gamma }^C\right)-v_{LR}t\right\}}.$$

$$\left|{⟨a,b⟩}_{\beta }^c\right|\le \chi \Vert a\Vert \Vert b\Vert e^{-\rho \mathrm{d}\left(\text{supp}(a),\text{supp}(b)\right)},$$

$${⟨a,b⟩}_{\beta }^c:={⟨ab⟩}_{\beta }-{⟨a⟩}_{\beta }{⟨b⟩}_{\beta }.$$

$$c(t)$$

$$c(t)\le \sum _{i,j,k\in \mathbb{Z}}\left|c_{ijk}(t)\right|,$$

$$c_{ijk}(t)$$

$$|{⟨a,b⟩}_{\beta }^c|\le |{⟨ab⟩}_{\beta }|+|{⟨a⟩}_{\beta }{⟨b⟩}_{\beta }|\le 2\Vert a\Vert \Vert b\Vert .$$

$$\Vert {a_x(t)|}_{\Gamma }\Vert$$

Full text

**Introduction.—**Quantum chaos was an active area of research in the 80’s and 90’s Gutzwiller (1990); Haake (2010); Stöckmann (2007). The main success of the field was a random matrix theory (RMT) classification of universal properties of quantum systems whose classical counterparts are chaotic. The classical limits of such systems have positive Lyapunov exponents, which characterise exponential sensitivity to initial conditions—the so-called butterfly effect. However, since the (classical) definition of the Lyapunov exponent is based on the concept of phase-space trajectories, one cannot unambiguously translate it to the quantum realm.

Nevertheless, it has been argued that a weaker property of dynamical mixing—a decay of almost all connected temporal correlators—is sufficient to establish universal quantum chaotic behaviour, such as random matrix statistics of energy spectra Andreev et al. (1996) or the universal exponential decay of Loschmidt echoes Gorin et al. (2006). In the theory of dynamical systems, complex (mixing) dynamics that displays no exponential butterly effect is referred to as weak chaos (see Ref. Klages and references therein). Examples of such dynamical systems include generic polygonal billiards in which nearby trajectories deviate only linearly with time, while correlation functions nevertheless exhibit mixing Casati and Prosen (1999, 2000).

The study of dynamical mixing (now called scrambling) and Lyapunov chaos in quantum mechanics was recently revived by the high-energy physics community, initially in the context of the propagation of information in black hole backgrounds Sekino and Susskind (2008). In 2014, Kitaev proposed to quantify chaos in interacting quantum many-body systems Kitaev (2014) in terms of the following out-of-time-ordered (four-point) correlation function (OTOC):

[TABLE]

where $w_{x}$, $v_{x}$ are local observables and $\langle\bullet\rangle_{\beta}$ denotes the thermal expectation value at inverse temperature $\beta$. The concept is based on a work by Larkin and Ovchinnikov Larkin and Ovchinnikov (1969) from 1969, where OTOC was connected to the instability of semi-classical trajectories of electrons scattered by impurities in a superconductor. Consequently, extended quantum systems were defined as chaotic if there exists a pair of local observables, $w$ and $v$, such that the OTOC (1) grows exponentially at early times Larkin and Ovchinnikov (1969); Maldacena et al. (2016):

[TABLE]

Motivated by the semi-classical picture, $\lambda_{L}$ is referred to as the Lyapunov exponent and $v_{B}$ the butterfly velocity.

A multitude of works examining the properties of quantum chaos have recently been written both from the high-energy perspective (typically in models with long-range interactions and in theories with holographic gravity duals) and from the condensed matter perspective (typically in experimentally more feasible models with local interaction) Shenker and Stanford (2014); Roberts et al. (2015); Maldacena et al. (2016); Polchinski and Rosenhaus (2016); Maldacena and Stanford (2016); Jensen (2016); Roberts and Swingle (2016); Blake (2016a, b); Blake and Donos (2016); Lucas and Steinberg (2016); Miyaji (2016); Stanford (2016); Swingle and Chowdhury (2016); Turiaci and Verlinde (2016); Aleiner et al. (2016); Patel and Sachdev (2016); Huang et al. (2016); He and Lu (2016); Fan et al. (2016); Chen et al. (2016); Hosur et al. (2016); Halpern (2016); Fu and Sachdev (2016); Swingle et al. (2016); Gaerttner et al. (2016); Shen et al. (2016); Chen (2016); Davison et al. (2016); Tsuji et al. (2017); Rozenbaum et al. (2016); Buijsman et al. (2016).

In this work, we investigate systems with local interactions with extensive number $N\to\infty$ of degrees of freedom, but with a finite local Hilbert space dimension $D$. In any model with a finite $D$ (including all fermionic and spin lattice models), in which local operators $u$, $v$ are bounded, the exponential growth in (2) can be bounded by operator norm inequalities (the triangular inequality, $\|ab\|\leq\|a\|\|b\|$ and ${\langle a\rangle}_{\beta}\leq\|a\|$):

[TABLE]

Thus, the OTOC can only grow exponentially up to a finite (scrambling) time $t^{*}$, after which it remains bounded by a constant. This is consistent with the observations made in other works on OTOC’s (of local observables) in fermionic systems where OTOC’s were always observed to reach a plateau Aleiner et al. (2016); Huang et al. (2016); He and Lu (2016); Fan et al. (2016); Chen et al. (2016); Hosur et al. (2016). As already noted in Patel and Sachdev (2016), the only way for the exponential time evolution to persist to late times is if there is a small prefactor multiplying the exponential function in (2). Even in the Sachdev-Ye-Kitaev (SYK) model with long-range interactions, this prefactor is $1/N$, which becomes small as $N\to\infty$ Polchinski and Rosenhaus (2016). Exponential growth (2) of the OTOC is therefore at best a transient effect in systems of interest to this work.

If interactions are local, $C(x,t)$ can be further bounded by the Lieb-Robinson theorem (LRT) Lieb and Robinson (1972) (see also Roberts and Swingle (2016)):

[TABLE]

In this case, for $t\ll t^{*}=\left|x\right|/v_{\!L\!R}$, the OTOC is even more suppressed. The interpretation of this effect is clear: namely, $t^{*}$ is the time in which $C(x,t)$ enters the causal cone. Before $t^{*}$, $C(x,t)$ is almost zero, while after $t^{*}$, it is bounded by (3) and saturates at a plateau. The dynamics can only be non-trivial near the edge of the causal-cone (or for $t\sim t^{*}$), where $C(x,t)$ can vary greatly. This is consistent with Aleiner et al. (2016); Huang et al. (2016).

Another important fact is that momentum operators—the observables that Ref. Larkin and Ovchinnikov (1969) originally used to compute the Lyapunov exponent of the semiclassical trajectories—are unbounded. Therefore, if we wanted to preserve the semiclassical justification of the OTOC, which is necessary to be able to speak about quantum chaos, the quantum observables under consideration must have unbounded spectra.

These observations can be summarised in the intuitive statement that if chaos is to fully develop over long time, the observables have to provide enough “space” for this to happen; they need to be unbounded. Indeed, this is the case with general observables in bosonic systems (usually studied in holography). However, this condition is not fulfilled by local observables in fermionic or spin systems, or more generally, in systems with a finite $D$. On the other hand, extensive observables in such theories do satisfy the unbounded spectrum criterium and therefore have the capacity to fully unveil the system’s dynamical properties and quantum chaos. Motivated by this fact, we propose a new measure of quantum chaos: the density of the OTOC (dOTOC) of (non-local) extensive operators $V\equiv\sum_{x\in\Lambda}v_{x}$, $W\equiv\sum_{x\in\Lambda}w_{x}$, with $w_{x},v_{x}$ local. It is defined on a $d-$dimensional lattice $\Lambda$ with $N$ sites as the centralised second moment of the commutator

[TABLE]

The disconnected part, which is just the square of the standard dynamical susceptibility (i.e. the response function), has been subtracted to make the dOTOC well defined in the thermodynamic limit (TL) for any temperature. Because of the cyclicity of the trace, this term vanishes at $\beta=0$ (this will occur in the model that we study below). Using the LRT and the clustering property of thermal states, which holds for any temperature in $d=1$ Araki (1969) and for sufficiently high temperature in $d>1$ Kliesch et al. (2014), in Appendix A, we rigorously prove that the dOTOC satisfies a uniform (in $N$) polynomial bound

[TABLE]

where $A$ is an $(N,t)-$independent constant. The same bound equally holds in the TL, $c(t):=\lim_{N\to\infty}c^{(N)}(t)$.

Moreover, we report below the results of extensive numerical and analytical calculations, which demonstrate that possibly the simplest non-trivial locally interacting quantum chaotic spin system: the kicked Ising (KI) quantum spin chain Prosen (2000, 2002), exhibits linear growth of the dOTOC of extensive magnetisation observables, $c(t)\propto t$. An exception is the integrable KI model (equivalent to a free fermion model), for which we show analytically that its dOTOC of extensive quadratic observables (in fermionic variables) saturates, $c(t\to\infty)={\rm const}$. Since the KI model seems to be generic, we further conjecture that the bound (6) is not optimal and that typical one-dimensional, non-integrable and locally interacting models exhibit linear growth of dOTOC’s.

As a consequence, theories under consideration in this work are not expected to exhibit any late-time butterfly effect, but as we know from results in the RMT, can still be chaotic. In reference to classical mixing systems without the butterfly effect, we term the phenomenon of infinite polynomial growth of dOTOC’s weak quantum chaos.

**Kicked quantum Ising model.—**The Hamiltonian of the one-dimensional KI model consists of the Ising-interaction term $H_{\text{Ising}}=\sum_{j}J\sigma_{j}^{x}\sigma_{j+1}^{x}$ and the kick term $H_{\text{kick}}=\sum_{j}h\left(\sigma_{j}^{z}\cos\varphi+\sigma_{j}^{x}\sin\varphi\right)$:

[TABLE]

where $\sigma^{\alpha}_{j}$ are local Pauli spin operators. The model has three parameters: the Ising coupling $J$, the magnitude of the external magnetic field $h$ and the inclination of the external magnetic field $\varphi$. KI is a periodic (in time) system with the Floquet propagator:

[TABLE]

Because of the temporal periodicity, KI dynamics can be viewed as discrete in time, or as a quantum cellular automaton. The effect of a perturbation on a single lattice site propagates in a causal-cone with speed $1$. Namely, information can spread only by one site, left or right, within one period (kick of the magnetic field). Random matrix analysis Pineda and Prosen (2007); Kukuljan and Prosen (2016) revealed that KI is chaotic.

The system has a further nice property of being integrable (quasi-free) for transverse magnetic field, $\varphi=0$, and non-integrable (and interacting) for $\varphi>0$. Thus, $\varphi$ serves as a handy parameter which allows us to study integrability breaking. See e.g. Prosen (2002, 2007) for a survey of elementary dynamical properties of the KI model.

Here, we study the KI chain with $N$ spins and evaluate the dOTOC (5) $c^{(N)}_{\alpha}(t)$ for a (non-local) extensive magnetisation, $W=V=M_{\alpha}=\sum_{j=1}^{N}\sigma_{j}^{\alpha}$, which can either be transverse ($\alpha=z$) or parallel ($\alpha=x$) to the direction of the Ising interaction. We take $\beta=0$ as an infinite-temperature Gibbs ensemble is the only meaningful equilibrium state for periodically driven systems, which generically heat up to infinite temperature. We use three different approaches, two numerical methods for the general inclination ($0\leq\varphi\leq\frac{\pi}{2}$) and an analytical solution for the transverse field case $\varphi=0$. In the first, appropriate for small system sizes (up to $N\sim 12$), we used the exact numerical Floquet operator (8). The second method, used for intermediate system sizes (up to $N\sim 22$), was a Monte-Carlo wave-function sampling based on typicality arguments (explained in Appendix B). The analytical solution in the TL for the integrable (transverse) case and transverse magnetisation $M_{z}$, was found using fermionisation. We outline the main steps for obtaining the analytical solution in what is to follow.

**Analytical solution.—**For the transverse field ($\varphi=0$), KI is a quasi-free model. If, furthermore, the (extensive) observable of interest is simple enough, the dOTOC allows for an analytic solution in terms of Jordan-Wigner transformation of Pauli spins into staggered Majorana fermion operators

[TABLE]

obeying the anti-commutation relations $\left\{w_{i},w_{j}\right\}=2\delta_{ij}$. The Floquet operator (8) then takes the following form:

[TABLE]

It is clear from (10) that the KI model is free for $\varphi=0$ 111We note that even free theories can exhibit complicated entangled collective behaviour when one considers the dynamics of composite operators. See e.g. Grozdanov and Polonyi (2015); Polonyi (2015).. Now, the transverse magnetisation can be expressed as a sum of quadratic Majorana operators:

[TABLE]

which enables the analytic computation of the dOTOC of $M_{z}$ 222Longitudinal magnetisation is instead a sum of infinite strings of fermions, $M_{x}=\sum_{j}\prod_{k<j}\left(-{\rm i}w_{2k}w_{2k+1}\right)w_{2j}$. Computation of the OTOC for $M_{x}$ is therefore significantly more involved.. Power-expanding the Floquet operator (10) and using $\left(w_{i}w_{j}\right)^{2}=-1$ for $i\neq j$, $U$ further simplifies to

[TABLE]

Since the transverse field model is free, it is convenient to work in the Fourier transformed Majorana basis:

[TABLE]

with shorthand notation $\underline{w}(\theta)=\begin{pmatrix}w(\theta)\\ w^{\prime}(\theta)\end{pmatrix}.$ One can show (Appendix C) that the Floquet propagator in the Heisenberg picture, $\mathscr{U}\underline{w}(\theta):=\begin{pmatrix}U^{\dagger}w(\theta)U\\ U^{\dagger}w^{\prime}(\theta)U\end{pmatrix}$, takes the following form in Fourier transformed Majorana basis:

[TABLE]

This $2\times 2$ unitary matrix valued symbol can be diagonalised as:

[TABLE]

where

[TABLE]

and $V(J,h,\theta)$ is given explicitly in Appendix D.

Knowing that the KI Majorana fermions in the Fourier basis time evolve as $\underline{w}(\theta,t)=\mathscr{U}(\theta)^{t}\underline{w}(\theta,0)$ allows us to define the real space propagator as:

[TABLE]

for $a,b\in\left\{1,2\right\}$. $K^{kj}$ can then be computed from $\mathscr{U}(\theta)$ (Appendix E):

[TABLE]

Using the propagator (18), we can compute the infinite temperature OTOC of the transverse magnetisation, $c^{(N)}_{z}(t)$. First, we express the terms in (5) using (11), e.g. ${\langle\sigma_{i}^{z}(t)\sigma_{j}^{z}\sigma_{k}^{z}(t)\sigma_{l}^{z}\rangle}$ as an eight-fermion expectation value ${\langle w_{2i}(t)w_{2i+1}(t)w_{2j}w_{2j+1}w_{2k}(t)w_{2k+1}(t)w_{2l}w_{2l+1}\rangle}$. Then, using $\eqref{eq:NotationSymbol}$, these are expressed as the product of four propagators (one for each time-dependent fermion) times an equal-time eight-fermion expectation value, with terms summed over four spatial and spin indices (see Appendix F for details). Simple algebraic manipulations then lead to the final expression for the dOTOC in the TL:

[TABLE]

where we used the following notation: $R_{1}:=(l_{1}-j,l_{1})$, $R_{3}:=(l_{3}-j,l_{3})$, $S:=(s_{j},s_{0}),$ together with the notation $v=(v(1),v(2))$ for vector components and $\tilde{1}:=2$, $\tilde{2}:=1$. We can use the formula (19) in two different ways. For intermediate times $t\sim 50$, we can perform the integral in (18) exactly and evaluate the sums in (19), which, because of the causal-cone spreading of information, now become finite sums (see Appendix G.1 for details).

Secondly, we can use the stationary phase approximation in combined (18), (15) and (19) to compute the large-$t$ asymptotic behaviour of the dOTOC. In this way, we prove that for large times, $c_{z}(t)$ is a constant (dependent only on $J$ and $h$). In other words, the dOTOC of quadratic extensive observables in the integrable KI model saturates to a plateau. Details are explained in Appendix G.2.

**Results and discussion.—**In summary, we observe two distinct behaviours of the OTOC density for extensive observables in a one-dimensional KI model. For a generic situation, unless the model is integrable and the observable quadratic, the extensive dOTOC grows linearly with time. In fact, numerical results for finite system sizes saturate to a plateau at $t\sim N/2$, but this is simply due to a finite size effect—a consequence of the causal cone coming around the periodic boundary. This plateau grows with an increasing system size $N$ and we expect that it disappears in the TL $N\rightarrow\infty$. In the regime where the model is integrable (free) and the observable is simple (quadratic in fermion operators), the dOTOC saturates to a genuine plateau despite the fact that the spectrum of the observable is unbounded. The latter statement was proven in this work by finding an explicit analytic solution for $c_{z}(t)$ from which the expression for the height of the plateau for a given set of parameters $J$ and $h$ could be found. The results of the time dependence of the extensive dOTOC for different scenarios are presented in Figure 1 and explained in the caption.

For the integrable case with $\varphi=0$, the quasiparticle spectral gap closes on the line of $J=h$ in the parameter space and the system exhibits a Floquet analogue of a quantum phase transition, i.e. $\kappa\left(J=h,\theta=0\right)=-\kappa\left(J=h,\theta=0\right)=0$ (cf. (15), (16)). It is interesting to ask whether the OTOC also reflects this transition in any way. What we find is that the plateau height ceases to be smooth for $J=h$. Beyond that, we also checked the slope of the OTOC for longitudinal magnetisation $M_{x}$ in the vicinity of this line. It turns out that the slope exhibits a peak, but not exactly on the line $J=h$. This could be the effect of a small system size, which was necessary for numerics. It is plausible that the peak may align with $J=h$ in the TL.

This work should be considered as a starting point for future investigations of quantum, weakly chaotic systems, which exhibit dynamical late-time mixing but do not display any exponential butterfly effect due locality of interactions and finiteness of the local Hilbert space. In such systems, the standard OTOC rapidly plateaus and is therefore not a good measure of chaos. This observation led us to propose of a new measure of chaos: density of the OTOC of non-local extensive operators. We have proven (Appendix A) that such correlators always exhibit a polynomial bound and can thus be widely used to diagnose and classify quantum chaos. In the case of the non-integrable KI model studied here, the growth is linear. Intuitively, it seems apparent that in locally interacting systems, information propagates slower than in an all-to-all interacting theory like the SYK model. The speed is limited by the Lieb-Robinson velocity. However, what is less apparent is that such systems can still be chaotic; a result established by an RMT analysis Pineda and Prosen (2007).

Lastly, we note that in order to study chaos in strongly coupled, large-$N$ theories (even in those that do exhibit the buttery effect), it would be interesting to extend holographic calculations to computations of OTOC’s of non-local, smeared operators. For detailed future analyses, we will likely need to utilise the full machinery of holographic $n$-point function calculations Barnes et al. (2010); Arnold and Vaman (2011); Grozdanov and Starinets (2015) that will extend beyond studying gravitational shock waves Shenker and Stanford (2014); Roberts et al. (2015).

Acknowledgements.— We thank M. Medenjak for fruitful discussions and E. Ilievski for useful comments. The work has been supported by the ERC grant OMNES and grants P1-0025, P1-0044 of Slovenian Research Agency (ARRS). S. G. is supported in part by a VICI grant of the Netherlands Organization for Scientific Research (NWO) and by the Netherlands Organization for Scientific Research/Ministry of Science and Education (NWO/OCW).

Appendix A Proof of the polynomial bound on the density of extensive OTOC

In this section, we first prove that the density of the extensive OTOC (dOTOC) for 1D locally interacting translationally invariant lattice systems with finite local Hilbert space dimension cannot grow faster than with the third power of time. Then, we directly extend our theorem, Eq. (6) of the paper, to $d-$dimensional regular lattices.

To derive the bound, we will take advantage of two important theorems that hold for locally interacting lattice systems. The first, the Lieb-Robinson theorem (LRT) Lieb and Robinson (1972) states that for any locally interacting lattice system there exist positive constants $\xi$, $\mu$ and $v_{LR}$, such that for any two operators $a$ and $b$:

[TABLE]

Here, $\text{supp}(a)\subset\mathbb{Z}$ denotes the support of a local operator $a$ and ${\rm d}\left(\bullet,\bullet\right)$ is the distance between two sets. Roughly speaking, the theorem says that the commutator of two local observables grows in a causal-cone, spreading with velocity $v_{\!L\!R}$.

The authors of Ref. Bravyi et al. (2006) have found an elegant and useful reformulation of the LRT. Let $\Gamma$ be a subset of the lattice of $N$ sites and define

[TABLE]

where $\Gamma^{C}$ denotes the set complement, to be a projection of the operator $a$ on the sublattice $\Gamma$. Note that $\text{supp}(a|_{\Gamma})=\Gamma$. Then, for a given locally interacting system, the LRT is equivalent to Bravyi et al. (2006); Ilievski and Prosen (2013)

[TABLE]

The second theorem that we will need is the exponential clustering property of thermal states Araki (1969); Matsui (2003). For a thermal state of a one-dimensional locally interacting system, there exist positive constants $\chi$ and $\rho$, such that the following inequality is satisfied by any two operators $a$ and $b$:

[TABLE]

where, in order to make the expressions in this section more compact, we have introduced the notation for the connected (bipartite) correlation function:

[TABLE]

An analogous result is true for locally interacting Hamiltonians on arbitrary $d-$dimensional lattices for sufficiently high temperatures Kliesch et al. (2014). As we will show, the three bounds, (20), (22) and (23), imply a polynomial bound for the dOTOC.

Our goal is to compute an upper bound on the dOTOC:

[TABLE]

Note that in the third line, we used the translational invariance of the system.

The bound can be established by first using the triangular inequality

[TABLE]

and then by finding the appropriate bounds for individual terms. To take advantage of the exponential clustering property, at every time $t$, we will separate the $i$–$k$ plane into two domains. The first, $\left|k-i\right|\leq 2v_{\!L\!R}t$, is the region where the causal-cones of the two commutators in (25) are overlapping. There, the exponential clustering cannot be used, but the region is bounded in $\left|k-i\right|$ which will yield a finite contribution to the upper bound. In the second region, $\left|k-i\right|>2v_{\!L\!R}t$, the causal-cones are well separated so we will be able to use the exponential clustering to produce a finite upper bound. The contributions from both regions will be summed up in the end to get the overall upper bound on the OTOC. We will treat $c_{ijk}(t)$ slightly differently in the two regions:

[TABLE]

In the proof, the following obvious bound on the connected correlation functions will be useful,

[TABLE]

We will also need a rigorous estimate for the norm of a projected operator $\left.a\right|_{\Gamma}=a-\left(a-\left.a\right|_{\Gamma}\right)$:

[TABLE]

In the first line, we used the triangular inequality. In the second line, we utilised the fact that unitary time evolution preserves the norm, followed by an application of the second version of the LRT, Eq. (22). In the last line, we used the fact that the exponential of a non-positive function can be bound by 1.

For $\left|k-i\right|\leq 2v_{\!L\!R}t$ :

[TABLE]

where we used Eq. (28) to bound the connected correlator and the LRT (20) to bound the norms of the commutators. We can now sum over the corresponding domain to get the contribution to the overall upper bound:

[TABLE]

In the second line, we used the fact that each individual sum is independent of all other coefficients.

For $\left|k-i\right|>2v_{\!L\!R}t$ :

We can now use the linearity of commutators and thermal expectation values to write $c_{ijk}(t)$ as a sum of four terms and then bound it from above, again using the triangular inequality:

[TABLE]

Let us now find the upper bounds for each of the terms individually.

Term $\left(I\right)$: Since the supports of the two commutators are well separated (the first commutator is different from identity only on $\Gamma_{i}$, the second one only on $\Gamma_{k}$), we can use the exponential clustering property of thermal states to bound the term. This is the only place in the proof where this property of thermal states is used. However, here, it is indeed crucial:

[TABLE]

In the first line, we used the exponential clustering property. In the second line, we used the LRT together with (29) and the fact that the commutator is non-zero only if $v_{0}$ and $v_{j}$ are located inside the supports of $w_{i}(t)|_{\Gamma_{i}}$ and $w_{k}(t)|_{\Gamma_{k}}$, respectively. In the third line, we used $\Theta\left(\bullet\right)\leq 1$, where $\Theta$ is defined as $\Theta(\text{true})=1$, $\Theta(\text{false})=0$.

Note that if one wanted to compute the density of the disconnected OTOC, this bound would still be valid, but only in the infinite temperature regime $\beta=0$, where the expectation values of the commutators vanish because of the cyclicity of the trace. At finite temperature, an estimate obtained using exponential clustering gives a divergent contribution upon summation over $\sum_{\left|k-i\right|>2v_{\!L\!R}t}\sum_{j}$, indicating that the disconnected dOTOC is generically not a well defined quantity in the thermodynamic limit.

Term $\left(II\right)$: Here, we will use (28) to bound the term. To obtain the bound, which will give a non-divergent contribution when summed over $\sum_{\left|k-i\right|>2v_{\!L\!R}t}\sum_{j}$, we will take advantage of a convenient fact that the first commutator can be bound in two different ways—using two different versions of the LRT. One version will give us exponential damping when $\left|k-i\right|$ grows to infinity, the other version when $\left|i-j\right|$ grows to infinity. For a given combination of $i,j,k$, we then take the minimum of the two bounds, which results in a convergent bound upon summation over the domain:

[TABLE]

Since we only want to prove that the term will contribute to the upper bound no more (no faster) than polynomially in time, we were allowed to make some of the terms in the third line larger by a constant factor. In the fourth line, we used the fact that the functions appearing in the exponent next to $-\mu$ are non-negative.

Term $\left(III\right)$: In analogy with the previous term:

[TABLE]

Term $\left(IV\right)$: Writing again all possible combinations of different versions of the LRT and taking the minimum:

[TABLE]

In the fourth line, we used the fact that $\min$ cannot decrease if we simply omit a couple of (non-negative) functions and if we multiply some of the remaining functions by a constant factor. We have taken the common term of the two remaining functions out of the minimum.

We now have the estimates for all of the four terms so we are ready to sum them over the domain to get the contribution to the bound for the dOTOC. Since the setting is reflection symmetric (upon exchanging $i$ and $k$), it is enough to compute (all the terms) for $k>i+2v_{\!L\!R}t$ and double the result. In the case of $k>i+2v_{\!L\!R}t$, the subsets are:

[TABLE]

and the distances:

[TABLE]

We can plug these into the expressions (33-78) and evaluate the sums (all in the form of geometric series). We find that the contribution to the upper bound on the dOTOC coming from the region $\left|k-i\right|>2v_{\!L\!R}t$ is of order of $t^{2}$:

[TABLE]

since each of the bounds gives an $\mathscr{O}\left(t^{2}\right)$ contribution upon the summation.

By adding this that to the result for $\left|k-i\right|\leq 2v_{\!L\!R}t$, we arrive to the end of the proof. Hence, we have established that the dOTOC cannot grow faster than with the third power of time:

[TABLE]

This result can be straightforwardly extended in two ways:

• •

Even without taking the thermodynamic limit in (25), we can still find a bound on $c^{(N)}(t)$ by using exactly the same formal steps. We only have to assume that the finite $N$ lattice is periodic, so translational invariance can be used.

• •

We may consider any regular $d$-dimensional lattice in the regime where the temperature is sufficiently high for the generalisation of the exponential clustering property to hold Kliesch et al. (2014). In this case, each summation over a positional index (with an appropriate constraint) yields a factor that scales as ${\cal O}(t^{d})$, rather than ${\cal O}(t)$. With these results in hand, we finally arrive at the general polynomial bound stated in Eq. (6) of the main text.

Appendix B Numerical evaluation of OTOC based on typicality

The approximative numerical method for evaluating the OTOC for intermediate system sizes, $N\sim 22$, that we used is based on Levy’s lemma (also referred to as the measure concentration, or typicality). The lemma states, roughly, that in a large enough Hilbert space, the expectation value of a well-behaved observable on a single randomly chosen quantum state will be exponentially close in probability to the ensemble average of the observable. That is, in a large enough Hilbert space, almost any state is typical. For a precise formulation, see for example Refs. Gerken (2013); Müller et al. (2011). Here, we will approximate ensemble averages by averaging over a set $\left\{\left|\Psi_{\text{rand}}\right\rangle\right\}$ of random states in the Hilbert space. In this case, for an observable A, typicality arguments lead to

[TABLE]

where $|{\cal S}|$ denotes the cardinality of the set ${\cal S}$ (i.e. the number of random states used in the calculation). Rather than estimating the error of such an approximation by analytical arguments, we will estimate it numerically by computing variances.

The numerical method for computing the OTOC is then constructed as follows. We generate two sets $\left\{\left|\Psi_{1}\right\rangle\right\}$ and $\left\{\left|\Psi_{2}\right\rangle\right\}$ of random (normalised) vectors in the $2^{N}$ dimensional Hilbert space. Then, we can compute

[TABLE]

In the second line, we inserted a partition of unity and approximated it by $\mathds{1}=\sum_{\left|\Psi\right\rangle\in\mathscr{H}}\left|\Psi\right\rangle\left\langle\Psi\right|\approx\frac{2^{N}}{\left|\left\{\left|\Psi\right\rangle\right\}\right|}\sum_{\left\{\left|\Psi\right\rangle\right\}}\left|\Psi\right\rangle\left\langle\Psi\right|$. In the third line, we defined $|\widetilde{\Psi}\rangle:=V\left|\Psi\right\rangle$ and switched from the Heisenberg to the Schrödinger picture. A similar expression can be derived for the other term appearing in the OTOC, namely $\left\langle VW(t)W(t)V\right\rangle_{\beta=0}=\left\langle V^{2}W^{2}(t)\right\rangle_{\beta=0}$.

In evaluating the dynamics of $|\Psi_{1,2}(t)\rangle$, the computation of the action of one-spin (or two-spin) unitary operators on a vector is numerically very efficient. This is because the operators act only on two (or four) among the composite spin indices of the vector at the time. Thus, the operation only requires $2\cdot 2^{N}$ (or $4\cdot 2^{N}$) computational steps. The evaluation of the expression (84) is composed only of one-spin and two-spin operations and can therefore be computed within $\mathscr{O}\left(N2^{N}\right)$ computational steps. For a detailed discussion of such an algorithm, see for example Prosen (2007).

Finally, the statistical error of such an approximation is estimated by

[TABLE]

where $\sigma_{c^{(N)}}$ is the variance of $\left\{c^{(N)}_{\left|\Psi_{1}\right\rangle,\left|\Psi_{2}\right\rangle}\right\}$ if $c^{(N)}_{\left|\Psi_{1}\right\rangle,\left|\Psi_{2}\right\rangle}$ is an expression like RHS of Eq. (84) for a fixed state pair $\left|\Psi_{1}\right\rangle,\left|\Psi_{2}\right\rangle$ (omitting the averaging).

Appendix C Kicked quantum Ising propagator in the Fourier basis

Propagators of translationally invariant free models can be simplified by writing them in the Fourier transformed basis of Majorana fermions. This is useful for computing the action of powers of the propagator and therefore obtaining real time dynamics.

The Fourier transformed Majorana fermions can be written as

[TABLE]

In most of the calculations here, we can safely assume that the chain is infinite, $N=\infty$, and hence, $\theta\in[-\pi,\pi)$ is a continuous quasi-momentum parameter. It is convenient to introduce the shorthand spinor notation

[TABLE]

Written in the Heisenberg picture, the propagator in this basis acts as a $2\times 2$ unitary matrix

[TABLE]

As an example, let us explicitly compute the expression for the propagator of the transverse KI field ($\varphi=0$) in the Fourier basis:

[TABLE]

The kick term acts as

[TABLE]

where in the first line, we used the definition of the propagator in the Heisenberg picture (88). In the second line, we employed the definition of the Fourier transformed Majorana fermions (86). Then, in the third line, we used the fact that $w_{i}$ has a non-trivial product only with terms containing $w_{i}$ (this follows from the Majorana anti-commutation relations, $\left\{w_{i},w_{j}\right\}=2\delta_{ij}$) and finally, in the last line, we again utilised the definition of the Fourier transform.

Similarly, for the Ising propagator,

[TABLE]

We note that in the second line, we shifted the summation index.

Then, for the action of the entire Floquet propagator (note the correct order), we recover Eq. (14) from the main text:

[TABLE]

Appendix D Spectrum of the kicked quantum Ising model

The Floquet propagator $\mathscr{U}$ of a general quadratic model in the Fourier basis is a $2\times 2$ unitary matrix. We can diagonalise it to the following form:

[TABLE]

where $V$ is a unitary eigenvector matrix with elements

[TABLE]

Diagonalising (94) and using $\arccos\left(z\right)=-{\rm i}\ln\left(z+\sqrt{z^{2}-1}\right)$, we find the eigenvalues to be

[TABLE]

We see that the Floquet quasiparticle dispersion relation $\kappa(J,h,\theta)$ has three extrema: two maxima at $\theta=\pm\pi$ and a minimum at $\theta=0$. For $J\neq h$, we have $\kappa\left(J\neq h,\theta\right)>0$ so the system has a spectral gap (see Fig. 3). For $J=h$, $\kappa\left(J=h,\theta=0\right)=-\kappa\left(J=h,\theta=0\right)=0$, and so the gap closes. This is a Floquet-type analogue of a quantum critical line in the $J$–$h$ plane. On the $J=h$ line, the eigenphase $\kappa\left(J=h,\theta\right)$ has only two extrema at $\theta=\pm\pi$.

For the eigenvectors, we find

[TABLE]

where

[TABLE]

At the extremal points, $\theta\in\left\{-\pi,0,\pi\right\}$, the eigenvectors simplify to

[TABLE]

where

[TABLE]

Note that $\Sigma_{0/\pm\pi}^{2}=1$.

Appendix E Real space propagator

The propagator of a general translationally invariant Floquet system, quadratic in fermionic operators (that is, a free model), can be expressed in the Fourier transformed basis as a unitary $2\times 2$ matrix that depends on quasi-momentum. The same formalism applies for general time-dependent situation through the application of the Trotter-Suzuki formula. For example, we could obtain the results for the time-independent transverse field Ising model by setting $J\to dt\,J$, $h\to dt\,h$ and then taking the limit of $dt\to 0$.

By using the expression for $\mathscr{U}(\theta)$ and

[TABLE]

we compute the time-evolution of the Majorana fermions in the spatial basis; that is, the real-space propagator $K$, defined by

[TABLE]

or equivalently,

[TABLE]

for $a,b\in\left\{1,2\right\}$.

Let us start from Eq. (107). Then,

[TABLE]

At the same time, by definition (and by using Eq. (86)), the above expression is equal to

[TABLE]

Using the inverse Fourier transform

[TABLE]

along with the definition (108), the equation (111) can then be further expressed as

[TABLE]

Now, taking into account the translational invariance of the system, i.e. $K^{kj}=:K^{j-k}$, and thereby introducing a new index $l:=j-k$, the expression becomes

[TABLE]

Summing over $k$ and using $\sum_{n}e^{{\rm i}nx}=2\pi\sum_{k}\delta\left(x-2\pi k\right)$, we have

[TABLE]

Finally, by integrating over $\varphi$, we get

[TABLE]

Comparing Eqs. (110) and (116), we see that

[TABLE]

or equivalently, by performing the inverse Fourier transform

[TABLE]

Appendix F The dOTOC of transverse magnetisation

We now want to use the general form of the propagator to compute the high temperature limit of the dOTOC for transverse magnetisation in free fermionic systems:

[TABLE]

with

[TABLE]

Plugging (120) into (119) and using the definition (109), we obtain

[TABLE]

Note that we changed the summation indices compared to those used in the main text ($i$, $k$ to $l_{1}$, $l_{3}$ and $l$ to $k$). Taking into account the translational invariance of the KI, we can fix one of the indices in the first sum and replace the summation over that index with an overall multiplication by the number of particles in the system $N$, with $N\rightarrow\infty$. It is convenient to fix the last index to $k=0$. The expression then simplifies to

[TABLE]

where the temporal dependence of $K$ is omitted for compactness of notation. The expression is formal and will be simplified in what is to follow. Furthermore, we also have

[TABLE]

The key to simplifying the expressions $\left(I\right)$ and $\left(II\right)$ are the anti-commutation relations

[TABLE]

or the equation for the pair correlation function that follows from them:

[TABLE]

We can use it together with the Wick’s theorem to compute the eight-fermion correlation functions appearing in (122) and (123). We see that the infinite temperature expectation values in (122) and (123) are only non-zero if all of the Majorana fermions in them appear in pairs. In particular, it is helpful to consider the cases of $j\neq 0$ and $j=0$ separately.

a) $j\neq 0$

First, consider the terms in the correlator for $\left(I\right)$ that have the form

[TABLE]

where the empty slots ($\underline{\,\,\,\,\,\,}$) have to be filled by $w_{2j}$, $w_{2j+1}$, $w_{0}$, $w_{1}$, each appearing exactly once. This gives $24$ possible permutations. Then, for the corresponding correlator in $\left(II\right)$ (with the first two pairs of terms interchanged), we want

[TABLE]

so that the combination results in a non-zero term in (119). Note that the only allowed combinations are those where we have only one among $w_{2j}$ and $w_{2j+1}$ in the first two slots. The remaining slot, among the first two slots, has to be filled by either $w_{0}$ or $w_{1}$. This is also a direct consequence of the anti-commutation relations. We are left with $16$ possible permutations, which we write out explicitly:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now, we introduce the following notation:

[TABLE]

the spinor “inversion”:

[TABLE]

and the following notation for the vector components

[TABLE]

This allows us to summarise the above $16$ lines in a compact formula:

[TABLE]

where the summation runs over the $3$ spatial indices: $j$, $l_{1}$, $l_{3}$ and $4$ permutations: $s_{j}$, $s_{0}$, $p_{1}$, $p_{3}$. Here, $s_{j}$ denotes the ”spin” of the $j$-type fermion in the first pair of slots ($w_{2j+s_{j}-1}$) and $s_{0}$ the ”spin” of the [math]-type fermion in the first two slots ($w_{s_{0}-1}$). Furthermore, $p_{1}$ denotes the permutation of fermions in the first pair of slots and $p_{3}$ the permutation of fermions in the third pair of slots. Note also that

[TABLE]

b) $j=0$

In this case, we have to fill in the following correlator:

[TABLE]

Since the present fermions are already contracted, we can fill the empty slots with arbitrary two pairs of fermions $w_{2\bar{j}+s_{j}-1}$, $w_{2\bar{j}+s_{j}-1}$ and $w_{2l+s_{0}-1}$,$w_{2l+s_{0}-1}$. Again, we want to have

[TABLE]

so that the terms in (119) do not end up cancelling out. It is easy to check that the only way to achieve this is to set either $\bar{j}$ or $l$ to zero, with the remaining index being non-zero. We choose $l=0$ and $\bar{j}\neq 0$. As before, we again have $16$ possible permutations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which we can express as

[TABLE]

where we used $R_{i}:=\left(\begin{matrix}l_{i}-\bar{j}\\ l_{i}\end{matrix}\right)$ and $S:=\left(\begin{matrix}s_{\bar{j}}\\ s_{0}\end{matrix}\right)$.

Renaming the dummy summation index in (137) from $\bar{j}$ to $j$, summing (137) and (132) and plugging it together with $\left(I\right)=-\left(II\right)$ into the expression (119) for $c_{z}(t)$, we finally recover Eq. (19) from the main text:

[TABLE]

Eq. (138) is a general expression for the dOTOC of transverse magnetisation and holds for any free fermion model.

Appendix G Application of the formula for the dOTOC of $M_{z}$

The equation (138) can be used in two ways: for exact numerical computation at intermediate times and for the analytical computation of the long-time asymptotics.

G.1 Intermediate times

For intermediate times, $t<50$, we proceed by first computing the power $\mathscr{U}^{t}$ of KI Floquet propagator (94) for given numerical values of the parameters $J$ and $h$:

[TABLE]

where $U_{n}$ are $2\times 2$ matrices whose elements depend only on $J$ and $h$. It is then easy to see from (118) that

[TABLE]

thereby enabling a direct computation of the real space propagator for given numerical values of $J$ and $h$.

From the above calculation, we also learn that $K^{l}(t)\neq 0$ only for $\left|l\right|\leq t$. This is a direct observation of the fact that the information in KI spreads in a sharp causal-cone with the speed of propagation equal to $1$. Hereon, it follows that the sums over $j$, $l_{1}$ and $l_{3}$ in (138) do not need to be taken over the entire $\mathbb{Z}$ but only over the finite intervals $j\in\left[-2t,2t\right]-\left\{0\right\}$, $l_{1},l_{3}\in\left[\max\left(-t,j-t\right),\min\left(t,j+t\right)\right]$. The number of terms in this sum is proportional to $t^{3}$ so the summation can be efficiently carried out for intermediate $t$. The result is numerically exact.

G.2 Long-time asymptotics

To find the long-time asymptotic behavior of $c_{z}(t)$, we express (138) using (118):

[TABLE]

Introducing the following notation

[TABLE]

we then have

[TABLE]

Here, we took the $j$ sum over the entire $\mathbb{Z}$ in the first term and then subtracted the $j=0$ case in the second term. Performing (formally) the $j$, $l_{1}$ and $l_{3}$ summations and taking into account $\sum_{n}e^{{\rm i}nx}=2\pi\sum_{k}\delta\left(x-k2\pi\right)$, we get

[TABLE]

Finally, integrating over $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ in the first term and $\theta_{1}$ and $\theta_{3}$ in the second, we get:

[TABLE]

It can be shown that the integrand $I_{2}$ vanishes:

[TABLE]

The remaining integration over $\theta$ in $I_{1}$ is in general difficult to perform but we can find the asymptotic behavior for large $t$.

To complete this task we can take advantage of the fact that $\mathscr{U}$ is a unitary matrix. Using the form (95), the powers of $\mathscr{U}$ are simply:

[TABLE]

We can then compute the integral (146) in the large $t$ regime by using the stationary phase approximation Širca and Horvat (2012):

[TABLE]

where $\xi_{j}$ denotes (all of) the local extrema of $\psi(\theta)$, i.e. $\psi^{\prime}(\xi_{j})=0$, on the interval $\left[a,b\right]$.

G.2.1 The case of the kicked quantum Ising model

The considerations so far have been general and can be applied to any quadratic fermion model. For the KI model, we can use the expressions for eigenvalues and eigenvectors from Section D.

Let us introduce the following notation:

[TABLE]

Plugging the results from Section (D) into (146), we see that integrand $I_{1}$ can be written the following form:

[TABLE]

Here, all $P$’s are polynomials in their arguments. The indices denote the power of the term $e^{{\rm i}t\kappa(\theta)}$ multiplying a particular polynomial. All eigenvactor components $\left\{V\right\}$ and the eigenvalue $\kappa$ also depend on parameters $J$ and $h$.

The large $t$ behavior of the integrals $\int_{-\pi}^{\pi}{\rm d}\theta$ of all the terms except for $P_{0}$ can be obtained using stationary phase approximation (149). Plugging in the elements of the eigenvectors at stationary points of $\kappa(\theta)$, that is (105), we see that $P_{-2}$, $P_{2}$, $P_{-4}$, $P_{4}$ vanish at these points. This means that in the large $t$ regime, $c_{z}(t)$ is constant for the transverse field KI.

The only remaining integral

[TABLE]

can be evaluated numerically to get the asymptotic (constant) value of the $c_{z}(t)$ for any given $J$ and $h$ (See the inset of Figure 1 of the main text).