Anomalous Subdiffusion from Subsystem Symmetries

TL;DR

This paper introduces automaton quantum circuits with subsystem symmetries that exhibit anomalous subdiffusion, revealing slower-than-diffusive charge spreading and non-Gaussian distributions, supported by numerical simulations and operator hydrodynamics.

Contribution

It characterizes automaton quantum circuits with subsystem symmetries and demonstrates their role in causing anomalous subdiffusion and non-Gaussian charge distributions in quantum dynamics.

Findings

Charge autocorrelator decays as log(t)/√t in 2D with line symmetries

Charge autocorrelator decays as 1/t^{3/4} in 3D with planar symmetries

Numerical evidence of non-Gaussian charge distributions and quantum chaos signatures

Abstract

We introduce quantum circuits in two and three spatial dimensions which are classically simulable, despite producing a high degree of operator entanglement. We provide a partial characterization of these "automaton" quantum circuits, and use them to study operator growth, information spreading, and local charge relaxation in quantum dynamics with subsystem symmetries, which we define as overlapping symmetries that act on lower-dimensional submanifolds. With these symmetries, we discover the anomalous subdiffusion of conserved charges; that is, the charges spread slower than diffusion in the dimension of the subsystem symmetry. By studying an effective operator hydrodynamics in the presence of these symmetries, we predict the charge autocorrelator to decay ($i$) as $\log(t)/\sqrt{t}$ in two dimensions with a conserved $U(1)$ charge along intersecting \emph{lines}, and ($ii$) as $1/t^{3/4}$ in three spatial dimensions with intersecting \emph{planar} $U(1)$ symmetries. Through large-scale studies of automaton dynamics with these symmetries, we numerically observe charge relaxation that is consistent with these predictions. In both cases, the spatial charge distribution is distinctly non-Gaussian, and reminiscent of the diffusion of charges along a fractal surface. We numerically study the onset of quantum chaos in the spreading of local operators under these automaton dynamics, and observe power-law broadening of the ballistically-propagating fronts of evolving operators in two and three dimensions, and the saturation of out-of-time-ordered correlations to values consistent with quantum chaotic behavior.

Tables

Table 1. Table 1: Subdiffusive Dynamics with U ​ ( 1 ) 𝑈 1 U(1) Subsystem Symmetries – Summary of analytical results for the subdiffusive decay of charges in two and three spatial dimensions, in the presence of overlapping line-like and planar U ​ ( 1 ) 𝑈 1 U(1) symmetries. These predictions for line symmetries in d = 2 𝑑 2 d=2 and planar symmetries in d = 3 𝑑 3 d=3 are consistent with numerical studies of automaton dynamics with these symmetries.

Subsystem Symmetry	Charge Decay
	$d=2$	$d=3$
Line Symmetry	$\log \left(t\right)/\sqrt{t}$	${\log }^2\left(t\right)/\sqrt{t}$
Plane Symmetry	—	$1/t^{3/4}$

Table 2. Table 2: Properties of a Generic Automaton Unitary Dynamics , as compared to Clifford dynamics, a generic integrable dynamics, and the dynamics of a quantum chaotic system, each with a D 𝐷 D -dimensional Hilbert space. The typical recurrence time for a state evolving under an integrable dynamics in one dimension grows polynomially in the system size, since each quasiparticle takes a time O ​ ( log ⁡ D ) 𝑂 𝐷 O(\log D) to return to its initial position; furthermore, ( † † \dagger ) assumes a translationally invariant, Floquet Clifford circuit von Keyserlingk et al. ( 2018 ) . Operator spreading, as quantified by out-of-time-ordered correlations, reveals a power-law broadening of the front of an evolving operator in time for automaton circuits in this work, which is close to what has been observed in random unitary dynamics in various spatial dimensions Nahum et al. ( 2018 ) , unlike what is observed for integrable dynamics Gopalakrishnan et al. ( 2018 ) . The computational cost of simulating one timestep of evolution, by applying O ​ ( log ⁡ D ) 𝑂 𝐷 O(\log D) gates is indicated, with results for integrable systems from Ref. Valiant ( 2002 ); Osborne ( 2006 ) . The memory overhead is the memory cost to store an evolving state under the indicated dynamics at long times; this is generally O ​ ( D ) 𝑂 𝐷 O(D) as a product initial state will generally become volume-law entangled under each of these dynamics. For Clifford dynamics, the memory cost is significantly reduced for initial states in the Pauli basis Gottesman ( 1998 ) . For integrable dynamics ( ‡ ‡ \ddagger ), the operator entanglement is believed to grow logarithmically in time Alba et al. ( 2019 ); Prosen and Znidaric ( 2007 ) , which reduces the cost of storing an evolving state, though the state may eventually become volume-law entangled.

Features	$\begin{array}{c}\text{Integrable}\\ \mathrm{𝐃𝐲𝐧𝐚𝐦𝐢𝐜𝐬}\end{array}$	$\begin{array}{c}\mathrm{𝐂𝐥𝐢𝐟𝐟𝐨𝐫𝐝}\\ \mathrm{𝐃𝐲𝐧𝐚𝐦𝐢𝐜𝐬}\end{array}$	$\begin{array}{c}\mathrm{𝐀𝐮𝐭𝐨𝐦𝐚𝐭𝐨𝐧}\\ \mathrm{𝐃𝐲𝐧𝐚𝐦𝐢𝐜𝐬}\end{array}$	$\begin{array}{c}\mathrm{𝐐𝐮𝐚𝐧𝐭𝐮𝐦}\mathrm{𝐂𝐡𝐚𝐨𝐭𝐢𝐜}\\ \mathrm{𝐃𝐲𝐧𝐚𝐦𝐢𝐜𝐬}\end{array}$
$\begin{array}{c}\mathrm{𝐺𝑒𝑛𝑒𝑟𝑎𝑡𝑒𝑠}\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\text{-}\mathrm{𝐿𝑎𝑤}\\ \mathrm{𝑆𝑡𝑎𝑡𝑒}\mathrm{𝐸𝑛𝑡𝑎𝑛𝑔𝑙𝑒𝑚𝑒𝑛𝑡}\end{array}$	Yes	Yes	$\begin{array}{c}\text{For states in}\\ \text{a specific basis}\end{array}$	Yes
$\begin{array}{c}\mathrm{𝐺𝑒𝑛𝑒𝑟𝑎𝑡𝑒𝑠}\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\text{-}\mathrm{𝐿𝑎𝑤}\\ \mathrm{𝑂𝑝𝑒𝑟𝑎𝑡𝑜𝑟}\mathrm{𝐸𝑛𝑡𝑎𝑛𝑔𝑙𝑒𝑚𝑒𝑛𝑡}\end{array}$	$\begin{array}{c}\text{Yes}\end{array}$	No	Yes	Yes
$\begin{array}{c}\text{<em>Recurrence Time </em>}\\ \text{(product initial state in 1D)}\end{array}$	$\mathrm{poly}\left(\log D\right)$	$\log D$ ($†$)	$D$	$\exp \left(\lambda D\right)$
Operator Spreading	$\begin{array}{c}\text{Ballistic Growth,}\\ \text{Diffusive Broadening}\end{array}$	$\begin{array}{c}\text{Ballistic Growth,}\\ \text{No Broadening}\end{array}$	$\begin{array}{c}\text{Ballistic Growth,}\\ \text{Power-Law Broadening}\end{array}$	$\begin{array}{c}\text{Ballistic Growth,}\\ \text{Power-Law Broadening}\end{array}$
$\begin{array}{c}\text{<em>Operations to simulate </em>}\\ \text{<em>one timestep of evolution</em>}\end{array}$	$\mathrm{poly}\left(\log D\right)$	$\mathrm{poly}\left(\log D\right)$	$D$	$D^2$
Memory Cost	$D$ ($‡$)	${\log }^{2}D$	$D$	$D$

Equations

$$U|m⟩=e^{i{\theta }_m}|\pi (m)⟩,$$

$$\mathcal{O}\to U^{†}\mathcal{O}U$$

$$\mathrm{CCNOT}(\theta {)}_{123}\equiv 1-{\Pi }_{12}+{\Pi }_{12}{e}^{i\theta }{X}_3.$$

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

$${\mathcal{S}}_n(t)=U(t){\mathcal{S}}_{n}U(t{)}^{†}=\sum _{m=0}^{D-1}{a}_{nm}(t){\mathcal{S}}_m.$$

$$a_{nm}(t)=\frac{1}{D}\sum _{k=0}^{D-1}(-1{)}^{\mathbf{\pi }(k)\cdot \mathbf{m}+\mathbf{k}\cdot \mathbf{n}}.$$

