Operator Spreading in Random Unitary Circuits

TL;DR

This paper provides exact analytical results and models for operator spreading in random quantum circuits, revealing universal diffusive and KPZ scaling behaviors in various dimensions, and connects these to classical droplet growth and Tracy-Widom distributions.

Contribution

It introduces exact solutions and coarse-grained models for operator spreading in Haar-random circuits across multiple dimensions, linking quantum dynamics to classical growth processes.

Findings

OTOC satisfies a biased diffusion equation in 1+1D

Front broadening scales as t^{1/2} in 1+1D

In higher dimensions, the front width scales as t^{1/3} in 2+1D and t^{0.24} in 3+1D

Abstract

Random quantum circuits yield minimally structured models for chaotic quantum dynamics, able to capture for example universal properties of entanglement growth. We provide exact results and coarse-grained models for the spreading of operators by quantum circuits made of Haar-random unitaries. We study both 1+1D and higher dimensions, and argue that the coarse-grained pictures carry over to operator spreading in generic many-body systems. In 1+1D, we demonstrate that the out-of-time-order correlator (OTOC) satisfies a biased diffusion equation, which gives exact results for the spatial profile of the OTOC, and the butterfly speed $v_{B}$. We find that in 1+1D the `front' of the OTOC broadens diffusively, with a width scaling in time as $t^{1/2}$. We address fluctuations in the OTOC between different realizations of the random circuit, arguing that they are negligible in comparison to the broadening of the front. Turning to higher D, we show that the averaged OTOC can be understood exactly via a remarkable correspondence with a classical droplet growth problem. This implies that the width of the front of the averaged OTOC scales as $t^{1/3}$ in 2+1D and $t^{0.24}$ in 3+1D (KPZ exponents). We support our analytic argument with simulations in 2+1D. We point out that, in a lattice model, the late time shape of the spreading operator is in general not spherical. However when full spatial rotational symmetry is present in 2+1D, our mapping implies an exact asymptotic form for the OTOC in terms of the Tracy-Widom distribution. For an alternative perspective on the OTOC in 1+1D, we map it to the partition function of an Ising-like model. As a result of special structure arising from unitarity, this partition function reduces to a random walk calculation which can be performed exactly. We also use this mapping to give exact results for entanglement growth in 1+1D circuits.

Equations

$$\mathcal{C}(x,t)\equiv \frac{1}{2}\mathrm{Tr}\rho [{\mathcal{O}}_0(t),Y_x{]}^{†}[{\mathcal{O}}_0(t),Y_x]$$

$$\mathcal{C}(x,t)=1-\mathrm{Tr}{\rho }_{\infty }{\mathcal{O}}_0(t)Y_x{\mathcal{O}}_0(t)Y_x.$$

$$U(t^{\prime },t^{\prime }-1)=\{\begin{array}{ll}{⨂}_{x\in 2\mathbb{Z}}{U}_{x,x+1}(t^{\prime },t^{\prime }-1)& \text{if t′ is even,}\\ {⨂}_{x\in 2\mathbb{Z}+1}{U}_{x,x+1}(t^{\prime },t^{\prime }-1)& \text{if t′ is odd.}\end{array}$$

$$\mathcal{C}(x,t)$$

$$=-\frac{1}{2}\mathrm{Tr}{\rho }_{\infty }[X_0(t),Y_x{]}^2.$$

$$v_B(q)=\frac{q^2-1}{q^2+1}.$$

$$\sigma (q,t)=\frac{2q}{q^2+1}\sqrt{t}.$$

$$\overline{\mathcal{C}(x,t)}\simeq \Phi \left(\frac{v_{B}t+x}{\sigma (t)}\right)\Phi \left(\frac{v_{B}t-x}{\sigma (t)}\right).$$

$$\overline{\mathcal{C}(x,t)}=$$

$$(1-\xi )$$

$$+\xi$$

$$p$$

$$g(n,a,p)$$

$$\overline{\mathcal{C}}(x,t)$$

$${\partial }_t\overline{\rho }(x,t)$$

$$X_0(t)=\sum _{\mathcal{S}}{a}_{\mathcal{S}}(t)\mathcal{S}$$

$$\mathrm{Tr}{\rho }_{\infty }\mathcal{S}{\mathcal{S}}^{\prime }={\delta }_{\mathcal{S}{\mathcal{S}}^{\prime }},$$

$$\sum _{\mathcal{S}}{a}_{\mathcal{S}}(t{)}^2=1.$$

$$\rho (x,t)$$

$$[X_0(t),Y_x{]}^2$$

$$\displaystyle=\Bigg{(}\sum_{\begin{subarray}{c}\mathcal{S}:\mathcal{S}_{x}=\sigma^{y},\sigma^{z}\end{subarray}}2a_{\mathcal{S}}(t)\mathcal{S}Y_{x}\Bigg{)}^{2}.$$

$$\mathcal{C}=-\frac{1}{2}\mathrm{Tr}{\rho }_{\infty }[X_0(t),Y_x{]}^2=\sum _{\begin{array}{c}\mathcal{S}:{\mathcal{S}}_x={\sigma }^y,{\sigma }^z\end{array}}2a_{\mathcal{S}}(t{)}^2.$$

$$\overline{a_{{\mathcal{S}}^{\prime }}(t+1{)}^2}$$

$$W_{{\mathcal{S}}^{\prime }\mathcal{S}}$$

$$\mu (x,t)=\sum _{\mathcal{S}:{\mathcal{S}}_x\ne I}{a}_{\mathcal{S}}(t{)}^2$$

$$\overline{\mathcal{C}(x,t)}=\frac{q^2}{q^2-1}\overline{\mu (x,t)}$$

$$\overline{\mu (x,t)}$$

$$\overline{\mathcal{C}(x,t)}=\sum _{x^{\prime }\ge x}\overline{\rho (x^{\prime },t)}.$$

$$\overline{\rho (t+1,x)}$$

$$\overline{\rho (t+1,x+1)}$$

$$p=\frac{q^2-1}{q^4-1}=\frac{1}{q^2+1}$$

$${\partial }_t\overline{\rho (x,t)}$$

$$v_B(q)$$

Full text

For an alternative perspective on the OTOC in 1+1D, we map it to the partition function of an Ising-like statistical mechanics model. As a result of special structure arising from unitarity, this partition function reduces to a random walk calculation which can be performed exactly. We also use this mapping to give exact results for entanglement growth in 1+1D circuits.

I Introduction

A key challenge for many-body physics is to identify universal properties of quantum dynamics and the approach to thermalization. Particularly important are universal results that hold for generic quantum systems. Examples of such universal properties include the existence of effective light cones for the propagation of quantum information Lieb and Robinson (1972) and the existence of universal scaling forms for the growth and saturation of the von Neumann entanglement entropy in 1+1D Calabrese and Cardy (2005, 2016); Kim and Huse (2013); Zhang et al. (2015); Liu and Suh (2014); Asplund et al. (2015); Casini et al. (2016); Ho and Abanin (2015); Nahum et al. (2016); Mezei and Stanford (2016); Serbyn et al. (2013); Huse et al. (2014); Vosk and Altman (2013); Kaufman et al. (2016).

By definition, generic systems lack the structures (for example large numbers of symmetries or conservation laws) that allow for exact results in typical solvable many-body systems. Surprisingly, insight into generic systems can come from studying dynamics with even less structure than a generic Hamiltonian system, such as the dynamics generated by a random quantum circuit. Random circuit dynamics provide a minimally structured model with which real Hamiltonian dynamics can be compared Hayden and Preskill (2007); Sekino and Susskind (2008); Brown and Fawzi (2012); Lashkari et al. (2013); Hosur et al. (2016); Oliveira et al. (2007); Dahlsten et al. (2007); Žnidarič (2008); Brandao et al. (2016). Despite its simplicity, this model is able to capture universal scaling forms for entanglement growth both in 1+1D and in higher dimensions Nahum et al. (2016). Random circuits are also toy models for information scrambling in black holes and other strongly coupled systems Hayden and Preskill (2007); Sekino and Susskind (2008); Brown and Fawzi (2012); Lashkari et al. (2013); Hosur et al. (2016); Oliveira et al. (2007); Dahlsten et al. (2007); Žnidarič (2008); Brandao et al. (2016); Shenker and Stanford (2015).

In this paper we provide both exact results and coarse-grained descriptions for the spreading of quantum operators under random circuit dynamics, as measured by the ‘out-of-time-order correlator’ (OTOC). The OTOC originally appeared in the study of quasi-classical approximations to superconductivity Larkin and Ovchinnikov (1969), and is closely related to the commutator norm that appears in Lieb-Robinson bounds Lieb and Robinson (1972), but it has been studied recently as a means of quantifying the scrambling of quantum information Kitaev (2014); Shenker and Stanford (2014a, b); Maldacena et al. (2016). It has been argued that early-time exponential growth of the OTOC is a characteristic feature of chaotic quantum systems, and such growth has been obtained within the AdS/CFT correspondence and in a range of physical systems Aleiner et al. (2016); Stanford (2016); Patel and Sachdev (2017); Banerjee and Altman (2017); Gu et al. (2016); Roberts and Swingle (2016); Chowdhury and Swingle (2017); Patel et al. (2017). The OTOC has also been applied to characterize slow dynamics in the presence of disorder and in the many-body localized phase Swingle and Chowdhury (2017); Nahum et al. ; Huang et al. (2016); Chen et al. (2016); Chen (2016); He and Lu (2017); Fan et al. (2017); Chapman and Miyake (2017). It has been calculated in 2D conformal field theories Roberts and Stanford (2015) and integrable chains Dóra and Moessner (2016), and studied numerically in nonintegrable 1D systems Bohrdt et al. (2016); Luitz and Lev (2017); Leviatan et al. (2017). Following theoretical proposals Zhu et al. (2016); Swingle et al. (2016); Yao et al. (2016), experiments addressing the OTOC have been conducted Gärttner et al. (2016); Li et al. (2016); Wei et al. (2016).

Random quantum circuits provide an ideal theoretical setting for the exact calculation of quantities such as the OTOC. While the behaviour of the OTOC in a random circuit is interesting in its own right, we conjecture that the long-distance properties of the OTOC that we derive will also be applicable to deterministic dynamics. Therefore we believe that these results will provide a useful starting point for understanding the generic spatial structure of spreading operators.

An operator $\mathcal{O}_{0}$ which is initially localized near the spatial origin (say, on a single site of a spin chain) will evolve under Heisenberg time evolution into a vastly more complicated operator $\mathcal{O}_{0}(t)=U^{\dagger}(t)\mathcal{O}_{0}U(t)$ that acts nontrivially on many sites. The ‘size’ of $\mathcal{O}_{0}(t)$ is the size of the region in which $\mathcal{O}(t)$ fails to commute with a typical local operator $Y_{x}$ at position $x$. This may be measured by

[TABLE]

where the expectation value has been taken in an appropriate Gibbs state. (For our purposes this will be taken to be the infinite temperature Gibbs state $\rho_{\infty}$, which is the state to which random circuit dynamics equilibrate.) To make the connection with the out-of-time-order correlator (OTOC), we may expand out the commutators in (1). For simplicity let us assume for the moment that the operators $\mathcal{O}_{0}$ and $Y_{x}$ are both Pauli-like operators squaring to the identity. We then have

[TABLE]

The second term, in which the operators are not time ordered, is the OTOC.

At a given time $t$, the range of $x$ where the commutator ${\mathcal{C}}(x,t)$ is significantly larger than zero gives a measure of the size of the operator. This region typically grows ballistically Roberts et al. (2015), even when local conserved quantities exhibit diffusive transport Aleiner et al. (2016); Patel et al. (2017); Bohrdt et al. (2016); Luitz and Lev (2017).111Strongly disordered Hamiltonians in 1+1D provide counterexamples to this ballistic spreading. The immediate natural questions about ${\mathcal{C}}(x,t)$ include: what is the ‘butterfly’ velocity $v_{B}$ associated with this ballistic growth? What is the spatial structure of ${\mathcal{C}}(x,t)$? Is there a ‘hydrodynamic’ equation for ${\mathcal{C}}(x,t)$ at large time and distance scales? Are there important differences between 1+1D and higher dimensions? We will answer all these questions for the case where the time evolution operator $U(t)$ is a circuit composed of Haar random unitaries, as in Fig. 1.

We demonstrate that, both in 1D and in higher dimensions, operator spreading and the growth of the OTOC can be mapped to classical stochastic growth models. We show via an exact calculation that operator spreading in 1+1D can be understood in terms of an equation involving diffusion and drift. The ‘front’ of the operator propagates at a finite velocity $v_{B}$. However the front also broadens diffusively, so that its width is proportional to $\sqrt{t}$ (Fig. 2). We conjecture that this physics also occurs in generic (nonintegrable) 1D systems undergoing deterministic Hamiltonian dynamics. For random circuit dynamics, we must also consider how the averaged correlator $\overline{{\mathcal{C}}}$ differs from the correlator ${\mathcal{C}}$ within a given realization of the random circuit. We argue that fluctuations between realizations are small: typical variations in the front position between different realizations are $O(t^{1/4})$, so negligible in comparison with the $\sqrt{t}$ broadening of the front.

Turning to higher dimensions, we show by an exact mapping that there is a remarkable relationship with a classical droplet growth problem in the Kardar–Parisi–Zhang universality class Kardar et al. (1986). (To avoid confusion, we note that this is not related to the connection between entanglement growth and KPZ introduced in Nahum et al. (2016).) We use this relationship to quantify the broadening of the ‘front’ of a growing operator in a higher dimensional random circuit. In 2+1D the front broadens like $t^{1/3}$ Kardar et al. (1986), and in 3+1D like $t^{0.240}$ Pagnani and Parisi (2015).222The phase diagram of the KPZ equation in higher dimensions Kardar et al. (1986) indicates that in 4+1D and above, two distinct universality classes may be possible for operator growth in a random circuit, one with a growing front width and one without. In the two-dimensional case, and in the absence of lattice anisotropies, recent breakthroughs in the theory of interface growth Johansson (2000); Prähofer and Spohn (2000); Calabrese et al. (2010); Dotsenko (2010); Sasamoto and Spohn (2010a, b, c); Amir et al. (2011); Calabrese and Le Doussal (2011); Prolhac and Spohn (2011); Doussal and Calabrese (2012); Imamura and Sasamoto (2012, 2013) also translate to an exact form for the OTOC, in terms of the celebrated Tracy-Widom distribution (Fig. 3). The broadening of the front of the OTOC is summarized in Fig. 2. Again, we conjecture that these universal scaling forms extend to nonintegrable models with time-independent Hamiltonians, although we note that a previous calculation in a different setting has instead found a front that does not broaden with time, and is governed by a traveling wave equation Aleiner et al. (2016) (see also Chowdhury and Swingle (2017); Patel et al. (2017)). (A traveling wave equation arises from our mappings if we make a certain mean field treatment, Appendix. I. But this mean field is not valid in physical dimensionalities.)

In higher dimensions the shape of the spreading operator333That is, the shape of the spatial region in which ${\mathcal{C}}$ has saturated. We can neglect here the broadening of the front, since at late times the length scale associated with this broadening is parametrically smaller than the size of the operator.

is also of interest. At first sight one might expect the shape of the operator to be asymptotically spherical at late times. Instead, we argue that in systems with an underlying lattice, which have only discrete spatial symmetries, the spreading operator will not become spherical. Its asymptotic shape is determined by a model-specific velocity function $v(\hat{n})$, the speed of the front depending on the local normal vector $\hat{n}$. We verify this for random circuits by simulation in 2+1D.

The results above are for random circuits composed of generic (Haar-random) unitary matrices. It is interesting to compare with random circuits composed of unitaries from the Clifford group, a discrete subgroup which leads to efficiently simulable dynamics Gottesman (1998); Gütschow et al. (2010). In this case the dynamics of the operator is much simpler Nahum et al. (2016), and randomness-induced fluctuations are much more severe. But remarkably the results for the averaged OTOC $\overline{{\mathcal{C}}}$ coincide with the results for generic unitaries. This is a consequence of the fact that the Clifford group is a unitary 2-design DiVincenzo et al. .

In one dimension we give a complementary exact calculation of the averaged OTOC, using a mapping to the partition function of a classical Ising model. These Ising degrees of freedom have a similar origin to those found in calculations in random tensor networks Hayden et al. (2016). We show that special structure arising from the unitarity of the quantum circuit means that this partition function is exactly calculable for any value of the local Hilbert space dimension.

Another important question is how the speed $v_{B}$ associated with operator spreading relates to the speed $v_{E}$ which can be associated with entanglement growth following a quench in 1D Calabrese and Cardy (2005, 2009, 2016); Kim and Huse (2013); Liu and Suh (2014); Asplund et al. (2015); Casini et al. (2016); Ho and Abanin (2017); Mezei and Stanford (2016); Nahum et al. (2016). Refs. Nahum et al. (2016); Mezei and Stanford (2016) pointed out that in general $v_{E}$ is smaller than $v_{B}$, unlike the situation in a 1D conformal field theory Calabrese and Cardy (2016). We extend this here, showing that arbitrarily small values of $v_{E}/v_{B}$ can be achieved without any fine-tuning. (In a system with quenched, i.e. time-independent, disorder it is even possible to have $v_{B}>0$ but $v_{E}=0$ Nahum et al. .)

We also use the Ising mapping described above to give an exact calculation of the average entanglement purity (the exponential of the second Renyi entropy) for a random circuit, complementing the scaling picture, in terms of a coarse-grained minimal cut, of Ref. Nahum et al. (2016).

Contents

1. I Introduction

2. II Operator dynamics in 1+1D

3. II.1 Hydrodynamic equation for averaged OTOC

4. II.2 Hydrodynamic description including fluctuations

5. II.3 Comparison with Clifford circuit dynamics

6. III Higher dimensions

7. III.1 Higher dimensions: setup and mapping to classical growth

8. III.2 Classical model in 2+1D: analytical and numerical results

9. III.3 Scaling of the OTOC in 2+1D

10. III.4 Scaling of the OTOC in $3+1$D and above

11. III.5 Shape of the operator at late times

12. III.6 Formal viewpoint

13. IV Exact calculation of OTOC in ‘spacetime’ picture

14. IV.1 Reduction to Ising spins

15. IV.2 Partition function for two directed paths

16. IV.3 Bounds on Fluctuations

17. V Entanglement growth

18. V.1 Nonuniversality of the ratio $v_{E}/v_{B}$

19. VI Outlook

20. A Evolving Distribution on Operator Strings

21. B Velocity and diffusion constant for lattice diffusion equation

22. C Noisy diffusion equation

23. D Random Clifford operators

24. E Anomalous behaviour of the front for $\phi=0$

25. F Shape of a spreading droplet for weakly varying $v(\phi)$

26. G Haar average formula

27. H Proof of Eq. (106)

28. I A mean field approximation

II Operator dynamics in 1+1D

We begin by defining the random circuit dynamics which we consider in 1+1D, and describing the ‘hydrodynamic’ continuum picture we propose for the OTOC in 1+1D. In Sec. IV we give an alternative exact calculation of the OTOC, confirming and extending the results below.

II.1 Hydrodynamic equation for averaged OTOC

We consider time evolution by a quantum circuit on an infinite 1D spin chain where each spin (qudit) has local Hilbert space dimension $q$. The structure of the quantum circuit is shown in Fig. 2a. Two-site unitaries are applied to ‘even’ bonds on even time steps and ‘odd’ bonds on odd time steps (a ‘running bond’ layout in the language of bricklaying). Each two-site unitary is drawn independently from the uniform distribution on the two-site unitary group $\mathbb{U}(q^{2})$. Formally, our time evolution operator is ${U(t)=U(t,t-1)U(t-1,t-2)\cdots U(1,0)}$, where a single layer of the circuit is given by

[TABLE]

Each two-site unitary $U_{x,x+1}(t^{\prime},t^{\prime}-1)$ is Haar random and independent of all of the others.

Given an operator $\mathcal{O}$, we write $\mathcal{O}(t)=U(t)^{\dagger}\mathcal{O}U(t)$. We will evaluate the following out-of-time order correlator with respect to this time evolution:

[TABLE]

Here, $\rho_{\infty}$ is the infinite temperature Gibbs state, i.e., the mixture of all possible spin configurations with equal weights. $X_{0}$ is a Hermitian operator located at the origin of the spin chain, and $Y_{x}$ is a Hermitian operator located at site $x$. We take both $X$ and $Y$ to be traceless, and normalized such that $\operatorname{Tr}X^{2}=\operatorname{Tr}Y^{2}=q$. For example if $q=2$ (the spin–1/2 chain) we can take $X$ and $Y$ to be Pauli matrices at sites [math] and $x$, respectively.

Since the unitaries in the circuit are random, we must distinguish between averaged quantities (denoted by $\mathbb{E}_{U}$, or whenever unambiguous by overline $\overline{[\cdots]}$) and quantities within a given realization of randomness. However, we will argue that fluctuations induced by the random circuit are small in a certain sense, meaning that the spatial profile of ${\mathcal{C}}(x,t)$ in a given realization of the circuit is, at large times, parametrically close to the average value $\overline{{\mathcal{C}}(x,t)}$.

Fig. 2 is a schematic of the spatial profile we will show for $\overline{{\mathcal{C}}(x,t)}$ at fixed large time. The ‘size’ of the operator is determined by a butterfly speed which is

[TABLE]

Within a region of size $\sim 2v_{B}(q)t$ the commutator ${\mathcal{C}}(x)$ has saturated to a value very close to unity. Note that for finite $q$ the butterfly velocity is smaller than the ‘naive’ speed limit of unity, which is set by the geometry of the quantum circuit, while in the limit $q\rightarrow\infty$ they coincide. The ‘front’ of the operator, i.e. the region in which ${\mathcal{C}}$ varies between 0 and 1, broadens diffusively. The width of the front is proportional to

[TABLE]

More precisely, letting $\Phi$ denote the cumulative density function of the Gaussian distribution, ${\Phi(y)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{y}e^{-x^{2}/2}\mathrm{d}x}$ (which tends to zero for $y\ll 0$ and to $1$ for $y\gg 0$), we have

[TABLE]

In Sec. IV we will see that Eq. (7) is the partition function of an Ising-like statistical mechanics problem, and will derive an exact formula on the lattice (without any continuum approximation):

[TABLE]

where

[TABLE]

and

[TABLE]

Here we show how Eq. (7) can be related to a continuum hydrodynamic equation which is asymptotically accurate at large times. For $x$ near the operator’s right hand front, ${\mathcal{C}}$ is related to a diffusing conserved density $\rho$:

[TABLE]

[TABLE]

We will explain the quantity $\rho$ below.

To begin with, focus on the spin–1/2 chain ($q=2$). At time $t$ we may write the operator in the basis of products of Pauli matrices Roberts et al. (2015); Ho and Abanin (2017); Mezei and Stanford (2016); Dahlsten et al. (2007); Žnidarič (2008),

[TABLE]

where $X_{0}(t=0)$ is a single site operator. Here the ‘string’ $\mathcal{S}$ can be any product of Pauli matrices on distinct sites. The number of strings in the sum generically grows exponentially with time (at the naive lightcone velocity, set by the geometry of the circuit). The $\mathcal{S}$ are normalized as

[TABLE]

and $X_{0}$ is also normalized so $\operatorname{Tr}\rho_{\infty}X_{0}^{2}=1$, implying

[TABLE]

It is useful also to introduce $\rho(x,t)$, the ‘fraction’ of strings ending at $x$:

[TABLE]

We observe that the OTOC is determined by $a_{\mathcal{S}}(t)^{2}$ as follows. Let $Y_{x}$ be the Pauli matrix $\sigma^{y}$ at site $x$. (This choice does not sacrifice generality due to Haar randomness of the circuit.) Since distinct Pauli matrices anti-commute, we see

[TABLE]

Due to the orthonormality in Eq. (12), we have

[TABLE]

This tells us that if we determine the evolution of $\overline{a_{\mathcal{S}}(t)^{2}}$, then the averaged OTOC is also determined. The dynamics of $\overline{a_{\mathcal{S}}(t)^{2}}$ turns out to be remarkably simple, as shown in Refs. Dahlsten et al. (2007); Žnidarič (2008). It is best understood if we first consider a system of just two sites (rather than an infinite chain) over which a Haar random unitary is applied at time $t$. It is straightforward to calculate (see Appendix A) that for arbitrary $q$

[TABLE]

where

[TABLE]

Note two features. First, the result is linear in $\overline{a_{\mathcal{S}}(t)^{2}}$. Second, $\mathcal{S}^{\prime}$ must be the identity if and only if $\mathcal{S}$ is, but otherwise $\overline{a_{\mathcal{S}^{\prime}}(t+1)}$ is a constant for all $\mathcal{S}^{\prime}\neq I$. In other words, the random unitary introduces a (fictitious) Markov process on the probabilistic ensemble $\{(\mathcal{S},\overline{a_{\mathcal{S}}(t)^{2}})\}$ of strings Dahlsten et al. (2007); Žnidarič (2008). This Markov process describes a single string $\mathcal{S}$ which is stochastically updated over time. If $\mathcal{S}$ is nontrivial, each update maps it to any nontrivial string, with uniform probability. The generalization to multiple spins is immediate: for each pair of spins that interact in a given timestep, the stochastic update is applied to the corresponding two-site substring of $\mathcal{S}$. This Markov process will also be used in higher dimensional setting below, as it is not specific to the 1+1D setting. Note that the fictitious stochastic dynamics, which involves a single evolving string, is entirely different from the stochastic dynamics of the operator $X_{0}(t)$ itself (which is a superposition of exponentially many strings).

Returning to the average of the OTOC, we realize that it only matters whether or not the string component of $X_{0}(t)$ at the site $x$ is the identity. In the ensemble $\{(\mathcal{S},a_{\mathcal{S}}(t)^{2})\}$, the fraction

[TABLE]

of strings that occupy the site $x$, may fail to commute with $Y_{x}$. There are $q^{2}-1$ possible nontrivial operators at the the site, which are all equally probable in the ensemble of string components of $X_{0}(t)$. In the present case of $q=2$, this yields444 For general $q$, one has to start with an operator basis that obeys our normalization condition in Eq. (12). It is easy to construct such a basis. Define $X=\sum_{k\in\mathbb{Z}_{q}}\ket{k+1}\bra{k}$ and $Z=\sum_{k\in\mathbb{Z}_{q}}e^{2\pi ik/q}\ket{k}\bra{k}$. Then, the discrete group generated by these two matrices contains exactly $q^{2}$ elements up to unimportant phase factors. These are not hermitian, but no problem arises if one considers $|a_{\mathcal{S}}|^{2}$. Over Haar random unitaries, one easily obtains $\overline{{\mathcal{C}}(x,t)}=\frac{q^{2}}{q^{2}-1}\overline{\mu(x,t)}$.

[TABLE]

In turn, the average occupation number ${\mu(x,t)}$ can be related to the endpoint density $\rho$, assuming that $x$ is far to the right of the left-hand front of the operator:

[TABLE]

The constant of proportionality $\mu_{0}$ has been determined by assuming local equilibration of the structure of the strings.555To find $\mu_{0}$, make the ansatz that each $\mu(x)$ is independent from $\mu(x^{\prime})$ for $x\neq x^{\prime}$. Under this ansatz, the probability that a pair of sites is partially or fully occupied is $1-(1-\mu_{0})^{2}$, and such an occupied pair evolves to fill one of the pair with probability $1-p$. Therefore, setting $\mu_{0}=(1-p)(2\mu_{0}-\mu_{0}^{2})=\frac{q^{2}-1}{q^{2}}$ yields the stationary state.

Therefore

[TABLE]

It is natural to conjecture that local equilibration of the strings, together with the exponentially large number of strings contributing to $\rho$, will make this identity valid asymptotically even without the average.

It remains to analyze the dynamics of $\overline{\rho}(x,t)$. The above Markov process implies a simple autonomous dynamics for $\overline{\rho}$:

[TABLE]

where

[TABLE]

is calculated by counting the non-identity two-site operators $S$ that have the identity at $x+1$, and the overline denotes averaging over unitaries applied up to a given time.

Recalling that unitaries are applied on even and odd bonds alternately, Eq. 23 gives a complete description of the dynamics of $\overline{\rho(x,t)}$. This is a lattice diffusion equation for the conserved density $\overline{\rho}$. Formally, $\overline{\rho}$ behaves like the probability density for a random walker who starts at the origin, and who prefers to travel to the right since $p<\frac{1}{2}$. In the continuum (i.e. at long timescales) $\overline{\rho}$ satisfies a simple diffusion equation,

[TABLE]

whose drift and diffusion constants are determined in Appendix. B:

[TABLE]

The peak in $\rho$ corresponds to the front of the spreading operator $X_{0}(t)$. It travels at speed $v_{B}(q)$ and broadens as $\sigma(q,t)$ (Eq. 6). We emphasize that this fictitious random walker should not be thought of as ‘the endpoint’ of the operator $X_{0}(t)$, which is a superposition of many strings with different endpoints.

From (22), or in the continuum

[TABLE]

we see that $\overline{C}$ obeys the same equation as $\overline{\rho(x,t)}$ but with different boundary conditions,

[TABLE]

Taking into account the similar behaviour at the left hand front gives (7).

Above we had to make two (very natural) assumptions. One was that we can ignore the interaction between the left end and the right end, and the other was that the occupation density $\overline{\mu(x,t)}$ reaches its equilibrium value. In Section IV we give an exact calculation of the averaged OTOC (including exact results for finite $t$ and $x$, not necessarily large) without making any approximation.

II.2 Hydrodynamic description including fluctuations

Having determined the averaged OTOC, the key question is about the fluctuations between different realizations of the random circuits. From the point of view of exact results this is a much harder question (it is possible to obtain bounds in regions far from the front: we return to this in Sec. IV.3). However, we conjecture that the universal physics of fluctuations in $\rho(x,t)$ can be obtained by upgrading Eq. 25 to a stochastic diffusion equation for the random quantity $\rho(x,t)$. This description indicates that fluctuations are strongly suppressed at late times. Since the diffusive broadening is present in a single realization (i.e. is not an artefact of disorder averaging) it is natural to conjecture that it will also be present in generic non-random 1D many-body systems.

Microscopically we expect noise in both the diffusion constant and the drift, but we restrict to noise in the latter since it is more relevant in the RG sense:

[TABLE]

Here $\eta(x,t)$ is white noise, uncorrelated in space and time.

The statistical properties of this equation are easy to analyze. In the absence of the noisy drift term, $\rho(x,t)$ forms a ‘wavepacket’ whose width grows like $\sqrt{t}$ and whose center of mass is at $x_{\text{cm}}=v_{B}t$. When the noisy drift is turned on, it induces statistical fluctuations in $x_{\text{cm}}$ whose magnitude scales with time as

[TABLE]

A quick way to see this is to ask what the drift velocity has been in a given realization, averaged over the spacetime region visited by the wavepacket. The wavepacket visits a spacetime volume of order $\int^{t}\mathrm{d}t^{\prime}\sqrt{t^{\prime}}\sim t^{3/2}$. Averaging the drift velocity $\eta(x,t)$ over this spacetime volume yields $v_{\text{av}}\sim t^{-3/4}$. The typical random displacement of the wavepacket is thus of order $\Delta x_{\text{cm}}\sim v_{\text{av}}t\sim t^{1/4}$. A standard perturbative calculation in Appendix C reproduces this exponent $1/4$, which also characterizes the spreading of directed waves in random media Saul et al. (1992).

The quantity $\Delta x_{\text{cm}}$ is parametrically smaller than $\sqrt{t}$, the width of the front of the averaged commutator. Therefore this heuristic argument indicates that the front profile of the averaged OTOC also applies to the OTOC within a given instance of the random circuit. This is somewhat surprising. To see why, let us contrast the above Haar random dynamics with Clifford dynamics for $q=2$.

II.3 Comparison with Clifford circuit dynamics

The Clifford group is a discrete subgroup of the unitary group, defined by the property that any Pauli matrix is mapped to a product of Pauli matrices. When the quantum circuit consists of Clifford operators, an initial Pauli matrix remains a single string (rather than evolving into a superposition of exponentially many strings as for dynamics with generic unitaries) and the endpoint density $\rho(x,t)$ is localized on a single site for all times.

However, uniformly random Clifford circuits have a crucial relationship with Haar random circuits. Under a uniformly random Clifford update on a pair of sites, a nontrivial operator is mapped with equal probability to any of the nontrivial operators, and thus the dynamics satisfies the master equation of the Markov process in Eq. (17) Nahum et al. (2016). As a result, the averaged quantities such as the average end point density $\overline{\rho(x,t)}$, the average occupation number $\overline{\mu(x,t)}$, and, most importantly, the average OTOC $\overline{{\mathcal{C}}(x,t)}$, obey exactly the same dynamics as the Haar random case. Formally, this is a consequence of the fact that random Clifford operators form a unitary 2-design DiVincenzo et al. ; see Appendix D for the definition of design and a proof for random Clifford. One may say that Clifford dynamics realizes the a priori-fictitious Markov process in a physical system.666 All the statements here hold for any prime power $q$ such as $q=2,3,4,5,7,8,9,11,13,16,\ldots$.

Despite the equivalence of averaged quantities, the quantities within a realization are entirely different. In the Clifford case the endpoint density $\rho$ and the OTOC ${\mathcal{C}}$ are strongly fluctuating, while we have argued that for generic unitaries they are self-averaging (fluctuations are parametrically small).

III Higher dimensions

We now address the structure of the out of time order correlator ${\mathcal{C}}(x,t)$ in spatial dimensions greater than one, by exploiting the relationship between the averaged OTOC and a fictitious classical Markov process (Sec. II.1). We show that this process is a classical droplet growth problem whose universal physics can be understood in terms of the Kardar–Parisi–Zhang equation Kardar et al. (1986). By taking appropriate averages, we then obtain exact universal exponents and scaling forms for the OTOC in a circuit composed of Haar random unitary matrices. We conjecture that these scaling forms also apply to more realistic Hamiltonian dynamics in non-integrable lattice models and field theories.

Somewhat surprisingly, we show that the ‘shape’ of the spreading operator at late times does not become spherical, unless the microscopic dynamics has symmetry under continuous spatial rotations. In a lattice model, the spreading operator remembers forever that the lattice has only discrete point group symmetries. Our argument for this is not specific to random circuit dynamics. The point is simply that ‘the’ butterfly velocity $v_{B}$, which sets the speed at which the operator’s front moves, generically depends on the front’s orientation, resulting in an anisotropic profile for the spreading operator at long times. Another surprising outcome, given previous work in the context of many-body perturbation theory including Ref. Aleiner et al. (2016), is that for the dynamics considered here the averaged OTOC $\overline{{\mathcal{C}}}$ does not satisfy a local differential equation.

In 2+1D, when lattice anisotropy is absent (e.g. in an appropriate continuum model) or negligible, recent results in KPZ theory Johansson (2000); Prähofer and Spohn (2000); Calabrese et al. (2010); Dotsenko (2010); Sasamoto and Spohn (2010a, b, c); Amir et al. (2011); Calabrese and Le Doussal (2011); Prolhac and Spohn (2011); Doussal and Calabrese (2012); Imamura and Sasamoto (2012, 2013) yield the the full functional form of $\overline{{\mathcal{C}}(x,t)}$ as a function of position and time. For an initially localized operator, this is expressed in terms of the GUE Tracy Widom distribution Sasamoto and Spohn (2010a); Amir et al. (2011) (which describes the extremal eigenvalue statistics for the Gaussian Unitary Ensemble of Hermitian matrices Tracy and Widom (1994); Forrester (1993)).

III.1 Higher dimensions: setup and mapping to classical growth

We now describe the unitary dynamics for which we wish to study operator spreading and the OTOC. We choose a circuit where in each timestep Haar-random two-site unitaries are applied to bonds of a $d$-dimensional cubic lattice in a manner that generalizes the 1+1D protocol. We describe the 2+1D case for concreteness; the generalization to higher dimensions is immediate. The periodicity of the circuit is 4 layers. Four successive layers cycle through the four columnar ‘dimer coverings’ of the square lattice as shown schematically in Fig. 4 and Fig. 5, so that the site at the origin interacts sequentially with its neighbours at $\boldsymbol{x}=(0,1)$, $(-1,0)$, $(0,-1)$, $(1,0)$.

The reduction to a classical stochastic process in terms of the fictitious occupation numbers

[TABLE]

proceeds just as in 1D (Sec. II and Appendix. A). Consider two adjacent sites $\boldsymbol{x}$ and $\boldsymbol{y}$ which undergo a joint update in a given timestep. If both sites are initially empty ($n(\boldsymbol{x})=n(\boldsymbol{y})=0$) they remain so after the update. If at least one of the sites is initially occupied ($n(\boldsymbol{x})=1$ or $n(\boldsymbol{y})=1$ or both) then the configuration after the update can be $n(\boldsymbol{x})=1$, $n(\boldsymbol{y})=0$ with probability $p$, or $n(\boldsymbol{x})=0$, $n(\boldsymbol{y})=1$ with the same probability, or $n(\boldsymbol{x})=n(\boldsymbol{y})=1$ with probability $1-2p$, where as before

[TABLE]

If we consider the OTOC for a spreading operator which is initially localized at a single site, then the corresponding classical model is initialized with $n=1$ at the origin and $n=0$ everywhere else. A possible evolution of $n(\boldsymbol{x})$ in a single timestep is shown in Fig. 4.

Recall that the Haar-avaraged OTOC is related to the mean occupation number for this Markov process at time $t$ by the relation

[TABLE]

as illustrated schematically in Fig. 6. The averages on the two sides of the above equation have different meanings. On the left, the bar denotes an average over realizations of a unitary circuit, and ${\mathcal{C}}$ is a correlator for this quantum dynamics. On the right, the angle brackets denote an average in a classical stochastic process. The real number ${\mathcal{C}}$ and the integer $n$ are only related after averaging. As we noted above, the fictitious Markov process can be physically realized by random Clifford dynamics, whenever $q$ is a prime power.

III.2 Classical model in 2+1D:

analytical and numerical results

The ‘seed’ at the origin grows to produce a cluster of linear size $\sim t$. In the interior of this cluster the state equilibrates rapidly to a state in which nearby sites are essentially uncorrelated, with average occupation $\left\langle{n(x,t)}\right\rangle_{\text{classical}}=(q^{2}-1)/q^{2}$. In a given realization there is an interface between the occupied and unoccupied regions which is sharp on length-scales of the order of the lattice spacing. The evolution of the droplet is very similar to well-studied growth models such as the Eden model Eden (1961), and reduces to the stochastic growth of this one-dimensional interface. The size of the occupied region grows linearly in time, with statistical fluctuations in the shape of the interface. (The average shape in our 2D model is not circular, but has only four-fold rotational symmetry; we discuss this in Sec. III.5.)

Typically such growth processes are in the universality class of the Kardar-Parisi-Zhang (KPZ) equation Kardar et al. (1986). Consider a section of the interface and let $\xi$ be a coordinate parallel to the interface and $h$ its height in the perpendicular direction. The KPZ equation is

[TABLE]

where $\zeta$ is uncorrelated spatiotemporal noise. The constant $c$ contributes to the average normal growth rate for the interface, while the $\nu$ term describes diffusive smoothing of sharp features. Finally, the non-linear $\lambda$ term encodes the dependence of average growth rate on the slope. This equation renormalizes to a nontrivial fixed point. One of its most basic properties is the fact that the fluctuations in the height at a given position $\xi$ grow with time as $t^{\beta}$, with $\beta=1/3$.

Let us write the shape of the droplet as a parameterized curve in polar coordinates, with $R(\theta)$ the radius at angle $\theta$ from the origin. (As mentioned above, the interface is sharp on an $O(1)$ lengthscale, and therefore $R(\theta)$ is well defined up to an $O(1)$ uncertainty; this is sufficient since the properties we discuss below are on parametrically larger lengthscales when $t$ is large.) From KPZ scaling we would expect

[TABLE]

with the exactly-known exponent $\beta=1/3$. We will discuss the nonuniversal function $r(\theta)$ below and in Sec. III.5, and we will discuss more detailed universal properties in the next section.

We have examined the growth of the droplet for spin-1/2 degrees of freedom ($q=2$) on the square lattice, by tracking the average, evolving support of a cluster over $M=2\times 10^{3}$ realizations of the classical dynamics up to time $t=1000$. We store only the density $\langle{n}(\boldsymbol{x},t)\rangle$, averaged over all $M$ realizations, as a function of position and time, as this is the quantity with a direct interpretation in the quantum setting. We have also investigated smaller values of $q$: these do not have an interpretation in the quantum circuit, but in the classical model decreasing $q$ simply corresponds to increasing the probability $p$ in the update. At each time slice, the form of $\langle{n}(\boldsymbol{x},t)\rangle$ is fitted, along cuts through lattice symmetry axes, to extract the cluster size and the the width of the front (where $\left\langle n(\boldsymbol{x},t)\right\rangle$ is appreciably different both from zero and from its $t\rightarrow\infty$ value). We observe linear growth of the size as expected. Note that the fluctuations in the second equation of (35) imply that the width of the front region is expected to scale like $t^{1/3}$.

Fig. 8 (top) shows the growing width of the front for cuts along the diagonal, $\theta=\pm\pi/4$. There, at the largest times we can access, the fitted exponent is $\beta=0.3305\pm 0.0269$, extracted from a fit to the blue data points in Fig. 8. As expected, this value is consistent with the KPZ value $\beta=1/3$.

A slight surprise is that the behaviour along the axis, e.g. at $\theta=0$ is rather different: see Fig. 8 (bottom), which does not show KPZ growth. Generically the only stable fixed point for the growth of a 1D interface is believed to be the KPZ fixed point. However anomalous growth is possible in this model, for $q$ greater than a critical value $q_{c}\lesssim 2$, when the direction of the front’s local normal vector is fine-tuned to coincide with one of the axes, as occurs at $\theta=0$. In this regime, a front with normal parallel to a lattice axis moves at a speed exactly equal to the naive light-cone speed, $v_{B}=2$, and does not roughen. This is a known phenomenon in various lattice growth models in discrete time which have synchronous parallel updates, and can be understood by a relationship with directed percolation Richardson (1973); Durrett and Liggett (1981); Savit and Ziff (1985); Kertész and Wolf (1989); Krug et al. (1990): see Appendix. E for an explanation. While interesting, this phenomenon is an artefact of the specific discrete spacetime geometry we have chosen, which could be eliminated by modifying this geometry,777The effect disappears for smaller $q$. For example for $q=1.4$ we see clear KPZ growth both at $\theta=0$ (fitted exponent value $\beta=0.3223\pm 0.0199$) and at $\theta=\pi/4$ ($\beta=0.3304\pm 0.0149$). and we certainly do not expect it to be relevant to continuous time dynamics. (It would be interesting to look for this effect in appropriate deterministic Floquet dynamics, however.) It has an effect on the shape of the droplet, which we discuss in Sec. III.5.

We now discuss the OTOC scaling that results from the KPZ mapping, neglecting effects of lattice anisotropy (which we will return to in Sec. III.5).

III.3 Scaling of the OTOC in 2+1D

We have already mentioned the basic consequence of KPZ growth, which is the $t^{1/3}$ broadening of the front. But, unusually for a nontrivial fixed point, not only the exact critical exponents but also certain exact scaling functions are known for the growth of an interface in 1D Johansson (2000); Prähofer and Spohn (2000); Calabrese et al. (2010); Dotsenko (2010); Sasamoto and Spohn (2010a, b, c); Amir et al. (2011); Calabrese and Le Doussal (2011); Prolhac and Spohn (2011); Doussal and Calabrese (2012); Imamura and Sasamoto (2012, 2013) (see Sasamoto and Spohn (2010d); Kriecherbauer and Krug (2010); Corwin (2012); Halpin-Healy and Takeuchi (2015) for reviews). We can now apply this information to the OTOC to obtain scaling functions which we propose are generic.

To simplify things let us consider a case where lattice anisotropy is absent or very weak, so that the spreading operator is circular and the OTOC depends only on a radial coordinate and time. Weak anisotropy could certainly be engineered in an appropriate random circuit. More importantly, we conjecture that the scaling form below captures universal scaling in realistic rotationally invariant many-body systems and field theories.

For the growth of a droplet, the probability distribution of the interface radius is given by the GUE Tracy Widom distribution Sasamoto and Spohn (2010a); Amir et al. (2011) (which has been been observed experimentally in striking experiments on the growth of a turbulent domain in liquid crystals Takeuchi et al. (2011); Takeuchi and Sano (2012)). Following convention we write

[TABLE]

where the non-universal constants $v_{B}$ and $c$ are of order one, and $\chi(\theta,t)$ is a random variable whose mean and variance are of order one at large times. Focussing on a fixed value of $\theta$, the cumulative probability distribution of $\chi$ at a fixed time is $t$–independent at large times and given by the Tracy Widom distribution $F_{2}$:

[TABLE]

Remarkably, this allows us to fix the full functional form of $\overline{{\mathcal{C}}(\boldsymbol{x},t)}$ in two dimensions, in the case where lattice anisotropy is absent. In polar coordinates $(r,\theta)$, and in the continuum, Eq. 33 is

[TABLE]

where $\Theta$ is the Heavyside step function. The right hand side is precisely the probability that $\chi$ is greater than $(r-v_{B}t)/ct^{1/3}$. We suppress the $\theta$ dependence since we are assuming rotational symmetry:

[TABLE]

The form of $\overline{{\mathcal{C}}(r,t)}$ is shown in Fig. 9. The asymptotic behaviour near the trailing edge, close to saturation, i.e. for $[v_{B}t-r]/ct^{1/3}\gg 1$, is Deift et al. (2008); Baik et al. (2008)

[TABLE]

where $b_{1}=2^{1/24}e^{\zeta^{\prime}(-1)}$ with $\zeta^{\prime}(-1)\approx-0.165$, the derivative of the Riemann zeta function. Near the leading edge, $[r-v_{B}t]/ct^{1/3}\gg 1$,

[TABLE]

The former asymptotic expansion of $F_{2}$ was achieved only recently Deift et al. (2008); Baik et al. (2008).

One can also consider operator spreading with other initial conditions. For example we can initialize an operator in a half-plane so that ${\mathcal{C}}(x,t)$ has a straight, rather than a circular, front. The scaling form for ${\mathcal{C}}(x,t)$ will then be given by the Tracy Widom distribution of the Gaussian orthogonal ensemble, denoted $F_{1}$. The objects $F_{1}$ and $F_{2}$ are of fundamental importance in a broad range of mathematical and physical problems and it would be very interesting to see whether any of these connections shed light on operator growth.

III.4 Scaling of the OTOC in $3+1$D and above

The basic features of the 3+1D case are very similar to those in 2+1D. The KPZ equation extends to an interface of arbitrary dimensionality Kardar et al. (1986). For the the 3+1D quantum problem, the dimensionality of the interface is two and the critical exponent $\beta$ relevant to the width is $\beta\simeq 0.240$ (Ref. Pagnani and Parisi (2015) and references therein). The analogue of $F_{2}$ which yields the universal form of $\overline{{\mathcal{C}}}$ is not known analytically, but has been determined numerically Halpin-Healy and Takeuchi (2015). Numerical simulations for the 3+1D random circuit, along the lines of those above, would be feasible.

Dimensions equal to or higher than 4+1 are of course inaccessible experimentally, but they are nonetheless interesting because in these high dimensionalities the KPZ equation yields a phase transition as a function of the strength of nonlinearity Kardar et al. (1986). Both a rough phase, in which fluctuations grow as $t^{\beta}$ with $\beta>0$, and a smooth phase, where fluctuations remain of order one as $t\rightarrow\infty$, exist. It would be interesting to know whether both phases are accessible in appropriate many-body systems.

III.5 Shape of the operator at late times

It is interesting to consider the shape of the spreading operator at late times — i.e. the shape of the growing spatial region in which the OTOC $\overline{{\mathcal{C}}(x,t)}$ has already saturated to its late time value (to within an exponentially small correction). Rescaling distances by a factor of $t^{-1}$ gives a ‘droplet’ of $O(1)$ size, which we expect to reach a fixed asymptotic shape. In this scaling limit the width of the front is negligible, so the front can be treated simply as a curve. What is its asymptotic shape?

At first glance, one might expect that the asymptotic shape is a circle in two spatial dimensions and a sphere in higher dimensions. For example, this would be expected if the OTOC satisfied a local nonlinear differential equation in which derivatives higher than 2 could be neglected, as discrete lattice symmetries would ensure that such an equation had symmetry under continuous spatial rotations. Instead, we argue that for many body systems on the lattice the shape of the operator is model-dependent and retains information about the discrete symmetries of the lattice, even at arbitrarily late times. For the random circuit model this follows immediately from the mapping to domain growth processes, for which anisotropy is a well-known feature Freche et al. (1985); Hirsch and Wolf (1986); Meakin (1988); Wolf (1987); Krug et al. (1991). Figure. 10 shows the shape of the droplet in the present model for various values of $q$.

For concreteness consider the 2D case (similar statements hold in higher dimensions). The asymptotic droplet shape is described by a radius $R(\theta)$ depending on the polar coordinate $\theta$. Since the size of the operator is large at large times, the curvature of the front is parametrically small, except possibly at isolated $\theta$ values where $R(\theta)$ is not smooth. Away from such isolated points, the local velocity of the front, in the direction of its normal vector, can depend only on the orientation of this local normal vector. This dependence is captured by a velocity function $v_{B}(\phi)$, where $\phi=\phi(\theta)$ is the angle of the local normal vector to the $x$-axis. A priori $v_{B}(\phi)$ is constrained only by lattice symmetry; for example on the square lattice

[TABLE]

It is evident that the asymptotic shape cannot be a circle except when $v(\phi)$ is a constant function. Since the front of the operator advances by $v(\phi)\mathrm{d}t$ in the normal direction $\hat{n}$, the distance between the front and the origin grows by $v(\hat{n})\mathrm{d}t/\hat{n}\cdot\hat{r}$, which must be equal to $h(\hat{r})\mathrm{d}t$. Expressing the normal vector in terms of $h$, one obtains

[TABLE]

where $\phi(\theta)$ is the angle of the normal at polar position $\theta$ on the interface. This equation is solved by a geometrical construction described in Ref. Wolf (1987); Krug et al. (1991): $h(\theta)=\min_{\phi}\frac{v_{B}(\phi)}{\cos(\phi-\theta)}$. When the effect of lattice anisotropy is weak (as is likely to be the case in many realistic situations when the relevant degrees of freedom are long-wavelength modes), we expect $v_{B}(\phi)$ to be a smooth, weakly varying function, and we may also solve for the shape perturbatively in ${w(\phi)=v_{B}^{\prime}(\phi)/v_{B}(\phi)}$, as described in Appendix. F. Restoring the time dependence,

[TABLE]

However when $v_{B}(\phi)$ varies sufficiently strongly, the asymptotic shape $R(\theta)$ can include sharp corners or straight segments on the boundary: in this regime the perturbative solution above is no longer appropriate.

For many-body systems in continuous time we expect $v(\phi)$ to be analytic. In the present lattice model, $v(\phi)$ is analytic for $q<q_{c}$ ($q_{c}\lesssim 2$) while for $q>q_{c}$ this function is nonanalytic near $\phi=0$ as a result of the anomalous behaviour of a lattice-aligned front: $v(\phi)\simeq 2+\text{const.}\,|\phi|$ Krug et al. (1990). This leads to flat facets near $\theta=0$ in the asymptotic shape Krug et al. (1990). This change in the surface morphology as $q$ is varied is shown in Fig. 10.

For the random circuit model, it is straightforward to determine $v_{B}(\theta)$ in the extreme limit $q=\infty$ where growth becomes deterministic. The propagation of the front in this limit is similar to that of the ‘next nearest neighbour’ deterministic Eden model introduced in Ref. Krug and Spohn (1988) and has the same nonanalytic angular dependence of the velocity888This is seen by following the evolution of a front through the four layers comprising the time step. Krug and Spohn (1988):

[TABLE]

In this limit, the growing operator is simply a square.

Fig. 11 shows the angular dependence of the radius for the 2+1D random circuit dynamics at $q=2$, for several values of the time, showing a clear anisotropy. Note also that $R(\theta=0)\rightarrow 2t$ at late times. In the light of our 1D results where, for finite $q$, $v_{B}$ is always less than the speed associated with the naive lightcone, it is remarkable that in a higher dimensional circuit it is possible for the OTOC front to propagate at the maximal speed in some directions. However we emphasize that this effect relies on the specific discrete spacetime geometry.

III.6 Formal viewpoint

Before returning to 1D, we restate the higher-dimensional results of Sec. III.1 in a more formal language which parallels our discussion in 1D. We introduce a density on clusters, $C$, where $C$ is a collection of sites:

[TABLE]

Here $\mathrm{supp}(\mathcal{S})$ is the support of $\mathcal{S}$. After coarse-graining, we can represent $C$ by a closed surface of spherical topology, namely the boundary of the coarse-grained cluster. Therefore $\rho(C)$ is the natural analogue of the ‘endpoint density’ $\rho(x)$ in 1D. The surface growth picture implies that the effective dynamics of ${\rho}(C)$ are the dynamics of the probability distribution of a growing interface. Therefore, when this is KPZ, $\rho(C)$ satisfies the Fokker-Planck equation corresponding to the KPZ equation. We will discuss this further elsewhere.

IV Exact calculation of OTOC in ‘spacetime’ picture

We now given an analytical treatment of the OTOC from a ‘spacetime’ point of view. This leads to connections with domain walls in an effective Ising model. Similar Ising degrees of freedom have appeared in work on random tensor networks Hayden et al. (2016). Here the effective Ising model looks complicated at first sight, but turns out to be much simpler than those encountered in random (non-unitary) tensor networks, due to special structure arising from unitarity.

This spacetime picture may be much more generalizable than the dynamical point of view above. In Sec. V we will use it to calculate an entanglement–related quantity. In the future, we hope that the tools introduced in this section will be generalizable to higher moments of the OTOC which capture fluctuations ($\overline{{\mathcal{C}}^{2}}$ etc.), or higher powers of the commutator, or to a direct calculation of the von Neumann entropy.

Our exact result for the OTOC for arbitrary $x$ and $t$ (not necessarily large) is

[TABLE]

where

[TABLE]

and

[TABLE]

There are no approximations in Eq. (48). Approximating $g$ by the cumulative density function ${\Phi(y)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{y}e^{-x^{2}/2}\mathrm{d}x}$ of the Gaussian distribution, we reproduce Eqs. 5, 6, 7 above. This approximation is valid when $t$ is large.

Although the spin chain is spatially infinite in both directions and so is our quantum circuit, the time evolved operator $U^{\dagger}(t)X_{0}U(t)$ is supported only on the interval $[-t,t-1]$ of length $2t$. Therefore, it suffices to consider an observable $Y$ inserted in this interval, and our correlator becomes the trace of a $q^{2t}\times q^{2t}$ matrix. The infinite temperature Gibbs state reduces to the identity matrix divided by the dimension $q^{2t}$. Expanding the commutator, we see

[TABLE]

The Haar average of the first term is easy to evaluate. The observable $X_{0}^{2}$ is conjugated by a unitary $U_{-1,0}(1,0)$ and after taking the Haar average becomes proportional to the identity. The constant of proportionality is fixed by the trace-preserving condition. By the normalization convention, $\operatorname{Tr}X_{0}^{2}=q=\operatorname{Tr}Y_{x}^{2}$, and therefore the Haar average of the first term is equal to $q^{-2t}\operatorname{Tr}I=1$. The second term $F$ contains all the complexity.

Observe that the local unitaries form a square lattice that is rotated by $45^{\circ}$. It is thus natural to introduce null coordinates as

[TABLE]

Due to the cylic property of the trace, the only unitaries in the circuit that could affect the correlator are those in the intersection (a rectangle) of the future light cone of $X_{0}$ and the past light cone of $Y_{x}$. From now on, let us use $\ell_{u}$ and $\ell_{v}$ to denote the linear sizes of this intersection along $u$- and $v$-direction, respectively. There are $\ell_{u}\ell_{v}$ local unitaries contained in the intersection of the lightcones.

IV.1 Reduction to Ising spins

For each local unitary $U$ the expression $F$ contains two $U$s and two $U^{\dagger}$s. We will see that averaging over the local unitaries allows us to express $F$ as a partition function for a set of classical Ising spins. To see why such Ising spins arise, consider the standard expression for the Haar average of a single unitary matrix in $\mathbb{U}(n)$:

[TABLE]

(See Appendix G for a self-contained derivation of this formula.) It is convenient to regard the above expression as a matrix whose rows are labelled by the multi-index $(a^{\prime},b^{\prime},c^{\prime},d^{\prime})$ and whose columns are labelled by $(a,b,c,d)$. Note that two types of contraction appear for the unprimed indices, namely $\delta_{ab}\delta_{cd}$ and $\delta_{ad}\delta_{bc}$, and similarly for the primed ones. Correspondingly, in bra–ket notation the above matrix can be written in terms of two vectors which we denote $\ket{I_{\uparrow}}$ and $\ket{I_{\downarrow}}$ (the reason for the notation will become clear below):

[TABLE]

In the natural basis, these vectors are

[TABLE]

With this definition, we may write (53) as

[TABLE]

with

[TABLE]

We see that the unitary may be associated with a pair of classical Ising ‘spins’ $s$ and $s^{\prime}$.

For the application of interest to us the unitaries are two-site unitaries acting on the $q^{2}$-dimensional space associated with spins $i$ and $i+1$. In this case it is easy to see that

[TABLE]

with

[TABLE]

and now $\alpha,\ldots,\delta$ run over the $q$ basis vectors associated with a given spin. The vectors $\ket{\uparrow}$ and $\ket{\downarrow}$ have norm 1 and satisfy

[TABLE]

When we consider $F$, the spins arising from each unitary in the circuit will form an interacting network. The interactions between spins from the same unitary will be given by $w(s,s^{\prime})$, while the interactions between spins from different unitaries will arise from the inner products of kets $\ket{\uparrow,\downarrow}_{i}$ associated with a given spin.

As a final piece of notation, we generalize the definition of $\ket{I_{\uparrow}}$ and $\ket{I_{\downarrow}}$. Given an operator $\mathcal{O}_{ab}$ on the $n$-dimensional space, we define $n^{4}$-dimensional vectors $\ket{\mathcal{O}_{\uparrow}}$ and $\ket{\mathcal{O}_{\downarrow}}$ via

[TABLE]

Choosing $\mathcal{O}$ to be the identity gives the vectors $\ket{I_{\uparrow,\downarrow}}$.

Before we evaluate $F$ for arbitrary $\ell_{u}$ and $\ell_{v}$, let us consider the simplest case where $\ell_{u}=\ell_{v}=1$.

[TABLE]

In the second line, $X^{\prime}=I\otimes X$ and $Y^{\prime}=I\otimes Y$. The third line is a trivial rearrangement of the second, and the fourth employs the formal correspondence between matrices and normalized vectors introduced above.

The Haar average of the tensor product of four unitaries is given by (53) with $n=q^{2}$. To complete the evaluation of $\overline{F}$ we observe that $\braket{(I\otimes X)_{\uparrow}}{(I_{q^{2}})_{\uparrow}}=q^{-4}(\operatorname{Tr}X)^{2}(\operatorname{Tr}I_{q})^{2}=0$, $\braket{(I\otimes X)_{\uparrow}}{(I_{q^{2}})_{\downarrow}}=q^{-4}(\operatorname{Tr}X^{2})(\operatorname{Tr}I_{q})=q^{-2}$, $\braket{(I\otimes Y)_{\downarrow}}{(I_{q^{2}})_{\uparrow}}=q^{-4}(\operatorname{Tr}Y^{2})(\operatorname{Tr}I_{q})=q^{-2}$, and $\braket{(I\otimes Y)_{\downarrow}}{(I_{q^{2}})_{\downarrow}}=q^{-4}(\operatorname{Tr}Y)^{2}(\operatorname{Tr}I_{q})^{2}=0$. This gives $\overline{F}=-1/(q^{4}-1)$.

When $\ell_{u},\ell_{v}>1$, we map the layout of local unitaries to a partition function for the spins $s$, $s^{\prime}$ in Eq. 55. To facilitate the mapping, we decompose the input bra $\bra{I_{s^{\prime}}}$ and output ket $\ket{I_{s}}$ into separate ‘legs’ corresponding to the two physical spins, as in Eq. (57),

[TABLE]

Similarly, the vectors encountered above for the case $\ell_{u}=\ell_{v}=1$ can be decomposed as $\ket{(IX)_{\uparrow}}=\ket{\uparrow}\ket{X_{\uparrow}}$ and $\ket{(IY)_{\downarrow}}=\ket{\downarrow}\ket{Y_{\downarrow}}$, which satisfy

[TABLE]

The expression Eq. (61) is now depicted as in Fig. 12.

It is now clear that for general $\ell_{u}$, $\ell_{v}$ we may regard the array of unitaries as a tensor network composed of tensors of the form (62). The boundaries of this tensor network — i.e. the external legs of the array of $\ell_{u}\times\ell_{v}$ unitaries — involve inner products with fixed vectors. Two of the boundary legs are dressed with $q\ket{X_{\uparrow}}$ and $q\bra{Y_{\downarrow}}$; see Fig. 12. Apart from these, the external legs on the top boundary are dressed by states $q\bra{\downarrow}$, while those on the bottom boundary are dressed with $q\ket{\uparrow}$. In addition $\overline{F}$ includes an overall dimension factor $q^{-2t}$ coming from the infinite temperature Gibbs state. For convenience, we absorb the overall dimension factor $q^{-2t}$ into the vectors on the lower boundaries; these vectors are taken to be normalized, whereas the boundary bras in the top boundaries have norm $q$.

We may now interpret $\overline{F}$ as a partition function for the Ising spins $s_{u,v}$ and $s^{\prime}_{u,v}$ which according to Eq. (62) are associated with the unitary at position $(u,v)$. These spins take the values $\uparrow$, $\downarrow$. The weight associated with the ‘bond’ between $s_{u,v}$ and $s^{\prime}_{u,v}$ comes from the single–unitary Haar average and is $\frac{q^{4}}{q^{4}-1}$ if $s_{u,v}=s^{\prime}_{u,v}$, and $\frac{-q^{2}}{q^{4}-1}$ if $s_{u,v}\neq s^{\prime}_{u,v}$. The leg of the tensor network connecting the unitary at $(u,v)$ to that on its lower right at $(u,v-1)$ yields an interaction between $s^{\prime}_{u,v}$ and $s_{u,v-1}$ which comes simply from the inner product $\braket{s^{\prime}_{u,v-1}}{s_{u,v}}$. This gives weight $1$ if $s^{\prime}_{u,v}=s_{u,v-1}$, and weight $\frac{1}{q}$ if $s^{\prime}_{u,v}\neq s_{u,v-1}$.

We have thus mapped the Haar average of the out-of-time correlator to a partition function for Ising degrees of freedom (with the $q$–dependence residing in the interactions on the bonds). At first sight, this may appear to be a formidable problem. Note in particular that some configurations have negative weight. However, a simplification is possible, as a result of the unitarity of the underlying dynamics. A hint that a simplification is possible comes from the fact that the expression for $F$, Eq. (49), becomes trivial if one of the operators $X_{0}$ and $Y_{x}$ is the identity operator. In the Ising language this corresponds only to a slight change of boundary conditions.

The simplification is effected by integrating out the ‘bra’ variable $s^{\prime}_{u,v}$ from each unitary. This generates a three-spin interaction among the ‘ket’ variables $s_{u,v}$, $s_{u-1,v}$, and $s_{u,v-1}$. The calculation is straightforward and yields the table of weights in Fig. 14. For example, if $s_{u,v}=s_{u,v-1}=s_{u-1,v}$, then the weight is $\frac{q^{4}}{q^{4}-1}\cdot 1\cdot 1+\frac{-q^{2}}{q^{4}-1}\cdot\frac{1}{q}\cdot\frac{1}{q}=1$.

The fact that the weight is zero for two of the configurations means that only a very restricted subset of Ising configurations are allowed. We will show that these can be summed exactly by viewing the configurations in terms of domain walls.

Let us specify the new boundary conditions. The above rules apply along the bottom boundaries due to our normalization convention for the boundary kets, except for the site where the observable ket $\ket{X_{\uparrow}}$ is dangling.

[TABLE]

The top boundary bras, which have norm $q$, follow the rule that

[TABLE]

IV.2 Partition function for two directed paths

Now the problem is reduced to a partition function of Ising variables with the three-body interaction and the boundary interaction. We first simplify the partition function by relating it to one with modified boundary conditions as follows. We denote the weight of a given configuration by $W_{X,Y}(s)$, where the subscripts indicate the dependence on the boundary conditions induced by the operators $X$ and $Y$. Because of the last rule in Eq. (90), the spin at the site where $Y_{x}$ is attached — null coordinate $(\ell_{u},\ell_{v})$ — has to be $s_{\ell_{u},\ell_{v}}=\,\uparrow$. As a result we can replace $Y_{\downarrow}$ with $\downarrow_{\text{fixed}}$, which according to Eq. (90) gives the same weight when $s_{\ell_{u},\ell_{v}}=\,\uparrow$. Let us denote the weight of a configuration $s$ under this modified top boundary condition by $W_{X}(s)$, dropping the subscript $Y$. We may then write the desired quantity $\overline{F}=\sum_{s}W_{X,Y}(s)$ as

[TABLE]

We claim that the first term is equal to 1, and thus

[TABLE]

The claim can be shown in two ways. First, the out-of-time correlator ${1-\overline{F}}$ must vanish if the operator $Y$ is replaced by the identity. The boundary vector $\ket{Y_{\downarrow}}$ then becomes $\ket{\downarrow}$, and the partition function for $\overline{F}$ becomes precisely $\sum_{s}W_{X}(s)$. Therefore ${1-\sum_{s}W_{X}(s)=0}$. The other way to show the claim is by directly integrating out the Ising variables inductively, starting from the top line with respect to the all–$\downarrow$ boundary condition along top boundary. This is a nontrivial consistency check on our reduction.

Now we focus on $W_{X}(s)$ with the variable at null-coordinate $(\ell_{u},\ell_{v})$ fixed to be $\downarrow$. If the bottom variable where $X$ is attached is $\uparrow$, then the second rule in Fig. 14 together with the boundary condition along the bottom boundary dictates that all the bulk variables be $\uparrow$. This cannot be fulfilled for the top-right variable, implying that the weight is zero.

Hence, we have fixed two Ising variables in the bulk to be $\downarrow$ where the observables $X_{0}$ and $Y_{x}$ are attached. Let us think of domain walls instead of spins. The key point is the first rule in Fig. 14, which leads to the domain walls being directed, drastically simplifying the partition function. (Since the top boundaries have a different interaction, let us speak of domain walls to mean disagreeing edges in the south and west of a square encompassing an Ising variable, with the two bottom edges of the square at null coordinate $(1,1)$ excluded.) If we follow a domain wall from the top to the bottom, it should always go down-left or down-right. This implies that there are two non-intersecting domain walls extending from the bottom to the top boundary. The starting vertices of the right and left domain walls have null coordinates $(1,0)$ and $(0,1)$, respectively. Each domain wall has length $t-1$, giving the weight factor

[TABLE]

from the weight table in Fig. 14.

A domain wall can be deformed without changing the weight to the partition function. There is essentially one local deformation of the domain wall. One can easily see from Fig. 14 that the weights of the two configurations in Fig. V are the same.

However, the end position of the domain wall at the top boundary does affect the weight, and we will need to count the number of domain walls for each ending position. For the right domain wall, the end position can be at a vertex $(t-v,v)$ in null coordinates for some $v=0,1,\ldots,\ell_{v}-1$. Likewise, the end position of the left domain wall can be $(u,t-u)$ for some $u=0,1,\ldots,\ell_{u}-1$. The weight from the top boundary interaction is then

[TABLE]

by Eq. (90).

It remains to count the number of domain walls given their end positions. The right domain wall connects $(1,0)$ to $(t-v,v)$ while the left domain wall connects $(0,1)$ to $(u,t-u)$, with the constraint that they must not intersect. To handle the constraint, we use a reflection trick. Regard the domain walls as random walkers. The right random walker randomly chooses between up-left and up-right edges, to proceed from $A=(1,0)$ to $B=(t-v,v)$. Similarly, the left random walker proceeds from $C=(0,1)$ to $D=(u,t-u)$. Any pair of paths $A\to B$ and $C\to D$ that have a common point can be viewed as a pair of paths $A\to D$ and $C\to B$. Conversely, any pair of paths $A\to D$ and $C\to B$, which must meet at at a point, can be viewed as a pair of paths $A\to B$ and $C\to D$ with a common point. Therefore, the number of pairs of paths from $A\to B$ and $C\to D$ without intersection is the number of all unrestricted pairs from $A\to B$ and $C\to D$, minus the number of all unrestricted pairs from $A\to D$ and $C\to B$. The number of our domain wall configurations is therefore

[TABLE]

where the second factor vanishes when $u=0$ or $v=0$.

Finally we combine the results above:

[TABLE]

This correctly reproduces the answer $q^{4}/(q^{4}-1)$ when $\ell_{u}=\ell_{v}=1$. This can be conveniently rewritten as

[TABLE]

where

[TABLE]

[TABLE]

Further simplification is possible since ${g(t-1,a)}\simeq{g(t-1,a-1)}$ for large $t$.

[TABLE]

where

[TABLE]

IV.3 Bounds on Fluctuations

Here we estimate the fluctuation of ${\mathcal{C}}(t,x)$ due to the randomness of the unitaries. One might wish to calculate this fluctuation directly, using a similar technique that we employ for the average of ${\mathcal{C}}(t,x)$, but we found the exact computation unwieldy as it involves high powers of unitaries. Nevertheless, we can argue that the fluctuations are negligible in two regimes.

Since the random variable ${\mathcal{C}}(t,x)$ takes values between 0 and 2, the variance is upper bounded by $2\bar{\mathcal{C}}$. Therefore, the standard deviation is upper bounded by

[TABLE]

This bound is valid for any $t,x$, but only meaningful when $|x|\gg v_{B}t$. This basically says that there is almost no “leakage” of operators beyond the lightcone defined by $v_{B}$. (In passing, we note that one can also use Markov inequality $\Pr[X\geq a]\leq a^{-1}\mathbb{E}X$ which holds for any positive random variable $X$ and a positive number $a$ to have a probability tail bound.)

In the opposite regime where $|x|\ll v_{B}t$, we have shown that the average $\overline{{\mathcal{C}}(t,x)}$ is almost 1; the discrepancy is upper bounded by $O(1)\exp(-(v_{B}t-|x|)^{2}/2\sigma^{2})$. Thus, in this regime the fluctuation is basically given by

[TABLE]

where $F$ is defined in Eq. (49) and $N$ is the dimension of the Hilbert space of spins where $UXU^{\dagger}Y$ is supported on. Here $U$ includes all the local unitaries in the evolution quantum circuit.

To estimate the fluctuation, we consider a slightly different system where $2ct$ spins form a ring, where $c$ is some absolute constant that depends on $q$ only. If $c>1$, this does not modify the dynamics at all, since the evolved operator $UX_{0}U^{\dagger}$ is supported on $2t$ spins. For $c<1$, while we do not insist that this allows us to compute the fluctuation rigorously, we anticipate that qualitative conclusions from this modified setting carry over to the original open chain system.

In Appendix H, we show that if $c=M_{q}^{-1}=\tilde{O}(q^{-2})$, then for all $|x|<ct$

[TABLE]

That is, deep in the lightcone, there is a region in spacetime bounded by a nonzero speed where the fluctuation of ${\mathcal{C}}(t,x)$ is suppressed exponentially in $t$. It is likely that this is only a bound, rather than a tight estimate of fluctuation. Eq. (106) is proved using previous results on approximate unitary designs Brandao et al. (2016), and estimates for $\mathbb{E}_{U}F^{2}(U)$ when $U$ is truly Haar random Hastings .

V Entanglement growth

Entanglement can by quantified in various ways, but perhaps the simplest measure is the entanglement purity $\mathcal{P}=\operatorname{Tr}\,(\operatorname{Tr}_{A^{c}}\ket{\psi}\bra{\psi})^{2}\leq 1$, where $A$ is some region. A pure state $\ket{\psi}$ on $A\cup A^{c}$ is entangled if and only if the purity is not equal to 1. The logarithm of the entanglement purity is the Renyi-2 entropy

[TABLE]

In this section, we calculate exactly the average purity of ‘half’ of the infinite chain, for arbitrary $t$, under the evolution protocol in Section II, with an initial product state. (Previous work has obtained a bound on the saturation time for $q=2$ Žnidarič (2008).) The calculation technique will be very similar to the OTOC calculation; the difference is only in the boundary conditions.

Let $A$ be the left half of our chain, and $B$ be the right half. The initial pure density matrix is $\rho(t=0)=\cdots\otimes P_{-1}\otimes P_{0}\otimes P_{1}\otimes\cdots$, where $P_{i}$ is a projector $\ket{i}\bra{i}$, the state at site $i$. If $U$ is the full time-evolution unitary consisting of local unitaries, then the entanglement purity $\mathcal{P}$ across the cut between $A$ and $B$ is

[TABLE]

The notation $\ket{\uparrow},\ket{\downarrow}$ is the same as in Sec. IV.1. For a one-dimensional projector $P$ on $q$-dimensional space, the $q^{4}$-dimensinal vector $\ket{P_{\uparrow}}$ satisfies

[TABLE]

The expression for the purity can be thought of as a partition function for classical Ising spins as in Sec. IV.1. There are two Ising spins associated with each local unitary; see Eq. (55). Due to Eq. (109), for any configuration of the Ising spins, the weight factor from the bottom boundary is $q^{-|A|-|B|}$, which cancels the factor $q^{|A|+|B|}$ in front of Eq. (108). Hence, the average purity is simply the sum of weights from the domain wall in the bulk (e.g. see Fig. ).

In Sec. IV.1, we first integrated out the ‘bra’ Ising variables $s^{\prime}$, but here we find it simpler to integrate out the ‘ket’ Ising variables $s$. The transition rules of Fig. 14 are now upside down, but otherwise the same. Then, we have a single domain wall starting from the top boundary to reach the bottom. Any domain wall has length exactly $t$, giving rise to weight $\left(\frac{q}{q^{2}+1}\right)^{t}$. The domain wall can choose between left-down or right-down moves as it proceeds from the top, and therefore there are $2^{t}$ domain walls. We conclude that

[TABLE]

We may define the ‘purity speed’

[TABLE]

This quantity gives a bound on the growth rate of the second Renyi entropy:

[TABLE]

The inequality is because the function $f(x)=-\log x$ is convex. Note that this expression bounds the growth rate of $S_{2}$ but does not fix it. The distribution of $S_{2}$ fluctuates in a window of small size compared to its mean Nahum et al. (2016),999 Ref. Nahum et al. (2016) argued that the width of the distribution scales as $t^{1/3}$. The mean value is of order $t$.

but since $S_{2}$ appears in the exponential in ${\overline{q^{-S_{2}}}}$, this does not rule out the possibility that this quantity is affected by rare anomalously small values of $S_{2}$, making it very different from $q^{-\overline{S_{2}}}$.

The von Neumann entropy $S_{vN}$ is always greater than or equal to $S_{2}$, so the growth rate $v_{E}$ of $S_{vN}$ is also bounded by $v_{P}$:

[TABLE]

where the expansion is for large $q$.

In Ref. Nahum et al. (2016) we argued that the universal fluctuations of the entanglement in random circuit dynamics may be understood in terms of a coarse-grained minimal cut, of random shape, through the random circuit. This picture may be contrasted with the domain wall calculation of the averaged purity, which reduces to a statistical mechanics problem without quenched randomness. This is reminiscent of the difference between a quenched and an annealed average in the statistical mechanics of disordered systems Cardy (1996). A direct exact calculation of $\overline{S_{2}}$ (not to mention $\overline{S_{\text{vN}}}$, or of the fluctuations in the entropy) for finite101010 In the limit $q\rightarrow\infty$ it is easy to show that $v_{E}=v_{P}$ in the present model.

$q$ would be much more difficult than the calculation above, as a replica-like limit Cardy (1996) would be required to handle the logarithm. However structure arising from unitarity might make this calculation tractable. This is an interesting task for the future.

The scaling limit of the representation obtained in this section, where we take length and time scales to be large and of the same order, yields a ‘deterministic’ domain wall configuration. This is simply a vertical line for the infinite geometry considered here.111111This is because the $\sqrt{t}$ fluctuations in the transverse position of the domain wall are negligible compared to $t$; compare Nahum et al. (2016) where the minimal cut configuration is also deterministic in the scaling limit. We expect that extending the calculation to higher dimensions will give, in the scaling limit, a formula for $-{\log\overline{\mathcal{P}}}$ as the ‘energy’ of a minimal surface (representing the Ising domain wall) which has a deterministic coarse-grained geometry, obtained from an effective elastic energy minimization problem. This is precisely the scaling picture proposed in Ref. Nahum et al. (2016) for the growth of entanglement in higher-dimensional systems.

V.1 Nonuniversality of the ratio $v_{E}/v_{B}$

In Ref. (Nahum et al., 2016), see also Ref. Mezei and Stanford (2016), we showed that the speed $v_{E}$ associated with entanglement growth is in general smaller than the operator growth speed $v_{B}$, and gave explicit models displaying a ratio $v_{E}/v_{B}<1$. In these models121212 This was determined analytically for a certain large $q$ model, distinct from that here, and numerically for various circuits composed of Clifford gates. this ratio happened to be $1/2$. Values close to 1/2 were also found numerically in Ref. Mezei and Stanford (2016) and Ref. Leviatan et al. (2017). These results might lead one to wonder whether this value is in some sense generic. Here we show that it is not, by constructing random circuit dynamics, involving interactions of large but finite range, which give arbitrarily small values of $v_{E}/v_{B}$ without any fine tuning. The construction uses random unitaries made up of ‘staircases’ of length $O(q^{2})$ which are made up of smaller random unitaries (Fig. 16). When $q$ is large, we obtain a ratio $v_{E}/v_{B}$ which is at most of order $1/q^{2}$. (In a deterministic spin chain with quenched spatial disorder it is even possible to have $v_{E}/v_{B}=0$ Nahum et al. , but here we insist on statistical translational invariance: i.e. we insist that the probability distribution for the circuit is invariant under translations.)

Consider quantum circuit dynamics in which ‘staircase unitaries’ are applied at random locations and at random times in a Poissonian fashion. A staircase is a collection of 2–site unitaries arranged as in Fig. 16. Left and right-oriented staircases are applied with equal probability. The staircase acts on $\ell$ bonds and we take $\ell$ large but finite, satisfying $\ell/q^{2}\gg 1$. Let $r$ be the rate at which staircases are dropped at a given location.

A single staircase can increase the entanglement across a given bond by at most 2 units, implying $v_{E}\leq 2r\ell$. On the other hand a single staircase can move the endpoint of an operator a long way when $\ell\gtrsim q^{2}\gg 1$. The random walk picture of Sec. II.1 shows that in the limit of large $\ell/q^{2}$, a staircase advances the front of the OTOC by an average distance $\sim q^{2}/2$. This involves an average over the two staircase orientations, only one of which is effective in advancing the front a long distance. The large value is because, when $q$ is large, the small value of $p=1/({q^{2}+1})$ (Eq. 24) means the random walker can ‘run’ a long way up a rightward-oriented staircase before falling off. The previous implies $v_{B}\simeq q^{2}r\ell/2$ at leading order in $\ell$. This yields a ratio $v_{E}/v_{B}\lesssim 4/q^{2}$ in this regime, which can be made arbitrarily small by taking $q$ (and hence $\ell$) to be large.

VI Outlook

We have argued that universal scaling forms for the out-of-time-order correlator can be obtained using mappings to paradigmatic problems in classical statistical mechanics. In one dimension we gave an extremely simple hydrodynamic picture in terms of diffusion. In higher dimensions we gave a mapping to classical surface growth and the KPZ equation.131313 Ref. Nahum et al. (2016) obtained a connection between entanglement growth in 1+1D random unitary circuits and the KPZ equation. To avoid confusion we emphasize that the connection with KPZ discussed here is physically entirely distinct from that one, and is not even in the same number of spatial dimensions.

These mappings were derived exactly for random unitary circuits, which are natural ‘least structured’ models for chaotic quantum dynamics in situations where conserved quantities are not playing an important role. We have conjectured that the universal scaling forms found here also apply to OTOCs at asymptotically late times in generic, nonintegrable many body systems and quantum field theories. It will be interesting to test this conjecture in other situations where calculations are possible.

This picture differs from that obtained in a number of previous calculations using many-body perturbation theory Aleiner et al. (2016); Patel and Sachdev (2017); Patel et al. (2017); Chowdhury and Swingle (2017), and it will be interesting to understand the reasons for these differences. Ref. Aleiner et al. (2016) found an operator front that did not broaden in time, whereas here we find a broadening front in all dimensions below $4+1$. Additionally, in Ref. Aleiner et al. (2016) the OTOC was found to obey a local, nonlinear traveling wave equation, which is unlike what we found for random circuits. In 1D we obtained a linear hydrodynamic equation, while in higher dimensions $\overline{{\mathcal{C}}}(x,t)$ in a random circuit is not governed by a local differential equation at all, contrary to standard lore about OTOCs.

Interestingly, a mean field approximation to the classical growth process would yield a local differential (or rather difference) equation for the OTOC, of traveling wave form. This is discussed in Appendix. I. However, the mean field approximation is not valid in physical dimensionalities.

Assuming that the results here do indeed have applications to realistic systems with Hamiltonians that are fixed in time, it will be interesting to consider extensions of the present coarse-grained pictures which take conserved quantities into account.

We have also given exact results for entanglement growth in 1+1D which support the scaling ideas put forward in Nahum et al. (2016), as discussed in Sec. V. In this picture (in any D) entanglement growth is determined by a minimal surface in spacetime, whose geometry becomes well-defined141414But model dependent above 1+1D in an appropriate scaling limit and is determined essentially by an elastic minimization problem. Furthermore, it was argued in that paper and in Ref. Mezei and Stanford (2016) that generically $v_{E}<v_{B}$, where $v_{E}$ is the speed characterizing the growth of entanglement. Here we have shown that it is possible to have arbitrarily small $v_{E}/v_{B}$ in a random quantum circuit.

The effective Ising partition functions for calculating the OTOC and the purity turned out to have interesting structure, making them drastically simpler than they appeared at first sight, and much simpler than the analogous partition function for a non-unitary tensor network Hayden et al. (2016). It would be very interesting to explore whether similar simplifications occur when the averaging involves higher powers of the unitary circuit. If so this would permit calculations of, say, modified versions of the OTOC involving higher powers of the commutator, or a direct calculation of the fluctuations. Even more interesting would be a direct calculation of the von Neumman entropy, which would have to use a replica limit to handle the logarithm.

OTOCs involving higher powers of the commutator are important for comparison with Lieb-Robinson bound. The OTOC considered here can be thought of as the squared Frobenius norm of the commutator divided by the Hilbert space dimension, whereas the Lieb-Robinson bound is on the operator norm of the commutator. The two norms are related as the operator norm is always upper bounded by the Frobenius norm, but our results do not put any nontrivial bound on the operator norm, due to the large dimension factor. The exact relation of the two quantities is yet to be determined.

In addition to exploring implications for realistic many-body systems, interesting questions remain that are specific to the random circuit context. (Note that, at the most basic level, our results show that operator growth saturates the naive causal lightcone of the quantum circuit as $q\rightarrow\infty$, but not for finite $q$.) The randomness in the circuit necessarily implies statistical fluctuations in all observables including ${\mathcal{C}}(x,t)$. We have argued that these statistical fluctuations are (perhaps counterintuitively) a subleading effect at late times. We have shown this in regimes far from the front of the OTOC by giving inequalities, and we have given a heuristic argument for it in the region near the front. This argument was based on a phenomenological extension of the hydrodynamic equation for $\overline{{\mathcal{C}}(x,t)}$ in the 1D case to allow for statistical fluctuations in ${\mathcal{C}}(x,t)$ (Eq. 29). It would be desirable to give a microscopic derivation of Eq. 29. (For the entanglement entropy, statistical fluctuations were investigated in Ref. Nahum et al. (2016).) It also remains to characterize the classical growth problem in Section III more fully, for example by obtaining the nonuniversal constants via an approximate analytic treatment.

The KPZ equation is connected to a remarkable array of topics in classical statistical mechanics Halpin-Healy and Zhang (1995), including the directed polymer in a random medium Huse and Henley (1985) and one-dimensional hydrodynamics Forster et al. (1977), and has beautiful experimental applications Takeuchi et al. (2011); Takeuchi and Sano (2012). Through the Tracy-Widom distribution Tracy and Widom (1994), it is also connected to random matrix theory and an array of combinatorial problems (for example the longest increasing subsequence problem and the statistics of random permutations Baik et al. (1999); Tracy and Widom (2001)). It will be interesting to explore which members of this array can shed light on operator growth.

Related work: While this manuscript was being finalized, we became aware of related work von Keyserlingk et al. , to appear in the same arXiv posting. We also alert the reader to forthcoming numerical work on operator spreading Jonay and Huse .

Acknowledgements.

We thank D. Huse, C. Jonay, D. Chowdhury, J. Ruhman and H. Spohn for useful discussions. We are grateful to T. Veness and J. Chalker for helpful discussions regarding droplet shape. JH thanks M. B. Hastings for useful discussions. AN was supported by the Gordon and Betty Moore Foundation under the EPiQS initiative (grant No. GBMF4303) and by EPSRC Grant No. EP/N028678/1. SV is supported by the DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award de-sc0010526 and was also supported for part of this work by the National Science Foundation under Grant No. NSF PHY11-25915.

Appendix A Evolving Distribution on

Operator Strings

In this Appendix we give a more detailed explanation of the relationship between the dynamics of the coefficients $a_{\mathcal{S}}$ and a Markov process Dahlsten et al. (2007); Žnidarič (2008) and of the derivation of the diffusion picture. We consider Haar-random, local unitary dynamics. In an $N$-site system with a $q$-dimensional Hilbert space at each site, a Hermitian operator that has evolved under the unitary circuit $\mathcal{O}(t)={U}(t)^{\dagger}\,\mathcal{O}\,{U}(t)$ may expanded in a basis of SU$(q^{N})$ generators $\{\mathcal{S}\}$ as

[TABLE]

Our normalization convention is $\mathrm{Tr}(\mathcal{S}\mathcal{S}^{\prime})=q^{N}\delta_{\mathcal{S}\mathcal{S}^{\prime}}$, so that $a_{\mathcal{S}}(t)=q^{-N}\mathrm{Tr}(\mathcal{O}(t)\mathcal{S})$. The squared coefficient $a_{\mathcal{S}}(t)^{2}$ evolves as

[TABLE]

where $\mathcal{U}$ is a layer of $m$-site unitaries that were applied at time $t-1$. In the second line, we have written $\mathcal{U}=\prod_{\boldsymbol{r}}U_{\boldsymbol{r}}$ where $\boldsymbol{r}$ is the coordinate of disjoint, $m$-site clusters on which the unitary $U_{\boldsymbol{r}}\in\mathrm{U}(q^{m})$ acts, and we have also decomposed $\mathcal{S}=\prod_{\boldsymbol{r}}S_{\boldsymbol{r}}$ as a product of basis elements acting on these $m$-site clusters. These operators are normalized according to $\mathrm{tr}[S_{\boldsymbol{r}}S_{\boldsymbol{r}}^{\prime}]=q^{m}\delta_{S_{\boldsymbol{r}},S_{\boldsymbol{r}}^{\prime}}$ and $\operatorname{tr}[S_{\boldsymbol{r}}]=q^{m}\delta_{S_{\boldsymbol{r}},1}$. The Haar average of the above expression is given by

[TABLE]

And so, the Haar-averaged $a_{\mathcal{S}}(t)^{2}$ evolves linearly

[TABLE]

with the real, symmetric matrix

[TABLE]

Averaging again over the unitaries applied in the previous timesteps gives an equation for $P_{\mathcal{S}}(t)\equiv\overline{a_{\mathcal{S}}^{2}(t)}$

[TABLE]

which is formally a master equation for a fictitious Markov process Dahlsten et al. (2007); Žnidarič (2008); at a given time there is a single string $\mathcal{S}$ which is updated stochastically in each time step, via local updates involving a cluster of $m$ sites. From the form of $W_{\mathcal{S}\mathcal{S}^{\prime}}$ we see that the local update on $m$ sites is performed by replacing a non-trivial generator on the cluster randomly by any one of the $q^{2m}-1$ non-trivial generators. We emphasize that this fictitious Markov process is not the true unitary dynamics of the operator $\mathcal{O}(t)$.

This fictitious classical stochastic process dramatically simplifies through the following observations. We focus here on one spatial dimension with updates on bonds. First, observe that the matrix elements $W_{\mathcal{S}\mathcal{S}^{\prime}}$ only depend on the support of the generators $\mathcal{S}$ and $\mathcal{S}^{\prime}$, so that (119) gives rise to a simpler Markov process for the binary occupation number $n(x)$, which is 1 if the corresponding generator has support at site $x$ and 0 otherwise (if $\mathcal{S}$ acts as the identity at $x$). Formally, the probability distribution of the occupation numbers is given by

[TABLE]

where the prime indicates that the sum is only over strings $\mathcal{S}$ that are compatible with the configuration $n(x)$. Further, the endpoint of the string observes an autonomous Markovian dynamics. Since $m=2$, updates involving the endpoint either include the site to the right of it which is empty or that to the left which may be empty or full. The dynamical rule above implies that the probabilities for the position of the endpoint after the update are independent of whether the leftward site was initially occupied or empty. Formally the probability distribution for the position of the endpoint in this fictitious dynamics is

[TABLE]

which is precisely $\overline{\rho(x,t)}$, as defined in Sec 2. Therefore, for an endpoint at $x$ or $x+1$, a single update applied to the sites $x$ and $x+1$ leaves the endpoint at $x$ with probability $p=(q^{2}-1)/(q^{4}-1)=1/q^{2}+1$ and at $x+1$ with probability $1-p$. This establishes the claim in Sec. 2 for the evolution of $\overline{\rho(x,t)}$ in a single timestep.

Appendix B Velocity and diffusion constant for lattice diffusion equation

Our layout of the evolution operator is such that the local unitaries alternate between even and odd bonds. In other words, a bond at time step $t$ is either at the left or right of the bond at $t-1$. Thus, it suffices to count the left and right moves to specify the position of the right end-bond of $X_{0}(t)$. As described in the main text the probability of a left move is $p$. Let $u\geq 0$ be the number of right moves, and $v\geq 0$ be the number of left moves. We have $u+v=t$, and $u-v$ (or $u-v\pm 1$) is the spatial coordinate of the right endpoint. Therefore, the probability distribution of the position of the right end bond is

[TABLE]

This is correctly normalized since $\sum_{u+v=t}f(u,v)=(1-p+p)^{t}=1$. Then, the probability that a site $x$ is left to the right end of $X_{0}(t)$ is

[TABLE]

where $\Phi$ is the cumulative density function of the normal distribution, and

[TABLE]

Appendix C Noisy diffusion equation

Starting with Eq. 29, WLOG rescale space so $D=1$ and set $v=0$ by going to the moving frame. Let $\rho_{0}$ be the solution without noise, $\rho_{0}=(4\pi t)^{-1/2}e^{-x^{2}/4t}$. In terms of the Green’s function

[TABLE]

The centre of mass position of the wavepacket within a given realisation is ${x_{\text{cm}}=\int_{x}x\rho(x,t)}$, so, if $\overline{x_{cm}}$ is the centre of mass position averaged over realisations,

[TABLE]

Averaging over the noise with ${\left\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\right\rangle}=\lambda{\delta(x-x^{\prime})\delta(t-t^{\prime})}$,

[TABLE]

For the leading order scaling, we replace $\rho$ with $\rho_{0}$ on the right hand side. Then dimensional analysis applied to the integral gives

[TABLE]

This is a statistical variation in $x_{\text{cm}}$ of order $t^{1/4}$, in agreement with the heuristic argument and with Ref. Saul et al. (1992). This variation is small compared to the width of $\rho$, indicating that $\rho-\rho_{0}\ll\rho$ at late times. The typical size of $\partial_{x}\rho$ near the peak is $O(1/t)$, so ${x_{\text{cm}}-\overline{x_{\text{cm}}}}\sim t^{1/4}$ corresponds to $\rho-\rho_{0}\sim t^{-3/4}$, as compared with $\rho\sim t^{-1/2}$. The approximation above is therefore self-consistent.

Appendix D Random Clifford operators

Here we review that the left- and right- invariant probability distribution over the Clifford group on $n$ $q$-dimensional qudits is a unitary 2-design when $q$ is a prime number. In other words, for a qudit of prime power dimension $q^{n}$, the unitary group $U(q^{n})$ has a finite subgroup that is a unitary 2-design, and there is a linear operator basis that remains closed under conjugations by this subgroup. This is a well-known result DiVincenzo2001, but we include it here for readers’ convenience.

To define the Clifford group, we first need the Pauli group. Define $X=\sum_{j=0}^{q-1}\ket{j+1\mod q}\bra{j}$ and $Z=\sum_{j=0}^{q-1}e^{2\pi ij/q}\ket{j}\bra{j}$. Then, the Pauli group is the subgroup of $U(q^{n})$ generated by matrices $X_{1},Z_{1},\ldots,X_{n},Z_{n}$ where

[TABLE]

The Clifford group is defined to be the normalizer of the Pauli group in $U(q^{n})$. The Pauli group quotiented out by its center $\langle\omega=e^{2\pi i/q}\rangle$ is abelian since $XZX^{\dagger}=\omega^{-1}Z$, and is isomorphic to the additive group $\mathbb{Z}_{q}^{2n}$. We define $P_{v}$ for $v\in\mathbb{Z}_{q}^{2n}$ to be an element of the Pauli group (Pauli operator) as

[TABLE]

The center of the Pauli group is also contained in the center of $U(q^{n})$, and therefore the conjugation action by the Clifford group on the Pauli group induces an action on $\mathbb{Z}_{q}^{2n}$. It turns out that this group $\mathcal{S}$ of action consists precisely of those that preserves the symplectic form

[TABLE]

over $\mathbb{Z}_{q}$.

A probability distribution $\nu$ of unitary matrices to form a 2-design means that

[TABLE]

where $\mu$ is the Haar probability distribution over $U(q^{n})$, and $U^{*}$ is the complex conjugate of $U$. Tautologically, the Haar distribution is a 2-design. This is equilvalent to having that

[TABLE]

for any $q^{n}\times q^{n}$ matrices $O$ and $O^{\prime}$. Since Pauli operators (the elements of the Pauli group defined above) generates over the complex numbers the full operator algebra, it is enough to have Eq. (133) with $O$ and $O^{\prime}$ being Pauli operators.

Let $\nu$ be the left-invariant (hence right-invariant) probability distribution over the Clifford group. This is the uniform distribution over the finite Clifford group. Consider a $\mathbb{C}$-linear map $\Pi_{\nu}$ on the set of operators defined by

[TABLE]

Since $\nu$ is a left-invariant distribution over a group of unitaries, $\Pi_{\nu}$ is a projector (which is hermitian under the Hilbert-Schmidt inner product). Since the Clifford group includes the Pauli group, we have for arbitrary $a,b\in\mathbb{Z}_{q}^{2n}$

[TABLE]

The use of inverse $P_{b}^{\dagger}$ here instead of $P_{b}$ is for notational convenience later.

Now, observe that for any nonzero $x,y\in\mathbb{Z}_{q}^{2n}$ there exists a symplectic transformation $S\in\mathcal{S}$ such that $y=Sx$. For this step, it is essential that $q$ is prime. By the right-invariance of $\nu$ by $S$, we see $\Pi_{\nu}(P_{x}\otimes P_{x}^{\dagger})=\Pi_{\nu}(P_{y}\otimes P_{y}^{\dagger})$. This implies that

[TABLE]

We claim that this is a linear combination of the identity operator and the swap operator $\mathcal{F}=\sum_{u,v=0}^{q-1}\ket{u}\bra{v}\otimes\ket{v}\bra{u}$. This is easily verified once we expand $\mathcal{F}$ in the Pauli operator basis; using $\operatorname{Tr}(\mathcal{F}O\otimes O^{\prime})=\operatorname{Tr}(OO^{\prime})$ for any $q^{n}\times q^{n}$ matrices $O$ and $O^{\prime}$, we see that $\mathcal{F}\propto\sum_{x\in\mathbb{Z}_{q}^{2n}}P_{x}\otimes P_{x}^{\dagger}$.

The identity operator and the swap operator commute with $U\otimes U$ where $U\in U(q^{n})$. This implies that $\Pi_{\nu}(P_{a}\otimes P_{b}^{\dagger})$ commutes with $U\otimes U$, and hence is equal to $\Pi_{\mu_{\text{Haar}}}\circ\Pi_{\nu}(P_{a}\otimes P_{b}^{\dagger})$. By the right-invariance of the Haar distribution $\mu_{\text{Haar}}$, we conclude that Eq. (133) is proved.

When $q$ is not prime, any probability distribution over the Clifford group fails to be a unitary 2-design. Let $n=1$. Since the image of $\Pi_{\mu}$ is a linear combination of the identity and the swap, we must have (see App. G below)

[TABLE]

When $q=6$, there are non-identity Pauli operators $P$ and $Q$ such that $P^{2}=I$ and $Q^{3}=I$. By Eq. (136), we have $\Pi_{\mu}(P\otimes P^{\dagger})=\Pi_{\mu}(Q\otimes Q^{\dagger})\neq 0$. However, $\Pi_{\nu}(P\otimes P^{\dagger})$ is a linear combination of Pauli operators, each of which squares to identity, whereas $\Pi_{\nu}(Q\otimes Q^{\dagger})$ is a linear combination of those that cube to identity, so they cannot be equal.

Appendix E Anomalous behaviour of the

front for $\phi=0$

Above we noted that for sufficiently large $p$, $p>p_{c}$, the lattice growth process which we consider has anomalous behaviour when the front is oriented parallel to a lattice plane. This is a known phenomenon in various lattice growth models in discrete time which have synchronous parallel updates and is well understood in terms of directed percolation Richardson (1973); Durrett and Liggett (1981); Savit and Ziff (1985); Kertész and Wolf (1989); Krug et al. (1990).

In the regime $p>p_{c}$ the lattice-aligned ($\phi=0$) front has a speed $v(\phi=0)=2$ which is precisely the maximum possible speed allowed by causality. In this regime the front is pinned to the ‘light front’ and is not rough (i.e. the width is of order one). (Exactly at $p_{c}$, the aligned front is logarithmically rough Kertész and Wolf (1989).) For our lattice model it appears that $p_{c}\lesssim 2$.

This phenomenon is easily understood via a correspondence with directed percolation Savit and Ziff (1985). First consider a straight, lattice-aligned front in the trivial deterministic limit $p=0$ ($q=\infty$). Apart from possibly on the first time step, this flat front advances by two lattice spacings every period: the front keeps pace with the ‘light cone’ which is the line $x=2t$. Let $\tilde{n}(y,t)=0,1$ denote the occupation numbers of the column of sites at the lightcone: $\tilde{n}(y,t)$ is the occupation number of the site at position $(2t,y)$ at time $t$. When $p=0$ we have $\tilde{n}(y,t)=1$. We are interested in the density $\left\langle\tilde{n}\right\rangle$ (averaged over $y$) at late times when $p$ is nonzero. If this density remains finite, that means the front has an $O(1)$ width, and is attached to the light cone. If it instead tends to zero, the front detaches from the light cone, and we expect to recover standard KPZ roughening. Note that, in order to determine $\tilde{n}$ at time $t+1$, it is sufficient to know only $\tilde{n}$ at time $t$. The dynamics of the occupation numbers $\tilde{n}(y,t)$ are as follows. Under a horizontal dimer update (which advances the lightfront) each occupied $y$ has a chance $(1-p)$ of becoming unoccupied. Under a vertical update pairs of adjacent $y$ undergo the pairwise update described in the main text. This allows occupied sites to ‘reproduce’. This is therefore a birth-death process of the directed percolation type Cardy (1996). When $p$ is large the death rate is small and the reproduction rate is large, and the process is in an ‘active’ phase with $\left\langle\tilde{n}\right\rangle>0$, while when $p$ is small the population of occupied sites dies out.

Appendix F Shape of a spreading droplet for weakly varying $v(\phi)$

Consider an asymptotic front shape described by the parameterized curve $(\theta,r_{t}(\theta))$ in polar coordinates, which grows simply by rescaling: $r_{t}(\theta)=t\times r(\theta)$. Let $\phi(\theta)$ be the angle of the front’s normal (to the $x$ axis) at polar position $\theta$. The radial growth rate is ${\dot{r}_{t}(\theta)=v(\phi(\theta))/\cos[\phi(\theta)-\theta]}$. Since the curve grows by rescaling we have $\partial_{\theta}[\dot{r}_{t}(\theta)/r_{t}(\theta)]=0$. Note that

[TABLE]

Combining these gives Wolf (1987)

[TABLE]

Therefore at a location where $r(\theta)$ is smooth we either have $\phi^{\prime}(\theta)=0$, i.e. a straight segment, or

[TABLE]

If the solution is everywhere smooth then the above equation must be satisfied everywhere. (Such solutions exist for sufficiently weakly varying $v$.) It is straightforward to solve this equation in powers of $w$:

[TABLE]

We find that the RHS involves only total derivatives of periodic functions. (Just from looking at Eq. 139 this is at first sight surprising since it emerges from various cancellations.) Therefore, integrating the right hand side according to (137) gives a periodic $r(\theta)$.

For a formal explanation for why the expansion of $\tan(\phi-\theta)$ contains only total derivatives of periodic functions, consider a flow in the space of functions $v(\phi)$ which interpolates between the function of interest and the trivial function $v(\phi)$=const. Let $v_{1}(\phi)$ and $v_{2}(\phi)$ be two functions that are infinitesimally close on this flow and let $\phi_{1}(\theta)$ and $\phi_{2}(\theta)$ be the corresponding solutions. Assuming that $\phi_{1}(\theta)$ is periodic and corresponds to a periodic $r(\theta)$ we show that this property is inherited by $\phi_{2}(\theta)$ to order $\phi_{2}-\phi_{1}$. Using (137), (139) we obtain

[TABLE]

As required, the RHS is indeed the total derivative of a periodic function (note that $\phi_{2}$ does not appear on the RHS). Integrating along the flow then establishes the property for general $v(\phi)$ at the formal level — i.e. assuming that the solution evolves smoothly during the flow.

Appendix G Haar average formula

Here we review a standard formula for the average of matrix elements of unitary matrix with respect to the Haar probability measure $\mu$ on $U(N)$. Let us abbreviate $\int_{U(N)}\mathrm{d}\mu(U)$ as $\mathbb{E}_{U}$. We are going to prove that

[TABLE]

where $F$ is the swap operator on $(\mathbb{C}^{N})^{\otimes 2}$. Evaluating a particular matrix element, we have

[TABLE]

Proof of Eq. (144).

The average is a matrix on $\mathcal{H}=(\mathbb{C}^{N})^{\otimes 2}$ that commutes with every $U^{\otimes 2}$. Hence, the average is block-diagonal in the basis where the representation of $U(d)$ is block-diagonal. The irreps appearing in $\mathcal{H}$ are the symmetric subspace and the anti-symmetric subspace. In each irrep, the average must be proportional to the identity $I_{\pm}$ by the Schur’s lemma, and we need to evaluate the trace in order to determine the constant of proportionality. The projection onto the (anti-)symmetric subspace is $(I\pm F)/2$ where $F$ is the swap operator: $F\ket{ac}=\ket{ca}$. So the trace is

[TABLE]

This must be equal to $C_{\pm}\operatorname{Tr}(I_{\pm})=C_{\pm}N(N\pm 1)/2$ Therefore, the average is equal to $\sum_{s=\pm}C_{s}(I+sF)/2.$ ∎

Appendix H Proof of Eq. (106)

Let $N$ be the Hilbert space dimension of $n$ $q$-dimensional qudits; $N=q^{n}$. For any $N\times N$ unitary $U$, denote by $U^{\otimes t,t}$ the tensor product $(U\otimes U^{*})^{\otimes t}$, where $U^{*}$ is the complex conjugate of $U$. Let $\mu$ be the Haar probability distribution on $U(N)$, and define for any probability distribution $\nu$ on $U(N)$, a real number

[TABLE]

Here, $\left\|{\cdot}\right\|_{\infty}$ denotes the maximum singular value. Due to left and right invariance of $\mu$, it follows that $\Pi_{\mu}^{2}=\Pi_{\mu}=\Pi_{\nu}\Pi_{\mu}=\Pi_{\mu}\Pi_{\nu}$. ($\Pi_{\nu}$ is not in general a projector.) Therefore,

[TABLE]

where $\nu^{*m}$ is the $m$-fold convolution of $\nu$, i.e., $\nu^{*m}$ is the distribution of the product $U_{1}U_{2}\cdots U_{m}$ when every $U_{i}$ obeys distribution $\nu$.

Now, let $\nu$ be the distribution on $U(N)$ obtained by applying one layer of even bond local Haar random unitaries ($U(q^{2})$) and then one layer of odd bond local Haar random unitaries. Brandao-Harrow-Horodecki’s result Brandao et al. (2016) implies that

[TABLE]

Their theorem does not directly cover this, but they have lemmas that are good enough for our purpose; eq. (48) of the CMP version is what we actually need.

Consider $f(U)=N^{-2}(\operatorname{Tr}UXU^{\dagger}YUXU^{\dagger}Y)^{2}\geq 0$ where all the matrices $U,X,Y$ are $N\times N$. $f(U)$ can be thought of as $\bra{\tilde{X}}U^{\otimes 4,4}\ket{\tilde{Y}}$ for some vectors $\ket{\tilde{X}}$ and $\ket{\tilde{Y}}$. Assume $\operatorname{Tr}(X)=\operatorname{Tr}(Y)=0$, but $\operatorname{Tr}(X^{2})=\operatorname{Tr}(Y^{2})=N$. Then, the Euclidean norms of $\ket{\tilde{X}}$ and $\ket{\tilde{Y}}$ are both $N^{2}$. Normalizing so that $\ket{X}:=\ket{\tilde{X}}/N^{2}$ and $\ket{Y}:=\ket{\tilde{Y}}/N^{2}$, we can write $f(U)=N^{2}\bra{X}U^{\otimes 4,4}\ket{Y}$.

By Eqs. (148) and (149), we have

[TABLE]

If $\left|{\mathbb{E}_{U\sim\mu}f(U)}\right|\leq q^{-cn}$, then $\mathbb{E}_{U\sim\nu^{*m}}f(U)\leq 2q^{-cn}$ whenever $m/n\geq(c+2)M_{4,q}\log q$.

Hastings’ Schwinger-Dyson trick Hastings gives

[TABLE]

Therefore, whenever $m/n\geq 4M_{4,q}\log q$, we have

[TABLE]

Appendix I A mean field approximation

Ref. Aleiner et al. (2016) argued, on the basis of Keldysh perturbation theory, that in various circumstances the out-of-time-order correlator would satisfy a traveling wave equation such as the Fisher-KPP equation (the details of this equation depending on the physical system). An example is the Fisher-KPP equation itself:

[TABLE]

The key feature is $\lambda$ term, which means that if ${\mathcal{C}}$ is ‘seeded’ with a small nonzero value, it will increase to a value close to one on a timescale of order $\lambda^{-1}$ (and then saturate). This equation has stable solutions describing a front propagating with a speed $v_{B}=2\sqrt{D\lambda}$. This front does not broaden.

This phenomenology is very different from the picture which we have obtained from the random circuit and the mapping to classical growth processes. Recall that in 1D we related ${\mathcal{C}}$ to a homogeneous (linear) equation, and in higher dimensions we found that ${\mathcal{C}}$ was not governed by a partial differential equation.

The purpose of this Appendix is to show that a traveling wave picture can emerge from our mappings if we make a certain mean field approximation. This mean field approximation is not valid in the systems we have studied — it is an uncontrolled approximation which does not capture the true behavour either at short or at large times. However in variant models a small parameter could be present which justified the mean field approximation up to some finite but large timescale. In this situation we expect that mean field will nevertheless break down at asymptotically long times, with the front eventually roughening in the manner discussed in the text.

Recall that for the random circuit we have

[TABLE]

where $n(x,t)$ is the occupation number in the fictitious classical cluster growth problem. Let us consider the joint probability distribution $P(\{n\};t)$ for this occupation number. This distribution involves nontrivial correlations between sites which are crucial for capturing the correct asymptotic behaviour. Nevertheless let us explore the mean field approximation in which we pretend all sites are independent, $P(\{n\};t)=\prod_{x}P_{x}(n(x);t)$, with

[TABLE]

For simplicity, consider a model on the hypercubic lattice in $d$ dimensions (with coordination number $z=2d$) in which unitaries (‘updates’) are applied to bonds in a Poissonian fashion at rate $\Gamma/2$ per bond. This continuous time protocol does not change the basic point but it simplifies the equations. Write $m=\left\langle n\right\rangle$. Note that if we update a bond which contains at least one fictitious particle, the subsequent (conditionally) averaged density on that bond is $1-p$. This implies

[TABLE]

where the first term is the probability that site $x$ does not receive an update in the interval $\Delta t$. Making the mean field approximation, $\left\langle\delta_{n(x),0}\delta_{n(y),0}\right\rangle$ factorizes into $(1-m(x))(1-m(y))$, so that

[TABLE]

The first term on the right is a ‘death rate’. The second term is spreading. The third term is a correction to overcounting in the second term. An analogous equation could be written down for the regular circuit considered in the main text, but we would have to use discrete time.

Eq. 158 is a lattice traveling wave equation. This is most apparent if we make a formal expansion in the lattice spacing $a$ to second order (valid, given the approximations already made, if the solution is slowly varying). Recalling $p=1/(q^{2}+1)$ and Eq. 155,

[TABLE]

This differs from the Fisher-KPP equation only in the ${\mathcal{C}}$–dependence of the diffusion constant, and we expect similar properties.

Above, the mean field limit was an unjustified formal approximation. We could of course construct random circuit models in which the (lattice) mean field approximation was quantitatively accurate up to a large time, for example by using long range interactions or a large coordination number to reduce the effect of correlations. However at long times, in physical dimensionalities, we expect the front to roughen so that the mean field traveling wave picture breaks down. In (unphysically) high dimensions, mean field may be valid even at late times (recall that the phase diagram of the KPZ equation allows for a non-roughening phase in high dimensions, as discussed in the text).