$$|n(t)⟩\equiv \sum _{n=0}^{D-1}{a}_{nm}(t)|m⟩,$$

$$W(t)=H^{\otimes N}U(t)H^{\otimes N},$$

$${\rho }_A(n)\equiv {\mathrm{Tr}}_B|\mathbf{n}(t)⟩⟨\mathbf{n}(t)|.$$

$$\displaystyle\overline{\,\mathrm{Tr}\,\rho_{A}(n)^{2}\,}=\left\{\begin{array}[]{cc}1&(n=0)\\ &\\ \displaystyle{D_{A}^{-1}+D_{B}^{-1}-D^{-1}}&(n\neq 0)\end{array}\right.,$$

$$U(t{)}^{†}\mathcal{O}U(t)=\sum _{m=0}^{D-1}{e}^{i({\theta }_m-{\theta }_{\tau \pi (m)})}|{\pi }^{-1}\tau \pi (m)⟩⟨m|$$

$$\frac{D!}{2^{D/2}(D/2)!}\stackrel{D\to \infty }{\sim }{D}^{(D/2)\left[1+O({\log }^{-1}D)\right]},$$

$$t_{\mathrm{prod}}\sim D/2$$

$$t_{\mathrm{rand}}\stackrel{D\to \infty }{\sim }\exp \left[\lambda \sqrt{D/\log D}\right],$$

$$F(\mathbf{r},t)\equiv \frac{1}{D}\mathrm{Tr}[{\mathcal{O}}_{\mathbf{r}}(t){\mathcal{O}}_0^{\prime }(0){\mathcal{O}}_{\mathbf{r}}(t){\mathcal{O}}_0^{\prime }(0)],$$

$$F(\mathbf{r},t)=1-2\sum _{\{\mathcal{S},{\mathcal{O}}_0^{\prime }\}=0}{a}_{\mathcal{S}}(t{)}^2.$$

$$\displaystyle F(\boldsymbol{r},t)\sim\left\{\begin{array}[]{cc}0&\,\,\,\,v_{B}t\gg|\boldsymbol{r}|\\ &\\ 1&\,\,\,\,v_{B}t\ll|\boldsymbol{r}|\end{array}\right..$$

$$\overline{\mathrm{Tr}[U^{†}{Z}_{0}U{Z}_0{U}^{†}{Z}_{0}U{Z}_0]}=0,$$

$$a_{\mathcal{S}}(t)$$

$$\frac{1}{D}\sum _{\mathbf{r}}\mathrm{Tr}[Z_{\mathbf{r}}\mathcal{O}(t)]=\sum _{\mathbf{r}}{a}_{Z_{\mathbf{r}}}(t)=const.$$

$$\frac{1}{D}\mathrm{Tr}[Z_{\mathbf{r}}(t)Z_0]\sim e^{-(r^2/D_{0}t)}/(D_{0}t{)}^{d/2},$$

$$C_x=\sum _y{Z}_{x,y}\text{and}{C}_y=\sum _x{Z}_{x,y},$$

$$U_{\mathrm{flip}}$$

$$P\equiv \frac{1}{8}(1-Z_1{Z}_2)(1-Z_2{Z}_3)(1-Z_3{Z}_4)$$

$$\displaystyle U_{\mathrm{flip}}^{\dagger}Z_{1}U_{\mathrm{flip}}=\frac{3}{4}Z_{1}+\frac{1}{4}\Big{[}Z_{2}+Z_{3}-Z_{4}\Big{]}$$

$$\displaystyle+\frac{1}{4}\Big{[}Z_{1}Z_{2}Z_{4}+Z_{1}Z_{3}Z_{4}-Z_{1}Z_{2}Z_{3}-Z_{2}Z_{3}Z_{4}\Big{]}.$$

$$G_{\mathbf{r},t}=\frac{1}{D}\mathrm{Tr}[Z_{\mathbf{r}}(t)Z_0(0)].$$

$$G_{0,t}\sim \frac{a\log (t)+b}{\sqrt{t}}.$$

$$\sum _{r_x}{G}_{\mathbf{r},t}={\delta }_{r_y,0},\sum _{r_y}{G}_{\mathbf{r},t}={\delta }_{r_x,0}.$$

$$G_{\mathbf{r},t+1}$$

$$+\lambda \left[G_{\mathbf{r}+x,t}+G_{\mathbf{r}-x,t}+G_{\mathbf{r}+y,t}+G_{\mathbf{r}-y,t}\right]$$

Full text

The dynamics of interacting, quantum many-body systems provide a rich source of open problems in theoretical physics. Recent developments include advances in our understanding of quantum thermalization Deutsch (1991); Srednicki (1994); Rigol et al. (2008) and many-body localization Gornyi et al. (2005); Basko et al. (2006); Nandkishore and Huse (2015); Abanin et al. (2018), which have provided new paradigms for non-equilibrium quantum matter. The study of quantum dynamics has recently been revolutionized by the study of random quantum circuits Nahum et al. (2017, 2018, 2018); Khemani et al. (2018); Rakovszky et al. (2017); von Keyserlingk et al. (2018); Jonay et al. (2018), which have been a source of theoretically tractable problems that shed light on the dynamics of more general quantum systems.

A richer set of quantum dynamical phenomena arise in the presence of conservation laws. A global $U(1)$ symmetry leads to the diffusion of the conserved charges Khemani et al. (2018); Rakovszky et al. (2017). Overlapping conservation laws lead to more striking features in quantum dynamics. For example, the conservation of both the total $U(1)$ charge and dipole moment in one spatial dimension can lead to the breaking of ergodicity and localization, and a “shattering” of the Hilbert space into exponentially many dynamically disconnected sectors Pai et al. (2018); Sala et al. (2019); Khemani and Nandkishore (2019). Most studies of quantum dynamics with or without conservation laws have thus far been limited to systems in one spatial dimension.

In this work, we extend the study of quantum dynamics in two and three spatial dimensions, by introducing a class of quantum dynamics – termed automaton dynamics – for which the evolution of various correlation functions of local operators is classically simulable, and which retain certain key features of the dynamics of a chaotic quantum system. These dynamics share the property that the Heisenberg evolution of a local operator is complex, in a manner that resembles the evolution under a more general, chaotic quantum dynamics; nevertheless, these quantum dynamics are simulable since they do not generate any entanglement when acting on a complete set of product states in an appropriate basis. We identify key features of the most general automaton dynamics – in the Heisenberg evolution of local operators, the recurrence times for an initial state, and in the generation of operator entanglement – that distinguish them from other kinds of classically simulable quantum dynamics, and conjecture that the evolution of local observables under these dynamics can be quantifiably similar to that of a more generic, chaotic quantum system with the same symmetries. Two examples of automaton circuits, which necessarily generate operator entanglement, have been previously studied in one spatial dimension, and one of these is known to be integrable Gopalakrishnan (2018); Gopalakrishnan and Zakirov (2018); Alba et al. (2019).

We apply our understanding of these automaton dynamics to the study of quantum dynamics with subsystem symmetries – symmetries that act along overlapping, subdimensional “manifolds” of the system – and find a rich and unexpected behavior involving anomalous subdiffusion of charge, slower not only than diffusion in that physical dimension, but also slower than diffusion in the dimension of the submanifold acted on by the symmetry. Our motivation to study these dynamics arises from the existence of exotic quantum phases of matter with immobile, fractionalized excitations (“fracton” phases) Chamon (2005); Vijay et al. (2015); Haah (2011); Nandkishore and Hermele (2018) in which these symmetries, such as the conservation of multipole moments of charge Pretko (2017), are emergent properties of these stable phases Vijay et al. (2016), that severely restrict the dynamics of these exotic, fractionalized excitations.

In this work, we specifically consider Floquet dynamics that conserve $U(1)$ charge along intersecting lines in two and three dimensions, and three dimensional dynamics that conserve $U(1)$ charge along planes, with no additional symmetries; effectively, these are infinite-temperature dynamics in which these subsystem symmetries are exact. For a two dimensional system with line-like symmetries, we find that the local charge relaxes as $\log(t)/\sqrt{t}$, whereas for a three dimensional system with planar symmetries it relaxes as $1/t^{3/4}$. These results originate from an analytical explanation in terms of the emergent operator hydrodynamics of the problem, and are consistent with extensive numerical studies of automaton dynamics with these symmetries. Our results are summarized in Table 1. This anomalous subdiffusion is tied to a distinctly non-Gaussian charge distribution, and reminiscent of classical fractional subdiffusion Metzler and Klafter (2000); Henry et al. (2010).

An efficient numerical simulation of dynamical (both time-ordered and out-of-time-ordered) correlation functions in automaton circuits is possible due to the fact that these correlations may be computed using a Monte Carlo sampling procedure. This provides a unique example of a quantum dynamics, where we can measure these correlations in two and three spatial dimensions for a highly entangled operator evolution. In two dimensions, we are able to simulate systems with up to $N=396\times 396$ sites and circuits with up to 32000 layers. In three dimensions, we present results for systems with $N=128\times 128\times 128$ sites and circuit depths up to 21000 layers. Observables other than the conserved charge display apparently “quantum chaotic” behavior; out-of-time-ordered correlations of these observables propagate ballistically, and saturate to values that suggest that the structure of these Heisenberg-evolved operators equilibrates in a manner that is consistent with quantum chaos. The widths of the ballistically-propagating fronts of these operators grow in time as a power-law $t^{\alpha}$ with exponent $\alpha=0.308(18)$ in two dimensions and $\alpha=0.220(5)$ in three dimensions. These should be compared with the predicted values $\alpha=1/3$ and $\alpha=0.24$ in two and three dimensions, respectively, for random unitary quantum circuits Nahum et al. (2018).

This paper is structured as follows. In Sec. I we define and present general features of automaton dynamics, and argue that a generic automaton evolution can lead to the complex Heisenberg evolution of local operators, as quantified by the generation of a large degree of “operator entanglement”. For certain operators, we are able to solve for the Heisenberg evolution explicitly. We then present the numerical algorithm employed for our simulations. Our comparison of automaton evolution to other quantum dynamics is presented in Table 2. In Sec. II we present our results on the anomalous subdiffusion of charges in the presence of various subsystem symmetries. In Sec. III we study the evolution of out-of-time-ordered correlation functions under these automaton dynamics, and show that the Heisenberg evolution of operators other than the conserved charges appears quantifiably similar to that of a quantum chaotic system. We conclude in Sec. IV with a discussion of the implications of our results.

I Automaton Dynamics

We begin by defining the “automaton” dynamics of a quantum system as a unitary evolution that ($i$) does not generate any entanglement in an appropriate basis of product states, but that ($ii$) leads to the non-trivial evolution of local operators under Heisenberg evolution.

More precisely, an automaton unitary operator $U$ acting on an appropriate set of orthonormal product states in a $D$-dimensional Hilbert space – labeled $|m\rangle$, with $m\in\{0,\,\ldots,D-1\}$ – simply permutes these states up to a phase factor, i.e.

[TABLE]

where $\pi\in S_{D}$ is an element of the permutation group on $D$ elements. An automaton unitary will generally create entanglement when acting on product states in a different basis. We refer to Eq. (1) as the “automaton constraint” for the remainder of this section.

While the evolution of the product states $\{|m\rangle\}$ is simple, the Heisenberg evolution of a local operator

[TABLE]

can be complex, in a manner that resembles the Heisenberg evolution of an operator under a more generic, chaotic quantum dynamics; we will quantify this similarity in later sections. For this reason, automaton unitary evolution is of interest as an example of classically simulable, quantum dynamics that may capture broader features of the dynamics of chaotic quantum systems.

Before studying the evolution of local operators under automaton dynamics, it is useful to contrast the simulability of these dynamics with integrable dynamics, Clifford unitary evolution, and the dynamics of a chaotic quantum system. Our comparison is summarized in Table 2. We note, in particular, that automaton dynamics are qualitatively the opposite of a Clifford quantum circuit. In a Clifford unitary dynamics, any Pauli operator (for a spin-(1/2) system, these are the Pauli $X$, $Y$, $Z$ operators, and their products) simply evolves into a product of Pauli operators under Heisenberg evolution. While operator evolution is not complex, Clifford unitary operators can generate a high degree of quantum entanglement when acting on any wavefunction.

An important example of a unitary operator that can generate an automaton quantum circuit – defined as a product of local, automaton unitary operators that each do not generate entanglement in the same basis – is a variant of the controllled-controllled NOT gate, a three-qubit gate which we define as

[TABLE]

Here $\Pi_{12}=|{\uparrow_{1}\uparrow_{2}}\rangle|\langle{\uparrow_{1}\uparrow_{2}}|$ is the projector onto the up-state of spins $1$ and $2$, while $X_{3}$ is the Pauli $X$ operator acting on spin $3$. Since this unitary transformation conditionally flips the third qubit based on the states of the first two qubits, it generates no entanglement in the Pauli $Z$ basis. However, the Heisenberg evolution of a local operator is more complex, e.g. $Z_{3}\rightarrow(1/2)\left[Z_{3}-Z_{1}Z_{3}-Z_{2}Z_{3}-Z_{1}Z_{2}Z_{3}\right]$. We note that ($i$) $\mathrm{CCNOT}(\theta=0)$ is the Toffoli gate, which is universal for classical, reversible computing Nielsen and Chuang (2010), as any invertible Boolean function can be constructed using this gate and ancilla bits, and that ($ii$) the Toffoli gate, and the single-qubit Hadamard gate, define a universal gate set Shi (2002).

I.1 Growth of the Operator Entanglement

We now study the growth of the complexity of an operator under the Heisenberg evolution by an automaton quantum circuit. For simplicity of presentation, and without loss of generality, we restrict our attention to a spin-(1/2) system. Under Heisenberg evolution by a chaotic unitary dynamics, a local Pauli operator $\mathcal{O}$ in a chaotic quantum system will evolve into sums of arbitrary products of Pauli operators as

[TABLE]

where the sum is over all operators $\{\mathcal{S}\}$ that are products of Pauli operators on distinct sites in the system. Since $\mathcal{O}(t)^{2}=1$, the real coefficients $a_{\mathcal{S}}(t)$ satisfy the condition $\sum_{\mathcal{S}}a_{\mathcal{S}}(t)^{2}=1$, for all times $t$.

We now consider an automaton unitary operator $U(t)$ that generates no entanglement when acting on states in the Pauli $Z$-basis. In other words, a projection operator $|m\rangle\langle m|$ in this basis evolves as $U(t)|m\rangle\langle m|U^{\dagger}(t)=|\pi(m)\rangle\langle\pi(m)|$. A simple consequence of this is that a Pauli $Z$ operator will only evolve into a sum of products of other Pauli $Z$ operators under Heisenberg evolution.

While this restricted evolution is a special feature of the automaton dynamics, we now show that under a generic automaton evolution, a Pauli $Z$ operator can grow to develop a weight on any possible operator within this sub-space. More precisely, we demonstrate that the “operator entanglement” – which quantifies the complexity of the growing operator $U^{\dagger}(t)\,Z\,U(t)$ – can develop volume-law scaling under an automaton evolution.

To proceed, we let $\mathcal{S}_{n}$ denote a product of $Z$ operators at sites corresponding to the binary representation of the integer $n\in\{0,\ldots,D-1\}$. For example, the operator $\mathcal{S}_{5}=Z_{1}Z_{3}$ (since the binary representation of $5$ is $10100\cdots$). These operators evolve as

[TABLE]

The “automaton constraint” permits us to determine in Appendix A that the coefficients $a_{nm}(t)$ are exactly

[TABLE]

Here, $\boldsymbol{\pi}(k)$ is the binary vector representation of the integer $\pi(k)$, and $\boldsymbol{\pi}(k)\cdot\boldsymbol{m}$ and $\boldsymbol{k}\cdot\boldsymbol{n}$ denote the dot product of the appropriate vectors. As we show in Appendix A, the fact that the matrix $a_{nm}(t)$ is unitary is a consequence of the unitarity of the original dynamics.

To quantify the growth of the complexity of the Pauli operator $\mathcal{S}_{n}(t)$, it is natural to study its “operator entanglement”, which we define as the entanglement of the evolving wavefunction

[TABLE]

where $|{m}\rangle$ denotes the state of a fictitious $N$-spin system, corresponding to the binary representation of the integer $m$, so that a “$1$” (“[math]”) in this binary representation corresponds to a down (up) spin; for example, for $m=5$, we have the state $|{m}\rangle=|\downarrow\uparrow\downarrow\uparrow\uparrow\uparrow\cdots\rangle$. We observe that the unitary operator $W(t)$, whose matrix elements are $\langle{m}\,|\,W(t)\,|\,{n}\rangle\equiv a_{nm}(t)$, can be written as

[TABLE]

where $H$ is the single qubit Hadamard gate, which acts on Pauli operators as $HXH=Z$, $HZH=X$. In other words, the operator entanglement evolves according to the same automaton circuit in a rotated basis; in this basis, a generic automaton unitary operator will eventually produce a high degree of entanglement.

Instead of demonstrating this explicitly by studying the growth of operator entanglement for a particular automaton evolution, we quantify the operator entanglement generated by a random automaton unitary operator. We bi-partition the state $|\boldsymbol{n}(t)\rangle$ into an $A$ and $B$ subsystem, of Hilbert space dimension $D_{A}$ and $D_{B}=D/D_{A}$, respectively. The reduced density matrix for the $A$ system is defined as

[TABLE]

We determine in Appendix B, that the purity of this density matrix for a random automaton unitary, corresponding to a random permutation $\pi\in S_{D}$, is

[TABLE]

where the line $\overline{\cdots}$ denotes an average over automaton unitary operators that do not generate entanglement in the Pauli $Z$ basis. For any finite subsystem $A$, taking the thermodynamic limit $D\rightarrow\infty$, $D_{A}/D\rightarrow 0$ yields $\overline{\,\mathrm{Tr}\,\rho_{A}(n)^{2}\,}=D_{A}^{-1}$ when $n\neq 0$, which is the maximally-entangled value of the density matrix. We observe that, $\mathrm{Tr}\,\rho_{A}(0)^{2}=1$ since the state $|\boldsymbol{n}=0\rangle$ is an eigenstate of $W(t)$, as this state corresponds to the trivial evolution of the identity operator.

We now argue that time evolution must generate a high degree of operator entanglement for an initial operator that is a product of Pauli $X$ operators. Instead of attempting to solve for their Heisenberg evolution, we observe that any product of these operators – which we label $\mathcal{O}$ – implements a permutation $\tau\in S_{D}$ when acting on states in the $Z$ basis. As a result,

[TABLE]

if the automaton unitary $U(t)$ generates no entanglement in the $Z$ basis. If we let $\theta_{n}=0$ for all $n$, then the Heisenberg-evolved operator $\mathcal{O}(t)=U(t)^{\dagger}\mathcal{O}U(t)$ can only act as an element of its conjugacy class, i.e. as an element of the form $\sigma^{-1}\tau\sigma$, for some $\sigma\in S_{D}$. Conjugacy classes of the permutation group are labeled by their cycle type, defined as the lengths of all cycles in an element of that class. Since $\mathcal{O}$ squares to the identity and is traceless, the permutation $\tau$ consists of $D/2$ independent transpositions. The size of this conjugacy class is

[TABLE]

which is the number of possible operators that $\mathcal{O}(t)$ can evolve into. As this is much larger than the $O(D^{2})$ operators which are simple products of Pauli operators in the system, we conclude that $\mathcal{O}(t)$ must exhibit a high degree of operator entanglement, for a randomly chosen automaton unitary $U(t)$. If the phases in this automaton unitary are non-zero ($\theta_{n}\neq 0$), then $\mathcal{O}$ no longer evolves into a countable set of operators. Nevertheless, the above argument for the growth of the operator entanglement remains valid.

I.2 The Recurrence Time and Quantum Chaos

We compare the recurrence times – when a given initial state returns to itself, so that $|\langle\psi|U(t)|\psi\rangle|\sim O(1)$ – for a typical automaton circuit, with that of a chaotic quantum system. For a chaotic quantum system, the Poincare recurrence time scales exponentially in the Hilbert space dimension of the system $t_{\mathrm{rec}}\sim\exp(\lambda D)$. In contrast, let $U$ be a random automaton operator that generates a single timestep of a Floquet unitary operator $U(t)=U^{t}$, so that $U$ is proportional to a random permutation of product states in a $D$-dimensional Hilbert space. The behavior of the recurrence time for $U(t)$ can be varied. First, the recurrence time for a random product state for which $U$ generates no entanglement will grow as

[TABLE]

as we show in Appendix C. In contrast, the recurrence time for a random state (not necessarily a product state) will be much larger, as this corresponds to the order of a random element of the permutation group $\pi\in S_{D}$ whose asymptotic form as $D\rightarrow\infty$ Goh and Schmutz (1991) gives

[TABLE]

where $\lambda$ is an $O(1)$ constant.

Finally, to compare the operator evolution under an automaton dynamics with the evolution in a chaotic quantum system, we consider the out-of-time-ordered correlation function

[TABLE]

where $\mathcal{O}_{\boldsymbol{r}}(0)$ and $\mathcal{O}^{\prime}_{0}(0)$ are two initially local operators, at the indicated positions. The out-of-time-ordered correlation function probes the structure of the evolving operator $\mathcal{O}_{\boldsymbol{r}}(t)$, which we may expand as a sum of products of Pauli operators as in Eq. (4). For $\mathcal{O}^{\prime}_{0}(0)$ a Pauli string, we observe that

[TABLE]

The quantity $\sum_{\{\mathcal{S},\mathcal{O}^{\prime}_{0}\}=0}a_{\mathcal{S}}(t)^{2}$ is just the probability that $\mathcal{O}_{\boldsymbol{r}}(t)$ has weight on a Pauli operator that anti-commutes with $\mathcal{O}^{\prime}$. For a chaotic quantum system, a typical operator grows ballistically under Heisenberg evolution, with a butterfly velocity $v_{B}$, and within this growing region, the operator is equally likely to develop a weight on Pauli operators that commute or anti-commute with $\mathcal{O}^{\prime}_{0}$ so that $\sum_{\{\mathcal{S},\mathcal{O}^{\prime}_{0}\}=0}a_{\mathcal{S}}(t)^{2}\rightarrow(1/2)$ and

[TABLE]

Therefore, OTOC ($i$) has a ballistically growing front, and ($ii$) within this region the OTOC asymptotically vanishes due to the equilibration of the structure of $\mathcal{O}_{\boldsymbol{r}}(t)$. We will confirm that these features, which hold for chaotic quantum dynamics, are also present for automaton unitary circuits in our numerical studies. We note that as before, one can show using Eq. (6) that for a random automaton dynamics $U$, that

[TABLE]

where $\overline{\cdots}$ again denotes an average over the choice of automaton circuits.

Finally, the front of the OTOC is believed to broaden as a power-law in time $\sim t^{\alpha}$ in a chaotic quantum system, where $\alpha$ is a dimension-dependent exponent. In $d=1$, the ends of an operator perform a biased random walk, leading to the exponent $\alpha=1/2$ Nahum et al. (2018); von Keyserlingk et al. (2018), while in higher dimensions, these exponents are related to probability distributions for classical stochastic growth processes Nahum et al. (2018). We verify that a power-law broadening of the operator front occurs for the OTOC in the automaton circuits we consider.

I.3 Simulating Automaton Evolution

We now discuss how dynamical correlation functions may be efficiently calculated for an automaton evolution, using classical Monte Carlo techniques. Let $U(t)$ be an automaton unitary evolution that generates no entanglement in the Pauli $Z$ basis. We wish to determine the weight $a_{\mathcal{S}}(t)=(1/D)\,\mathrm{Tr}[\mathcal{O}_{\boldsymbol{r}}(t)\mathcal{S}]$ of an evolving Pauli operator $\mathcal{O}_{\boldsymbol{r}}(t)\equiv U(t)^{\dagger}\mathcal{O}_{\boldsymbol{r}}U(t)$ on each basis string $\mathcal{S}$. We observe that

[TABLE]

This quantity can be easily calculated using classical methods since $|n(t)\rangle=U(t)|n\rangle$ remains a product state for any product state $|n\rangle$ in the computational basis. A further speedup is obtained because for strings diagonal in the computational basis (e.g. products of Pauli $Z$ operators), $\langle m|\mathcal{S}|n\rangle=\pm\delta_{mn}$, so one only needs to sample over $|n\rangle$, whereas for a string off-diagonal in the computational basis (e.g. $X$) one only needs to sample over $|n\rangle$ and $|m\rangle=X_{i}|n\rangle$, instead of needing to sample over $|n\rangle$ and $|m\rangle$ independently. This idea of classically sampling a quantum wave function is standard in variational Monte Carlo methods. Here, we extend this idea to sampling a class of entangled quantum operators. The coefficients of the basis strings in an expansion of the highly entangled ‘variational quantum operator’ are generated by the action of the unitary circuit on an initially local operator. Using this method we can study the unitary time evolution of very large quantum systems. In this paper, we look at 2D circuits with up to $396^{2}$ sites and circuit depths of up to 32000 layers, as well as 3D circuits with up to $128^{3}$ sites and circuit depths of up to 21000 layers.

II Results on Automaton Circuits with Subsystem Symmetries

In this section, we study automaton quantum dynamics that possess an extensive set of intersecting, global symmetries. In this case, the conserved charges cannot move in isolation, but engage in a complex, correlated motion. We demonstrate through analytical arguments and numerical studies of these quantum dynamics, that the conserved charges in systems with overlapping ($i$) line-like $U(1)$ symmetries in two spatial dimensions, and ($ii$) planar $U(1)$ symmetries in three spatial dimensions, evolve subdiffusively.

Before proceeding to a discussion of subsystem symmetries, we note that if the unitary circuit of interest possesses a continuous, global symmetry, then there exists a set of conserved charges, which constrain the Heisenberg evolution of operators. For example, for a unitary evolution that preserves the total $Z$ spin – so that $\sum_{\boldsymbol{r}}Z_{\boldsymbol{r}}$ is a conserved operator – the coefficients appearing in the expansion of any operator $\mathcal{O}(t)$ satisfy the constraint

[TABLE]

That is, the existence of a symmetry implies there exists a set of conserved strings which consist of a single Pauli operator. These are known as the conserved charges. By dividing an operator into a conserved and non-conserved part, one can gain a full understanding of the operator dynamics. It will be important to compare to the known results of unitary circuits with a global $U(1)$ symmetry. For such a circuit, the support of the non-conserved operators grows ballistically with time, while the conserved operators essentially perform a random walk. This leads to diffusive spreading of the conserved $U(1)$ charge. In turn, the weights of the conserved operator strings $a_{\mathcal{S}}$ take the form of a Gaussian

[TABLE]

where $D_{0}$ is the diffusion constant, and $d$ is the spatial dimension.

II.1 Dynamics with Line-like Subsystem Symmetries in Two Dimensions

We now consider quantum dynamics with subsystem $U(1)$ symmetries. In two spatial dimensions, the simplest such symmetry corresponds to overlapping $U(1)$ charges along intersecting lines. On the square lattice, we may consider dynamics that preserve the line charges

[TABLE]

which are the total $Z$ spin on each row and column of the square lattice, respectively.

Consider the case where a single charge is placed at the origin $(0,0)$ on the lattice, so that $C_{x}=\delta_{x,0}$ and $C_{y}=\delta_{y,0}$. This charge can only move while preserving the subsystem symmetries if we also allow for the creation of new charges. In particular, the 2D charge can move along a row or column if it emits a charge dipole in the direction perpendicular to its motion. We can implement such motion via an automaton unitary gate which preserves the subsystem symmetries.

In two dimensions on the square lattice, there are no two-site gates which preserve the line-like subsystem symmetries. The smallest such gate must act on four sites that lie at the corners of a rectangle. One such gate which preserves the subsystem symmetries is given by the “plaquette flip” unitary operator

[TABLE]

where $P$ is the projection operator

[TABLE]

with the indices $\{1,2,3,4\}$, as labeled as in Fig. 1. The projection operator $P$ projects into the subspace spanned by the two states shown in Fig. 1. The gate $U$ then flips between the two states in this subspace and acts trivially on all other states in the four site region. This automaton unitary operator generates operator entanglement; for example, we observe that

[TABLE]

We note that the most general unitary operator acting on an elementary plaquette that preserves the line symmetries, is one that performs an arbitrary rotation within the two-dimensional sub-space spanned by the two states shown in Fig. 1 – consisting of the states where each pair of neighboring spins is anti-aligned – plus any phase gate that is diagonal in the Pauli $Z$ basis.

We can therefore construct a circuit model which is composed only of these four-site unitary gates. In two dimensions, in addition to acting such gates on the fundamental plaquettes of the square lattice, we also include next-nearest-neighbor gates. In this case, we apply the same four-site unitary to the four corner sites of the rectangular plaquettes shown in Fig. 2. This is done to improve the ergodicity of the model. We then act the three gate types sequentially on all plaquette coverings of the square lattice in the pattern shown in Fig. 2. Therefore one periodic timestep of our circuit consists of 16 circuit ‘layers’.

II.1.1 Numerical Study

We start by numerically studying the 2D automaton circuit with linear subsystem symmetries, where we apply the 4-site unitary gates defined in Eq. (30) on all plaquettes of the square lattice, in the configuration shown in Fig. 2. We wish to determine the evolution of the correlation function

[TABLE]

We begin by numerically studying the auto-correlation function when $\boldsymbol{r}=0$. The results are shown in Fig. 3. While one may naively expect the same form for this correlation function as in the 1D case, since the conserved charges of the $Z(t)$ operator are constrained to move along the 1D rows and columns of the lattice, we find numerically that

[TABLE]

This is in sharp contrast to the 1D and 2D case with a global symmetry where we expect $G(t)\sim t^{-d/2}$. We note that it is only by going to very large system sizes and large circuit depths through the automaton unitary evolution, that we can resolve this logarithmic term in the correlation function.

Furthermore, we numerically determine the full space-time correlation function $G(\boldsymbol{r},t)$. Again, in systems with a global conserved charge, we expect this quantity to show Gaussian behavior. Our results, shown in Fig. 4, vary dramatically from this expectation. In particular, we see a distinct cusp shape in the spatial correlation function near the $t=0$ location of the charge.

This anomalous motion of the charge can be understood from the simple time evolution of the $Z$ operator given in Eq. 31. Notice that a charge cannot move freely in this model, but instead can only move by also emitting a dipole in the direction perpendicular to the motion. This means that motion of conserved charges in this model necessarily requires the creation of negative charges. Evidently, the creation of these dipoles results in the creation of a medium which slows the diffusion of the conserved charge. We note that this type of anomalous diffusion can be considered a purely quantum phenomenon, since these negative charges are associated with the sign of the operator wave function.

II.1.2 Analytical Study of the Charge Subdiffusion

We can gain some analytic understanding of these dynamics by observing that the overlap of $Z_{0}$ with the conserved line charge is $\mathrm{Tr}[Z_{0}C_{x}]=D\delta_{x,0}$ and $\mathrm{Tr}[Z_{0}C_{y}]=D\delta_{y,0}$. Therefore, for any dynamics that conserves these charges, the dynamical correlation function $G(\boldsymbol{r},t)$ satisfies the constraint that

[TABLE]

From this, we may construct the simplest evolution equation for $G_{\boldsymbol{r},t}$, after a single timestep of a discrete-time unitary evolution that preserves the line charges $C_{x}$, and $C_{y}$. Assuming that the dynamics are invariant under four-fold rotations on the square lattice, the simplest difference equation for the evolution of the dynamical correlations is given by

[TABLE]

where $\lambda$ is a free parameter of the evolution. We note that this is, in fact, an accurate description of the first timestep of the automaton unitary evolution considered in our numerical study with nearest-neighbor gates only (if we take $\lambda=1/8$); however, Eq. (35) neglects “backflow” effects – that non-conserved operators appearing in the Heisenberg evolution of $Z_{\boldsymbol{r}}(t)$ can develop a non-negligible weight on the operator $Z_{0}$ at sufficiently long times – which are present in the automaton dynamics that we numerically simulate.

Taking the continuum limit in time, and going to momentum space by defining $G(\boldsymbol{k},t)=\frac{1}{\sqrt{N}}\sum_{\boldsymbol{r}}\exp(i\boldsymbol{k}\cdot\boldsymbol{r})G(\boldsymbol{r},t)$, we find that Eq. (35) becomes

[TABLE]

where $f(\boldsymbol{k})=8\lambda\sin^{2}(k_{x}/2)\sin^{2}(k_{y}/2)$. We may analytically determine $G(|\boldsymbol{r}|=0,t)$ by taking the inverse Fourier transform, which may be performed to give $G(0,t)=\,_{\frac{1}{2}}F_{\frac{1}{2}}(1;1;-8\lambda t)$, where $\,{}_{a}F_{b}(p;q;z)$ is the generalized hypergeometric function. The long time expansion of this expression yields

[TABLE]

as we show in Appendix D; this agrees with our numerical study of the automaton dynamics. We note that the same result may be obtained by taking the continuum limit of Eq. (35) in both space and time, to obtain the partial differential equation

[TABLE]

By performing a Fourier transform, we may solve for $G(\boldsymbol{r},t)$ in terms of special functions, whose expansion at long times yields the same result as in Eq. (37).

We can also numerically simulate the finite difference equation Eq. (35). The solution $G(\boldsymbol{r},t)$ along the line $\boldsymbol{r}=(x,0)$ is plotted in Fig. 4. The insets show the behavior in the two limits $x<<\sqrt{t}$ and $x>>\sqrt{t}$, at fixed times $t_{0}$. We find that this is consistent with

[TABLE]

where $x_{0}$ is a constant which we can fit numerically.

Finally, it is interesting to note that, by virtue of conserving charge along lines, this circuit also preserves the dipole moment in both the $x$ and $y$ directions. In Ref.Pai et al. (2018) it was suggested that dipole conserving circuits in two dimensions should display localization of charge. In contrast, our results here establish that charge does spread, albeit subdiffusively, with the equation for charge spreading containing a term that involves the square of the Laplace operator. This suggests that the hydrodynamics for dipole conserving circuits proposed in Pai et al. (2018) cannot be complete, and furthermore suggests that the missing ingredient may be a Laplacian squared term.

Finally, we may consider the dynamics of a conserved charge in a three-dimensional system with intersecting line-like symmetries along three orthogonal directions, which are given by the total $Z$ spin along each of these directions. As we demonstrate in Appendix E, the three-dimensional generalization of Eq. (38) – given by $\partial_{t}G(\boldsymbol{r},t)=\lambda\,\partial_{x}^{2}\partial_{y}^{2}\partial_{z}^{2}G(\boldsymbol{r},t)$ – may be solved in terms of a special function, whose asymptotic form when $\lambda t\gg x^{2}y^{2}z^{2}$ yields the result that

[TABLE]

We have not confirmed this through numerical simulations of an automaton circuit, and we leave a detailed study of automaton dynamics with line-like symmetries in three spatial dimensions to future work.

II.2 Dynamics with Planar Symmetries in Three Dimensions

We can also consider the case of intersecting planar symmetries on the 3D cubic lattice. In this case, a $U(1)$ charge must be conserved on all $xy$, $xz$ and $yz$ planes. For spin-(1/2) degrees of freedom on the sites of a cubic lattice, we may define these charges as the total $Z$ spin on each plane,

[TABLE]

Automaton dynamics that respect these symmetries, along with the symmetries of the cubic lattice may be implemented using four-site unitary gates. In particular, the same $U_{\mathrm{flip}}$ unitary gate in Eq. (29) can be applied on the planes of the cube shown in Fig. 1, plus all symmetric rotations of these planes. Notice that these gates conserve the more strict line symmetries in the plane which they are applied, but only conserve the $U(1)$ charge on the plane for the perpendicular planes which they intersect.

We also study 3D circuit by applying the 4 sites gates sequentially on plaquettes lying in each of the planes shown in Fig. 1 (plus the symmetry allowed cubic rotations). In this case, each time step consists of 36 layers of gates (4 layers per plane times nine planes), as opposed to the 16 layers required for the 2D case.

We study the same quantities in the 3D model where charges are constrained to move only along two dimensional planes. The numerics in this case find that the autocorrelation function scales like

[TABLE]

That is, in this 3D case we find a clear violation of the diffusion law.

We also study the spatial distribution of $G(\boldsymbol{r},t)$ at fixed time for a 2D slice of the system, as shown in Fig. 6. We see that again the spatial distribution is clearly not Gaussian, and again shows a cusp near $x=0$. The charge density is greater along the lines which are shared between pairs of planes which have nonzero net charge.

Analytically we may understand this subddiffusive behavior by employing similar techniques as we employed in Sec. II.1. The conservation of planar $U(1)$ charges requires that the auto-correlation function $G_{\boldsymbol{r},t}$ is constrained, after summing over any $xy$, $yz$, or $xz$ plane, analogous to the constraint derived for line-like subsystem symmetries in Sec. II.1. The simplest, discrete-time evolution of this correlation function that is consistent with these constraints, and with the symmetries of the cubic lattice is

[TABLE]

where $\lambda$ is, again, a free parameter of the dynamics. Taking the continuum limit of this equation in space and time yields the partial differential equation for the coarse-grained correlation function

[TABLE]

which yields the result that $G(\boldsymbol{k},t)=\exp[-(\lambda t/4)\left[k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2}]\right]$ in momentum space. Performing the inverse Fourier transform, and re-scaling $k\rightarrow k^{\prime}=kt^{1/4}$ yields the desired result that $G(0,t)\overset{t\rightarrow\infty}{\sim}t^{-3/4}$.

As before, a more careful analysis of the difference equation (46), where we only take the continuum limit in time, while Fourier transforming in space, yields

[TABLE]

While we are unable to analytically perform the inverse Fourier transform to obtain the exact autocorrelation function from these expressions, we observe that $f(\boldsymbol{k})$ vanishes along the three lines $\boldsymbol{k}=(k_{x},0,0)$, $(0,k_{y},0)$, $(0,0,k_{z})$ that meet at the origin. As a result, the spatial continuum limit considered previously – which corresponds to an expansion of (48) near $\boldsymbol{k}\rightarrow 0$, where $f(\boldsymbol{k})$ vanishes most rapidly – is justified in understanding the long-time asymptotics of $G(0,t)$.

III Operator Spreading in Automaton Circuits

We now numerically study the behavior of the out-of-time-ordered correlation function (OTOC) – defined in Eq. (18) as $F(\boldsymbol{r},t)=({1}/{D})\text{Tr}\left[\mathcal{O}_{\boldsymbol{r}}(t)\mathcal{O}^{\prime}_{0}\mathcal{O}_{\boldsymbol{r}}(t)\mathcal{O}^{\prime}_{0}\right]$ – in order to quantify how local operators evolve under Heisenberg evolution by the automaton quantum circuits that we have considered. The primary purpose of our study of the OTOC is to argue that the “non-conserved part” of Heisenberg evolution of a local operator under automaton unitary dynamics – defined as the portion of the evolution of an operator that has no overlap with the conserved charges in the evolution – is quite generic, and resembles the evolution under the most general, symmetry-preserving unitary gates.

We study the following out-of-time-ordered correlation functions

[TABLE]

which we refer to as the “ZX” and “XX” OTOC’s for the remainder of this section.

As detailed in Sec. I.2, the ballistic growth of a non-conserved operator in a chaotic quantum system, and the equilibration of the local structure of these operators implies ($i$) that the OTOC has a ballistically-propagating front, and ($ii$) that $F(\boldsymbol{r},t)\rightarrow 0$ as $r/(v_{B}t)\rightarrow 0$. The saturation of the OTOC also coincides with the development of volume-law “operator entanglement” Hosur et al. (2016). Finally, we observe that ($iii$) the front of the OTOC is believed to broaden as $\sim t^{\alpha}$, where $\alpha$ is a dimension-dependent exponent; in one spatial dimension, the ends of an operator perform a biased random walk, leading to $\alpha=1/2$, while in higher dimensions, these exponents are related to the probability distributions for classical stochastic growth processes Nahum et al. (2018). We will verify these three features, which are characteristic of operator spreading in chaotic systems, in the out-of-time-ordered correlations of the automaton circuits considered previously.

In addition, we expect the ZX OTOC to have “tails” connecting the ballistically propagating front of the OTOC, to its value at position $\boldsymbol{r}$, due to the fact that the slow-moving conserved operators appearing in the Heisenberg evolution of $Z_{\boldsymbol{r}}$ can “emit” non-conserved operators, which then propagate ballistically Khemani et al. (2018). In contrast, no such tails are expected for the XX OTOC, which should saturate (exponentially) rapidly to its asymptotic value at times $t>r/v_{B}$. As a result, we consider a scaling form for the XX OTOC,

[TABLE]

which is motivated by similar scaling forms for the OTOC in a random unitary circuit Nahum et al. (2018). Here, $v_{B}(\Omega)$ and $\alpha(\Omega)$ are direction-dependent butterfly velocities and broadening exponents, respectively, and $\Omega$ denotes a $d$-dimensional angle.

We measure these out-of-time-ordered quantities using classical Monte Carlo techniques. For example, we may write the ZX OTOC as

[TABLE]

where $|n^{\prime}\rangle=X_{0}|n\rangle$. Therefore, measuring the ZX OTOC is as simple as keeping track of the classical evolution of the two states $|n\rangle$, $|n^{\prime}\rangle$. The XX OTOC may be measured in a similar manner, by evolving the two states $|n\rangle$ and $|m\rangle$, applying a flip operator at site $\boldsymbol{r}$ to each state, and evolving the resulting states backwards in time. The average overlap between the resulting states gives the XX OTOC. In this case, every space-time measurement at point $(\boldsymbol{r},t)$ requires an independent simulation, making the XX OTOC measurements somewhat slower to simulate in practice.

III.1 Numerical Study of the OTOC

We now present the results of our measurements of the out-of-time-ordered correlator.

In Fig. 7, we see the growth of the XX OTOC along the line $\boldsymbol{r}=(x,0)$, for our 2D automaton circuit. We observe a ballistically spreading light cone, and that the OTOC is very nearly zero inside the light cone and one outside. From similar measurements of the XX OTOC along other directions, we observe a nearly isotropic butterfly velocity. We also see in Fig. 7c, that for a fixed position, the OTOC appears to decay exponentially in time, within the times that we are able to perform numerical simulations. This is consistent with the expectation that the OTOC should decay to zero at sufficiently long times. Finally, shown in Fig. 8 is a scaling collapse of the XX OTOC, to the functional form in Eq. 52; the optimal scaling collapse is to the butterfly velocity $v_{B}$ and front broadening exponent $\alpha$, which are given by

[TABLE]

Here, $v_{\mathrm{max}}$ is the lightcone velocity along the line $\boldsymbol{r}=(x,0)$ for our automaton circuit. This should be compared with $\alpha=1/3$, which is the observed broadening for a two-dimensional, random unitary circuit, for which the growth of the OTOC is related to the stochastic growth of a two-dimensional cluster. The fluctuations in the growing edge of the cluster are believed to be the same as the fluctuations in the height of a stochastically growing, one-dimensional interface, which is known to grow in time as $t^{1/3}$ Nahum et al. (2018).

In Fig. 9, we look at the ZX OTOC, and we observe a “tail” behind the ballistically-propagating front of the OTOC, due to the fact that the slow-moving operators that have an overlap on the subsystem symmetry charges continue to “emit” non-conserved operators that propagate ballistically. We fit the tail of the OTOC data in two and three dimensions to the functional form

[TABLE]

and find that $\beta=0.50(1)$ for our 2D circuit with line-like subsystem symmetries, while $\beta=1.02(7)$ for the 3D circuit with planar symmetries. Trendlines with this power-law behavior are shown in Fig. 9.

For the 3D case, we are unable to calculate the XX OTOC to long enough times that we are able to perform a satisfactory scaling collapse to quantify the broadening of the front of a non-conserved operator. Instead, we restrict our attention to the ZX OTOC along the line $\boldsymbol{r}=(x,0,0)$. We obtain a convincing scaling collapse of the ballistically propagating front of this OTOC – with a front broadening exponent $\alpha=0.220(5)$ – in Fig. 8b; this is very close to the predicted value of the broadening exponent obtained from studies of operator spreading in Haar-random unitary circuits Nahum et al. (2018), which predict $\alpha=0.24$. We also observe that the tail follows a power law form as in Eq. (55), except with exponent $\beta=1.02(7)$, as shown in Fig. 9.

III.2 Recurrence times

Finally, we study the recurrence time for product states in the Pauli $Z$ basis, evolving under the automaton dynamics; these product states only evolve to other product states in the same basis. The recurrence time $t_{\mathrm{prod}}(n)$ is then the circuit depth for which the state returns to its initial value, i.e. $U^{t_{\mathrm{prod}}(n)}|n\rangle=|n\rangle$. In Fig. 10, we show the distribution of these recurrence times, taken over the full set of states $|n\rangle$ in the computational basis, for different system sizes. From the peaks of these distributions, we obtain the typical recurrence time as a function of system size, as shown in Fig. 10b). We see that $t_{\mathrm{prod}}$ grows linearly in the Hilbert space dimension $D$ of the system, as predicted for a random automaton dynamics. This is considerably slower than for fully chaotic systems, where the recurrence time would be doubly exponential in the number of sites, but considerably faster than translation invariant, Floquet Clifford circuits, where one can prove that $t_{rec}=\mathcal{O}(\log D)$ von Keyserlingk et al. (2018).

IV Conclusions

We have discussed a new class of circuits in two and three dimensions, which produce a quantum dynamics that is classically simulable. In particular, the evolution of certain correlation functions under the action of these circuits is classically simulable even though the circuits generate operator entanglement, since these dynamics avoid generating state entanglement when acting on product states in the computational basis. We have introduced a set of circuits that respect subsystem symmetries, conserving a $U(1)$ charge along either lines in two space dimensions, or along planes in three space dimensions. These constraints lead to dynamics that lies in a new universality class, with charge spreading subdiffusively both with respect to the physical dimension and with respect to the dimension of the submanifold in which the symmetries act. In particular, for a two dimensional system with line like symmetries, charge spreads as $\log(t)/\sqrt{t}$, whereas for a three dimensional system with planar symmetries, charge spreads as $1/t^{3/4}$. This behavior is captured by a finite difference equation, the continuum limit of which yields an ‘anomalous subdiffusion’ equation in which the Laplace operator is replaced by something akin to a square of the Laplace operator.

While the dynamics we have considered are structured (and have to be, in order to be classically simulable), we have nevertheless shown that they do exhibit several features commonly associated with chaos. In particular, OTOC’s involving non-conserved operators exhibit a ballistic lightcone with a front that broadens as a dimension-dependent power-law. Furthermore these OTOC’s appear to saturate to zero at long times, which suggests that the bulk structure of a non-conserved operator equilibrates in a manner that is consistent iwth quantum chaos. OTOCs between a conserved and non-conserved operator exhibit in addition a ‘universal tail’ behind the propagating front, as has also been observed in a random circuit with a global $U(1)$ symmetryKhemani et al. (2018). As such, we conjecture that our results may hold in more generic quantum dynamics with the same symmetries.

This work opens up a new direction for the exploration of quantum dynamics, rendering accessible questions regarding chaos and operator spreading in higher dimensions in a class of circuits that are much less structured and more generic than any previously known classically-simulable higher dimensional circuit of which we are aware. It also opens the door to a new class of dynamical phenomena produced by subsystem conservation laws.

A number of extensions of this work immediately present themselves for consideration. One possibility is to move away from the square/cubic lattices considered herein, and to consider instead subsystem symmetries on more general crystal lattices. Symmetries involving fractal subsystems, as in Type II fracton models Vijay et al. (2016); Haah (2011) are important to consider. Either of these extensions holds the promise of qualitatively new universality classes for dynamical behavior. Extensions of this work to Hamiltonian systems (rather than quantum circuits) are of interest. We expect analogous behavior in the Hamiltonian setting. Other potential extensions include considering subsystem symmetry groups other than $U(1)$. A final possibility is to move away from gates with strictly local support to e.g. quasilocal gates with exponential tails in real space. We leave such questions for future work.

Acknowledgements.

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-17-1-0183 (JI and RMN). SV is supported by the Harvard Society of Fellows. SV and RMN also acknowledge the hospitality of the KITP, where this work was initiated, during a visit to the program Dynamics of Quantum Information. The KITP is supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

Appendix A Operator Spreading under Automaton Dynamics

It is convenient to study the growth of the complexity of an operator evolving under the automaton dynamics as follows. We consider $N$ spins with Hilbert space dimension $D=2^{N}$. Let $n\in\{0,1,\ldots,D-1\}$ be an integer, and let $\boldsymbol{n}=(n_{1},n_{2},\ldots,n_{N})$ be an $N$-component vector that corresponds to the binary representation of integer $n$. We further define $\mathcal{S}_{n}$ to be a product of $Z$ operators that corresponds to the binary representation $n$. Specifically,

[TABLE]

For example, $\mathcal{S}_{0}$ is simply the identity operator, while $\mathcal{S}_{7}=Z_{1}Z_{2}Z_{3}$ since seven is represented as $111000\cdots$ in binary. We note that $\frac{1}{D}\mathrm{Tr}(\mathcal{S}_{n}\mathcal{S}_{m})=\delta_{nm}$ where $D$ is the total Hilbert space dimension.

Under the automaton dynamics, these operators evolve as

[TABLE]

In contrast, for a more general unitary dynamics, each of the operators $\{\mathcal{S}_{n}\}$ will evolve into sums of arbitrary products of Pauli operators in the system.

It is now convenient to define the following wavefunction in a spin-(1/2) system with the same total Hilbert space dimension

[TABLE]

where $|m\rangle$ is the representation of integer $m$ as the state of a spin-(1/2) system in the Pauli-$Z$ basis, with an up (down) spin appearing whenever a 0 (1) occurs in the binary representation of $m$; as a concrete example, the state

[TABLE]

To summarize, Eq. (57) describes the Heisenberg evolution of operators under the automaton dynamics, while Eq. (58) describes a fictitious dynamics of a many-body wavefunction in the same Hilbert space. While the entanglement of a wavefunction in the $Z$ basis does not grow under the automaton dynamics, the states of the fictitious spin system do become highly entangled as they evolve, reflecting the fact that typical operators become quite entangled under the automaton evolution.

We now identify a few key points about the fictitious dynamics of the state $|n(t)\rangle$.

1. (I)

The evolution of this wavefunction is unitary, so that

[TABLE]

where $W(t)W(t)^{\dagger}=W(t)^{\dagger}W(t)=1$. 2. (II)

$W(t)$ is exactly

[TABLE]

where $H$ is the single-qubit Hadamard gate that acts on Pauli matrices as $HXH=Z$, $HZH=X$. That is $W(t)$ is the same automaton unitary operator acting in a rotated basis.

Proof of (I): Since $\langle{m|W(t)|n}\rangle=a_{nm}(t)$, we observe that

[TABLE]

so that $W(t)$ is indeed unitary.

Proof of (II): The automaton unitary $U(t)$ acts as a permutation $\pi\in S_{D}$ on the $D$-dimensional set of product states in the Pauli $Z$ basis. We define the projection operator

[TABLE]

The automaton unitary acts as

[TABLE]

This relation places a restriction on the coefficients $a_{mn}(t)$. We first rewrite the projection operator as

[TABLE]

where $\boldsymbol{n}=(n_{1},\ldots,n_{N})$ is the $N$-component binary vector that corresponds to the binary representation of $n$, and where $\boldsymbol{n}\cdot\boldsymbol{m}$ is the dot product of the two vectors. The Heisenberg evolution of this projector may be written as

[TABLE]

so that

[TABLE]

Using the fact that $(1/D)\sum_{m}(-1)^{\boldsymbol{m}\cdot(\boldsymbol{n}+\boldsymbol{k})}=\delta_{n,k}$, we have found that

[TABLE]

This is equivalent to the statement that we wish to prove. We observe that if $|m\rangle$ is a product state in the Pauli $Z$ basis, where the spin configuration corresponds to the binary representation of the integer $m$ as in Eq. 59, then applying a Hadamard gate to every spin yields

[TABLE]

Now, we observe that

[TABLE]

We conclude that

[TABLE]

Appendix B Operator Entanglement Entropy

We calculate the purity of the wavefunction in Eq. (58), which quantifies the complexity of the Pauli $Z$ operator spreading from the automaton dynamics. Consider a bi-partitioning of the state $|n(t)\rangle$ in to an $A$ and $B$ subsystem, with Hilbert space dimensions $D_{A}$ and $D_{B}=D/D_{A}$, respectively. Furthermore, let $\{|m_{A,B}\rangle\}$ denote a complete set of states in $A$ and $B$, respectively; as before, we let $\boldsymbol{m}_{A,B}$ denote the binary vector representation of these states, so that the dot product $\boldsymbol{n}\cdot\boldsymbol{m}=\boldsymbol{n}_{A}\cdot\boldsymbol{m}_{A}+\boldsymbol{n}_{B}\cdot\boldsymbol{m}_{B}$.

Using this notation, we may write down the reduced density matrix as

[TABLE]

where the matrix elements of the reduced density matrix are given by

[TABLE]

Substituting the expression for $a_{nm}(t)$ and performed the sum over $m_{B}$, $k_{B}$ in the above expression yields

[TABLE]

Here, $\pi^{-1}$ is the inverse of the permutation $\pi\in S_{D}$, while $\boldsymbol{\pi^{-1}}(q)$ is the vector representation of the element $\pi^{-1}(q)$, as before. Using this expression, we find that the purity (exponential of the second Rényi entropy is given by

[TABLE]

We now calculate the purity, averaged over all permutations $\pi\in S_{D}$. In this case, we observe that for non-zero $\boldsymbol{n}$, the quantity $\boldsymbol{n}\cdot\boldsymbol{\pi^{-1}}(q)$ can be an even or odd integer with equal probability, for a randomly chosen permutation $\pi$. Therefore, for non-zero $\boldsymbol{n}$, we observe that

[TABLE]

Similarly, the average of the expression in Eq. (B) only contributes if pairs of elements are equal. Therefore for non-zero $\boldsymbol{n}$,

[TABLE]

We may use this expression to compute the purity averaged over all choices of permutations, which we denote as $\overline{\,\mathrm{Tr}[\rho_{A}^{2}]\,}$. Substituting Eq. (B) into Eq. (B) yields the result that

[TABLE]

Appendix C Recurrence Time Distribution

Consider the permutation group $S_{D}$. If $\pi\in S_{D}$ has a length-$\ell$ cycle that includes a particular element (say the element $1$), then we may write $\pi$ as the product of two commuting permutations, i.e.

[TABLE]

with $\sigma\in S_{D-\ell}$ acting on the $D-\ell$ elements that do not appear in the length-$\ell$ cycle. For a fixed set of elements $k_{2},\ldots,k_{\ell}$, the number of such permutations $\pi\in S_{D}$ is precisely $(\ell-1)!(D-\ell)!$. Furthermore, the number of ways to choose these elements out of $D$ elements is \left(\begin{array}[]{c}D-1\\ \ell-1\end{array}\right). Therefore, the number of permutations $\pi\in S_{D}$ with a length-$\ell$ cycle that contains the element $1$ is precisely

[TABLE]

Alternatively, the probability that a random permutation in $S_{D}$ has a length-$\ell$ cycle containing a particular, fixed element $m$ is

[TABLE]

For a permutation $\pi\in S_{D}$, we define $\ell_{n}(\pi)$ to be the length of the cycle that element $n$ is in. Using this, we may write that the average number of cycles $\langle{N}\rangle$ in a permutation is

[TABLE]

Similarly, the average return time for a random element $\langle t\rangle$ is given by

[TABLE]

Therefore, we conclude that

[TABLE]

Appendix D Dynamics with 1D Subsystem Symmetries in Two Dimensions

We assume that the correlation function $G_{\boldsymbol{r},t}\equiv\mathrm{Tr}\left[Z_{\boldsymbol{r}}(t)Z_{0}(t)\right]/D$, where $D=2^{N}$ is the Hilbert space dimension of the $N$-site system evolves according to the difference equation,

[TABLE]

Taking the continuum limit in time, so that $G_{\boldsymbol{r},t+1}\approx G(\boldsymbol{r},t)+dG(\boldsymbol{r},t)/dt$, and going to momentum space,

[TABLE]

we find that

[TABLE]

where

[TABLE]

which then gives the solution

[TABLE]

given the initial condition $G(\boldsymbol{r},0)=\delta_{\boldsymbol{r},0}$.

We now evaluate $G(0,t)$: in the thermodynamic limit, we make the replacement $N^{-1}\sum_{\boldsymbol{k}}\rightarrow\int d^{2}\boldsymbol{k}/(2\pi)^{2}$, so that

[TABLE]

Performing the integral over $k_{x}$ yields

[TABLE]

where $I_{0}(z)$ is the modified Bessel function of the first kind. Finally, the integral over $k_{y}$ yields

[TABLE]

where $\,{}_{a}F_{b}(p;q;z)$ is the generalized hypergeometric function. The asymptotic expansion of $G(0,t)$ at long times is then

[TABLE]

More precisely, the asymptotic expansion is

[TABLE]

Resolving the $\ln(t)/\sqrt{t}$ behavior would require going to very long times (note the logarithm becomes approximately twice as large as the constant correction only when $t\sim 10^{4}$).

We now determine the late-time scaling of the correlation function $G(\boldsymbol{r},t)$

[TABLE]

along the line $\boldsymbol{r}=(x,0)$. Hereafter, we refer to this correlation function as $G(x,t)$. Performing the integral over $k_{y}$, we obtain

[TABLE]

, where

[TABLE]

At long times $t\gg 1$, the quantity $(t/2)\sin^{2}(k/2)$ is large over the interval $k\gtrsim 1/\sqrt{t}$. The scaling of the correlation function in this regime may be approximated by replacing $g(k,t)$ by its asymptotic form at long times, for small $k$,

[TABLE]

This integral is logarithmically divergent at small $k$. When $t\gg x^{2}\gg 1$, we may evaluate the integral perturbatively in $x/\sqrt{t}$ to find that

[TABLE]

So that the correlation function “peaks” when $t\sim x^{2}$, and the maximum peak height scales as $\sim 1/x$.

Appendix E Line-Like Symmetries in 3D

For a three-dimensional system with intersecting, line-like symmetries, we consider the natural generalization of Eq. (38) in the space-time continuum limit:

[TABLE]

where $\lambda$ is a free parameter of the evolution, so that

[TABLE]

where $\mathcal{G}^{m,n}_{p,q}$ is the Meijer G-function, and the vectors $\boldsymbol{v}$ and $\boldsymbol{w}$ are defined as

[TABLE]

The asymptotic form $t\gg x^{2}y^{2}z^{2}$ of this integral yields the leading, long-time behavior

[TABLE]