Quantum Entanglement Growth under Random Unitary Dynamics

Quantum Entanglement Growth under Random Unitary Dynamics PHYSICAL REVIEW X 7, 031016 (2017) 1,3 1 1,2 1 Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Kavli Institute for Theoretical Physics, University of California Santa Barbara, Santa Barbara, California 93106, USA Theoretical Physics, Oxford University, 1 Keble Road, Oxford OX1 3NP, United Kingdom (Received 1 September 2016; revised manuscript received 16 April 2017; published 24 July 2017) Characterizing how entanglement grows with time in a many-body system, for example, after a quantum quench, is a key problem in nonequilibrium quantum physics. We study this problem for the case of random unitary dynamics, representing either Hamiltonian evolution with time-dependent noise or evolution by a random quantum circuit. Our results reveal a universal structure behind noisy entanglement growth, and also provide simple new heuristics for the “entanglement tsunami” in Hamiltonian systems without noise. In 1D, we show that noise causes the entanglement entropy across a cut to grow according to the celebrated Kardar- Parisi-Zhang (KPZ) equation. The mean entanglement grows linearly in time, while fluctuations grow like 1=3 2=3 ðtimeÞ and are spatially correlated over a distance ∝ ðtimeÞ . We derive KPZ universal behavior in three complementary ways, by mapping random entanglement growth to (i) a stochastic model of a growing surface, (ii) a “minimal cut” picture, reminiscent of the Ryu-Takayanagi formula in holography, and (iii) a hydrodynamic problem involving the dynamical spreading of operators. We demonstrate KPZ universality in 1D numerically using simulations of random unitary circuits. Importantly, the leading-order time dependence of the entropy is deterministic even in the presence of noise, allowing us to propose a simple coarse grained minimal cut picture for the entanglement growth of generic Hamiltonians, even without noise, in arbitrary dimensionality. We clarify the meaning of the “velocity” of entanglement growth in the 1D entanglement tsunami. We show that in higher dimensions, noisy entanglement evolution maps to the well- studied problem of pinning of a membrane or domain wall by disorder. DOI: 10.1103/PhysRevX.7.031016 Subject Areas: Condensed Matter Physics, Quantum Information, Statistical Physics I. INTRODUCTION goes on. This irreversible growth of entanglement— quantified by the growth of the von Neumman entropy— The language of quantum entanglement ties together is important for several reasons. It is an essential part condensed matter physics, quantum information, and high- of thermalization, and as a result has been addressed in energy theory. The von Neumann entanglement entropy is diverse contexts ranging from conformal field theory known to encode universal properties of quantum ground [1–4] and holography [5–12] to integrable [13–19], non- states and has led to new perpectives on the AdS=CFT integrable [20–23], and strongly disordered spin chains correspondence. But the dynamics of the entanglement are [24–30]. Entanglement growth is also of practical impor- far less understood. The entanglement entropy is a highly tance as the crucial obstacle to simulating quantum nonlocal quantity, with very different dynamics to energy dynamics numerically, for example, using matrix product or charge or other local densities. Traditional many-body states or the density matrix renormalization group [31]. tools therefore do not provide much intuition about how The entanglement entropy, and even its time dependence, entanglement spreads with time, for example, after a is also beginning to be experimentally measurable in quantum quench (a sudden change to the Hamiltonian). cold atom systems [32–34]. In a very different context, We need to develop simple heuristic pictures, and simple black holes have motivated studies of how fast quantum long-wavelength descriptions, for entanglement dynamics. systems can scramble information by dynamically gen- If a many-body system is initialized in a state with erating entanglement [35–38]. Simple quantum circuits— low entanglement, the dynamics will typically generate quantum evolutions in discrete time—serve as useful toy entanglement between increasingly distant regions as time models for entanglement growth and scrambling [39–43]. For integrable 1D systems and rational CFTs, there is an appealing heuristic picture for entanglement growth fol- Published by the American Physical Society under the terms of lowing a quench in terms of spreading quasiparticles [1,3]. the Creative Commons Attribution 3.0 License. Further distri- However, this picture does not apply to general interacting bution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. systems [3,4,10,44]. 2160-3308=17=7(3)=031016(30) 031016-1 Published by the American Physical Society NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) In this paper, we propose new heuristic pictures for KPZ universality class entanglement growth in generic, nonintegrable systems, both in 1D and in higher dimensions. We arrive at these pictures by studying “minimally structured” models for quantum Classical Directed polymer in Stochastic dynamics: dynamics that are spatially local, and unitary, but surface growth random medium particle dynamics random both in time and space (“noisy”). Concretely, we focus on quantum circuit dynamics with randomly chosen quantum gates. Entanglement growth in these systems exhibits a remarkable universal structure in its own right, related to paradigmatic problems in classical statistical Growth of Hydrodynamics mechanics. But in addition, random circuits provide a entanglement through circuit of operator spreading theoretical laboratory that allows us to derive scaling pictures FIG. 1. The KPZ “triumvirate” is made up of three very for entanglement growth and the so-called “entanglement different problems in classical statistical mechanics which all tsunami” [8], which, we conjecture, generalize to quenches map to the KPZ universality class. As we discuss, each of them in many-body systems without noise. For example, we can be usefully related to entanglement in 1 þ 1D. propose a simple “minimal membrane” picture that can be used to derive scaling forms for the growth of the entangle- ment. We also argue that generically there is a well-defined summarized in Fig. 1. We show that entanglement growth “entanglement speed” v , but this is generically smaller than E can usefully be related to all three of the classical problems in the “butterfly speed” v governing operator growth, and we B in Fig. 1. give a physical explanation for this phenomenon. In the quantum setting, the directed polymer is related to We show that noisy entanglement growth allows a long- the “minimal cut,” a curve in space-time that bisects the wavelength description with an emergent universal structure. unitary circuit representing the time evolution. This picture Physically, the class of noisy dynamics includes closed, is more general than the surface growth picture, as it allows many-body systems whose Hamiltonian HðtÞ contains one to consider the entropy for any bipartition of the system. fluctuating noise terms, and also quantum circuits in which It also allows us to generalize from 1D to higher dimensions. qubits are acted on by randomly chosen unitary gates. In this The picture is reminiscent of the Ryu-Takayanagi prescrip- setting, we pin down both the leading-order deterministic tion for calculating the entanglement entropy of conformal behavior of the entanglement and the subleading fluctuations field theories in the AdS=CFT correspondence, which associated with noise. We argue that fluctuations and spatial makes use of a minimal surface in the bulk space [50], correlations in the entanglement entropy are characterized and analogous results for certain tensor network states by universal scaling exponents, expected to be independent [51–53]. Here, however, the cut lives in space-time rather of the details of the microscopic model. than in space, and in a noisy system its shape is random For noisy systems in one spatial dimension, we argue that rather than deterministic. (For a different use of the idea of a the critical exponents for entanglement growth are those of minimal cut in space-time, see Ref. [10].) In d þ 1 space- the Kardar-Parisi-Zhang (KPZ) equation, originally intro- time dimensions the minimal cut becomes a d-dimensional duced to describe the stochastic growth of a surface with time membrane pinned by disorder. This picture allows us to t [45]. In the simplest setting, we find that the height of this obtain approximate critical exponents for noisy entangle- ment growth in any number of dimensions. surface at a point x in space is simply the von Neumann This picture also leads us to a conjecture for entangle- entanglement entropy Sðx; tÞ for a bipartition which splits the system in two at x. The average entanglement grows linearly ment growth in systems without noise, both in 1D and in time, while fluctuations are characterized by nontrivial higher dimensions, as we discuss below. According to this conjecture, the calculation of the entanglement in higher exponents. We support this identification with analytical dimensions reduces to a deterministic elastic problem for arguments and numerical results for discrete time quantum the minimal membrane in space-time. In 1D, it results in evolution (unitary circuits). particularly simple universal scaling functions, which agree A remarkable feature of the KPZ universality class is that with scaling forms in holographic 1 þ 1D CFTs [8,10,44], it also embraces two classical problems that at first sight are very different from surface growth [45,46]. These connec- and which we suggest are universal for generic, non- tions lead us to powerful heuristic pictures for entanglement integrable, translationally invariant 1D systems. growth, in both 1D and higher dimensions. The KPZ The third member of the triumvirate in Fig. 1 is a noisy universality class embraces the statistical mechanics of a hydrodynamic equation describing the diffusion of inter- directed polymer in a disordered potential landscape [47] and acting (classical) particles in 1D. We show that this can be 1D hydrodynamics with noise (the noisy Burgers equation related to the spreading of quantum operators under the [48]). These problems, together with surface growth, are unitary evolution, giving a detailed treatment of the special sometimes known as the “KPZ triumvirate” [49]. They are case of stabilizer circuits. Note that while the minimal cut 031016-2 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) picture generalizes to higher dimensions, the KPZ and hydrodynamic pictures are special to 1D. A. Organization of the paper We propose that noisy dynamics are a useful toy model for quantum quenches in generic (nonintegrable, noncon- FIG. 2. Spin chain with open boundary conditions. SðxÞ formally–invariant) systems, even without noise. The logic denotes the entanglement entropy (von Neumann or Rényi of our approach is to pin down the universal behavior of depending on context) between the part of the chain to the left noisy systems (Secs. II–VI), to establish simple heuristics of bond x, indicated by the box, and the part to the right of bond x. capturing this behavior (Secs. III and IV), and then to extend these heuristic pictures to dynamics without noise by splitting the chain into two halves at x and tracing out (Secs. V and VIII). the left-hand side (Fig. 2). The nth Rényi entropy for a cut The detailed physics of the entanglement fluctuations at x is defined as (including KPZ exponents) certainly relies on noise. However, the coarser features of the dynamics are in fact S ðxÞ¼ log ðTrρ Þ: ð1Þ deterministic. These include the leading-order time depend- n x 1 − n ence of the entanglement entropy and mutual information when the length and time scales are large. We conjecture Logarithms are taken base q. In the limit n → 1, the Rényi that this leading-order behavior, as captured by the directed entropy becomes the von Neumann entropy: polymer and hydrodynamic pictures, carries over to Hamiltonian dynamics without noise. On the basis of this S ðxÞ¼ −Trρ log ρ : ð2Þ vN x x we address (Sec. V) features of entanglement growth that have previously been unclear. We argue that in generic 1D A basic constraint on the von Neumann entropy is that systems the entanglement growth rate can be interpreted as neighboring bonds can differ by at most one (this follows a well-defined speed v , but that this speed is generically from subadditivity of the von Neumann entropy): smaller than another characteristic speed, which is the speed v at which quantum operators spread out under the jS ðx þ 1Þ − S ðxÞj ≤ 1: ð3Þ vN vN dynamics (the butterfly speed). Section V also addresses universal scaling forms for the entanglement entropy in 1D. In this section, we focus on the growth of the bipartite In Sec. VIII, we discuss the geometry dependence of the entropies Sðx; tÞ with time, starting from a state with low dynamical entanglement in higher-dimensional systems: entanglement. [Here Sðx; tÞ, without a subscript, can denote we argue that there is again a scaling picture in terms of a any of the Rényi entropies with n> 0.] For simplicity, we minimal surface, but that more nonuniversal parameters take the initial state to be a product state, but we expect the enter into the time dependence than in 1+1D. same long-time behavior for any initial state with area-law entanglement. (The setup with area-law entanglement in the initial state is analogous to a quantum quench that starts II. SURFACE GROWTH IN 1D in the ground state of a noncritical Hamiltonian. We briefly We begin by studying entanglement growth under consider initial states with non-area-law entanglement in random unitary dynamics in one dimension. After sum- Sec. IX.) We argue that for noisy unitary dynamics, the marizing the KPZ universal behavior, we derive this universal properties of Sðx; tÞ are those of the Kardar- behavior analytically in a solvable model, using a mapping Parisi-Zhang equation: to a classical surface growth problem. In the following sections we provide alternative derivations of this universal ∂S λ 2 2 ¼ ν∂ S − ð∂ SÞ þ ηðx; tÞþ c: ð4Þ x x behavior by relating the minimal cut bound on the ∂t 2 entanglement to the classical problem of a directed polymer in a random environment, and by relating the spreading of This equation was introduced to describe the stochastic quantum operators to a 1D hydrodynamics problem. growth of a 1D surface or interface with height profile SðxÞ Consider a chain of quantum spins with local Hilbert [45]. It captures an important universality class which has space dimension q (for example, spin-1=2’s with q ¼ 2). found a wealth of applications in classical nonequilibrium We take open boundary conditions and label the bonds of physics, and its scaling properties have been verified the lattice by x ¼ 1; …;L. We consider only unitary in high-precision experiments [54,55]. The constant c in dynamics, so the full density matrix ρ ¼jΨihΨj represents Eq. (4) contributes to the positive average growth rate, a pure state. For now we consider the entanglement across a while ηðx; tÞ is noise which is uncorrelated in space and single cut at position x; we generalize to other geometries time. The ν term describes diffusive smoothing of sharp later in the paper. The reduced density matrix ρ is defined features. The nonlinear term, with coefficient λ, describes 031016-3 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) how the average growth rate depends on the slope; the t ≃ ; ð10Þ saturation negative sign is natural here, as we discuss in Sec. II A 3 2v [and implies that B in Eq. (6) below is positive]. with bonds closer to the boundary saturating sooner (Secs. III KPZ scaling is characterized by an exponent β governing and V). Saturation is also captured in the surface growth the size of fluctuations in the interface height, an exponent description, once we note that there are Dirichlet boundary α governing spatial correlations, and a dynamical exponent conditions on the entropy: Sð0;tÞ¼ SðL þ 1;tÞ¼ 0. z determining the rate of growth of the correlation length Note that the scaling described in Eqs. (6) and (8) implies (z ¼ α=β by a scaling relation). These are known exactly the existence of two distinct diverging length scales during [45]: entanglement growth. The fact that hSðx; tÞi is of order t implies that spins are entangled over distances of order t. β ¼ 1=3; α ¼ 1=2;z ¼ 3=2: ð5Þ In fact, we show in Sec. V that v t is a sharply defined length scale. But prior to saturation, the relevant length In our context the height of the surface is the bipartite scale for spatial variations in Sðx; tÞ is parametrically entanglement Sðx; tÞ. This is a random quantity that 2=3 smaller than v t; namely, ξðtÞ ∼ t . depends on the realization of the noise in the quantum Before deriving KPZ for entanglement, let us briefly dynamics. The mean height (entanglement) grows linearly consider the status of this equation. At first sight, we might in time, with a universal subleading correction: try to justify this description of Sðx; tÞ simply on grounds of symmetry and coarse graining. If we were describing hðx; tÞ ≡ hSðx; tÞi ¼ v t þ Bt : ð6Þ classical surface growth, we would appeal to translational symmetry in the growth direction (S → S þ const) in order Angle brackets denote an average over noise. The linearity to restrict the allowed terms, and would note that the right- of the leading t dependence is expected from rigorous hand side includes the lowest-order terms in ∂ and ∂ S. x x bounds for various 1 þ 1D random circuits [39–41]. Linear But for entanglement we cannot rely on this simple growth is also generic for quenches in translationally reasoning. First, the transformation S → S þ const is not invariant 1D systems [3,21]. The fluctuations grow as a symmetry (or even a well-defined transformation) of the quantum system. More importantly, it is not clear a priori 2 1=2 β wðx; tÞ ≡ ⟪Sðx; tÞ ⟫ ¼ Ct : ð7Þ that we can write a stochastic differential equation for Sðx; tÞ alone, since the full quantum state contains vastly We refer to w as the width of the surface. The ratio C=B is more information than Sðx; tÞ. Despite these differences universal (the constants v and B are not). The KPZ from simple surface growth, we show below that the fluctuations are non-Gaussian: remarkably, their universal above equation does capture the universal aspects of the probability distribution has been determined analytically entanglement dynamics. [56–66]. In the next section, we exhibit a solvable quantum model The correlation length governing spatial correlations in that maps to a classical surface growth problem that is the fluctuations grows with time as manifestly in the KPZ universality class. Then in the two following sections, we give heuristic arguments for more 1=z ξðtÞ ∼ t ; ð8Þ general systems by making connections with the other members of the KPZ triumvirate. Together with the results for the solvable model, these arguments suggest that KPZ and the equal time correlation function has the scaling form exponents should be generic for entanglement growth in 2 1=2 α any quantum system whose dynamics involves time- GðrÞ ≡ h½Sðx; tÞ − Sðx þ r; tÞ i ¼ r g(r=ξðtÞ): ð9Þ dependent randomness. In Sec. VI, we perform numerical checks on KPZ universality for quantum dynamics in On length scales 1 ≪ r ≪ ξðtÞ, the surface profile SðxÞ discrete time. resembles the trace of a 1D random walk: this is consistent with the exponent α ¼ 1=2. On scales r ≫ ξðtÞ, the fluctua- A. Solvable 1D model tions in Sðx; tÞ and Sðx þ r; tÞ are essentially uncorrelated. At short times the entanglement growth is affected by We now focus on a specific quantum circuit model for the initial conditions, while on very long time scales, of the dynamics of a spin chain with strong noise. We take random order of the system size, the entanglement saturates. unitaries to act on pairs of adjacent spins (i.e., on bonds) at Equations (6)–(9) apply prior to this saturation. In a finite random locations and at random times, as illustrated in Fig. 3. system the asymptotic hSðxÞi profile is that of a pyramid, For simplicity, we discretize time and apply one unitary per with a maximum at height x ¼ L=2, whose height is L=2, time step. (Dynamics in continuous time, with unitaries minus an Oð1Þ correction [67,68]. This profile is reached at applied to the links at a fixed rate in a Poissonian fashion, are a time equivalent.) We choose the initial state to be a product state, 031016-4 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) U For the random dynamics we are considering, the dynami- time step cal rule in Eq. (13) leads to a simple but nontrivial stochastic process. Before discussing its properties, we use Eq. (13) as a FIG. 3. Dynamical update in the solvable model. Application of starting point to show that, in the limit of large Hilbert space a random unitary U to a randomly chosen pair of adjacent spins. dimension, the von Neumann entropy (and in fact all the higher Rényi entropies) obeys the same dynamical rule. This requires an explicit calculation in the limit q → ∞.The von with S ðxÞ¼ 0 for all n and x. We choose the unitaries from Neumann entropy is of more interest than S , since the latter the uniform (Haar) probability distribution on the unitary behaves pathologically in many circumstances. [This is group for a pair of spins, Uðq Þ. This model is solvable in the because it simply counts up all the (nonzero) eigenvalues limit of large local Hilbert space dimension q. in the spectrum of ρ , regardless of how small they are. For example, Hamiltonian dynamics in continuous time—as 1. Dynamics of Hartley entropy opposed to unitary circuits like the above—will generally give an infinite growth rate for S , in contrast to the finite A useful starting point is to consider the n → 0 limit of 0 growth rate for S and the higher Rényi entropies.] the Rényi entropy, S . This is known as the Hartley entropy, vN and quantifies (the logarithm of) the number of nonzero eigenvalues of the reduced density matrix. Equivalently, the 2. Limit of large Hilbert space dimension Hartley entropy determines the minimal necessary value of The present quantum circuit dynamics lead to a solvable the local bond dimension in an exact matrix product model in the limit of large local Hilbert space dimension, representation [31,69] of the state: q → ∞. In this limit, all the Rényi entropies obey the dynamical rule in Eq. (13). S ðxÞ¼ logðbond dimension at xÞ: ð11Þ To show this, we consider the reduced density matrix for a cut at x, where x is the bond to which we are applying the Like the von Neumann entropy, the Hartley entropy of unitary in a given time step. We may write ρ ðt þ 1Þ in neighboring bonds can differ by at most one: terms of ρ ðtÞ and the applied unitary matrix. Averaging x−1 Trρ over the choice of this unitary, we then obtain: jS ðx þ 1Þ − S ðxÞj ≤ 1: ð12Þ 0 0 2 2 2 hTrρ ðt þ 1Þ i ¼ ½Trρ ðtÞ þ Trρ ðtÞ : x Haar x−1 xþ1 Recall that logarithms are base q. For the present, we keep q þ 1 q finite. For the random dynamics we describe above (Sec. II A), See Appendix B for details. In terms of the second Rényi the Hartley entropy obeys an extremely simple dynamical entropy S , this is rule. In a given time step, a unitary is applied at a random −S ðx−1;tÞ−1 −S ðxþ1;tÞ−1 2 2 bond, say, at x. Applying this unitary may change the q þ q −S ðx;tþ1Þ hq i ¼ : ð14Þ Haar Hartley entropy across the bond x; the entropy remains 1 þ 1=q unchanged for all other bonds. The rule for the change in S ðxÞ is that, with probability one, it increases to the The general constraint S ≤ S allows us to write 2 0 maximal value allowed by the general constraint Eq. (12): S ðx; tÞ¼ S ðx; tÞ − Δðx; tÞ; ð15Þ 2 0 S ðx; t þ 1Þ¼ minfS ðx − 1;tÞ;S ðx þ 1;tÞg þ 1: ð13Þ 0 0 0 with Δ ≥ 0. We now use Eqs. (13) and (14) to show that Δ is infinitesimal at large q. Rewriting Eq. (14) in terms of Δ, This “maximal growth” of S occurs with probability one and substituting Eq. (13), immediately shows when all unitaries are chosen randomly. Fine-tuned unitaries (e.g., the identity) may give a smaller value, but these choices Δðx;tþ1Þ Δðx−1;tÞ Δðxþ1;tÞ hq i <q þ q : ð16Þ Haar are measure zero with respect to the Haar distribution. We present a rigorous proof of Eq. (13) in Appendix A. For a simple bound (for a large system, this bound on The appendix also gives a heuristic parameter-counting max hq i is far from the tightest possible since we have not argument that suggests the same result, but as we explain, exploited the large size of the system), define Δ ðtÞ to max the more rigorous argument is necessary as the heuristic be the maximal value of Δðx; tÞ in the entire system. argument can be misleading. We note that Ref. [70] The equation above implies observed that the growth in bond dimension, when a Δ ðtþ1Þ Δ ðtÞ unitary is applied to a matrix product state, is upper max max hq i < 2q : ð17Þ Haar bounded by the right-hand side of Eq. (13) and used this to obtain an upper bound on bond dimension growth during We may iterate this by averaging over successively earlier a quantum computation. unitaries: 031016-5 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) negative sign, meaning that entanglement growth is slower when the coarse-grained surface has a nonzero slope. This is natural given the microscopic dynamics: if the slope is maximal in some region, local dynamics cannot increase the entropy there. 4. Entanglement speed in the solvable model FIG. 4. Surface growth model for entanglement Sðx; tÞ across a In the present model the difference in height between cut at x, in the large-q limit. Applying a unitary to bond x can two adjacent bonds is either ΔS ¼1 or ΔS ¼ 0. At early increase the height of the surface locally [Eq. (19)], correspond- stages of the evolution both possibilities occur. However, ing to dropping a “block” of height ΔS ¼ 1 or ΔS ¼ 2. one may argue that the “flat points,” where ΔS ¼ 0, ðln qÞΔ ðtÞ t become rarer and rarer at late times. [Flat points can max he i < 2 : ð18Þ Haar disappear by “pair annihilation” (Fig. 5, top left), and This shows that as q → ∞ at fixed time t, the probability can diffuse left or right (Fig. 5, top right), but cannot be distribution for Δ concentrates on Δ ¼ 0, so that S and S 2 0 created. As a result, their density decreases with time.] become equal. At late times the model therefore becomes equivalent to the This implies that the entanglement spectrum is flat, well-known “single step” surface growth model [71],in so, in fact, all the Rényi entropies obey Eq. (13) for the which ΔS ¼1 only. An appealing feature of this model is application of a unitary across bond x. that, for a certain choice of boundary conditions (BCs), the late-time probability distribution of the growing interface can be determined exactly [71]. [The solvable case corre- 3. Properties of the solvable model sponds to choosing periodic BCs in the classical problem. The dynamical rule we arrive at for the bipartite von (These BCs are useful for understanding the classical Neumann and Rényi entropies at large q, model, but they do not have an interpretation in terms of entanglement.) In this setting, the mean height grows Sðx; t þ 1Þ¼ minfSðx − 1;tÞ;Sðx þ 1;tÞg þ 1; ð19Þ indefinitely, but the probability distribution for the height fluctuations reaches a well-defined steady state.] This defines a stochastic surface growth model in which Sðx; tÞ shows that on scales smaller than the correlation length is always an integer-valued height profile (Fig. 4). The (and prior to saturation), the interface looks like a 1D remaining randomness is in the choice of which bond is random walk with uncorrelated ΔS ¼1 steps. This updated in a given time step. At each time step, a bond x is confirms the expected KPZ exponent α ¼ 1=2. It also chosen at random, and the height SðxÞ is increased to the allows the mean growth rate of the surface to be calculated maximal value allowed by the neighbors. Figure 5 gives [71]: the mean height increase in a given time step can be examples of local configurations before and after the calculated by averaging over the four possible initial central bond is updated. configurations for a bond and its two neighbors. After This model is almost identical to standard models for rescaling time so that one unit of time corresponds to an surface growth [71,72]. It is in the KPZ universality class (it is average of one unitary per bond, this gives an entanglement straightforward to simulate the model and confirm the growth rate [Eq. (6)]of expected KPZ exponents), and some nonuniversal properties can also be determined exactly (see below). Note that the v ¼ 1=2: ð20Þ boundary conditions S ¼ 0 on the right and the left, and the restriction jSðx þ 1Þ − SðxÞj ≤ 1, imply that the entangle- As we discuss below (Secs. IV and V), one can also ment eventually saturates in the expected pyramid profile. associate a speed v with the growth of quantum operators When we move to the continuum (KPZ) description of under the random dynamics; in the present large-q model, the interface [Eq. (4)], the nonlinear λ term appears with a this speed is (The result v ¼ 1 arises because in the large- q limit, the growth of a typical operator is limited only by the structure of the circuit. In Ref. [73], we give explicit derivations of v in random circuits for arbitrary q.) v ¼ 1: ð21Þ It is interesting to note that here v <v , contrary to 1D E B CFTs (where v ¼ v [1]) and contrary to previous con- E B jectures about generic systems [23].InSec. IV,we givean FIG. 5. Entanglement growth in the large-q model. Effect of applying a random unitary to the central bond, for four choices of appealing intuitive picture for why v can be smaller than v E B the initial local entropy configuration of three adjacent bonds. for the case of Clifford circuits. 031016-6 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) The mapping to surface growth gives us a clean derivation composite tensors L and R gives a representation of the P S P cut i i of universal entanglement dynamics in a solvable model. state as jψðtÞi ¼ L R jai jbi , where jai and a A A A i¼1 a;b b However, this surface growth picture is restricted to the jbi are basis states in A and A, respectively. This implies that entropy for a single cut (as opposed to the entropy of a region cut the Schmidt rank for a bipartition into A and A is at most q , with multiple end points) and to one dimension. It will be so that S , which is the logarithm of the Schmidt rank, is at useful to find a more general language which extends the most S . In turn, S ≤ S for any n ≥ 0.] cut n 0 above results. To do this, we now make a connection with the The best bound of this type is given by the minimal cut, second member of the KPZ triumvirate (Fig. 2), the statistical which passes through the smallest number of legs. We denote mechanics of a polymer in a random environment. the corresponding estimate for the entropy S ðxÞ.Ifthe min−cut geometry of the circuit is random, S ðxÞ and the min−cut III. DIRECTED POLYMERS AND MINIMAL CUT corresponding curve are also random. Finding S ðxÞ min−cut amounts to an optimization problem in a classical disordered In this section, we relate the dynamics of Sðx; tÞ to the geometry of a minimal cut through the quantum circuit system. which prepares the state (Fig. 6). This provides an alternative In the solvable large-q model, S ðxÞ in fact gives min−cut perspective on the exact result Eq. (19) for the solvable the von Neumann entropy exactly. This follows straight- model, and also a useful heuristic picture for noisy quantum forwardly from the results of the previous section (see dynamics in general. This line of thinking reproduces KPZ below). In a typical microscopic model, on the other hand, behavior in 1D. Importantly, it also allows us to generalize to S is only a bound on the true entropy. Nevertheless, min−cut higher dimensions and to more complex geometries. we conjecture that the following picture based on the Our starting point is the minimal cut bound for tensor minimal cut is generally valid as a coarse-grained picture: networks. This very general bound has been related to the i.e., that it correctly captures the universal properties of the entanglement dynamics in noisy systems. This conjecture is Ryu-Takayanagi formula for entanglement in holographic equivalent to the applicability of the KPZ description to conformal field theories [50–53,74], and has also been applied to unitary networks as a heuristic picture for generic noisy systems; further evidence for the latter is in entanglement growth [10]. Secs. IV and VI. The problem of finding the minimal curve is a version of a Consider again a random quantum circuit in 1 þ 1D, and a well-studied problem in classical statistical mechanics, curve like that in Fig. 6, which bisects the circuit and divides known as the directed polymer in a random environment the physical degrees of freedom into two at position x. Any (DPRE) [47,75]. Here, the “polymer” is the curve which such curve gives an upper bound on the entanglement: all the bisects the circuit, and its energy EðxÞ is equal to S ðxÞ,the Rényi entropies satisfy SðxÞ ≤ S , where S is the number cut cut cut number of legs it bisects. The spatial coordinate of the of “legs” that the curve passes through. (This relies only on polymer’s upper end point is fixed at x, while the lower end linear algebra: the rank of the reduced density matrix ρ is at S point is free. Finding S ðxÞ is equivalent to finding the cut min−cut most q .) [The cut divides the tensor network into two parts minimal value of the polymer’s energy. This corresponds to connected by S bonds. One part contains the physical cut the polymer problem at zero temperature; however, the legs for subsystem A and the other part contains those for the universal behavior of the DPRE is the same at zero and at complement. Regarding the two parts of the circuit as nonzero temperature. (For any finite temperature, the DPRE flows under renormalization to a zero-temperature fixed point at which temperature is an irrelevant perturbation.) DPRE models with short-range-correlated disorder are in the same universality class as the KPZ equation [45]. Let Eðx; tÞ be the minimal energy of the polymer in a sample of height t. We may increase t by adding an additional layer to the top of the sample. Eðx; t þ δtÞ can then be expressed recursively in terms of Eðy; tÞ for the various possible values of y. In the continuum limit, this leads to an equation for Eðx; tÞ which is precisely the KPZ equation (see Refs. [45,76] for more details of the mapping between DPRE and KPZ). The KPZ exponents we give in Sec. II may therefore be applied to the energy of the FIG. 6. Any cut through the unitary circuit that separates the polymer. The exponent z ¼ 3=2 also determines the length legs to the left and right of x (on the top boundary) gives an upper scale for transverse fluctuations of the polymer on large bound on Sðx; tÞ. The best such bound is given by the minimal cut length and time scales: (note that the cut shown here is not the minimal one). Finding the minimal cut in a random network is akin to finding the lowest- 2=3 Δx ∼ ðΔtÞ : ð22Þ energy state of a polymer in a random potential landscape. 031016-7 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) Since in our case the minimal E is simply S ,we For a finite interval of length l in an infinite system there is min−cut find that the latter executes KPZ growth. In light of the a crossover between a configuration with two vertical cuts previous section, this is not surprising. In fact, in our and one with a single horizontal cut, giving instead SðtÞ¼ solvable model, S is exactly equal to the true entan- 2v tfðl=2v tÞ. min−cut E E glement entropy (in the large-q limit). This follows from the These scaling forms are our first confirmation that v is fact that the recursive construction of Eðx; tÞ described really a speed, as well as a growth rate for the entanglement. above (on the lattice, rather than in the continuum) precisely We give an independent derivation of this fact for Clifford matches the large-q dynamics of Eq. (19). Examples of circuits in the following section, and test the above scaling nonunitary tensor networks in which the minimal cut bound form numerically in Sec. VI B. We discuss the interpreta- becomes exact are also known [52], including a large-bond- tion of v further in Sec. V. dimension limit similar to that we discuss here [53]. Note that fluctuations have dropped out of Eq. (23) as a The utility of the DPRE picture is that it is far more result of considering only the leading-order behavior of generalizable than the surface growth picture, which is Sðx; tÞ. These scaling forms agree with the results for restricted to the entropy across a single cut in 1D. As we holographic CFTs [8] and with a heuristic application of note above, the value of S in a given microscopic the minimal cut formula to a regular tensor network [10]. min−cut model is typically not equal to any of the physical entropies Here, we see them emerging from a simple and well-defined S with n> 0. Nevertheless, we conjecture that the DPRE coarse-grained picture, suggesting that they are universal for and KPZ pictures are valid universal descriptions for all all generic 1D systems, including, for example, translation- noisy models, so long as they are not fine-tuned or nonlocal. ally invariant but nonintegrable spin chains. [Reference [77] This includes noisy Hamiltonian dynamics in continuous includes numerical tests of scaling forms derived from the time (we discuss this case further in Sec. IX). If we restrict to directed polymer picture in deterministic systems, including the leading-order deterministic behavior, we can also make extensions to inhomogeneous systems (a chain with a weak link).] It is also worth noting that Eq. (24) is capable of conjectures about Hamiltonian systems without noise. distinguishing generic systems from (nonrelativistic) inte- grable systems. In the latter case the quasiparticle picture A. Scaling form for entanglement saturation applies and yields different profiles for SðtÞ [3,19].For At leading order in time, the growth of the height Sðx; tÞ relativistic systems in which the quasiparticle picture holds is deterministic: fluctuations are a subleading effect when t (rational CFTs), all quasiparticles travel at the same speed, is large. Similarly, Eq. (22) shows that the coarse-grained and as a result Eq. (24) does apply [1,3] (however, the minimal cut is essentially vertical (prior to saturation of the entanglement of more complex regions will differ between entropy): the length scale for its transverse fluctuations is the quasiparticle picture, on the one hand, and the results negligible in comparison with t. These pictures therefore from holographic systems and the minimal cut picture, on have well-defined and simple deterministic limits. They the other hand [3,4,44]). lead directly to deterministic scaling forms for the leading- In Secs. V and VIII we propose that the above picture in order behavior of the entanglement, which we discuss in terms of a coarse-grained minimal cut is the simplest way more detail in Sec. V. Here, we consider the simplest case, to understand the basic features of the entanglement the saturation of the entanglement entropy Sðx; tÞ across a tsunami for generic many-body systems (with or without single cut (or for a single interval). We reproduce a simple noise) in both 1D and higher dimensions. scaling function known from other contexts [8,10,21]. The definition of the entanglement growth rate implies IV. HYDRODYNAMICS OF OPERATOR that the energy E of such a vertical cut is v t to leading SPREADING order. The entanglement in a finite system grows at this rate until time t ¼ x=v , when it reaches its saturation An alternative way to think about the quantum dynamics is saturation E value S ¼ x. (Here, we are neglecting subleading terms and in terms of the evolution of local operators O . For example, a assuming x<L − x.) After this time a vertical cut is no Pauli operator initially acting on a single spin (e.g., O ≡ Y ; i i longer favorable: instead, the minimal cut exits the circuit we denote the Pauli matrices by X, Y, Z) will evolve with time via the left-hand side. Its shape is no longer unique, but it into an operator O ðtÞ which acts on many spins. Operators can be taken to be horizontal, and it has energy E ¼ x. typically grow ballistically [38], in the sense that the number This picture corresponds to a simple scaling form (again, of spins in the support of O ðtÞ grows linearly with t. In this neglecting subleading terms): section, we relate the growth of the bipartite entanglement to the spreading of operators. We focus on the special case of Sðx; tÞ¼ v tfðx=v tÞ; ð23Þ E E unitary evolution with Clifford circuits (defined below), but with we expect the basic outcomes to hold for more general unitary dynamics. We find that the entanglement growth rate u for u< 1 is not given by the rate at which a single operator grows, but is fðuÞ¼ ð24Þ 1 for u ≥ 1. instead determined by collective dynamics involving many 031016-8 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) number of sites to be L). These operators, denoted O (i ¼ 1; …;L), satisfy O jΨ i¼jΨ i: ð25Þ i 0 0 FIG. 7. Spreading of stabilizer operators defining the quantum state (Sec. IV). Each blue particle marks the right end point of For example, if the spins are initially polarized in the y some stabilizer (the rightmost spin on which it acts). Blue direction, we may take O ¼ Y . At a later time, the above particles hop predominantly to the right. Whenever a particle i i equation still holds, with each stabilizer O replaced with enters the right-hand region (A) the entanglement S increases by one bit. The particle density is described by the noisy Burgers the time-evolved stabilizer O ðtÞ¼ UðtÞO U ðtÞ, where i i equation, which maps to KPZ. A “hole” (empty circle) marks the UðtÞ is the unitary operator that evolves the initial state to left-hand end point of some stabilizer. the state at time t. [Note that O ðtÞ is not the standard Heisenberg picture operator, which would have U and U in the other order.] operators. Remarkably, in 1D these collective dynamics have In the following, we focus on evolution of the initial state a long wavelength hydrodynamic description. with unitary gates in the Clifford group [80]. Such gates have This hydrodynamic description turns out to be the noisy recently been used in toy models for many-body localization Burgers equation, which is related to the KPZ equation by a [29]. Entanglement generation in non-random Clifford simple change of variable and is the final member of the KPZ circuits has also been studied [43]. The defining feature of triumvirate shown in Fig. 1. In the present case the hydro- Clifford unitaries is that they have a simple action on Pauli dynamic mode is the density of certain fictitious “particles,” operators: single-spin Pauli operators are mapped to products shown in blue in Fig. 7. The quantum state is defined by a set of Pauli operators. of operators (Sec. IVA) which spread out over time, and the Any product of Pauli matrices can be written as a product particles are markers which show how far these operators of X and Z matrices, so to follow the dynamics of a given have spread. We derive their coarse-grained dynamics in stabilizer O ðtÞ, we need only keep track of which X Sec. IV B after introducing the necessary operator language. i i and Z operators appear in this product. Furthermore, the In generic many-body systems (with local interactions) overall sign of the stabilizer O ðtÞ does not affect the this process of operator growth is characterized by a speed known as the butterfly speed v . This speed defines an entanglement properties of a system undergoing Clifford evolution, so we do not keep track of it. By writing O ðtÞ as effective light cone within which the commutator between the spreading operator OðtÞ and a typical local operator v v v v is appreciable. The quantity v is a characteristic speed for 1x 1z Lx Lz O ðtÞ ∝ X Z …X Z ; ð26Þ i L L 1 1 the spreading of quantum information in a given model, and can be extracted from appropriate correlation functions. In we may specify any stabilizer by a binary vector v⃗ with 2L deterministic systems (time-independent evolution) v can components: depend on temperature, and typically does not saturate the well-known Lieb-Robinson bound [78,79]. Generic noisy v⃗ ¼ðv ;v ; …;v ;v Þ: ð27Þ systems equilibrate to infinite temperature, so in the present 1x 2x Lx Lz models there is no notion of temperature dependence—v is a constant defined entirely by the dynamics. For example, the first component of the vector v ¼ 1 if X 1 1 The scaling forms we discuss in the previous section appears in the product, and v ¼ 0 otherwise. The binary show that in 1D there is a well-defined speed v associated vector corresponding to a stabilizer O ¼ Y is i i with entanglement spreading. The following picture gives a physical interpretation of this speed, in terms of a certain set v⃗ ≡ ð0; …; 0; 1; 1; 0; …; 0Þ; ð28Þ of growing operators. However, it also shows that in general the speed v is smaller than the speed v .Thisisperhaps E B where the locations of the nonzero elements correspond to surprising: in 1D CFTs the two speeds are equal, and it X and Z . i i has been conjectured that they are equal in general [23]. We consider the dynamics in two stages. First, we (Note that we already encountered a solvable model with consider the evolution of a single operator. Then, we v ¼ v =2 in Sec. II A.) E B generalize this to understand the dynamics of the state. How does a single stabilizer O ðtÞ evolve? Applying a A. Stabilizer operators one- or two-site Clifford unitary to O ðtÞ corresponds to It is convenient to use the language of “stabilizer” applying simple local updates to the string v⃗ . Although the operators to describe the entanglement dynamics. We precise details of these updates are not crucial, we now give some explicit examples of gates that we encounter again in may define the initial state jΨ i by specifying L stabilizers under which it is invariant (in this section, we take the the numerical simulations. 031016-9 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) As single-site examples, consider the Hadamard and on v⃗ is invariant under the dynamics, regardless of the phase gates. The Hadamard is a rotation on the Bloch probabilities with which the gates are applied. By standard sphere [the rotation is by π around the (1,0,1) axis] that properties of Markov processes, this is the unique asymp- exchanges the X and Z axes, totic distribution to which the system tends, so long as the choice of Clifford gates is not fine-tuned to make the process nonergodic. (If the gate set includes each gate and R ¼ pffiffiffi ðX þ ZÞ; ð29Þ its inverse with the same probability, detailed balance is also obeyed, but this is not necessary.) We expect v⃗ to so applying a Hadamard to site i updates the string by equilibrate locally to this structureless state on an Oð1Þ time v ↔ v . The phase gate is a rotation around the Z axis ix iz scale, and similarly for the internal structure of operators which maps X to Y ¼ iX Z : i i i i smaller than L.] Now consider the dynamics of a quantum state. Once the pffiffiffiffi R ¼ Z: ð30Þ sign information in Eq. (26) is dropped, the relevant information in the state jΨðtÞi is contained in binary 1 L This means that an additional Z is generated whenever X i i vectors v⃗ ; …; v⃗ corresponding to the L stabilizers. We is present in the string, or equivalently v → v þ v iz iz ix may package this information in a 2L × L matrix: ðmod 2Þ. For a two-site example, consider the left and right 1⊤ L⊤ controlled- NOT (CNOT) gates acting on the leftmost spins ΨðtÞ¼ðv⃗ ; …; v⃗ Þ: ð32Þ in the chain. In the Z basis, the action of these operators is to flip the “target” spin if and only if the “control: spin is Each column corresponds to a stabilizer, and each row to a down: spin operator X or Z . The dynamical updates correspond i i to row operations (with arithmetic modulo two) on this ðLÞ matrix. For example, a Hadamard gate exchanges the rows CNOT ¼ ½ð1 þ Z Þþð1 − Z ÞX ; 1 1 2 corresponding to X and Z . i i A crucial point is that there is a large gauge freedom in ðRÞ CNOT ¼ ½ð1 þ Z Þþð1 − Z ÞX : ð31Þ 2 2 1 this definition of the state. This gauge freedom arises because we can redefine stabilizers by multiplying them ðLÞ Conjugating the Pauli matrices by CNOT yields: together. For example, if a state is stabilized by fX ;Z g, 1 2 then it is also stabilized by fX Z ;Z g, and vice versa. This 1 2 2 X → X X;Z → Z;X → X;Z → Z Z : freedom to redefine the stabilizers corresponds to the 1 1 2 1 1 2 2 2 1 2 freedom to make column operations on Ψ, or equivalently We see that the operator X is added to the string if X is the freedom to add the vectors v⃗ modulo two. Note that by 2 1 ðLÞ making such a “gauge transformation” we may be able to present (and similarly for Z and Z ). Applying CNOT 1 2 reduce the size of one of the stabilizers, giving a more therefore updates v⃗ by compact representation of the state. The final fact we need is an expression for the entropy v → v þ v ðmod 2Þ;v → v þ v ðmod 2Þ: 2x 2x 1x 1z 1z 2z S ðtÞ of a region A in terms of the stabilizers. Heuristically, ðRÞ this is given by the number of stabilizers that have spread CNOT acts similarly with the roles of the spins reversed. into region A from outside. More precisely, define I as the It is simple to argue that random application of such number of stabilizers that are independent when restricted operations causes the region of space in which v⃗ is nonzero to region A [81]. (Independence of the stabilizers corre- to grow ballistically. This corresponds to the operator spreading itself over a region of average size 2v t, where sponds to linear independence of the vectors v⃗ , with v is the operator spreading (butterfly) velocity for this arithmetic modulo two, once they are truncated to region system [79]. (For the present system, this velocity is also A.) The entropy is equal to [82,83] the analogue of the Lieb-Robinson velocity.) The value of v depends on the precise choice of dynamics, but it is the S ðtÞ¼ I − jAj; ð33Þ B A A same for all initial operators so long as the dynamics (the probability distribution on gates) is not fine-tuned. Further, where jAj is the number of sites in A. See Appendix C for a one may argue that the interior of the region where the simple derivation of Eq. (33). For Clifford dynamics all string v⃗ is nonzero is “structureless.” Within the interior, v⃗ Rényi entropies are equal, so we omit the Rényi index on S. rapidly “equilibrates” to become a completely random The maximal value of I is 2jAj,so S is bounded by jAj as A A binary string. [Consider the late-time dynamics of an expected. operator, or equivalently its string v⃗ ,inan L-site system. This formula has a simple interpretation. In the initial Random application of Clifford gates gives random dynam- product state we may take one stabilizer to be localized at ics to v. It is easy to see that the flat probability distribution each site, so I ¼jAj and the entanglement is zero. As time 031016-10 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) goes on, stabilizers that were initially localized outside of A commute, and this is preserved by the unitary dynamics grow and enter A. Each time a new independent operator and the redefinitions of the stabilizers). [Consider the case appears in A, the entanglement S ðtÞ increases by one bit. where ρ ðiÞ¼ 1: for example, let the corresponding sta- A r The linear independence requirement in the definition of I bilizer read O ¼ …X . Any stabilizer contributing to ρ ðiÞ A i l is crucial, as it leads to effective interactions between the must be of the form X … in order to commute with O.By stabilizers, which we discuss in the following section. the rule imposed in the text, this means that ρ ðiÞ ≤ 1.] From now on, we take A to consist of the spins to the Since there are a total of 2L particles which all have to live right of the bond x, and revert to the notation S ¼ S used somewhere, we have Eq. (34). A x in the rest of the text for the entanglement across a cut at x. With this convention, the dynamics of the bipartite entropy SðxÞ is simply related to the hopping dynamics of the particles. By Eq. (34) it suffices to consider only the r B. Coarse-grained operator dynamics particles: an l particle is just an r “hole.” We write the Each stabilizer O ðtÞ (labeled i ¼ 1; …;L) has a left and density ρ of r particles as ρ. See Fig. 7 for a typical a right end point l and r , marking the extremal spins i i configuration in a subregion of the system. included in the stabilizer. We view l and r as the positions i i The utility of the canonical form Eq. (34) is that the of two fictitious particles of type l and r, represented in independence requirement becomes trivial. One can easily white and blue, respectively, in Figs. 7 and 8. There are L of check that all the operators that have spread into A (the each type of particle in total. region to the right of x) are independent. (Consider the In the initial product state, O ð0Þ is a single Pauli stabilizers that act in region A, i.e., the stabilizers with operator on site i, say, Y . This means that each site has r >x. We may argue by contradiction that they remain one l particle and one r particle (since l ¼ r ¼ i), as i i independent after truncation to subsystem A. If not, this shown in Fig. 8 (left). As time increases, the r particles will means there is some product of the truncated stabilizers that typically move to the right and the l particles to the left. equals one. Let the rightmost spin appearing in any of these The nature of this motion depends on how we define the stabilizers be j. But by our convention for “clipping” the stabilizers. At first sight, the obvious choice is to define stabilizers, it is impossible for the Pauli matrices acting on O ðtÞ as the unitary time evolution of the initial stabilizer, spin j to cancel out when they are multiplied together. O ðtÞ¼ UðtÞY U ðtÞ. But in fact, it is useful to exploit the i i Therefore, the operators must, in fact, be independent.) gauge freedom in the choice of stabilizers to impose a Therefore, to find SðxÞ we need only count the number different “canonical” form. One result of this is that the of r particles to the right of the cut and subtract the number stabilizers effectively grow more slowly than the butterfly of sites: velocity v (discussed in the previous section) for the SðxÞ¼ ðρ − 1Þ: ð35Þ spreading of an operator considered in isolation. i>x Let ρ ðiÞ and ρ ðiÞ be the number of particles of each l r type at site i. The constraint that we impose is To reiterate, the entanglement increases by one every time an r particle drifts rightward across bond x (and decreases ρ ðiÞþ ρ ðiÞ¼ 2: ð34Þ l r by one if it drifts across in the other direction). To see that we can impose this constraint, consider the Now consider the dynamics of the particles. situation ρ ðiÞ¼ 3, so that there are three stabilizers that Microscopically, a dynamical time step involves (1) appli- start at i. The initial element of each string can be either X, cation of a unitary gate and (2) potentially a “clipping” of Y,or Z.If ρ ðiÞ¼ 3, it is impossible for all three initial stabilizers to enforce the canonical form. Effectively, the elements to be independent. We can then redefine one of particles perform biased diffusion, with the restriction that the stabilizers, by multiplying it by one or both of the more than two particles cannot share a site: others, in such a way that its length decreases by one. (By choosing the longer stabilizer we avoid adding length at the ρ ≤ 2: ð36Þ right-hand side.) Making reductions of this kind wherever possible guarantees that ρ ðiÞ ≤ 2, and also that if ρ ¼ 2, l l This constraint leads to “traffic jam” phenomena familiar the initial elements of the two stabilizers are distinct. from the so-called asymmetric exclusion process [84], and (And similarly for ρ .) With this convention it also follows to the same continuum description. Our essential approxi- that ρ ðiÞþ ρ ðiÞ ≤ 2: otherwise, the operators involved l r mation is to neglect the detailed internal structure of the could not commute, which they must (the initial stabilizers stabilizers and to treat the dynamics of the end points as effectively Markovian. We expect this to be valid at long length and time scales for the reason we mention in the previous section: the internal structure of the operators is essentially featureless and characterized by finite time FIG. 8. Left: The initial product state represented in terms of the fictitious particles. Right: A state with maximal SðxÞ. scales. 031016-11 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) We now move to a continuum description. The coarse- grained density obeys a continuity equation: ∂ ρ ¼ −∂ J; ð37Þ t x FIG. 9. Fluctuations at late times, after saturation of hSðxÞi,in with J the particle current. Further, there is a symmetry the Clifford case. When x ≪ L=2 it requires a rare fluctuation under spatial reflections, which exchange left and right end (fighting against the net drift) to remove a particle from region A, points (ρ ↔ ρ ). Writing leading to an exponentially small S ðxÞ − hSðxÞi. l r max ρ ¼ 1 þ Δρ; ð38Þ Let us revert to our previous notation, where the system has L þ 1 sites and bonds are labeled x ¼ 1; …;L. Without where Δρ is the deviation from the mean density, the loss of generality we take x ≤ L=2. When fluctuations are reflection symmetry is neglected, the region to the left of x is empty of r particles, x → −x; Δρ → −Δρ: ð39Þ and the entropy is maximal, S ðxÞ¼ x. Fluctuations will max reduce the average. But in order for SðxÞ to fluctuate To obtain a long wavelength description, we write the downward, a blue r particle must diffuse leftward from the current as a power series in Δρ and ∂ . Keeping the lowest- right half of the system in order to enter the region to the order terms that respect the symmetry, left of x, as in Fig. 9. This is a fluctuation by a distance ∼ðL=2 − xÞ. Such fluctuations are exponentially rare J ¼ c − ν∂ Δρ − ðΔρÞ þ η: ð40Þ x events, because they fight against the net rightward drift for the r particles. Thus, when L=2 − x is large, we expect These terms have a transparent physical meaning. The drift −αðL=2−xÞ constant c> 0 reflects the fact that the average motion is to S − hSðxÞi ∼ e : ð43Þ max the right (i.e., operators grow over time). The ν term is simple diffusion. The noise η reflects the randomness in the Our coarse-grained picture does not determine the numeri- dynamics. Most importantly, the nonlinear λ term is the cal constants. effect of the constraint Eq. (34). It reflects the fact that The detailed nature of these exponentially small correc- the current is maximal when the density is close to one. The tions will differ between Clifford circuits and more general current evidently vanishes when ρ ¼ 0, since there are no unitary circuits. (For example, in the Clifford case S is particles, but also when ρ ¼ 2 (the particles cannot move if independent of n, while in general the corrections will the density is everywhere maximal). Therefore, we depend on n [68].) Nevertheless, the functional form above expect λ > 0. agrees with the late-time result for generic gate sets, which From the above formulas, the density obeys is simply the mean entanglement in a fully randomized pure state [67,68]: 2 2 ∂ ρ ¼ ν∂ ρ þ ∂ ðρ − 1Þ − ∂ η; ð41Þ t x x x jAj−jAj −ðL=2−xÞ 2 4 S − hSðxÞi ≃ ¼ : ð44Þ known as the noisy Burgers equation [84]. The entangle- max 2 ln 2 4 ln 2 ment S ¼ ðρ − 1Þ obeys ∂ S ¼ J, leading to the KPZ equation: jAj¼ x and jAj¼ L − x þ 1 are the numbers of sites in A and its complement. For (generic) Clifford dynamics, the 2 2 ∂ S ¼ c þ ν∂ S − ð∂ SÞ þ η: ð42Þ probability distribution of the entanglement at asymptoti- t x x cally late times will be that of a random stabilizer state. This The sign of λ is in agreement with that obtained from the has been calculated in Ref. [85]. surface growth picture in Sec. II and from the directed polymer picture in Sec. III. While we focus here on V. ENTANGLEMENT TSUNAMI: SPEEDS dynamics of a restricted type (Clifford), this derivation AND SCALING FORMS of KPZ for entanglement provides independent support for It is not a priori obvious that the rate v governing the arguments in the previous sections. In the language of the particles, the initial state corre- entanglement growth can be viewed as a speed in generic sponds to uniform density ρ ¼ 1. Saturation of the entan- systems (see Ref. [79] for a recent discussion), although glement corresponds (neglecting fluctuations) to all of the r this is known to be the case in holographic CFTs [8].Our particles accumulating on the right-hand side and all of the l results in the directed polymer picture and in the operator particles on the left (Fig. 8), i.e., to a step function density. spreading picture suggest that v is indeed a well-defined As an aside, it is interesting to consider fluctuations in speed in generic systems. (We see in the previous section SðxÞ at late times, i.e., long after the saturation of hSðxÞi. that there is a simple visual interpretation of this speed in 031016-12 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) the stabilizer formalism.) However, this speed is in general In the 1D case, the deterministic scaling form for the smaller than the speed v which governs the spreading of entanglement (of an arbitrary region) which results from the an operator considered in isolation: “thermalization is leading-order directed polymer picture is rather simple, and slower than operator spreading.” is not new—it agrees with holographic results [8,44], and In the stabilizer context the difference between v and v as noted in Ref. [10], can also be obtained from a more E B arises because in enforcing Eq. (34) we “clip” the stabi- microscopic minimal cut picture in which the geometry of lizers, reducing their rate of growth. We believe the the minimal cut is highly nonunique. We propose that phenomenon of v being smaller than v to be general coarse graining fixes the geometry of the minimal cut. The E B (see also the result in Sec. II A 4) and relevant also to non- derivation of these scaling forms from a simple coarse- noisy dynamics. This picture is contrary to that of. e.g., grained picture suggests that they are universal in non- Ref. [23], where the operator spreading velocity is assumed integrable, translationally invariant systems. (These scaling to determine the entanglement growth rate. In the presence forms are generally not the same as those obtained from the of noise, one may also argue that a picture of independently quasiparticle picture for rational CFTs [3,4].) Our deriva- spreading operators underestimates the exponent governing tion also opens the door to generalizations to higher the growth of fluctuations. [Considering the unitary evo- dimensions (Sec. VIII B) and to 1D systems with quenched lution of a single operator in isolation, its right end point disorder [77]. executes a biased random walk, traveling an average We now consider some examples of the scaling of the 1=2 mutual information. This will help clarify the operational distance v t with fluctuations Oðt Þ. If we were to meaning of the speed v . neglect the independence requirement in Eq. (33), then To calculate the entanglement S of a region A,we must the entanglement would be estimated (incorrectly) as the take a cut, or multiple cuts, with end points on the boundary number of independently spreading operators which have points of A at the top of the space-time slice. These cuts can reached A. The mean of this quantity is v t and the 1=4 either be vertical, in which case they cost an energy’v t fluctuations are of order t . This is related to the differ- (to use the language of Sec. III), or they can connect two end ence between the KPZ universality class of surface growth, points, in which case we take them to be horizontal and to which is generic, and the Edwards-Wilkinson universality have an energy equal to their length. The entanglement S ðtÞ class, which applies when the strength of interactions is A is given by minimizing the energy of the cut configuration. fine-tuned to zero [45].] It is a continuous piecewise linear function, with slope The language of a tsunami is often used in discussing discontinuities when the geometry of the minimal cut entanglement spreading, so it is nice to see that–at least in configuration changes. To generalize the conjecture to 1D–entanglement spreading can be related to a hydro- systems without noise, we must allow for the fact that the dynamic problem. (The motivation for the tsunami termi- asymptotic value of the entanglement depends on the energy nology is the idea that for a region A, the entanglement S density of the initial state. We therefore replace the entan- is dominated by a subregion close to the boundary which glement S in the formulas with S=s , where s is the grows ballistically, like the advancing front of a tsunami.) eq eq equilibrium entropy density corresponding to the initial In higher dimensions the boundary of an operator has a more complicated geometry, so the hydrodynamic corre- energy density [2,8]. This ensures that the entanglement spondence we describe above does not generalize. entropy of an l-sized region matches the equilibrium thermal In order to understand the entanglement tsunami better, entropy when v t ≫ l=2, as required for thermalization. we now return briefly to the directed–polymer–in–a– Heuristically, s defines the density of “active” degrees of eq random–medium picture developed for noisy systems in freedom at a given temperature [79]. Sec. III. To clarify the meaning of v , consider the mutual information between two semi-infinite regions that are separated by a distance l (Fig. 10). With the labeling of A. Scaling forms for the entanglement tsunami the regions as in the figure, this is given by When all length and time scales are large, fluctuations in the entanglement are subleading. Neglecting them is equiv- I ¼ S þ S − S ¼ S þ S − S : ð45Þ AC A C A∪C A C B alent to saying that the coarse-grained minimal cut (prior to saturation) is a straight vertical line. This deterministic picture generalizes to the entanglement or mutual informa- tion of arbitrary regions, and also to higher dimensions (Sec. VIII). We conjecture that these pictures are valid for the long-time behavior of entanglement quite generally. The FIG. 10. Infinite chain with regions A, B, C marked. B is of setup relevant to us in the non-noisy case is a quench, in length l while A and C are semi-infinite. The mutual information which the initial state is a ground state of one Hamiltonian, between A and C is nonzero so long as l< 2v t: correlations and a different Hamiltonian is used for the evolution. exist over distances up to 2v t, not v t. E E 031016-13 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) We have S ¼ S ¼ v t for all times, since the appropriate A C E minimal cuts are vertical. If l> 2v t, S is given by two E B vertical cuts, so I vanishes. When l< 2v t, S is instead AC E B dominated by a horizontal cut, so that I ¼ 2v t − l. AC E The entanglement tsunami is sometimes taken to mean that at time t,a “boundary layer” of width v t inside a given region is entangled with the exterior. If this region were maximally entangled with the exterior, this would repro- duce the correct value of the entanglement across a cut (S ¼ v t). However, this picture is not correct: the result for the mutual information shows that correlations exist over FIG. 13. Bottom: Semi-infinite chain with regions A, B (length distances up to 2v t, not v t. So although v is a speed, it E E E l and l , respectively), and C adjacent to the boundary. Top: The A B should not be thought of as the speed at which the boundary mutual information between B and C for this geometry, for the of the entangled region moves. two regimes indicated. Although the rule for calculating the entanglement is almost trivial, the consequences are not always intuitively [For a simpler example of exponentially small values of obvious. First consider the case where the regions A and C the mutual information, consider hI i at infinite times in a above are finite rather than infinite (and embedded in an AC finite system. If the system contains L qubits and A ∪ C infinite chain); see Fig. 11. When the length d of the regions contains N qubits, the mutual information is exponentially A and C exceeds their separation l, the time dependence of small whenever N< L=2,and given by Eq. (44) as the mutual information is as shown in Fig. 11. [The sequence −1 −ðL−2NÞ hI of minimal cut configurations required for calculating S i ∼ ð2 ln 2Þ 2 .] [The sequence of cuts for S AC B in this case is simply (a), (c)]. in this case is (a), (b), (c) shown in Fig. 12.] By contrast, when Finally, consider the effect of a boundary. Take a semi- the separation l exceeds the length d, the mutual information infinite chain with regions A, B, C adjacent to the boundary is always zero (or more precisely, exponentially small). as in Fig. 13 (C is semi-infinite). Consider the mutual information between B and C, I ¼ S þ S − S .We BC B C A must distinguish the case l <l =2 from the case A B l >l =2. (In the former case, the first event is that the A B minimal cut at the boundary of A goes from being vertical to being horizontal; in the latter, the first event is that the two vertical cuts at the boundary of B are replaced by a horizontal one.) The resulting expressions for I are BC plotted in Fig. 13. VI. NUMERICAL EVIDENCE FOR KPZ GROWTH We now give numerical evidence that noisy entangle- FIG. 11. Bottom: Infinite chain with finite regions A and C each ment growth in 1D is in the KPZ universality class. We of length d, separated by distance l. Top: The mutual information study the time evolution of spin-1=2 chains with open between A and C in the case d> l. In the opposite regime the boundary conditions, taking the initial state to be a pro- mutual information vanishes. duct state with all spins pointing in the same direction (either the positive y or z direction) and keeping track of the entanglement entropy across each bond during the evolu- tion. The discrete time evolution is a circuit of one- and two-site unitaries. Figure 14 shows the structure of a single time step: two layers of two-site unitaries are applied, one layer on odd and one on even bonds, together with single- site unitaries. Each unitary is chosen independently and randomly (from a certain set specified below). We use the symbol R to denote a generic single-site unitary and U to denote a two-site unitary. We consider three kinds of dynamics, distinguished by FIG. 12. Sequence of minimal cut configurations (red lines) the choice of unitaries. To begin with, we study “Clifford determining the entropy of region B in Fig. 11. (a) gives way to (b) when 2v t ¼ l and (b) gives way to (c) when 2v t þ l ¼ 2d. evolution” in which the unitaries are restricted to the set of E E 031016-14 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) In all our simulations, each two-site unitary U in the circuit is chosen with equal probability from three pos- sibilities: the two types of CNOT gate [Eq. (31)] and the identity matrix: ðLÞ ðRÞ U ∈ f1; CNOT ; CNOT g: ð46Þ In this section, we discuss the simplest Clifford dynamics, which includes only these gates, and no one-site unitaries (R ¼ 1). When the initial state is polarized in the y– direction, this set of gates is sufficient to give nontrivial entanglement evolution, with universal properties that turn FIG. 14. Schematic structure of a layer in the quantum circuits out to be the same as those for more generic gate sets. We used for simulations. also study the “full” Clifford dynamics in which all the Clifford generators are used, choosing the single-site so-called Clifford gates (Sec. IV). Clifford evolution can be unitaries randomly from the three options simulated efficiently (in polynomial time) using the stabi- lizer representation we discuss in Sec. IV. This allows us to R ∈ f1;R ;R g: ð47Þ H P access very long times and to pin down KPZ exponents accurately. Next, we study more general dynamics for Results for this case are similar and are given in which polynomial-time classical simulation is impossible, Appendix. D. giving evidence that KPZ behavior holds more generally. To begin with, Fig. 15 shows the evolution of the The two types of non-Clifford dynamics we study here are bipartite von Neumann entropy SðxÞ (in units of log 2) referred to as the phase evolution and the universal for a single realization of the noise (i.e., a particular random evolution: we give details below. For these dynamics we circuit) in a system of L ¼ 459 sites. The curves show use a matrix product representation of the state imple- successively later times. Note that the entropy saturates mented via ITensor [86]. more rapidly closer to the boundary, because the maximum The fingerprints of KPZ behavior that we search for are entanglement across a bond is proportional to its distance the two independent critical exponents β and α (Sec. II). We from the boundary. At very late times Sðx; tÞ saturates to a extract β both from the fluctuations in the von Neumann pyramidlike profile representing close-to-maximal entan- entropy and from the corrections to the mean value glement. Our interest is in the stochastic growth prior to [Eqs. (6) and (7)], and we extract α from the spatial saturation, which we show is KPZ–like. All observables in correlations in the entanglement at distances shorter than the correlation length ξðtÞ [Eq. (9)]. For Clifford circuits, we also touch on the entanglement probability distribution. A. Clifford evolution Clifford circuits, or “stabilizer circuits,” are a special class of quantum circuits that play an important role in quantum information theory. As shown by Gottesman and Knill, they can be simulated efficiently, even when the entanglement entropy grows rapidly, by representing the quantum state in terms of stabilizers [87]: see Sec. IV. The time evolution operator for a Clifford circuit belongs to the Clifford group, a subgroup of the unitary group on the full Hilbert space. This group may be generated by a small set of local Clifford gates: the two-site controlled NOT gates [Eq. (31)] and the single-site Hadamard and phase gates R and R [Eqs. (29) and (30)]. For circuits H P FIG. 15. The von Neumann entropy Sðx; tÞ for a system of built from these gates, time evolving the state on L spins up length L ¼ 459, as a function of x, for several successive times to a time t takes a computational time of order Lt, and (t ¼ 340, 690, 1024, 1365, 1707, 2048, and 4096), in the Clifford measuring the entanglement across a given bond in the final evolution. This shows that the state evolves from a product state state takes a time of order L at most. This is in sharp to a near-maximally-entangled one. Prior to saturation the contrast to the exponential scaling that is inevitable for entanglement displays KPZ-like stochastic growth. Sðx; tÞ is in more general circuits. units of log 2. 031016-15 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) FIG. 16. The von Neumann entropy Sðx; tÞ in units of log 2, far from the boundaries, in a system of length L ¼ 1025 at various times (from bottom to top t ¼ 170, 340, 512, and 682) evolved with the Clifford evolution scheme. ξ schematically shows the 1=z typical correlation length Eq. (8), which grows in time like t . the following are measured far from the boundary, in order to avoid finite-size effects associated with saturation. FIG. 17. Top: Growth of the mean entanglement with time for Figure 16 shows successive snapshots for a subregion of the Clifford evolution with only CNOT gates (in units of log 2). a larger system of L ¼ 1025 bonds (times t ¼ 170, 340, The solid red curve is a fit using Eq. (49). The exponent β is found 512, 682, from bottom to top). The maximal slope that can to be β ¼ 0.33  0.01, in agreement with the KPZ prediction appear is 1, in accord with Eq. (3). Note the gradual β ¼ 1=3. Dashed line shows asymptotic linear behavior. Bottom: roughening of the surface and the growing correlation Growth in the fluctuations in the entanglement with time. The length. dashed line shows the expected asymptotic behavior, wðtÞ ∼ t , Figure 17 shows the height and width of the growing with β ¼ 1=3. The fit includes a subleading correction: Eq. (49), surface, with β ¼ 0.32  0.02. Error bars denote the 1σ uncertainty. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðtÞ¼hS ðx; tÞi;wðtÞ¼ ⟪S ðx; tÞ⟫; ð48Þ vN vN β ¼ 0.33  0.01; β ¼ 0.32  0.02: ð50Þ h w for very long times. These quantities are averaged over at Both estimates of β are in excellent agreement with the least 10 realizations. In each realization only the entan- KPZ value β ¼ 1=3. The solid lines in Fig. 17 show the fits glement across the center bond is used (therefore all data (the fit parameters are in Table I). The dashed lines show points are uncorrelated) and the system size is L ¼ at, the slopes corresponding to the expected asymptotic power 1=3 where a is chosen to avoid finite-size effects (see Sec VI B). laws, hðtÞ ∼ t and wðtÞ ∼ t . We obtain estimates β and β of the exponent β by fitting The analysis in Sec. IV implies that v is a well-defined h w the data to the expected forms [cf. Eqs. (6) and (7)]: velocity, and v t is a sharply defined length scale character- izing the range of entanglement in the state. We may β β η h w hðtÞ¼ v t þ Bt ;wðtÞ¼ Ct þ Dt : ð49Þ confirm this by measuring this length scale directly; see the section below. Here, η (with η < β ) captures subleading corrections. Note the small value of the subleading exponent η obtained We find from the fit. This implies that finite time corrections are TABLE I. Summary of all fitting parameters to Eq. (49) used in this section. The errors set the estimated 2σ uncertainty. Evolution type v β B β C η D E h w Clifford 0.1006  0.0001 0.33  0.01 0.66  0.04 0.32  0.02 0.28  0.05 0.08  0.10.16  0.04 Phase 0.133  0.03   0.54  0.04  0.223  0.004   0.168  0.03 Universal 0.202  0.001   0.09  0.005  0.14  0.003   0.36  0.01 031016-16 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) providing further support for KPZ universality in this system. B. Numerics on speeds and scaling forms We argue in Sec. V that in addition to determining the entanglement growth rate, v can also be viewed as a speed. This is the speed of the fictitious particles in Sec. IV. Operationally, the simplest manifestation of this speed is in the saturation behavior of the entanglement. The analytical arguments imply that to leading order (at large t and l) the entanglement across a cut at position l (l ≤ L=2) has the simple scaling form given above in Eq. (24): FIG. 18. The logarithmic derivative of the width dw=d log t u for u< 1 versus time for the Clifford evolution. The universal behavior S ¼ v tfðl=v tÞ;fðuÞ¼ ð51Þ A E E 1=3 1 for u ≥ 1. with exponent t is observed at shorter time scales compared with Fig. 17. This gives simply S ¼ v t for t < l=v , and S ¼ l for A E E A t > l=v . This means that there is no influence of the reduced if we plot the numerical derivative dw=d log t rather boundary at times t < l=v . (See also the numerical results 1=3 than w itself (both quantities scale as t at long times). This in Refs. [39,40], indicating sharp saturation in circuits 1=3 is done in Fig. 18. The data fit well to the t law even at short where interactions between any pair of spins are allowed.) times. This will be useful for the more general dynamics In Fig. 20, we test this result numerically for the Clifford where long times are not available. evolution. We set l ¼ L=2 and plot To complete the check of the two independent KPZ SðL=2;tÞ L=2 exponents, Fig. 19 shows the spatial correlator GðrÞ vs ð52Þ defined in Eq. (9), as a function of separation r, for three v t v t E E successive times. For small r, the correlation grows like a α 9 10 as a function of L, for several values of the time (t ¼ 2 , 2 , r with α ≃ 1=2, in agreement with the KPZ prediction for 11 12 this exponent. For distances r ≫ ξðtÞ, the correlator satu- 2 ,and 2 ). Here, v ¼ 0.1 is taken from the fits to Fig. 17. rates to a value proportional to wðtÞ. The figure gives an According to Eq. (51), this plot should converge for large t to idea of the size of the correlation length ξðtÞ for these times. aplotof fðuÞ against u. The results are in excellent Finally, recent advances in KPZ theory have yielded an agreement with the scaling form, confirming, for the case analytical expression for the full probability distribution of of Clifford circuits, that v is a meaningful velocity. the KPZ height field [56–66]. In Appendix F we show that It is also interesting to compare the entanglement velocity this analytical result compares well with Clifford numerics, v with the butterfly velocity v .Weobtain v from the E B B average spatial extent W of a growing Pauli string (see Sec. IVA) under the unitary Clifford dynamics at time t,as v ¼ W=2t. Remarkably, we find that v ¼ v =2 within B E B 2 1=2 FIG. 19. Correlation function GðrÞ¼h½SðrÞ − Sð0Þ i at time t ¼ 512, 1024, and 2048 for the Clifford evolution, showing FIG. 20. The entropy across the center of the chain (in units of excellent agreement with the KPZ prediction GðrÞ ∼ r with log 2) divided by v t versus L=2v t for various fixed values of t. E E χ ¼ 1=2 in the regime r ≪ ξðtÞ. This plot converges nicely to the scaling form in Eq. (51). 031016-17 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) “ ” FIG. 21. The average size W of a growing Pauli string as a function of time for two protocols, CNOT evolution (subscript “CNOT”) and the full Clifford evolution (subscript “Cliff”). The correspondence with the dashed lines, showing the average entanglement entropy multiplied by four, is consistent with v ¼ v =2. (Taken from a system of size L ¼ 1024.) E B FIG. 22. Top: Growth of the mean entanglement as a function numerical precision for both the CNOT-only Clifford dynam- of time for the universal and phase gate set fitted to Eq. (49), with ics and the “full” Clifford dynamics defined above. This is β set to 1=3. The dashed line shows the expected asymptotic shown in Fig. 21, wherewe plot W starting versus time for the behavior for comparison. Error bars indicate 1 standard deviation two protocols. The initial Pauli strings we consider in this (1σ) uncertainty. simulation are single-site Y operators. [The CNOT dynamics is not ergodic on the space of Pauli strings (unlike the full Universal evolution.—This set of gates, unlike the Clifford dynamics). Nevertheless, any operator grows in size others, is “universal” in the quantum information sense at the same rate v .] We compare W with 4 times the average (any unitary acting on the full Hilbert space of the spin entanglement entropy, 4SðtÞ. The two curves lie on top of chain can be approximated, arbitrarily closely, by a product each other, consistent with v =v ¼ 1=2. E B of gates from this set). The single-site gates include the We also find the same ratio for v =v in the exactly E B eightfold rotations mentioned above, together with the solvable large-q model (Sec. II A 2). However, it is possible Hadamard gate R [. (29)]. R is applied with probability to construct non-fine-tuned random circuits, involving Haar- H H 1=2 and the rotations with probability 1=16 each. random unitaries at finite q, in which the ratio is less than 1=2 Figure 22 shows the height and width hðtÞ and wðtÞ for the [73], so this value is not universal. A natural question is two protocols (averaged over 380 realizations for the phase whether it is generic for random Clifford circuits. evolution and 200 realizations for the universal evolution, and over bonds x with 20 <x< 480). The figure shows fits C. Universal and phase evolution to the forms in Eq. (49) with β and β fixed to the KPZ value h w The phase and universal dynamics take us outside the and η fixed to zero (fit parameters are in Table I). The fits with Clifford realm, and cannot be simulated efficiently on a Eq. (49) are consistent with the data. It is not possible to classical computer (in polynomial time). We give evidence extract precise estimates for β from the slope of the log-log that the correspondence with KPZ continues to hold in this plot of wðtÞ, although for the phase evolution the slope at late more generic situation. However, our results are not as times is in reasonable agreement with the expected KPZ conclusive as in the Clifford evolution as we do not have value, shown by the dashed gray trend line. access to such long times. For an alternative attack on β, we plot the numerical The simulations are performed on spin-1=2 chains of derivative dwðtÞ=d ln t. Recall that in the Clifford case the length L ¼ 500 bonds (501 spins) using the ITensor slope of this quantity (when plotted against time on a log- package [86]. The two types of dynamics are defined as log plot) has smaller finite-size corrections than the slope follows. [The two-site unitaries are always chosen from the for wðtÞ itself. The corresponding plot is shown in Fig. 23, set in Eq. (46); the initial state is taken polarized in the y for times up to t ¼ 25 (averaging over more than 5000 1=3 direction.] realizations). The dashed gray lines are the t trend lines. Phase evolution.—Each single-site unitary is chosen Results for both types of dynamics are in good agreement randomly and uniformly from the set of eightfold rotations with the expected slope β ¼ 1=3. about the z axis in spin space: R ¼ exp ðπinσ =8Þ, with Next, we examine the spatial correlator Eq. (9) in Fig. 24. n ∈ 1; …; 8. For both types of dynamics, the behavior for r ≪ ξðtÞ 031016-18 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) VII. FREE FERMIONS ARE NONGENERIC The growth of entanglement in systems of free particles is highly nongeneric. In the presence of noise, the entan- glement of a system of free particles on the lattice grows pffiffi only as S ∼ t, in contrast to the behavior S ∼ t of generic systems. The case of spatially homogeneous noise has been discussed recently [88]. The basic point is the same when the noise varies in space: the fact that the single-particle wave functions spread diffusively in the presence of noise implies that the entanglement cannot be larger than pffiffi Oð tÞ [88]. FIG. 23. The logarithmic derivative of the width dw=d log t As a concrete example, consider a short-range hopping versus time for the phase and universal evolution protocols. For Hamiltonian for free fermions, 1=3 comparison, we plot the universal behavior with exponent t in gray (dashed line). (The derivative is calculated using three data HðtÞ¼ H ðtÞc c ; ð53Þ ij j points. Errors are estimated from maximal and minimal slopes ij obtained within 1 standard deviation from the averaged data points.) with noisy matrix elements H ðtÞ. For simplicity, take the ij initial state to consist of particles localized at sites i ∈ S for some set S; for example, we could take S to consist of all the even-numbered sites: jΨð0i¼ c j0i: ð54Þ i∈S Under the evolution, each creation operator evolves into a superposition of creation operators, † † ðiÞ c → ψ ðj; tÞc ; ð55Þ i j ðiÞ where ψ ðj; tÞ is the solution of the time-dependent Schrödinger equation for a particle initially localized at ðjÞ i. In the absence of noise, ψ spreads ballistically, but in the presence of noise. it spreads only diffusively. The fact pffiffi that each creation operator is spread out over only Oð tÞ sites after a time t immediately implies that the mean pffiffi entanglement is at most of order t. (See also Ref. [88].) Note, however, that this argument does not tell us how large the fluctuations are. (Random unitary evolution of a single wave packet is discussed in Ref. [89].However,wemust 2 1=2 FIG. 24. Correlation function GðrÞ¼h½SðrÞ − Sð0Þ i at consider the full many-body wave function, since the three values of the time for the phase (top) and universal (bottom) formalism of Ref. [90] for the free-fermion density matrix gate sets, showing good agreement with the KPZ exponent value shows that the initially occupied orbitals do not simply α ¼ 1=2. contribute additively to the entanglement.) pffiffi We have confirmed numerically that hSi ∝ t for a agrees well with the KPZ exponent value α ¼ 1=2 at the noisy 1D hopping model, using the formalism of Ref. [90] largest available time. to construct the reduced density matrix. This is much The very long times accessible in the Clifford simulation slower than the linear-in-time growth of generic interacting allow us to establish KPZ exponents with high accuracy pffiffi models. The t scaling should apply for free fermions in there. For the more generic dynamical rules we cannot any number of dimensions. In 1D it also applies to certain reach the same level of precision, but, nevertheless, the noisy spin models via the Jordan-Wigner transformation: KPZ exponent values are compatible with the data. for example, the transverse field XY model: Next, we briefly discuss a fine-tuned situation in which entanglement dynamics are not KPZ-like: namely, when y y x x HðtÞ¼ fJ ðtÞ½σ σ þ σ σ þ h ðtÞσ g: ð56Þ the system is made up of free particles. Then, in Sec. VIII, i i i iþ1 i iþ1 i we move to higher dimensions. 031016-19 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) However, any generic perturbation to the spin chain spoils the free-fermion correspondence. We then expect the generic KPZ behavior to reassert itself. VIII. HIGHER DIMENSIONS We discuss several ways of thinking about entanglement growth in 1D. One of these, the directed polymer picture, generalizes naturally to higher dimensions: the polymer is simply replaced by a d-dimensional membrane embedded in (d þ 1)-dimensional space-time. As in 1D, we think of this membrane as a coarse-grained version of a minimal FIG. 26. Minimal membrane for a disk-shaped region in d ¼ 2. cut bisecting a unitary circuit. The membrane is subject to pinning by “disorder” in space-time arising from the dynamical noise. See Fig. 25 for the two-dimensional case. A. Universal fluctuations of SðtÞ in noisy systems We can explore two kinds of questions using this picture. Consider the entanglement SðtÞ for a region A whose First, we can examine universal properties that are specific to boundary ∂A has length or area j∂Aj. In the d ¼ 2 case, the noisy scenario: as in 1D, fluctuations are governed by shown in Fig. 25, j∂Aj is the length of the spatial boundary. universal exponents. Second, we can calculate leading-order Neglecting fluctuations, the “world volume” of the minimal properties of SðtÞ that do not involve fluctuations and that membrane scales as j∂Aj × t. This gives the leading scaling are, therefore, likely to be valid even in the absence of noise, of the membrane’s energy and, hence, of the entanglement. i.e., for dynamics with a time-independent Hamiltonian. As in 1D, subleading terms encode universal data. We now In higher dimensions the behavior of SðtÞ has nontrivial consider these terms. dependence on the geometry of the region for which we The pinning of a membrane or domain wall by disorder calculate the entanglement. We suggest the “minimal mem- is well studied [47,91–94] (a brief summary is in brane in space-time” as a simple and general heuristic for Appendix E). Translating standard results into the language such calculations. Below, we discuss the case of a spherical of the entanglement in a d-dimensional noisy quantum region (Sec. VIII B) and contrast our results with an alter- system, we find that in both d ¼ 1 and d ¼ 2 there is a native simple conjecture. For other toy models for entangle- unique dynamical phase with nontrivial critical exponents. ment spreading, see Refs. [10,23]. The same is true for continuum systems [more precisely, for Denoting the region for which we wish to calculate the systems with continuous (statistical) spatial translational entropy by A, and its boundary by ∂A, the membrane lives symmetry] in d ¼ 3. However, if a lattice is present, two in a space-time slice of temporal thickness t, and terminates stable phases (and thus a dynamical phase transition) are possible in d ¼ 3; one with nontrivial exponents and one at ∂A on the upper boundary of this time slice; see Fig. 25. with trivial ones. In the trivial phase, the membrane is For simple shapes and for times shorter than the saturation “smooth” and is pinned by the lattice. In the nontrivial time, the membrane also has a boundary on the lower slice, phases, the membrane is instead pinned by disorder in a as shown in Figs. 25 and 26. In this section, we focus on “rough” configuration. We discuss the nontrivial phases entanglement growth prior to saturation. (which are the only ones possible in d< 3 and for continuum systems in d ¼ 3). Generally, fluctuations have a weaker effect in higher dimensions than in 1D. For simplicity, take a quantum system which is infinite in one direction and of size L in the other d − 1 directions, and consider the entanglement for a cut perpendicular to the infinite direction. Since A and its complement are both infinite, SðtÞ will grow indefinitely for this geometry. However, there are two regimes, t ≲ L and t ≫ L (here, we drop a dimensionful prefactor). For times t ≲ L (see Appendix E for details): d−1 θþ1−d hSðtÞi ¼ L ðv t þ Bt þÞ; ð57Þ 2 1=2 ðd−1Þ=2 θ−ðd−1Þ=2 ⟪SðtÞ ⟫ ∝ L t ; ð58Þ where the exponent θ is defined below. This reproduces the FIG. 25. Minimal membrane picture for the entanglement of two regions in d ¼ 2. 1D result with θ ¼ β. Note that when d> 1, fluctuations 031016-20 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) are suppressed with respect to the mean by a factor of However, this simplifies in the limit of interest. We first 1=2 send t; R → ∞ with t=R fixed. In this limit r_ðt Þ remains j∂Aj : distant regions of the boundary give rise to finite, but the curvature terms become negligible (see, for essentially independent fluctuations which add up incoher- example, the explicit solution below), so we can write ently. In the opposite regime t ≫ L, the temporal dimen- E ¼ Eðr_Þ. Now we make the second approximation that sion of the membrane is much larger than its spatial dimensions, so there is a crossover to the 1D directed t=R is small, meaning that we can keep only the Oðr_ Þ polymer problem. However, the exponent of the higher- correction. dimensional problem appears in the universal L depend- The boundary condition at the top of the space-time slice is rðtÞ¼ R. We consider times prior to saturation, so the ence of the growth rate: membrane also has a free boundary on t ¼ 0. In the d−1 θ−1 SðtÞ¼ðvL þ wL Þt þ  : ð59Þ relevant limit, its energy is 1=3 (The higher corrections will include the t term associated 0 0 0 2 E ¼ 2πv dt rðt Þ½1 þ ar_ðt Þ þ : ð62Þ with the 1D universality class.) Numerically, the exponent is θ ¼ 0.84ð3Þ in d ¼ 2 and θ ¼ 1.45ð4Þ in d ¼ 3 [94]. The Minimal energy requires the boundary condition r_ð0Þ¼ 0. subleading exponent in Eq. (57) is negative for d> 1,so When t=R is small, we may expand in 1=R. This gives this correction may be hard to observe numerically. 0 2 02 rðt Þ ≃ R − ðt − t Þ=ð4aRÞ, as illustrated in Fig. 26. The corresponding entropy is B. Minimal membrane picture for dynamics without noise SðtÞ¼ 2πv Rt 1 − þ  : ð63Þ In higher dimensions we can ask how SðtÞ depends on 2 12aR the geometry of region A when this geometry is nontrivial. This calculation generalizes trivially to higher dimensions, Interestingly, the membrane picture makes predictions where the correction is of the same order. Corrections due about this that do not involve the noise-induced fluctua- to fluctuations come in with negative powers of t, and are tions, and that are likely also to be valid for Hamiltonian negligible in the limit we are discussing. dynamics without noise (with the replacement S → S=s eq Note that the first correction in the brackets in Eq. (63) is we discuss in Sec. V). of order ðt=RÞ , and not of order t=R. This result differs As an instructive special case, take A to be a disk-shaped from what one might naively have expected if one guessed region of radius R in d ¼ 2. (A ball in higher dimensions is that at time t an annulus of width v~ × t inside the disk is precisely analogous.) We assume continuous rotational entangled with the outside, where v~ is a tsunami velocity. symmetry, at least on average. At short times, the leading This picture gives an entropy proportional to the area of the scaling of the entanglement is SðtÞ ≃ 2πv Rt, since the annulus, world-sheet area of the membrane is approximately 2πR × t. However, there are corrections to this arising vt ~ 2 2 SðtÞ∼πR − πðR − vt ~ Þ ¼ 2πRvt ~ 1 − ; ð64Þ from the curvature of ∂A. We consider the limit of large t and large R with a fixed ratio t=R. In this regime, the effects of fluctuations may be leading to a negative correction of order t=R. The difference neglected, and instead the energetics of the membrane are between Eqs. (63) and (64) also indicates that a picture in terms of independently spreading operators is misleading, determined by deterministic elastic effects. We write the in agreement with what we find in 1D. energy of the membrane as It is interesting to note that in the regime where t=R is of E ¼ d SE; ð60Þ order 1, the full r_ dependence of Eðr_Þ plays a role. This suggests that an infinite number of nonuniversal parameters where d S is the membrane’s area element and E is its enter the expression for SðtÞ in this regime, and that there is “energy” density. We get SðtÞ by minimizing E with no general, universal scaling form for the entanglement of a appropriate boundary conditions. sphere in d> 1. However, we do expect saturation to Next, we Taylor expand E in terms of local properties remain discontinuous, as in 1D [Eq. (24)], occurring via a of the membrane. For a flat “vertical” membrane (i.e., with transition between an optimal membrane configuration that normal perpendicular to the t axis), E ¼ v . In general, reaches the bottom of the space-time slice and one (with however, E depends on the angle φ by which the surface E ¼ πR ) that does not. locally deviates from verticality, as well as, for example, the IX. OUTLOOK local curvatures κ and κ in the spatial and temporal s t directions. Using rotational symmetry to parametrize the Quantum quenches generate complex, highly entangled membrane by the radius rðt Þ, states whose dynamics cannot usually be tracked explicitly. For this reason, analytical approaches to quenches have 2 2 2 2 4 d SE ¼v rdθdtð1þar_ þbκ þcκ þcr_ þÞ: ð61Þ E t s typically relied on additional structure in the quantum 031016-21 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) dynamics: for example, integrability, or absence of inter- As a more speculative question in the domain of actions, or conformal invariance. This paper instead studies noisy dynamics, we may ask whether there exist time- dynamics that are as unstructured as possible. We propose independent Hamiltonians that show KPZ entanglement that noisy dynamics are a useful toy model for quantum fluctuations, despite the absence of explicit noise, in some quenches in generic (nonintegrable, non-conformally- dynamical regimes. We emphasize that this seems unlikely invariant) systems. on asymptotically long time scales for a generic system Many of our results are, of course, specific to noisy (since local reduced density matrices and observables dynamics: in particular, the emergence of KPZ behavior at will eventually thermalize), but it may hold on intermediate long wavelengths in 1D, and the detailed pictures for time scales in certain systems in which some degrees of entanglement growth afforded by the “KPZ triumvirate.” freedom act effectively as chaotic classical variables and But we suggest that some of our heuristic pictures apply to provide effective noise. non-noisy entanglement growth as well (with the replace- In the text, we discuss only initial states with area-law ment S → S=s mentioned above). We propose a general entanglement. A natural extension is to initial states with, for eq directed polymer picture or minimal membrane picture for example, submaximal volume-law entanglement. The natu- the scaling of the entanglement and mutual information ral expectation, say, in 1D, is that the directed polymer picture (Secs. V, VIII B) and we use the operator spreading picture extends to this case if we glue the unitary circuit to a tensor to clarify the meaning of the “entanglement velocity” and network representation of the initial state. Then the entropy its distinction from the operator spreading velocity Sðx; tÞ would include a fluctuating part with KPZ exponents (Sec. V). “Thermalization is slower than operator spread- together with a contribution from the initial state. Another ing” in generic 1D systems (i.e., in general ,v is smaller natural direction to explore is the role of conserved quantities. than v ): by contrast, this is not true in 1+1D CFTs [2],or Turning to higher dimensions, it would also be useful to in certain toy models [23]. It would be interesting to make test the higher-dimensional membrane pictures of Sec. VIII, more detailed comparisons with holographic models [8]. perhaps exploiting Clifford circuits to reduce the numerical Many interesting questions remain. First, within the difficulty of higher-dimensional dynamics. realm of noisy systems, an analytical treatment for the There are also many further questions regarding deter- regime with weak noise would be desirable, i.e., for ministic systems for which the tools we introduce here may dynamics of the form give insight. For example, the forthcoming Ref. [77] will discuss entanglement growth in disordered or inhomo- HðtÞ¼ H þ λH ðtÞ; ð65Þ 0 1 geneous spin chains from the point of view of surface where H is a time-independent many-body Hamiltonian, growth, while Ref. [73] will give results for the spreading of H ðtÞ represents noise, and λ is small. Our conjecture is that quantum operators. KPZ exponents apply for any nonzero value of λ [unless HðtÞ is fine-tuned]—i.e., that there is no universal distinction ACKNOWLEDGMENTS between continuous time dynamics and quantum circuits. (Note that there is no distinction between these two cases at We thank M. Kardar for useful discussions. A. N. and J. R. the level of conservation laws: once noise is added, energy is acknowledge the support of fellowships from the Gordon not conserved even in the continuous time case.) However, and Betty Moore Foundation under the EPiQS initiative our derivations and numerics correspond, roughly speaking, (Grant No. GBMF4303). A. N. also acknowledges support to the large λ regime. Perhaps the opposite regime could be from EPSRC Grant No. EP/N028678/1. S. V. is supported addressed using a more explicit renormalization group (RG) partly by the KITP Graduate Fellows Program and the U.S. treatment, although it is not obvious how to set this up. DOE Office of Basic Energy Sciences, Division of Materials Such a RG treatment might also shed light on the nature of Sciences and Engineering under Award No. DE-SC0010526. the entanglement spectrum or, equivalently, the dependence J. H. is supported by the Pappalardo Fellowship in Physics of S ðtÞ on the index n. While we believe that all the Rényi at MIT. entropies execute KPZ growth in the presence of noise, we have not pinned down the n dependence of the various APPENDIX A: GROWTH OF HARTLEY constants. The solvable models suggest that the leading- ENTROPY S IN 1D order behavior may be independent of n at large times. What is the appropriate scaling form for the spectrum? Limited Consider a one-dimensional quantum spin chain of local time scales prevent us from addressing this numerically Hilbert space dimension q, prepared initially in a product (except for Clifford circuits, where all the S are trivially state, and apply a sequence of random unitaries that couple two neighboring spins. The location of the local unitary at a equal). (The entanglement spectrum is one window on the structure of the quantum states generated by the random given time step is arbitrary. In the following we fix the dynamics. We can also ask in what ways these states differ location of the unitary, but take it to be Haar random. We prove that in this situation the Hartley entropy S from ground states of random Hamiltonians, when the amount of entanglement is similar.) generically (i.e., with probability 1) obeys 031016-22 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) S ðx; t þ 1Þ¼ min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1; ðA1Þ reduces to seven different cases excluding invalid ones 0 0 0 and those related by reflection, one easily checks that there if a unitary is applied at the bond x. The logarithm is of is always a choice among I, W, E that makes Eq. (A1) true. base q. Let us treat three exemplary cases here. If s and s are L R This formula can be interpreted in matrix-product-state entangled at time t, then Sðx − 1;tÞ¼ Sðx þ 1;tÞ and language. If d is the minimal value of the local bond Sðx; tÞ¼ Sðx − 1;tÞþ 1, so one chooses the identity I. dimension required for an exact matrix-product-state rep- If s is entangled with a spin on the left of s and s L L R resentation of the state, then S ðxÞ¼ log d . A heuristic is entangled with a spin on the right of s , then 0 x R parameter-counting argument for the local bond dimension, Sðx − 1;tÞ¼ Sðx þ 1;tÞ¼ 1 þ Sðx; tÞ. One chooses the given in Appendix A3, suggests Eq. (A1). swap W to obtain Sðx; t þ 1Þ¼ Sðx; tÞþ 2.If s and s L R However, a more rigorous proof is necessary as such are both unentangled, then one applies the entangling heuristic arguments can fail. In particular, one might unitary E to obtain Sðx;t þ 1Þ− 1 ¼ Sðx;tÞ¼ Sðx− 1;tÞ¼ naively conjecture a stronger statement: namely, that for Sðx þ 1;tÞ. any state at time t, if the unitary at bond x is Haar random, For the second part, recall that for any bipartite state then Eq. (A1) is true with probability 1. This conjecture is jψi¼ M jiijji; ðA3Þ i;j false; a counterexample is given in Appendix A2.We now i;j give a proof of Eq. (A1). the number of nonzero Schmidt coefficients is equal to the number of nonzero singular values of the matrix M, which 1. Proof of Eq. (A1) is nothing but the rank of M. For any positive integer r, the Our genericity proof consists of two parts. First, we show rank of M is smaller than r if and only if every r × r that given locations of unitaries, there exist certain unitaries submatrix has determinant zero; i.e., all r × r minors such that at each time step Eq. (A1) is true. Second, we vanish. Thus, a bipartite state jψi having Hartley entropy show the negation of Eq. (A1), (log of rank of M) strictly smaller than log r is expressed by a system of polynomial equations on the coefficients of jψi. S ðx; t þ 1Þ < min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1; ðA2Þ 0 0 0 If jψi is given by U …U U j0i, where j0i is a fixed t 2 1 happens if and only if a system of polynomial equations in product state, then the coefficients are some polynomials of the entries of the unitaries is satisfied. (The inequality “>” the entries of the unitaries U , and hence the equations that never holds, as we note in the main text.) By the first part of expresses vanishing determinants are polynomial equations the proof, the zero locus of these polynomial equations in the entries of the unitaries. does not cover the entire set of unitaries. Therefore, it is Our claim Eq. (A1) completely determines the Hartley only a submanifold of strictly smaller dimension, which entropy based on the location of unitaries, and therefore the implies it has measure zero. spatial configuration of the unitaries tells us which minors we For the first part, it is sufficient to consider only three should check. Namely, the size r of the minors we turn into types of local unitaries: the identity I, the swap gate W, the polynomial equations is given by (the exponential of) and a unitary E with the property that it turns a pair of the right-hand side of Eq. (A2). In other words, given a spatial −1=2 unentangled polarized spins, j11i, into q jiii,a i¼1 configuration of unitaries, the polynomial equations that maximally entangled state. Without loss of generality we express Eq. (A2) are determined. The polynomial equations may take the initial product state to be the polarized are over tL variables, and the actual number of equations is state j…1111…i. much larger yet finite. We do not need explicit expressions We show that using these three types of unitaries at the for these polynomials, only the fact of their existence. These given locations, one can construct a state whose entangle- polynomials might a priori read 0 ¼ 0; i.e., they could be ment entropy is given by Eq. (A1). Since Eq. (A1) defines the trivially satisfied. In that case, the solution to the polynomial entropy inductively, we only have to show it inductively, too. equation would be the entire set of unitaries, and Eq. (A1) At t ¼ 0, all the spins are unentangled, so we can simply could never be satisfied. However, we just showed in the first choose E for every designated location. Clearly, Eq. (A1) is part that this cannot happen because there exists a choice of satisfied. At later times, if we do not apply E except on unitaries for which Eq. (A1) is satisfied. This implies that the an unentangled pair of spins, then a spin can be either polynomial equations are nontrivial and define a measure unentangled or maximally entangled with a single other zero subset of the entire set of unitaries. This completes the spin. Therefore, at time t> 0, the spin s that is immedi- genericity proof. ately left to the bond x can be (i) unentangled, (ii) entangled with a spin to the left of s , (iii) entangled with the spin s L R 2. Counterexample to the stronger conjecture that is immediately to the right of the bond y,or (iv) entangled with a spin to the right of s . We show above that Eq. (A1) holds when all unitaries are chosen generically and the initial state is a product state. These are exclusive possibilities, and similarly s has four options. Enumerating all 16 cases, which in fact Naively one might make the stronger conjecture: that the 031016-23 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) update rule Eq. (A1) holds whenever a generic unitary U is U and U , we find that all (5 × 5) minors vanish, implying 0 1 applied to an arbitrary—possibly fine-tuned—state jΨi. that V has rank at most 4. Therefore, for this nongeneric We construct an explicit jΨi, which is a counterexample to initial state, this stronger conjecture. Consider four degrees of freedom ABCD. The spins B S ðx; t þ 1Þ ≠ min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1: ðA9Þ 0 0 0 and C have dimension 2 each, and A and D have dimension 3 each. (To conform with our consideration of spin chains, 3. Parameter-counting argument the subsystems A and D should be regarded as subspaces of two or more spin-1=2’s.) The most general form of a Consider a 1D state jΨi in a matrix product representa- quantum state on ABCD is tion. Labeling the states of the qubits (spins) by σ; σ ; … running from 1 to q, 2 1 X X XX T jaijbijcijdi: ðA4Þ abcd σ 0σ 0 jΨi¼ ð…A A …Þj…σσ …i: ðA10Þ a ;a a ;a x−1 x x xþ1 a;d¼0 b;c¼0 fσg fag We consider T ¼ T δ , i.e., C is in j0i, where abcd abd c;0 Since the state is not translationally invariant, we allow 0 1 0 1 the bond dimension d to vary from bond to bond 01 0 00 1 (a ¼ 1; …;d ). In an efficient representation, d is equal x x x B C B C 0 0 T ¼ @10 0A ;T ¼ @00 0A : a0d a1d to the rank of the reduced density matrix for a cut at x: 00 0 10 0 ad ad S ðxÞ d ¼ q : ðA11Þ ðA5Þ We ask how S ðxÞ changes when we apply a unitary U to (This does not give a normalized state, but we are only 0 the two spins, σ and σ , either side of bond x. This effects concerned about ranks.) the change (repeated indices are summed): The Hartley entropy for the cut A=BCD is simple to compute. As we remark in the previous section, it is the 0 0 σ 0σ τ 0τ 0 0 A A → U A A : ðA12Þ rank of the coefficient matrix. Interpreting this matrix as a a ;a a ;a σσ ;ττ a ;a a ;a x−1 x x xþ1 x−1 x x xþ1 linear map, the rank is the dimension of the image of the map from BCD to A. The image is precisely the linear span To update the matrix product representation, we must find 0 0 of columns of T and T . They have three linearly ~ ~ new matrices A and A which satisfy a0d a1d independent columns, implying that the Hartley entropy for σ 0σ A=BCD is log 3. Similarly, the rank of the coefficient τ 0τ 2 ~ ~ A A ¼ U 0 0A A : ðA13Þ a ;a a ;a σσ ;ττ a ;a a ;a x−1 x x xþ1 x−1 x x xþ1 matrix for ABC=D is the dimension of the linear span of the 0 0 rows of T and T , which reads 3. That is, the Hartley a0d a1d 0 ~ ~ In order to solve this equation for A and A , it will generally entropy for ABC=D is log 3. be necessary to increase the bond dimension at x to a new If Eq. (A1) were to be true for the generic choice of Haar 0 0 value d . Naively, the necessary value of d will generically x x random unitary on BC, then we should be able to find a be determined by equating the number of independent unitary on BC such that equations in Eq. (A13) with the number of degrees of ~ ~ freedom in A and A . (However, the previous section shows S ðAB=CDÞ¼ log 3 þ 1 ¼ log 6: ðA6Þ 0 2 2 that this expectation can fail for certain choices of A and A .) We show this cannot hold. Applying the unitary U on BC The number of equations is q d d , since this is the x−1 xþ1 the state, we obtain number of possible values for the external indices in 0 0 Eq. (A13). The numbers of degrees of freedom in A and 0 0 0 0 0 0 U T ¼ U T þ U T ; ðA7Þ b c ;bc abcd b c ;00 b c ;10 a0d a1d |fflfflffl{zfflfflffl} |fflfflffl{zfflfflffl} 0 0 0 2 b;c A are qd d and qd d , respectively. However, d of x−1 x xþ1 x x U U 0 1 these are redundant, because the state is unchanged by the 0 0 σ σ 0σ 0σ −1 ~ ~ ~ ~ transformation A → A M, A → M A , with M an where U and U are 2 × 2 matrices. The coefficient matrix 0 1 0 0 arbitrary d × d matrix. Equating the number of equations for the cut AB=CD is then x x with the number of independent degrees of freedom gives 0 0 V ¼ U ⊗ T þ U ⊗ T ; ðA8Þ 0 1 0 1 0 0 ðd − qd Þðd − qd Þ¼ 0: ðA14Þ x x−1 x xþ1 whose rank should be 6 if S ðAB=CDÞ¼ log 6. Computing all the minors of the 6 × 6 matrix V for arbitrary matrices Choosing the smallest solution, 031016-24 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) d ¼ q minfd ;d g: ðA15Þ Z jψi¼þjψi; ðC2Þ x x−1 xþ1 i This agrees with Eq. (A1) since S ðxÞ¼ log d . for all i ¼ 1; …;n. The condition that jψi is nonzero and 0 x unique is equivalent to the condition that the operator APPENDIX B: HAAR AVERAGE FOR Trρ g ðC3Þ Let ρ ðtÞ be the reduced density matrix for a cut at x, jGj g∈G obtained by tracing out the spins to the left of the cut. Each index on this matrix labels a configuration of the spins to is a projector of rank 1 [96,97]. Since jψi is in the image of the right of the cut. Let us temporarily label these spins this projector, we see 1; 2; …, and let the spin immediately to the left of the cut be denoted 0. The indices on the reduced density matrices are then jψihψj¼ g: ðC4Þ jGj g∈G σ ;σ ;σ ;… 0 1 2 σ ;σ ;… σ ;… 1 2 2 ρ ðtÞ ; ρ ðtÞ ; ρ ðtÞ : ðB1Þ μ ;μ ;μ ;… x−1 0 1 2 x μ ;μ ;… xþ1 μ ;… 1 2 2 Since this is a normalized pure density matrix, its trace is In the following we assume that repeated indices are equal to 1. But a Pauli matrix has the property that it is summed. After applying a unitary on bond x, traceless. Therefore, only the identity element on the right has nonzero trace: 0 0 σ ;σ ;σ ;… σ ;σ ;…  0 1 1 2 ρ ðt þ 1Þ ¼ U 0 0 U ρ ðtÞ : 0 0 0 0 x μ ;μ ;… τσ ;σ σ x−1 1 τμ ;μ μ μ ;μ ;μ ;… 1 2 0 1 1 0 1 0 1 1 1 2 ⊗n n 1 ¼ dimðC Þ ¼ 2 : ðC5Þ Let us average Trρ ðt þ 1Þ over the choice of unitary, for a jGj jGj fixed initial state: From this expression, it is straightforward to obtain 0 0 00 00 σ ;σ ;σ ;… μ ;μ ;μ ;… 2 2 expressions for reduced density matrices. Suppose the 2 0 1 0 1 hTrρ ðt þ 1Þ i¼ ρ ðtÞ ρ ðtÞ 0 0 00 00 x x−1 x−1 μ ;μ ;μ ;… σ ;σ ;σ ;… 2 2 0 1 0 1 n-qubit system is partitioned into two complementary × hU 0 0 U U 00 00U i: 0 0 00 00 regions A and B. Tracing out B,wehave τσ ;σ σ νμ ;μ μ 1 τμ ;μ μ 1 νσ ;σ σ 0 1 1 0 1 1 0 1 0 1 The Haar average for four elements of a UðdÞ matrix (here, ρ ¼ Tr ðgÞ: ðC6Þ A B d ¼ q , and each index on U represents a pair of spin g∈G indices) is Tr ðgÞ is nonzero if and only if the tensor component hU U 0 0U U i 0 0 a;b a ;b c;d Haar corresponding to B is identity, in which case c ;d ¼ fδ δ 0 0δ δ 0 0 þ δ 0δ 0 δ 0δ 0 g jBj a;c a ;c b;d b ;d a;c a ;c b;d b ;d Tr ðgÞ¼ 2 gj ; ðC7Þ d − 1 where gj denotes the tensor components of g correspond- − fδ δ 0 0δ 0δ 0 þ δ 0δ 0 δ δ 0 0g : ðB2Þ a;c a ;c b;d b ;d a;c a ;c b;d b ;d ing to A. The set of all gj such that Tr ðgÞ ≠ 0 can be regarded as a subgroup of G, which we denote by G . The index contractions give the result in the text: The formula for ρ now reads 2 2 −1 2 2 hTrρ ðt þ 1Þ i ¼ qðq þ 1Þ ðTrρ þ Trρ Þ: x Haar jBj x−1 xþ1 X X 2 jG j 1 ρ ¼ g ¼ g: ðC8Þ jAj ðB3Þ 2 jG j g∈G g∈G A A APPENDIX C: ENTANGLEMENT ENTROPY It is immediate that ρ is proportional to a projector since it OF STABILIZER STATES is a sum over a group. It follows that the rank of ρ is equal jAj A stabilizer state is a state of an n-qubit system defined to 2 =jG j. In particular, the (Rényi or von Neumann) by a complete set fg ; …;g g of commuting tensor entropy of ρ with base-2 logarithm is 1 n A products of Pauli matrices through equations g jψi¼þjψi: ðC1Þ Sðρ Þ¼jAj − log jG j: ðC9Þ i A 2 A The group generated by fg ; …;g g is naturally called a 1 n The subgroup G has period 2, and therefore log jG j is an stabilizer group, and denoted by G [87,95]. A trivial A A example is the all-spin-up state, defined as integer, which is equal to the number of independent 031016-25 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) stabilizers supported only on A. This expression for the and (30). respectively. (The Hadamard gate corresponds to entanglement entropy has also appeared in Refs. [82,83]. swapping the X and Z vectors, while the phase gate Now, regard the stabilizer group G as a binary vector space corresponds to adding the X vector to the Z vector.) V by ignoring the overall phase (sign) factors. Let Π be the The von Neumann entropy in units of log 2 and the truncation map retaining the components corresponding to corresponding width averaged over ∼2 × 10 realizations the region A, and similarly let Π be the truncation map for (except for the last data point, where ∼2 × 10 realizations B ¼ A. It is routine to check that V decomposes as V ⊕ are used for the average) are plotted in Fig. 27. The fit to 0 0 V ⊕ V for some subspace V ⊆V, where V and V are the the KPZ universal form Eq. (49) gives β ¼ 0.2  0.15 B A B h spans of stabilizers supported only on A and B, respectively. and β ¼ 0.3  0.04. We also obtain v ¼ 0.194  0.001, w E Both the truncation maps are injective on V . It follows that B ¼ 0.4  0.2, C ¼ 0.4  0.1, D ¼ 0.4  0.6, and η ¼ S ¼jBj − dim V ¼ dim ðΠ VÞ − jAj. This completes −0.4  0.8. These results are consistent with the KPZ A F B F A 2 2 universality and with the data presented in Fig. 17. the proof of Eq. (33). APPENDIX E: DETAILS OF STATISTICS APPENDIX D: NUMERICS FOR FULL OF MEMBRANES CLIFFORD EVOLUTION The exponents governing the membrane problem are In Sec. VI A, we present numerical results for random traditionally denoted θ and ζ, and are related by 2ζ − θ ¼ unitary evolution using only the CNOT gates Eq. (31). Here, 2 − d [47]. Consider a patch of the membrane with linear we present similar analysis using the full set of generators for dimensions scaling as l. This includes both its temporal the Clifford group, showing that the additional gates do not dimension and its internal spatial dimensions: after a modify the universal behavior. The additional single-site rescaling of time, the membrane is statistically isotropic gates are the Hadamard and phase gates defined in Eqs. (29) on large scales. The mean “energy” of this patch of d θ membrane scales as l þ const × l , with fluctuations of order l . The length scale for wandering of the membrane in the transverse direction is of order l . The numerical results we quote in the main text are in good agreement with an epsilon expansion about d ¼ 4, which gives ζ ≃ 0.208ð4 − dÞ [93] (see also Ref. [98]). The scaling forms for the entanglement we discuss in the text are easily found by regarding the membrane as made up of patches of appropriate linear size: size t for Eqs. (57) and (58) and size L for Eq. (59). Note that the geometry of the membrane, including the transverse length scale (which is Δx ∼ t for the regime t ≲ L), determines the dimensions of the space-time region around ∂A for which the final entanglement is sensitive to small changes in HðtÞ, i.e., in the history of the noise. APPENDIX F: ENTANGLEMENT PROBABILITY DISTRIBUTION As we mention in the main text, a remarkable recent advance in KPZ theory has been the derivation of the full universal probability distribution for the height of the surface at fixed position and fixed large time [56–66]; see Refs. [99–101] for reviews. In our case, this height corresponds to the entanglement S across a cut in a system undergoing noisy unitary dynamics. One may separate out FIG. 27. Top: Growth of the mean entanglement in units of log the nonuniversal growth rate v , and the nonuniversal 2 as a function of time for the random Clifford evolution (only constant D governing the scale of fluctuations, by writing CNOT gates). The red solid curve is a fit using the form Eq. (49). Dashed line shows asymptotic linear behavior. Bottom: Growth S ¼ v t þ Dt χ: ðF1Þ in the fluctuations in the entanglement with time. The exponent β E is found to be β ¼ 0.3  0.04, in agreement with the KPZ prediction β ¼ 1=3. The dashed line shows the expected asymp- The rescaled random variable χ is then expected to have totic behavior, wðtÞ ∼ t with β ¼ 1=3. a universal probability distribution PðχÞ at late times. 031016-26 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) states (Sec. IX) indicates that the lower end point of the polymer is then no longer free, and instead favors x ¼ 0.] In Fig. 28, we fit numerical data for the probability distribution of S to the expected Tracy-Widom form and, for comparison, to a Gaussian distribution. The data are for the “full” Clifford dynamics (defined in Sec. VI A) at time t ¼ 2048. Each fit involves two parameters, corresponding to the mean and the variance. The Tracy-Widom distribu- tion used is the theoretically expected one with β ¼ 1,but, in fact, the present data do not allow us to discriminate between TW and TW . The Tracy-Widom distribu- β¼1 β¼2 tion is a much better fit to the data than the Gaussian, as quantified in the caption of Fig. 28. This is further confirmation of KPZ universality in the Clifford case. A more detailed investigation of the probability distribution is beyond the scope of this paper, in view of finite-time effects at the accessible time scales. [1] P. Calabrese and J. Cardy, Evolution of Entanglement Entropy in One-Dimensional Systems, J. Stat. Mech. (2005) P04010. [2] P. Calabrese and J. Cardy, Entanglement Entropy and Conformal Field Theory, J. Phys. A 42, 504005 (2009). [3] P. Calabrese and J. Cardy, Quantum Quenches in 1 þ 1- Dimensional Conformal Field Theories, J. Stat. Mech. (2016) 064003. [4] C. T. Asplund, A. Bernamonti, F. Galli, and T. Hartman, FIG. 28. Observed probability distribution for the entanglement Entanglement Scrambling in 2D Conformal Field Theory, entropy across the central bond of chain of length L ¼ 2048 at J. High Energy Phys. 09 (2015) 110. time t ¼ 2048 under the full Clifford dynamics, fitted to two [5] V. Hubeny, M. Rangamani, and T. Takayanagi, A Covar- probability distributions. Top: Best fit to the Tracy-Widom (TW) iant Holographic Entanglement Entropy Proposal, J. High distribution with β ¼ 1. A fit to TW with β ¼ 2 is not shown, but Energy Phys. 07 (2007) 062. is indistinguishable at the scale of the figure. Bottom: Best fit to [6] J. Abajo-Arrastia, J. Aparício, and E. López, Holographic the Gaussian. Clearly, the Tracy-Widom distribution fits the data Evolution of Entanglement Entropy, J. High Energy Phys. better than Gaussian, as the latter shows systematic deviation. 11 (2010) 149. 2 −4 The 1 − R values for the fits are 2.1 × 10 for TW , [7] T. Hartman and J. Maldacena, Time Evolution of Entan- β¼1 −4 −3 glement Entropy from Black Hole Interiors, J. High Energy 2.0 × 10 for TW , and 1.6 × 10 for Gaussian. β¼2 Phys. 05 (2013) 014. [8] H. Liu and S. J. Suh, Entanglement Tsunami: Universal This probability distribution depends on the initial con- Scaling in Holographic Thermalization, Phys. Rev. Lett. dition for the surface. For a surface that is initially flat, PðχÞ 112, 011601 (2014). is the Tracy-Widom distribution with β ¼ 1. This is the [9] H. Liu and S. J. Suh, Entanglement Growth during case relevant to our setup, where Sðx; t ¼ 0Þ¼ 0. (In the Thermalization in Holographic Systems, Phys. Rev. D directed polymer interpretation, this corresponds to a setup 89, 066012 (2014). where the x–coordinate of the upper end point of the [10] H. Casini, H. Liu, and M. Mezei, Spread of Entanglement polymer is fixed but that of the lower end point is free; and Causality, J. High Energy Phys. 07 (2016) 077. [11] M. Alishahiha, A. F. Astaneh, and M. R. Mohammadi again, this is the setup relevant to our minimal cut picture.) Mozaffar, Thermalization in Backgrounds with Hyper- Other initial conditions for a growing surface can give scaling Violating Factor, Phys. Rev. D 90, 046004 (2014). different universal forms for PðχÞ—for example, the so- [12] S. Kundu and J. F. Pedraza, Spread of Entanglement for called “narrow wedge” initial condition gives the Tracy- Small Subsystems in Holographic CFTs, Phys. Rev. D 95, Widom distribution with β ¼ 2. (The latter distribution is 086008 (2017). likely to be relevant to noisy growth of entanglement [13] M. Fagotti and P. Calabrese, Evolution of Entanglement between two subsystems that are initially unentangled with Entropy Following a Quantum Quench: Analytic Results each other, but separately highly entangled.) [The gener- for the XY Chain in a Transverse Magnetic Field, Phys. alization of the directed polymer picture to entangled initial Rev. A 78, 010306 (2008). 031016-27 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) [14] I. Peschel and V. Eisler, Reduced Density Matrices and in a Quantum Many-Body System, Nature (London) 528, Entanglement Entropy in Free Lattice Models, J. Phys. A 77 (2015). [33] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. 42, 504003 (2009). Schittko, P. M. Preiss, and M. Greiner, Quantum Thermal- [15] J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. ization through Entanglement in an Isolated Many-Body Daley, Entanglement Growth in Quench Dynamics with System, Science 353, 794 (2016). Variable Range Interactions, Phys. Rev. X 3, 031015 [34] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, (2013). Measuring Entanglement Growth in Quench Dynamics of [16] P. Hauke and L. Tagliacozzo, Spread of Correlations in Bosons in an Optical Lattice, Phys. Rev. Lett. 109, 020505 Long-Range Interacting Quantum Systems, Phys. Rev. (2012). Lett. 111, 207202 (2013). [35] P. Hayden and J. Preskill, Black Holes as Mirrors: [17] A. Coser, E. Tonni, and P. Calabrese, Entanglement Quantum Information in Random Subsystems, J. High Negativity after a Global Quantum Quench, J. Stat. Mech. Energy Phys. 09 (2007) 120. (2014) P12017. [36] Y. Sekino and L. Susskind, Fast Scramblers, J. High [18] A. S. Buyskikh, M. Fagotti, J. Schachenmayer, F. Essler, Energy Phys. 10 (2008) 065. and A. J. Daley, Entanglement Growth and Correlation [37] C. Dankert, R. Cleve, J. Emerson, and E. Livine, Exact and Spreading with Variable-Range Interactions in Spin and Approximate Unitary 2-Designs and Their Application to Fermionic Tunneling Models, Phys. Rev. A 93, 053620 Fidelity Estimation, Phys. Rev. A 80, 012304 (2009). (2016). [38] D. A. Roberts, D. Stanford, and L. Susskind, Localized [19] V. Alba and P. Calabrese, Entanglement and Thermody- Shocks, J. High Energy Phys. 03 (2015) 051. namics after a Quantum Quench in Integrable Systems, [39] R. Oliveira, O. C. O. Dahlsten, and M. B. Plenio, Generic arXiv:1608.00614. Entanglement Can be Generated Efficiently, Phys. Rev. [20] A. M. Lauchli and C. Kollath, Spreading of Correlations Lett. 98, 130502 (2007). and Entanglement after a Quench in the One-Dimensional [40] O. C. O. Dahlsten, R. Oliveira, and M. B. Plenio, The Bose Hubbard Model, J. Stat. Mech. (2008) P05018. Emergence of Typical Entanglement in Two-Party Ran- [21] H. Kim and D. A. Huse, Ballistic Spreading of Entangle- dom Processes, J. Phys. A 40, 8081 (2007). ment in a Diffusive Nonintegrable System, Phys. Rev. Lett. [41] M. Žnidarič, Exact Convergence Times for Generation of 111, 127205 (2013). Random Bipartite Entanglement, Phys. Rev. A 78, 032324 [22] L. Zhang, H. Kim, and D. A. Huse, Thermalization of (2008). Entanglement, Phys. Rev. E 91, 062128 (2015). [42] F. G. S. L. Brandao, A. W. Harrow, and M. Horodecki, [23] W. W. Ho and D. A. Abanin, Entanglement Dynamics in Local Random Quantum Circuits Are Approximate Poly- Quantum Many-Body Systems, Phys. Rev. B 95, 094302 nomial-Designs, Commun. Math. Phys. 346, 397 (2016). (2017). [43] J. Gütschow, S. Uphoff, R. F. Werner, and Z. Zimborás, [24] M. Žnidarič,T. ž. Prosen, and P. Prelovšek, Many-Body Time Asymptotics and Entanglement Generation of Clif- Localization in the Heisenberg XXZ Magnet in a Random ford Quantum Cellular Automata, J. Math. Phys. 51, Field, Phys. Rev. B 77, 064426 (2008). 015203 (2010). [25] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded [44] S. Leichenauer and M. Moosa, Entanglement Tsunami in Growth of Entanglement in Models of Many-Body (1 þ 1)-Dimensions, Phys. Rev. D 92, 126004 (2015). Localization, Phys. Rev. Lett. 109, 017202 (2012). [45] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic Scaling of [26] M. Serbyn, Z. Papić, and D. A. Abanin, Universal Slow Growing Interfaces, Phys. Rev. Lett. 56, 889 (1986). Growth of Entanglement in Interacting Strongly Disor- [46] D. A. Huse, C. L. Henley, and D. S. Fisher, Respond, Phys. dered Systems, Phys. Rev. Lett. 110, 260601 (2013). Rev. Lett. 55, 2924 (1985). [27] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenom- [47] D. A. Huse and C. L. Henley, Pinning and Roughening of enology of Fully Many-Body-Localized Systems, Phys. Domain Walls in Ising Systems due to Random Impurities, Rev. B 90, 174202 (2014). Phys. Rev. Lett. 54, 2708 (1985). [28] R. Vosk and E. Altman, Many-Body Localization in One [48] D. Forster, D. R. Nelson, and M. J. Stephen, Large- Dimension as a Dynamical Renormalization Group Fixed Distance and Long-Time Properties of a Randomly Stirred Point, Phys. Rev. Lett. 110, 067204 (2013). Fluid, Phys. Rev. A 16, 732 (1977). [29] A. Chandran and C. R. Laumann, Semiclassical Limit for [49] T. Halpin-Healy, Directed Polymers versus Directed the Many-Body Localization Transition, Phys. Rev. B 92, Percolation, Phys. Rev. E 58, R4096 (1998). 024301 (2015). [50] S. Ryu and T. Takayanagi, Holographic Derivation of [30] D. J. Luitz, N. Laflorencie, and F. Alet, Extended Slow Entanglement Entropy from the Anti–de Sitter Space/ Dynamical Regime Close to the Many-Body Localization Conformal Field Theory Correspondence, Phys. Rev. Lett. Transition, Phys. Rev. B 93, 060201 (2016). 96, 181602 (2006). [31] F. Verstraete, V. Murg, and J. I. Cirac, Matrix Product [51] B. Swingle, Entanglement Renormalization and Hologra- States, Projected Entangled Pair States, and Variational phy, Phys. Rev. D 86, 065007 (2012). Renormalization Group Methods for Quantum Spin Sys- [52] F. Pastawski, B. Yoshida, D. Harlow, and J. Preskill, tems, Adv. Phys. 57, 143 (2008). Holographic Quantum Error-Correcting Codes: Toy [32] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Models for the Bulk/Boundary Correspondence, J. High Rispoli, and M. Greiner, Measuring Entanglement Entropy Energy Phys. 06 (2015) 149. 031016-28 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) [53] P. Hayden, S. Nezami, X.-L. Qi, N. Thomas, M. Walter, [73] A. Nahum, S. Vijay, and J. Haah “Operator Spreading in and Z. Yang, Holographic Duality from Random Tensor Random Unitary Circuits” (to be published). Networks, J. High Energy Phys. 11 (2016) 009. [74] G. Evenbly and G. Vidal, Tensor Network States and [54] K. A. Takeuchi and M. Sano, Evidence for Geometry- Geometry, J. Stat. Phys. 145, 891 (2011). Dependent Universal Fluctuations of the Kardar-Parisi- [75] M. Kardar, Roughening by Impurities at Finite Temper- Zhang Interfaces in Liquid-Crystal Turbulence, J. Stat. atures, Phys. Rev. Lett. 55, 2923 (1985). [76] M Kardar, Statistical Physics of Fields (Cambridge Phys. 147, 853 (2012). University Press, Cambridge, England, 2007). [55] K. A. Takeuchi, M. Sano, T. Sasamoto, and H. Spohn, [77] A. Nahum, J. Ruhman, and D. Huse, “Dynamics of Growing Interfaces Uncover Universal Fluctuations Behind Scale Invariance, Sci. Rep. 1, 34 (2011). entanglement and transport in 1D systems with quenched [56] P. Calabrese, P. Le Doussal, and A. Rosso, Free-Energy randomness” (to be published). Distribution of the Directed Polymer at High Temperature, [78] E. H. Lieb and D. W. Robinson, in Statistical Mechanics: Selecta of Elliott H. Lieb, edited by B. Nachtergaele, Europhys. Lett. 90, 20002 (2010). J. P. Solovej, and J. Yngvason (Springer, Berlin, 2004), [57] V. Dotsenko, Bethe Ansatz Derivation of the Tracy-Widom Distribution for One-Dimensional Directed Polymers, pp. 425–431. Europhys. Lett. 90, 20003 (2010). [79] D. A. Roberts and B. Swingle, Lieb-Robinson Bound and [58] T. Sasamoto and H. Spohn, One-Dimensional Kardar- the Butterfly Effect in Quantum Field Theories, Phys. Rev. Parisi-Zhang Equation: An Exact Solution and Its Uni- Lett. 117, 091602 (2016). versality, Phys. Rev. Lett. 104, 230602 (2010). [80] D. Gottesman, The Heisenberg Representation of [59] T. Sasamoto and H. Spohn, Exact Height Distributions for Quantum Computers, arXiv:quant-ph/9807006. [81] We truncate all stabilizers to region A by throwing away all the KPZ Equation with Narrow Wedge Initial Condition, Nucl. Phys. B834, 523 (2010). the spin operators acting outsideA. In this process some of the [60] T. Sasamoto and H. Spohn, The Crossover Regime for the stabilizers become trivial, and some become redundant, i.e., Weakly Asymmetric Simple Exclusion Process, J. Stat. equal to products of the others. I is the number of stabilizers Phys. 140, 209 (2010). that are independent when truncated to A. Equivalently, I is [61] G. Amir, I. Corwin, and J. Quastel, Probability Distribu- the rank of the matrix Ψ after the rows corresponding to the tion of the Free Energy of the Continuum Directed complement of A have been deleted; this is how we calculate Random Polymer in 1 þ 1 Dimensions, Commun. Pure the entropy numerically for Sec. VI A. Appl. Math. 64, 466 (2011). [82] A. Hamma, R. Ionicioiu, and P. Zanardi, Bipartite [62] P. Calabrese and P. Le Doussal, Exact Solution for the Entanglement and Entropic Boundary Law in Lattice Spin Kardar-Parisi-Zhang Equation with Flat Initial Condi- Systems, Phys. Rev. A 71, 022315 (2005). tions, Phys. Rev. Lett. 106, 250603 (2011). [83] A. Hamma, R. Ionicioiua, and P. Zanardi, Ground State [63] S. Prolhac and H. Spohn, Height Distribution of the Entanglement and Geometric Entropy in the Kitaev Model, Kardar-Parisi-Zhang Equation with Sharp-Wedge Initial Phys. Lett. A 337, 22 (2005). Condition: Numerical Evaluations, Phys. Rev. E 84, [84] T. Halpin-Healy and Y.-C. Zhang, Kinetic Roughening 011119 (2011). Phenomena, Stochastic Growth, Directed Polymers and [64] P. Le Doussal and P. Calabrese, The KPZ Equation with All That. Aspects of Multidisciplinary Statistical Mechan- Flat Initial Condition and the Directed Polymer with One ics, Phys. Rep. 254, 215 (1995). Free End, J. Stat. Mech. (2012) P06001. [85] O. C. O. Dahlsten and M. B. Plenio, Exact Entanglement [65] T. Imamura and T. Sasamoto, Exact Solution for the Probability Distribution of Bi-Partite Randomised Stabi- Stationary Kardar-Parisi-Zhang Equation, Phys. Rev. lizer States, Quantum Inf. Comput. 6, 527 (2006). Lett. 108, 190603 (2012). [86] ITensor, http://itensor.org. [66] T. Imamura and T. Sasamoto, Stationary Correlations for [87] D. Gottesman, A Class of Quantum Error-Correcting the 1D KPZ Equation, J. Stat. Phys. 150, 908 (2013). Codes Saturating the Quantum Hamming Bound, Phys. [67] D. N. Page, Average Entropy of a Subsystem, Phys. Rev. Rev. A 54, 1862 (1996). [88] G. Roósz, R. Juhász, and F. Iglói, Nonequilibrium Dy- Lett. 71, 1291 (1993). [68] C. Nadal, S. N. Majumdar, and M. Vergassola, Statistical namics of the Ising Chain in a Fluctuating Transverse Distribution of Quantum Entanglement for a Random Field, Phys. Rev. B 93, 134305 (2016). Bipartite State, J. Stat. Phys. 142, 403 (2011). [89] L. Saul, M. Kardar, and N. Read, Directed Waves in [69] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Random Media, Phys. Rev. A 45, 8859 (1992). Matrix Product State Representations, Quantum Inf. [90] I. Peschel, Calculation of Reduced Density Matrices Comput. 7, 401 (2007). from Correlation Functions, J. Phys. A 36, L205 [70] C. Chamon and E. R. Mucciolo, Virtual Parallel Comput- (2003). ing and a Search Algorithm Using Matrix Product States, [91] T. Nattermann, Ising Domain Wall in a Random Pinning Phys. Rev. Lett. 109, 030503 (2012). Potential, J. Phys. C 18, 6661 (1985). [71] P. Meakin, P. Ramanlal, L. M. Sander, and R. C. Ball, [92] M. Kardar, Domain Walls Subject to Quenched Impurities, Ballistic Deposition on Surfaces, Phys. Rev. A 34, 5091 J. Appl. Phys. 61, 3601 (1987). (1986). [93] D. S. Fisher, Interface Fluctuations in Disordered Systems: [72] J. M. Kim and J. M. Kosterlitz, Growth in a Restricted 5 − ϵ Expansion and Failure of Dimensional Reduction, Solid-on-Solid Model, Phys. Rev. Lett. 62, 2289 (1989). Phys. Rev. Lett. 56, 1964 (1986). 031016-29 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) [94] A. A. Middleton, Numerical Results for the Ground-State (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Interface in a Random Medium, Phys. Rev. E 52, R3337 Dagstuhl, (1995). Germany, 2013), Vol. 22, pp. 270–284. [95] A. R. Calderbank, E. M Rains, P. W. Shor, and N. J. A. [98] T. Halpin-Healy, Disorder-Induced Roughening of Diverse Sloane, Quantum Error Correction and Orthogonal Manifolds, Phys. Rev. A 42, 711 (1990). Geometry, Phys. Rev. Lett. 78, 405 (1997). [99] T. Kriecherbauer and J. Krug, A Pedestrian’s View on [96] A. Klappenecker and M. Rotteler, Beyond Stabilizer Interacting Particle Systems, KPZ Universality and Ran- Codes II: Clifford Codes, IEEE Trans. Inf. Theory 48, dom Matrices, J. Phys. A 43, 403001 (2010). 2396 (2002). [100] I. Corwin, The Kardar-Parisi-Zhang Equation and [97] N. Linden, F. Matus, M. B. Ruskai, and A. Winter, in Universality Class, Random Matrices Theory Appl. 01, Proceedings of the 8th Conference on the Theory of 1130001 (2012). Quantum Computation, Communication and Cryptogra- [101] T. Halpin-Healy and K. A. Takeuchi, A KPZ Cocktail- phy (TQC 2013), edited by S. Severini and F. Brandao Shaken, Not Stirred…, J. Stat. Phys. 160, 794 (2015). 031016-30 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Physical Review X American Physical Society (APS)

Quantum Entanglement Growth under Random Unitary Dynamics

Free
30 pages

Loading next page...
 
/lp/aps_physical/quantum-entanglement-growth-under-random-unitary-dynamics-TBawes03nl
Publisher
The American Physical Society
Copyright
Copyright © Published by the American Physical Society
eISSN
2160-3308
D.O.I.
10.1103/PhysRevX.7.031016
Publisher site
See Article on Publisher Site

Abstract

PHYSICAL REVIEW X 7, 031016 (2017) 1,3 1 1,2 1 Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Kavli Institute for Theoretical Physics, University of California Santa Barbara, Santa Barbara, California 93106, USA Theoretical Physics, Oxford University, 1 Keble Road, Oxford OX1 3NP, United Kingdom (Received 1 September 2016; revised manuscript received 16 April 2017; published 24 July 2017) Characterizing how entanglement grows with time in a many-body system, for example, after a quantum quench, is a key problem in nonequilibrium quantum physics. We study this problem for the case of random unitary dynamics, representing either Hamiltonian evolution with time-dependent noise or evolution by a random quantum circuit. Our results reveal a universal structure behind noisy entanglement growth, and also provide simple new heuristics for the “entanglement tsunami” in Hamiltonian systems without noise. In 1D, we show that noise causes the entanglement entropy across a cut to grow according to the celebrated Kardar- Parisi-Zhang (KPZ) equation. The mean entanglement grows linearly in time, while fluctuations grow like 1=3 2=3 ðtimeÞ and are spatially correlated over a distance ∝ ðtimeÞ . We derive KPZ universal behavior in three complementary ways, by mapping random entanglement growth to (i) a stochastic model of a growing surface, (ii) a “minimal cut” picture, reminiscent of the Ryu-Takayanagi formula in holography, and (iii) a hydrodynamic problem involving the dynamical spreading of operators. We demonstrate KPZ universality in 1D numerically using simulations of random unitary circuits. Importantly, the leading-order time dependence of the entropy is deterministic even in the presence of noise, allowing us to propose a simple coarse grained minimal cut picture for the entanglement growth of generic Hamiltonians, even without noise, in arbitrary dimensionality. We clarify the meaning of the “velocity” of entanglement growth in the 1D entanglement tsunami. We show that in higher dimensions, noisy entanglement evolution maps to the well- studied problem of pinning of a membrane or domain wall by disorder. DOI: 10.1103/PhysRevX.7.031016 Subject Areas: Condensed Matter Physics, Quantum Information, Statistical Physics I. INTRODUCTION goes on. This irreversible growth of entanglement— quantified by the growth of the von Neumman entropy— The language of quantum entanglement ties together is important for several reasons. It is an essential part condensed matter physics, quantum information, and high- of thermalization, and as a result has been addressed in energy theory. The von Neumann entanglement entropy is diverse contexts ranging from conformal field theory known to encode universal properties of quantum ground [1–4] and holography [5–12] to integrable [13–19], non- states and has led to new perpectives on the AdS=CFT integrable [20–23], and strongly disordered spin chains correspondence. But the dynamics of the entanglement are [24–30]. Entanglement growth is also of practical impor- far less understood. The entanglement entropy is a highly tance as the crucial obstacle to simulating quantum nonlocal quantity, with very different dynamics to energy dynamics numerically, for example, using matrix product or charge or other local densities. Traditional many-body states or the density matrix renormalization group [31]. tools therefore do not provide much intuition about how The entanglement entropy, and even its time dependence, entanglement spreads with time, for example, after a is also beginning to be experimentally measurable in quantum quench (a sudden change to the Hamiltonian). cold atom systems [32–34]. In a very different context, We need to develop simple heuristic pictures, and simple black holes have motivated studies of how fast quantum long-wavelength descriptions, for entanglement dynamics. systems can scramble information by dynamically gen- If a many-body system is initialized in a state with erating entanglement [35–38]. Simple quantum circuits— low entanglement, the dynamics will typically generate quantum evolutions in discrete time—serve as useful toy entanglement between increasingly distant regions as time models for entanglement growth and scrambling [39–43]. For integrable 1D systems and rational CFTs, there is an appealing heuristic picture for entanglement growth fol- Published by the American Physical Society under the terms of lowing a quench in terms of spreading quasiparticles [1,3]. the Creative Commons Attribution 3.0 License. Further distri- However, this picture does not apply to general interacting bution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. systems [3,4,10,44]. 2160-3308=17=7(3)=031016(30) 031016-1 Published by the American Physical Society NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) In this paper, we propose new heuristic pictures for KPZ universality class entanglement growth in generic, nonintegrable systems, both in 1D and in higher dimensions. We arrive at these pictures by studying “minimally structured” models for quantum Classical Directed polymer in Stochastic dynamics: dynamics that are spatially local, and unitary, but surface growth random medium particle dynamics random both in time and space (“noisy”). Concretely, we focus on quantum circuit dynamics with randomly chosen quantum gates. Entanglement growth in these systems exhibits a remarkable universal structure in its own right, related to paradigmatic problems in classical statistical Growth of Hydrodynamics mechanics. But in addition, random circuits provide a entanglement through circuit of operator spreading theoretical laboratory that allows us to derive scaling pictures FIG. 1. The KPZ “triumvirate” is made up of three very for entanglement growth and the so-called “entanglement different problems in classical statistical mechanics which all tsunami” [8], which, we conjecture, generalize to quenches map to the KPZ universality class. As we discuss, each of them in many-body systems without noise. For example, we can be usefully related to entanglement in 1 þ 1D. propose a simple “minimal membrane” picture that can be used to derive scaling forms for the growth of the entangle- ment. We also argue that generically there is a well-defined summarized in Fig. 1. We show that entanglement growth “entanglement speed” v , but this is generically smaller than E can usefully be related to all three of the classical problems in the “butterfly speed” v governing operator growth, and we B in Fig. 1. give a physical explanation for this phenomenon. In the quantum setting, the directed polymer is related to We show that noisy entanglement growth allows a long- the “minimal cut,” a curve in space-time that bisects the wavelength description with an emergent universal structure. unitary circuit representing the time evolution. This picture Physically, the class of noisy dynamics includes closed, is more general than the surface growth picture, as it allows many-body systems whose Hamiltonian HðtÞ contains one to consider the entropy for any bipartition of the system. fluctuating noise terms, and also quantum circuits in which It also allows us to generalize from 1D to higher dimensions. qubits are acted on by randomly chosen unitary gates. In this The picture is reminiscent of the Ryu-Takayanagi prescrip- setting, we pin down both the leading-order deterministic tion for calculating the entanglement entropy of conformal behavior of the entanglement and the subleading fluctuations field theories in the AdS=CFT correspondence, which associated with noise. We argue that fluctuations and spatial makes use of a minimal surface in the bulk space [50], correlations in the entanglement entropy are characterized and analogous results for certain tensor network states by universal scaling exponents, expected to be independent [51–53]. Here, however, the cut lives in space-time rather of the details of the microscopic model. than in space, and in a noisy system its shape is random For noisy systems in one spatial dimension, we argue that rather than deterministic. (For a different use of the idea of a the critical exponents for entanglement growth are those of minimal cut in space-time, see Ref. [10].) In d þ 1 space- the Kardar-Parisi-Zhang (KPZ) equation, originally intro- time dimensions the minimal cut becomes a d-dimensional duced to describe the stochastic growth of a surface with time membrane pinned by disorder. This picture allows us to t [45]. In the simplest setting, we find that the height of this obtain approximate critical exponents for noisy entangle- ment growth in any number of dimensions. surface at a point x in space is simply the von Neumann This picture also leads us to a conjecture for entangle- entanglement entropy Sðx; tÞ for a bipartition which splits the system in two at x. The average entanglement grows linearly ment growth in systems without noise, both in 1D and in time, while fluctuations are characterized by nontrivial higher dimensions, as we discuss below. According to this conjecture, the calculation of the entanglement in higher exponents. We support this identification with analytical dimensions reduces to a deterministic elastic problem for arguments and numerical results for discrete time quantum the minimal membrane in space-time. In 1D, it results in evolution (unitary circuits). particularly simple universal scaling functions, which agree A remarkable feature of the KPZ universality class is that with scaling forms in holographic 1 þ 1D CFTs [8,10,44], it also embraces two classical problems that at first sight are very different from surface growth [45,46]. These connec- and which we suggest are universal for generic, non- tions lead us to powerful heuristic pictures for entanglement integrable, translationally invariant 1D systems. growth, in both 1D and higher dimensions. The KPZ The third member of the triumvirate in Fig. 1 is a noisy universality class embraces the statistical mechanics of a hydrodynamic equation describing the diffusion of inter- directed polymer in a disordered potential landscape [47] and acting (classical) particles in 1D. We show that this can be 1D hydrodynamics with noise (the noisy Burgers equation related to the spreading of quantum operators under the [48]). These problems, together with surface growth, are unitary evolution, giving a detailed treatment of the special sometimes known as the “KPZ triumvirate” [49]. They are case of stabilizer circuits. Note that while the minimal cut 031016-2 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) picture generalizes to higher dimensions, the KPZ and hydrodynamic pictures are special to 1D. A. Organization of the paper We propose that noisy dynamics are a useful toy model for quantum quenches in generic (nonintegrable, noncon- FIG. 2. Spin chain with open boundary conditions. SðxÞ formally–invariant) systems, even without noise. The logic denotes the entanglement entropy (von Neumann or Rényi of our approach is to pin down the universal behavior of depending on context) between the part of the chain to the left noisy systems (Secs. II–VI), to establish simple heuristics of bond x, indicated by the box, and the part to the right of bond x. capturing this behavior (Secs. III and IV), and then to extend these heuristic pictures to dynamics without noise by splitting the chain into two halves at x and tracing out (Secs. V and VIII). the left-hand side (Fig. 2). The nth Rényi entropy for a cut The detailed physics of the entanglement fluctuations at x is defined as (including KPZ exponents) certainly relies on noise. However, the coarser features of the dynamics are in fact S ðxÞ¼ log ðTrρ Þ: ð1Þ deterministic. These include the leading-order time depend- n x 1 − n ence of the entanglement entropy and mutual information when the length and time scales are large. We conjecture Logarithms are taken base q. In the limit n → 1, the Rényi that this leading-order behavior, as captured by the directed entropy becomes the von Neumann entropy: polymer and hydrodynamic pictures, carries over to Hamiltonian dynamics without noise. On the basis of this S ðxÞ¼ −Trρ log ρ : ð2Þ vN x x we address (Sec. V) features of entanglement growth that have previously been unclear. We argue that in generic 1D A basic constraint on the von Neumann entropy is that systems the entanglement growth rate can be interpreted as neighboring bonds can differ by at most one (this follows a well-defined speed v , but that this speed is generically from subadditivity of the von Neumann entropy): smaller than another characteristic speed, which is the speed v at which quantum operators spread out under the jS ðx þ 1Þ − S ðxÞj ≤ 1: ð3Þ vN vN dynamics (the butterfly speed). Section V also addresses universal scaling forms for the entanglement entropy in 1D. In this section, we focus on the growth of the bipartite In Sec. VIII, we discuss the geometry dependence of the entropies Sðx; tÞ with time, starting from a state with low dynamical entanglement in higher-dimensional systems: entanglement. [Here Sðx; tÞ, without a subscript, can denote we argue that there is again a scaling picture in terms of a any of the Rényi entropies with n> 0.] For simplicity, we minimal surface, but that more nonuniversal parameters take the initial state to be a product state, but we expect the enter into the time dependence than in 1+1D. same long-time behavior for any initial state with area-law entanglement. (The setup with area-law entanglement in the initial state is analogous to a quantum quench that starts II. SURFACE GROWTH IN 1D in the ground state of a noncritical Hamiltonian. We briefly We begin by studying entanglement growth under consider initial states with non-area-law entanglement in random unitary dynamics in one dimension. After sum- Sec. IX.) We argue that for noisy unitary dynamics, the marizing the KPZ universal behavior, we derive this universal properties of Sðx; tÞ are those of the Kardar- behavior analytically in a solvable model, using a mapping Parisi-Zhang equation: to a classical surface growth problem. In the following sections we provide alternative derivations of this universal ∂S λ 2 2 ¼ ν∂ S − ð∂ SÞ þ ηðx; tÞþ c: ð4Þ x x behavior by relating the minimal cut bound on the ∂t 2 entanglement to the classical problem of a directed polymer in a random environment, and by relating the spreading of This equation was introduced to describe the stochastic quantum operators to a 1D hydrodynamics problem. growth of a 1D surface or interface with height profile SðxÞ Consider a chain of quantum spins with local Hilbert [45]. It captures an important universality class which has space dimension q (for example, spin-1=2’s with q ¼ 2). found a wealth of applications in classical nonequilibrium We take open boundary conditions and label the bonds of physics, and its scaling properties have been verified the lattice by x ¼ 1; …;L. We consider only unitary in high-precision experiments [54,55]. The constant c in dynamics, so the full density matrix ρ ¼jΨihΨj represents Eq. (4) contributes to the positive average growth rate, a pure state. For now we consider the entanglement across a while ηðx; tÞ is noise which is uncorrelated in space and single cut at position x; we generalize to other geometries time. The ν term describes diffusive smoothing of sharp later in the paper. The reduced density matrix ρ is defined features. The nonlinear term, with coefficient λ, describes 031016-3 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) how the average growth rate depends on the slope; the t ≃ ; ð10Þ saturation negative sign is natural here, as we discuss in Sec. II A 3 2v [and implies that B in Eq. (6) below is positive]. with bonds closer to the boundary saturating sooner (Secs. III KPZ scaling is characterized by an exponent β governing and V). Saturation is also captured in the surface growth the size of fluctuations in the interface height, an exponent description, once we note that there are Dirichlet boundary α governing spatial correlations, and a dynamical exponent conditions on the entropy: Sð0;tÞ¼ SðL þ 1;tÞ¼ 0. z determining the rate of growth of the correlation length Note that the scaling described in Eqs. (6) and (8) implies (z ¼ α=β by a scaling relation). These are known exactly the existence of two distinct diverging length scales during [45]: entanglement growth. The fact that hSðx; tÞi is of order t implies that spins are entangled over distances of order t. β ¼ 1=3; α ¼ 1=2;z ¼ 3=2: ð5Þ In fact, we show in Sec. V that v t is a sharply defined length scale. But prior to saturation, the relevant length In our context the height of the surface is the bipartite scale for spatial variations in Sðx; tÞ is parametrically entanglement Sðx; tÞ. This is a random quantity that 2=3 smaller than v t; namely, ξðtÞ ∼ t . depends on the realization of the noise in the quantum Before deriving KPZ for entanglement, let us briefly dynamics. The mean height (entanglement) grows linearly consider the status of this equation. At first sight, we might in time, with a universal subleading correction: try to justify this description of Sðx; tÞ simply on grounds of symmetry and coarse graining. If we were describing hðx; tÞ ≡ hSðx; tÞi ¼ v t þ Bt : ð6Þ classical surface growth, we would appeal to translational symmetry in the growth direction (S → S þ const) in order Angle brackets denote an average over noise. The linearity to restrict the allowed terms, and would note that the right- of the leading t dependence is expected from rigorous hand side includes the lowest-order terms in ∂ and ∂ S. x x bounds for various 1 þ 1D random circuits [39–41]. Linear But for entanglement we cannot rely on this simple growth is also generic for quenches in translationally reasoning. First, the transformation S → S þ const is not invariant 1D systems [3,21]. The fluctuations grow as a symmetry (or even a well-defined transformation) of the quantum system. More importantly, it is not clear a priori 2 1=2 β wðx; tÞ ≡ ⟪Sðx; tÞ ⟫ ¼ Ct : ð7Þ that we can write a stochastic differential equation for Sðx; tÞ alone, since the full quantum state contains vastly We refer to w as the width of the surface. The ratio C=B is more information than Sðx; tÞ. Despite these differences universal (the constants v and B are not). The KPZ from simple surface growth, we show below that the fluctuations are non-Gaussian: remarkably, their universal above equation does capture the universal aspects of the probability distribution has been determined analytically entanglement dynamics. [56–66]. In the next section, we exhibit a solvable quantum model The correlation length governing spatial correlations in that maps to a classical surface growth problem that is the fluctuations grows with time as manifestly in the KPZ universality class. Then in the two following sections, we give heuristic arguments for more 1=z ξðtÞ ∼ t ; ð8Þ general systems by making connections with the other members of the KPZ triumvirate. Together with the results for the solvable model, these arguments suggest that KPZ and the equal time correlation function has the scaling form exponents should be generic for entanglement growth in 2 1=2 α any quantum system whose dynamics involves time- GðrÞ ≡ h½Sðx; tÞ − Sðx þ r; tÞ i ¼ r g(r=ξðtÞ): ð9Þ dependent randomness. In Sec. VI, we perform numerical checks on KPZ universality for quantum dynamics in On length scales 1 ≪ r ≪ ξðtÞ, the surface profile SðxÞ discrete time. resembles the trace of a 1D random walk: this is consistent with the exponent α ¼ 1=2. On scales r ≫ ξðtÞ, the fluctua- A. Solvable 1D model tions in Sðx; tÞ and Sðx þ r; tÞ are essentially uncorrelated. At short times the entanglement growth is affected by We now focus on a specific quantum circuit model for the initial conditions, while on very long time scales, of the dynamics of a spin chain with strong noise. We take random order of the system size, the entanglement saturates. unitaries to act on pairs of adjacent spins (i.e., on bonds) at Equations (6)–(9) apply prior to this saturation. In a finite random locations and at random times, as illustrated in Fig. 3. system the asymptotic hSðxÞi profile is that of a pyramid, For simplicity, we discretize time and apply one unitary per with a maximum at height x ¼ L=2, whose height is L=2, time step. (Dynamics in continuous time, with unitaries minus an Oð1Þ correction [67,68]. This profile is reached at applied to the links at a fixed rate in a Poissonian fashion, are a time equivalent.) We choose the initial state to be a product state, 031016-4 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) U For the random dynamics we are considering, the dynami- time step cal rule in Eq. (13) leads to a simple but nontrivial stochastic process. Before discussing its properties, we use Eq. (13) as a FIG. 3. Dynamical update in the solvable model. Application of starting point to show that, in the limit of large Hilbert space a random unitary U to a randomly chosen pair of adjacent spins. dimension, the von Neumann entropy (and in fact all the higher Rényi entropies) obeys the same dynamical rule. This requires an explicit calculation in the limit q → ∞.The von with S ðxÞ¼ 0 for all n and x. We choose the unitaries from Neumann entropy is of more interest than S , since the latter the uniform (Haar) probability distribution on the unitary behaves pathologically in many circumstances. [This is group for a pair of spins, Uðq Þ. This model is solvable in the because it simply counts up all the (nonzero) eigenvalues limit of large local Hilbert space dimension q. in the spectrum of ρ , regardless of how small they are. For example, Hamiltonian dynamics in continuous time—as 1. Dynamics of Hartley entropy opposed to unitary circuits like the above—will generally give an infinite growth rate for S , in contrast to the finite A useful starting point is to consider the n → 0 limit of 0 growth rate for S and the higher Rényi entropies.] the Rényi entropy, S . This is known as the Hartley entropy, vN and quantifies (the logarithm of) the number of nonzero eigenvalues of the reduced density matrix. Equivalently, the 2. Limit of large Hilbert space dimension Hartley entropy determines the minimal necessary value of The present quantum circuit dynamics lead to a solvable the local bond dimension in an exact matrix product model in the limit of large local Hilbert space dimension, representation [31,69] of the state: q → ∞. In this limit, all the Rényi entropies obey the dynamical rule in Eq. (13). S ðxÞ¼ logðbond dimension at xÞ: ð11Þ To show this, we consider the reduced density matrix for a cut at x, where x is the bond to which we are applying the Like the von Neumann entropy, the Hartley entropy of unitary in a given time step. We may write ρ ðt þ 1Þ in neighboring bonds can differ by at most one: terms of ρ ðtÞ and the applied unitary matrix. Averaging x−1 Trρ over the choice of this unitary, we then obtain: jS ðx þ 1Þ − S ðxÞj ≤ 1: ð12Þ 0 0 2 2 2 hTrρ ðt þ 1Þ i ¼ ½Trρ ðtÞ þ Trρ ðtÞ : x Haar x−1 xþ1 Recall that logarithms are base q. For the present, we keep q þ 1 q finite. For the random dynamics we describe above (Sec. II A), See Appendix B for details. In terms of the second Rényi the Hartley entropy obeys an extremely simple dynamical entropy S , this is rule. In a given time step, a unitary is applied at a random −S ðx−1;tÞ−1 −S ðxþ1;tÞ−1 2 2 bond, say, at x. Applying this unitary may change the q þ q −S ðx;tþ1Þ hq i ¼ : ð14Þ Haar Hartley entropy across the bond x; the entropy remains 1 þ 1=q unchanged for all other bonds. The rule for the change in S ðxÞ is that, with probability one, it increases to the The general constraint S ≤ S allows us to write 2 0 maximal value allowed by the general constraint Eq. (12): S ðx; tÞ¼ S ðx; tÞ − Δðx; tÞ; ð15Þ 2 0 S ðx; t þ 1Þ¼ minfS ðx − 1;tÞ;S ðx þ 1;tÞg þ 1: ð13Þ 0 0 0 with Δ ≥ 0. We now use Eqs. (13) and (14) to show that Δ is infinitesimal at large q. Rewriting Eq. (14) in terms of Δ, This “maximal growth” of S occurs with probability one and substituting Eq. (13), immediately shows when all unitaries are chosen randomly. Fine-tuned unitaries (e.g., the identity) may give a smaller value, but these choices Δðx;tþ1Þ Δðx−1;tÞ Δðxþ1;tÞ hq i <q þ q : ð16Þ Haar are measure zero with respect to the Haar distribution. We present a rigorous proof of Eq. (13) in Appendix A. For a simple bound (for a large system, this bound on The appendix also gives a heuristic parameter-counting max hq i is far from the tightest possible since we have not argument that suggests the same result, but as we explain, exploited the large size of the system), define Δ ðtÞ to max the more rigorous argument is necessary as the heuristic be the maximal value of Δðx; tÞ in the entire system. argument can be misleading. We note that Ref. [70] The equation above implies observed that the growth in bond dimension, when a Δ ðtþ1Þ Δ ðtÞ unitary is applied to a matrix product state, is upper max max hq i < 2q : ð17Þ Haar bounded by the right-hand side of Eq. (13) and used this to obtain an upper bound on bond dimension growth during We may iterate this by averaging over successively earlier a quantum computation. unitaries: 031016-5 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) negative sign, meaning that entanglement growth is slower when the coarse-grained surface has a nonzero slope. This is natural given the microscopic dynamics: if the slope is maximal in some region, local dynamics cannot increase the entropy there. 4. Entanglement speed in the solvable model FIG. 4. Surface growth model for entanglement Sðx; tÞ across a In the present model the difference in height between cut at x, in the large-q limit. Applying a unitary to bond x can two adjacent bonds is either ΔS ¼1 or ΔS ¼ 0. At early increase the height of the surface locally [Eq. (19)], correspond- stages of the evolution both possibilities occur. However, ing to dropping a “block” of height ΔS ¼ 1 or ΔS ¼ 2. one may argue that the “flat points,” where ΔS ¼ 0, ðln qÞΔ ðtÞ t become rarer and rarer at late times. [Flat points can max he i < 2 : ð18Þ Haar disappear by “pair annihilation” (Fig. 5, top left), and This shows that as q → ∞ at fixed time t, the probability can diffuse left or right (Fig. 5, top right), but cannot be distribution for Δ concentrates on Δ ¼ 0, so that S and S 2 0 created. As a result, their density decreases with time.] become equal. At late times the model therefore becomes equivalent to the This implies that the entanglement spectrum is flat, well-known “single step” surface growth model [71],in so, in fact, all the Rényi entropies obey Eq. (13) for the which ΔS ¼1 only. An appealing feature of this model is application of a unitary across bond x. that, for a certain choice of boundary conditions (BCs), the late-time probability distribution of the growing interface can be determined exactly [71]. [The solvable case corre- 3. Properties of the solvable model sponds to choosing periodic BCs in the classical problem. The dynamical rule we arrive at for the bipartite von (These BCs are useful for understanding the classical Neumann and Rényi entropies at large q, model, but they do not have an interpretation in terms of entanglement.) In this setting, the mean height grows Sðx; t þ 1Þ¼ minfSðx − 1;tÞ;Sðx þ 1;tÞg þ 1; ð19Þ indefinitely, but the probability distribution for the height fluctuations reaches a well-defined steady state.] This defines a stochastic surface growth model in which Sðx; tÞ shows that on scales smaller than the correlation length is always an integer-valued height profile (Fig. 4). The (and prior to saturation), the interface looks like a 1D remaining randomness is in the choice of which bond is random walk with uncorrelated ΔS ¼1 steps. This updated in a given time step. At each time step, a bond x is confirms the expected KPZ exponent α ¼ 1=2. It also chosen at random, and the height SðxÞ is increased to the allows the mean growth rate of the surface to be calculated maximal value allowed by the neighbors. Figure 5 gives [71]: the mean height increase in a given time step can be examples of local configurations before and after the calculated by averaging over the four possible initial central bond is updated. configurations for a bond and its two neighbors. After This model is almost identical to standard models for rescaling time so that one unit of time corresponds to an surface growth [71,72]. It is in the KPZ universality class (it is average of one unitary per bond, this gives an entanglement straightforward to simulate the model and confirm the growth rate [Eq. (6)]of expected KPZ exponents), and some nonuniversal properties can also be determined exactly (see below). Note that the v ¼ 1=2: ð20Þ boundary conditions S ¼ 0 on the right and the left, and the restriction jSðx þ 1Þ − SðxÞj ≤ 1, imply that the entangle- As we discuss below (Secs. IV and V), one can also ment eventually saturates in the expected pyramid profile. associate a speed v with the growth of quantum operators When we move to the continuum (KPZ) description of under the random dynamics; in the present large-q model, the interface [Eq. (4)], the nonlinear λ term appears with a this speed is (The result v ¼ 1 arises because in the large- q limit, the growth of a typical operator is limited only by the structure of the circuit. In Ref. [73], we give explicit derivations of v in random circuits for arbitrary q.) v ¼ 1: ð21Þ It is interesting to note that here v <v , contrary to 1D E B CFTs (where v ¼ v [1]) and contrary to previous con- E B jectures about generic systems [23].InSec. IV,we givean FIG. 5. Entanglement growth in the large-q model. Effect of applying a random unitary to the central bond, for four choices of appealing intuitive picture for why v can be smaller than v E B the initial local entropy configuration of three adjacent bonds. for the case of Clifford circuits. 031016-6 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) The mapping to surface growth gives us a clean derivation composite tensors L and R gives a representation of the P S P cut i i of universal entanglement dynamics in a solvable model. state as jψðtÞi ¼ L R jai jbi , where jai and a A A A i¼1 a;b b However, this surface growth picture is restricted to the jbi are basis states in A and A, respectively. This implies that entropy for a single cut (as opposed to the entropy of a region cut the Schmidt rank for a bipartition into A and A is at most q , with multiple end points) and to one dimension. It will be so that S , which is the logarithm of the Schmidt rank, is at useful to find a more general language which extends the most S . In turn, S ≤ S for any n ≥ 0.] cut n 0 above results. To do this, we now make a connection with the The best bound of this type is given by the minimal cut, second member of the KPZ triumvirate (Fig. 2), the statistical which passes through the smallest number of legs. We denote mechanics of a polymer in a random environment. the corresponding estimate for the entropy S ðxÞ.Ifthe min−cut geometry of the circuit is random, S ðxÞ and the min−cut III. DIRECTED POLYMERS AND MINIMAL CUT corresponding curve are also random. Finding S ðxÞ min−cut amounts to an optimization problem in a classical disordered In this section, we relate the dynamics of Sðx; tÞ to the geometry of a minimal cut through the quantum circuit system. which prepares the state (Fig. 6). This provides an alternative In the solvable large-q model, S ðxÞ in fact gives min−cut perspective on the exact result Eq. (19) for the solvable the von Neumann entropy exactly. This follows straight- model, and also a useful heuristic picture for noisy quantum forwardly from the results of the previous section (see dynamics in general. This line of thinking reproduces KPZ below). In a typical microscopic model, on the other hand, behavior in 1D. Importantly, it also allows us to generalize to S is only a bound on the true entropy. Nevertheless, min−cut higher dimensions and to more complex geometries. we conjecture that the following picture based on the Our starting point is the minimal cut bound for tensor minimal cut is generally valid as a coarse-grained picture: networks. This very general bound has been related to the i.e., that it correctly captures the universal properties of the entanglement dynamics in noisy systems. This conjecture is Ryu-Takayanagi formula for entanglement in holographic equivalent to the applicability of the KPZ description to conformal field theories [50–53,74], and has also been applied to unitary networks as a heuristic picture for generic noisy systems; further evidence for the latter is in entanglement growth [10]. Secs. IV and VI. The problem of finding the minimal curve is a version of a Consider again a random quantum circuit in 1 þ 1D, and a well-studied problem in classical statistical mechanics, curve like that in Fig. 6, which bisects the circuit and divides known as the directed polymer in a random environment the physical degrees of freedom into two at position x. Any (DPRE) [47,75]. Here, the “polymer” is the curve which such curve gives an upper bound on the entanglement: all the bisects the circuit, and its energy EðxÞ is equal to S ðxÞ,the Rényi entropies satisfy SðxÞ ≤ S , where S is the number cut cut cut number of legs it bisects. The spatial coordinate of the of “legs” that the curve passes through. (This relies only on polymer’s upper end point is fixed at x, while the lower end linear algebra: the rank of the reduced density matrix ρ is at S point is free. Finding S ðxÞ is equivalent to finding the cut min−cut most q .) [The cut divides the tensor network into two parts minimal value of the polymer’s energy. This corresponds to connected by S bonds. One part contains the physical cut the polymer problem at zero temperature; however, the legs for subsystem A and the other part contains those for the universal behavior of the DPRE is the same at zero and at complement. Regarding the two parts of the circuit as nonzero temperature. (For any finite temperature, the DPRE flows under renormalization to a zero-temperature fixed point at which temperature is an irrelevant perturbation.) DPRE models with short-range-correlated disorder are in the same universality class as the KPZ equation [45]. Let Eðx; tÞ be the minimal energy of the polymer in a sample of height t. We may increase t by adding an additional layer to the top of the sample. Eðx; t þ δtÞ can then be expressed recursively in terms of Eðy; tÞ for the various possible values of y. In the continuum limit, this leads to an equation for Eðx; tÞ which is precisely the KPZ equation (see Refs. [45,76] for more details of the mapping between DPRE and KPZ). The KPZ exponents we give in Sec. II may therefore be applied to the energy of the FIG. 6. Any cut through the unitary circuit that separates the polymer. The exponent z ¼ 3=2 also determines the length legs to the left and right of x (on the top boundary) gives an upper scale for transverse fluctuations of the polymer on large bound on Sðx; tÞ. The best such bound is given by the minimal cut length and time scales: (note that the cut shown here is not the minimal one). Finding the minimal cut in a random network is akin to finding the lowest- 2=3 Δx ∼ ðΔtÞ : ð22Þ energy state of a polymer in a random potential landscape. 031016-7 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) Since in our case the minimal E is simply S ,we For a finite interval of length l in an infinite system there is min−cut find that the latter executes KPZ growth. In light of the a crossover between a configuration with two vertical cuts previous section, this is not surprising. In fact, in our and one with a single horizontal cut, giving instead SðtÞ¼ solvable model, S is exactly equal to the true entan- 2v tfðl=2v tÞ. min−cut E E glement entropy (in the large-q limit). This follows from the These scaling forms are our first confirmation that v is fact that the recursive construction of Eðx; tÞ described really a speed, as well as a growth rate for the entanglement. above (on the lattice, rather than in the continuum) precisely We give an independent derivation of this fact for Clifford matches the large-q dynamics of Eq. (19). Examples of circuits in the following section, and test the above scaling nonunitary tensor networks in which the minimal cut bound form numerically in Sec. VI B. We discuss the interpreta- becomes exact are also known [52], including a large-bond- tion of v further in Sec. V. dimension limit similar to that we discuss here [53]. Note that fluctuations have dropped out of Eq. (23) as a The utility of the DPRE picture is that it is far more result of considering only the leading-order behavior of generalizable than the surface growth picture, which is Sðx; tÞ. These scaling forms agree with the results for restricted to the entropy across a single cut in 1D. As we holographic CFTs [8] and with a heuristic application of note above, the value of S in a given microscopic the minimal cut formula to a regular tensor network [10]. min−cut model is typically not equal to any of the physical entropies Here, we see them emerging from a simple and well-defined S with n> 0. Nevertheless, we conjecture that the DPRE coarse-grained picture, suggesting that they are universal for and KPZ pictures are valid universal descriptions for all all generic 1D systems, including, for example, translation- noisy models, so long as they are not fine-tuned or nonlocal. ally invariant but nonintegrable spin chains. [Reference [77] This includes noisy Hamiltonian dynamics in continuous includes numerical tests of scaling forms derived from the time (we discuss this case further in Sec. IX). If we restrict to directed polymer picture in deterministic systems, including the leading-order deterministic behavior, we can also make extensions to inhomogeneous systems (a chain with a weak link).] It is also worth noting that Eq. (24) is capable of conjectures about Hamiltonian systems without noise. distinguishing generic systems from (nonrelativistic) inte- grable systems. In the latter case the quasiparticle picture A. Scaling form for entanglement saturation applies and yields different profiles for SðtÞ [3,19].For At leading order in time, the growth of the height Sðx; tÞ relativistic systems in which the quasiparticle picture holds is deterministic: fluctuations are a subleading effect when t (rational CFTs), all quasiparticles travel at the same speed, is large. Similarly, Eq. (22) shows that the coarse-grained and as a result Eq. (24) does apply [1,3] (however, the minimal cut is essentially vertical (prior to saturation of the entanglement of more complex regions will differ between entropy): the length scale for its transverse fluctuations is the quasiparticle picture, on the one hand, and the results negligible in comparison with t. These pictures therefore from holographic systems and the minimal cut picture, on have well-defined and simple deterministic limits. They the other hand [3,4,44]). lead directly to deterministic scaling forms for the leading- In Secs. V and VIII we propose that the above picture in order behavior of the entanglement, which we discuss in terms of a coarse-grained minimal cut is the simplest way more detail in Sec. V. Here, we consider the simplest case, to understand the basic features of the entanglement the saturation of the entanglement entropy Sðx; tÞ across a tsunami for generic many-body systems (with or without single cut (or for a single interval). We reproduce a simple noise) in both 1D and higher dimensions. scaling function known from other contexts [8,10,21]. The definition of the entanglement growth rate implies IV. HYDRODYNAMICS OF OPERATOR that the energy E of such a vertical cut is v t to leading SPREADING order. The entanglement in a finite system grows at this rate until time t ¼ x=v , when it reaches its saturation An alternative way to think about the quantum dynamics is saturation E value S ¼ x. (Here, we are neglecting subleading terms and in terms of the evolution of local operators O . For example, a assuming x<L − x.) After this time a vertical cut is no Pauli operator initially acting on a single spin (e.g., O ≡ Y ; i i longer favorable: instead, the minimal cut exits the circuit we denote the Pauli matrices by X, Y, Z) will evolve with time via the left-hand side. Its shape is no longer unique, but it into an operator O ðtÞ which acts on many spins. Operators can be taken to be horizontal, and it has energy E ¼ x. typically grow ballistically [38], in the sense that the number This picture corresponds to a simple scaling form (again, of spins in the support of O ðtÞ grows linearly with t. In this neglecting subleading terms): section, we relate the growth of the bipartite entanglement to the spreading of operators. We focus on the special case of Sðx; tÞ¼ v tfðx=v tÞ; ð23Þ E E unitary evolution with Clifford circuits (defined below), but with we expect the basic outcomes to hold for more general unitary dynamics. We find that the entanglement growth rate u for u< 1 is not given by the rate at which a single operator grows, but is fðuÞ¼ ð24Þ 1 for u ≥ 1. instead determined by collective dynamics involving many 031016-8 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) number of sites to be L). These operators, denoted O (i ¼ 1; …;L), satisfy O jΨ i¼jΨ i: ð25Þ i 0 0 FIG. 7. Spreading of stabilizer operators defining the quantum state (Sec. IV). Each blue particle marks the right end point of For example, if the spins are initially polarized in the y some stabilizer (the rightmost spin on which it acts). Blue direction, we may take O ¼ Y . At a later time, the above particles hop predominantly to the right. Whenever a particle i i equation still holds, with each stabilizer O replaced with enters the right-hand region (A) the entanglement S increases by one bit. The particle density is described by the noisy Burgers the time-evolved stabilizer O ðtÞ¼ UðtÞO U ðtÞ, where i i equation, which maps to KPZ. A “hole” (empty circle) marks the UðtÞ is the unitary operator that evolves the initial state to left-hand end point of some stabilizer. the state at time t. [Note that O ðtÞ is not the standard Heisenberg picture operator, which would have U and U in the other order.] operators. Remarkably, in 1D these collective dynamics have In the following, we focus on evolution of the initial state a long wavelength hydrodynamic description. with unitary gates in the Clifford group [80]. Such gates have This hydrodynamic description turns out to be the noisy recently been used in toy models for many-body localization Burgers equation, which is related to the KPZ equation by a [29]. Entanglement generation in non-random Clifford simple change of variable and is the final member of the KPZ circuits has also been studied [43]. The defining feature of triumvirate shown in Fig. 1. In the present case the hydro- Clifford unitaries is that they have a simple action on Pauli dynamic mode is the density of certain fictitious “particles,” operators: single-spin Pauli operators are mapped to products shown in blue in Fig. 7. The quantum state is defined by a set of Pauli operators. of operators (Sec. IVA) which spread out over time, and the Any product of Pauli matrices can be written as a product particles are markers which show how far these operators of X and Z matrices, so to follow the dynamics of a given have spread. We derive their coarse-grained dynamics in stabilizer O ðtÞ, we need only keep track of which X Sec. IV B after introducing the necessary operator language. i i and Z operators appear in this product. Furthermore, the In generic many-body systems (with local interactions) overall sign of the stabilizer O ðtÞ does not affect the this process of operator growth is characterized by a speed known as the butterfly speed v . This speed defines an entanglement properties of a system undergoing Clifford evolution, so we do not keep track of it. By writing O ðtÞ as effective light cone within which the commutator between the spreading operator OðtÞ and a typical local operator v v v v is appreciable. The quantity v is a characteristic speed for 1x 1z Lx Lz O ðtÞ ∝ X Z …X Z ; ð26Þ i L L 1 1 the spreading of quantum information in a given model, and can be extracted from appropriate correlation functions. In we may specify any stabilizer by a binary vector v⃗ with 2L deterministic systems (time-independent evolution) v can components: depend on temperature, and typically does not saturate the well-known Lieb-Robinson bound [78,79]. Generic noisy v⃗ ¼ðv ;v ; …;v ;v Þ: ð27Þ systems equilibrate to infinite temperature, so in the present 1x 2x Lx Lz models there is no notion of temperature dependence—v is a constant defined entirely by the dynamics. For example, the first component of the vector v ¼ 1 if X 1 1 The scaling forms we discuss in the previous section appears in the product, and v ¼ 0 otherwise. The binary show that in 1D there is a well-defined speed v associated vector corresponding to a stabilizer O ¼ Y is i i with entanglement spreading. The following picture gives a physical interpretation of this speed, in terms of a certain set v⃗ ≡ ð0; …; 0; 1; 1; 0; …; 0Þ; ð28Þ of growing operators. However, it also shows that in general the speed v is smaller than the speed v .Thisisperhaps E B where the locations of the nonzero elements correspond to surprising: in 1D CFTs the two speeds are equal, and it X and Z . i i has been conjectured that they are equal in general [23]. We consider the dynamics in two stages. First, we (Note that we already encountered a solvable model with consider the evolution of a single operator. Then, we v ¼ v =2 in Sec. II A.) E B generalize this to understand the dynamics of the state. How does a single stabilizer O ðtÞ evolve? Applying a A. Stabilizer operators one- or two-site Clifford unitary to O ðtÞ corresponds to It is convenient to use the language of “stabilizer” applying simple local updates to the string v⃗ . Although the operators to describe the entanglement dynamics. We precise details of these updates are not crucial, we now give some explicit examples of gates that we encounter again in may define the initial state jΨ i by specifying L stabilizers under which it is invariant (in this section, we take the the numerical simulations. 031016-9 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) As single-site examples, consider the Hadamard and on v⃗ is invariant under the dynamics, regardless of the phase gates. The Hadamard is a rotation on the Bloch probabilities with which the gates are applied. By standard sphere [the rotation is by π around the (1,0,1) axis] that properties of Markov processes, this is the unique asymp- exchanges the X and Z axes, totic distribution to which the system tends, so long as the choice of Clifford gates is not fine-tuned to make the process nonergodic. (If the gate set includes each gate and R ¼ pffiffiffi ðX þ ZÞ; ð29Þ its inverse with the same probability, detailed balance is also obeyed, but this is not necessary.) We expect v⃗ to so applying a Hadamard to site i updates the string by equilibrate locally to this structureless state on an Oð1Þ time v ↔ v . The phase gate is a rotation around the Z axis ix iz scale, and similarly for the internal structure of operators which maps X to Y ¼ iX Z : i i i i smaller than L.] Now consider the dynamics of a quantum state. Once the pffiffiffiffi R ¼ Z: ð30Þ sign information in Eq. (26) is dropped, the relevant information in the state jΨðtÞi is contained in binary 1 L This means that an additional Z is generated whenever X i i vectors v⃗ ; …; v⃗ corresponding to the L stabilizers. We is present in the string, or equivalently v → v þ v iz iz ix may package this information in a 2L × L matrix: ðmod 2Þ. For a two-site example, consider the left and right 1⊤ L⊤ controlled- NOT (CNOT) gates acting on the leftmost spins ΨðtÞ¼ðv⃗ ; …; v⃗ Þ: ð32Þ in the chain. In the Z basis, the action of these operators is to flip the “target” spin if and only if the “control: spin is Each column corresponds to a stabilizer, and each row to a down: spin operator X or Z . The dynamical updates correspond i i to row operations (with arithmetic modulo two) on this ðLÞ matrix. For example, a Hadamard gate exchanges the rows CNOT ¼ ½ð1 þ Z Þþð1 − Z ÞX ; 1 1 2 corresponding to X and Z . i i A crucial point is that there is a large gauge freedom in ðRÞ CNOT ¼ ½ð1 þ Z Þþð1 − Z ÞX : ð31Þ 2 2 1 this definition of the state. This gauge freedom arises because we can redefine stabilizers by multiplying them ðLÞ Conjugating the Pauli matrices by CNOT yields: together. For example, if a state is stabilized by fX ;Z g, 1 2 then it is also stabilized by fX Z ;Z g, and vice versa. This 1 2 2 X → X X;Z → Z;X → X;Z → Z Z : freedom to redefine the stabilizers corresponds to the 1 1 2 1 1 2 2 2 1 2 freedom to make column operations on Ψ, or equivalently We see that the operator X is added to the string if X is the freedom to add the vectors v⃗ modulo two. Note that by 2 1 ðLÞ making such a “gauge transformation” we may be able to present (and similarly for Z and Z ). Applying CNOT 1 2 reduce the size of one of the stabilizers, giving a more therefore updates v⃗ by compact representation of the state. The final fact we need is an expression for the entropy v → v þ v ðmod 2Þ;v → v þ v ðmod 2Þ: 2x 2x 1x 1z 1z 2z S ðtÞ of a region A in terms of the stabilizers. Heuristically, ðRÞ this is given by the number of stabilizers that have spread CNOT acts similarly with the roles of the spins reversed. into region A from outside. More precisely, define I as the It is simple to argue that random application of such number of stabilizers that are independent when restricted operations causes the region of space in which v⃗ is nonzero to region A [81]. (Independence of the stabilizers corre- to grow ballistically. This corresponds to the operator spreading itself over a region of average size 2v t, where sponds to linear independence of the vectors v⃗ , with v is the operator spreading (butterfly) velocity for this arithmetic modulo two, once they are truncated to region system [79]. (For the present system, this velocity is also A.) The entropy is equal to [82,83] the analogue of the Lieb-Robinson velocity.) The value of v depends on the precise choice of dynamics, but it is the S ðtÞ¼ I − jAj; ð33Þ B A A same for all initial operators so long as the dynamics (the probability distribution on gates) is not fine-tuned. Further, where jAj is the number of sites in A. See Appendix C for a one may argue that the interior of the region where the simple derivation of Eq. (33). For Clifford dynamics all string v⃗ is nonzero is “structureless.” Within the interior, v⃗ Rényi entropies are equal, so we omit the Rényi index on S. rapidly “equilibrates” to become a completely random The maximal value of I is 2jAj,so S is bounded by jAj as A A binary string. [Consider the late-time dynamics of an expected. operator, or equivalently its string v⃗ ,inan L-site system. This formula has a simple interpretation. In the initial Random application of Clifford gates gives random dynam- product state we may take one stabilizer to be localized at ics to v. It is easy to see that the flat probability distribution each site, so I ¼jAj and the entanglement is zero. As time 031016-10 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) goes on, stabilizers that were initially localized outside of A commute, and this is preserved by the unitary dynamics grow and enter A. Each time a new independent operator and the redefinitions of the stabilizers). [Consider the case appears in A, the entanglement S ðtÞ increases by one bit. where ρ ðiÞ¼ 1: for example, let the corresponding sta- A r The linear independence requirement in the definition of I bilizer read O ¼ …X . Any stabilizer contributing to ρ ðiÞ A i l is crucial, as it leads to effective interactions between the must be of the form X … in order to commute with O.By stabilizers, which we discuss in the following section. the rule imposed in the text, this means that ρ ðiÞ ≤ 1.] From now on, we take A to consist of the spins to the Since there are a total of 2L particles which all have to live right of the bond x, and revert to the notation S ¼ S used somewhere, we have Eq. (34). A x in the rest of the text for the entanglement across a cut at x. With this convention, the dynamics of the bipartite entropy SðxÞ is simply related to the hopping dynamics of the particles. By Eq. (34) it suffices to consider only the r B. Coarse-grained operator dynamics particles: an l particle is just an r “hole.” We write the Each stabilizer O ðtÞ (labeled i ¼ 1; …;L) has a left and density ρ of r particles as ρ. See Fig. 7 for a typical a right end point l and r , marking the extremal spins i i configuration in a subregion of the system. included in the stabilizer. We view l and r as the positions i i The utility of the canonical form Eq. (34) is that the of two fictitious particles of type l and r, represented in independence requirement becomes trivial. One can easily white and blue, respectively, in Figs. 7 and 8. There are L of check that all the operators that have spread into A (the each type of particle in total. region to the right of x) are independent. (Consider the In the initial product state, O ð0Þ is a single Pauli stabilizers that act in region A, i.e., the stabilizers with operator on site i, say, Y . This means that each site has r >x. We may argue by contradiction that they remain one l particle and one r particle (since l ¼ r ¼ i), as i i independent after truncation to subsystem A. If not, this shown in Fig. 8 (left). As time increases, the r particles will means there is some product of the truncated stabilizers that typically move to the right and the l particles to the left. equals one. Let the rightmost spin appearing in any of these The nature of this motion depends on how we define the stabilizers be j. But by our convention for “clipping” the stabilizers. At first sight, the obvious choice is to define stabilizers, it is impossible for the Pauli matrices acting on O ðtÞ as the unitary time evolution of the initial stabilizer, spin j to cancel out when they are multiplied together. O ðtÞ¼ UðtÞY U ðtÞ. But in fact, it is useful to exploit the i i Therefore, the operators must, in fact, be independent.) gauge freedom in the choice of stabilizers to impose a Therefore, to find SðxÞ we need only count the number different “canonical” form. One result of this is that the of r particles to the right of the cut and subtract the number stabilizers effectively grow more slowly than the butterfly of sites: velocity v (discussed in the previous section) for the SðxÞ¼ ðρ − 1Þ: ð35Þ spreading of an operator considered in isolation. i>x Let ρ ðiÞ and ρ ðiÞ be the number of particles of each l r type at site i. The constraint that we impose is To reiterate, the entanglement increases by one every time an r particle drifts rightward across bond x (and decreases ρ ðiÞþ ρ ðiÞ¼ 2: ð34Þ l r by one if it drifts across in the other direction). To see that we can impose this constraint, consider the Now consider the dynamics of the particles. situation ρ ðiÞ¼ 3, so that there are three stabilizers that Microscopically, a dynamical time step involves (1) appli- start at i. The initial element of each string can be either X, cation of a unitary gate and (2) potentially a “clipping” of Y,or Z.If ρ ðiÞ¼ 3, it is impossible for all three initial stabilizers to enforce the canonical form. Effectively, the elements to be independent. We can then redefine one of particles perform biased diffusion, with the restriction that the stabilizers, by multiplying it by one or both of the more than two particles cannot share a site: others, in such a way that its length decreases by one. (By choosing the longer stabilizer we avoid adding length at the ρ ≤ 2: ð36Þ right-hand side.) Making reductions of this kind wherever possible guarantees that ρ ðiÞ ≤ 2, and also that if ρ ¼ 2, l l This constraint leads to “traffic jam” phenomena familiar the initial elements of the two stabilizers are distinct. from the so-called asymmetric exclusion process [84], and (And similarly for ρ .) With this convention it also follows to the same continuum description. Our essential approxi- that ρ ðiÞþ ρ ðiÞ ≤ 2: otherwise, the operators involved l r mation is to neglect the detailed internal structure of the could not commute, which they must (the initial stabilizers stabilizers and to treat the dynamics of the end points as effectively Markovian. We expect this to be valid at long length and time scales for the reason we mention in the previous section: the internal structure of the operators is essentially featureless and characterized by finite time FIG. 8. Left: The initial product state represented in terms of the fictitious particles. Right: A state with maximal SðxÞ. scales. 031016-11 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) We now move to a continuum description. The coarse- grained density obeys a continuity equation: ∂ ρ ¼ −∂ J; ð37Þ t x FIG. 9. Fluctuations at late times, after saturation of hSðxÞi,in with J the particle current. Further, there is a symmetry the Clifford case. When x ≪ L=2 it requires a rare fluctuation under spatial reflections, which exchange left and right end (fighting against the net drift) to remove a particle from region A, points (ρ ↔ ρ ). Writing leading to an exponentially small S ðxÞ − hSðxÞi. l r max ρ ¼ 1 þ Δρ; ð38Þ Let us revert to our previous notation, where the system has L þ 1 sites and bonds are labeled x ¼ 1; …;L. Without where Δρ is the deviation from the mean density, the loss of generality we take x ≤ L=2. When fluctuations are reflection symmetry is neglected, the region to the left of x is empty of r particles, x → −x; Δρ → −Δρ: ð39Þ and the entropy is maximal, S ðxÞ¼ x. Fluctuations will max reduce the average. But in order for SðxÞ to fluctuate To obtain a long wavelength description, we write the downward, a blue r particle must diffuse leftward from the current as a power series in Δρ and ∂ . Keeping the lowest- right half of the system in order to enter the region to the order terms that respect the symmetry, left of x, as in Fig. 9. This is a fluctuation by a distance ∼ðL=2 − xÞ. Such fluctuations are exponentially rare J ¼ c − ν∂ Δρ − ðΔρÞ þ η: ð40Þ x events, because they fight against the net rightward drift for the r particles. Thus, when L=2 − x is large, we expect These terms have a transparent physical meaning. The drift −αðL=2−xÞ constant c> 0 reflects the fact that the average motion is to S − hSðxÞi ∼ e : ð43Þ max the right (i.e., operators grow over time). The ν term is simple diffusion. The noise η reflects the randomness in the Our coarse-grained picture does not determine the numeri- dynamics. Most importantly, the nonlinear λ term is the cal constants. effect of the constraint Eq. (34). It reflects the fact that The detailed nature of these exponentially small correc- the current is maximal when the density is close to one. The tions will differ between Clifford circuits and more general current evidently vanishes when ρ ¼ 0, since there are no unitary circuits. (For example, in the Clifford case S is particles, but also when ρ ¼ 2 (the particles cannot move if independent of n, while in general the corrections will the density is everywhere maximal). Therefore, we depend on n [68].) Nevertheless, the functional form above expect λ > 0. agrees with the late-time result for generic gate sets, which From the above formulas, the density obeys is simply the mean entanglement in a fully randomized pure state [67,68]: 2 2 ∂ ρ ¼ ν∂ ρ þ ∂ ðρ − 1Þ − ∂ η; ð41Þ t x x x jAj−jAj −ðL=2−xÞ 2 4 S − hSðxÞi ≃ ¼ : ð44Þ known as the noisy Burgers equation [84]. The entangle- max 2 ln 2 4 ln 2 ment S ¼ ðρ − 1Þ obeys ∂ S ¼ J, leading to the KPZ equation: jAj¼ x and jAj¼ L − x þ 1 are the numbers of sites in A and its complement. For (generic) Clifford dynamics, the 2 2 ∂ S ¼ c þ ν∂ S − ð∂ SÞ þ η: ð42Þ probability distribution of the entanglement at asymptoti- t x x cally late times will be that of a random stabilizer state. This The sign of λ is in agreement with that obtained from the has been calculated in Ref. [85]. surface growth picture in Sec. II and from the directed polymer picture in Sec. III. While we focus here on V. ENTANGLEMENT TSUNAMI: SPEEDS dynamics of a restricted type (Clifford), this derivation AND SCALING FORMS of KPZ for entanglement provides independent support for It is not a priori obvious that the rate v governing the arguments in the previous sections. In the language of the particles, the initial state corre- entanglement growth can be viewed as a speed in generic sponds to uniform density ρ ¼ 1. Saturation of the entan- systems (see Ref. [79] for a recent discussion), although glement corresponds (neglecting fluctuations) to all of the r this is known to be the case in holographic CFTs [8].Our particles accumulating on the right-hand side and all of the l results in the directed polymer picture and in the operator particles on the left (Fig. 8), i.e., to a step function density. spreading picture suggest that v is indeed a well-defined As an aside, it is interesting to consider fluctuations in speed in generic systems. (We see in the previous section SðxÞ at late times, i.e., long after the saturation of hSðxÞi. that there is a simple visual interpretation of this speed in 031016-12 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) the stabilizer formalism.) However, this speed is in general In the 1D case, the deterministic scaling form for the smaller than the speed v which governs the spreading of entanglement (of an arbitrary region) which results from the an operator considered in isolation: “thermalization is leading-order directed polymer picture is rather simple, and slower than operator spreading.” is not new—it agrees with holographic results [8,44], and In the stabilizer context the difference between v and v as noted in Ref. [10], can also be obtained from a more E B arises because in enforcing Eq. (34) we “clip” the stabi- microscopic minimal cut picture in which the geometry of lizers, reducing their rate of growth. We believe the the minimal cut is highly nonunique. We propose that phenomenon of v being smaller than v to be general coarse graining fixes the geometry of the minimal cut. The E B (see also the result in Sec. II A 4) and relevant also to non- derivation of these scaling forms from a simple coarse- noisy dynamics. This picture is contrary to that of. e.g., grained picture suggests that they are universal in non- Ref. [23], where the operator spreading velocity is assumed integrable, translationally invariant systems. (These scaling to determine the entanglement growth rate. In the presence forms are generally not the same as those obtained from the of noise, one may also argue that a picture of independently quasiparticle picture for rational CFTs [3,4].) Our deriva- spreading operators underestimates the exponent governing tion also opens the door to generalizations to higher the growth of fluctuations. [Considering the unitary evo- dimensions (Sec. VIII B) and to 1D systems with quenched lution of a single operator in isolation, its right end point disorder [77]. executes a biased random walk, traveling an average We now consider some examples of the scaling of the 1=2 mutual information. This will help clarify the operational distance v t with fluctuations Oðt Þ. If we were to meaning of the speed v . neglect the independence requirement in Eq. (33), then To calculate the entanglement S of a region A,we must the entanglement would be estimated (incorrectly) as the take a cut, or multiple cuts, with end points on the boundary number of independently spreading operators which have points of A at the top of the space-time slice. These cuts can reached A. The mean of this quantity is v t and the 1=4 either be vertical, in which case they cost an energy’v t fluctuations are of order t . This is related to the differ- (to use the language of Sec. III), or they can connect two end ence between the KPZ universality class of surface growth, points, in which case we take them to be horizontal and to which is generic, and the Edwards-Wilkinson universality have an energy equal to their length. The entanglement S ðtÞ class, which applies when the strength of interactions is A is given by minimizing the energy of the cut configuration. fine-tuned to zero [45].] It is a continuous piecewise linear function, with slope The language of a tsunami is often used in discussing discontinuities when the geometry of the minimal cut entanglement spreading, so it is nice to see that–at least in configuration changes. To generalize the conjecture to 1D–entanglement spreading can be related to a hydro- systems without noise, we must allow for the fact that the dynamic problem. (The motivation for the tsunami termi- asymptotic value of the entanglement depends on the energy nology is the idea that for a region A, the entanglement S density of the initial state. We therefore replace the entan- is dominated by a subregion close to the boundary which glement S in the formulas with S=s , where s is the grows ballistically, like the advancing front of a tsunami.) eq eq equilibrium entropy density corresponding to the initial In higher dimensions the boundary of an operator has a more complicated geometry, so the hydrodynamic corre- energy density [2,8]. This ensures that the entanglement spondence we describe above does not generalize. entropy of an l-sized region matches the equilibrium thermal In order to understand the entanglement tsunami better, entropy when v t ≫ l=2, as required for thermalization. we now return briefly to the directed–polymer–in–a– Heuristically, s defines the density of “active” degrees of eq random–medium picture developed for noisy systems in freedom at a given temperature [79]. Sec. III. To clarify the meaning of v , consider the mutual information between two semi-infinite regions that are separated by a distance l (Fig. 10). With the labeling of A. Scaling forms for the entanglement tsunami the regions as in the figure, this is given by When all length and time scales are large, fluctuations in the entanglement are subleading. Neglecting them is equiv- I ¼ S þ S − S ¼ S þ S − S : ð45Þ AC A C A∪C A C B alent to saying that the coarse-grained minimal cut (prior to saturation) is a straight vertical line. This deterministic picture generalizes to the entanglement or mutual informa- tion of arbitrary regions, and also to higher dimensions (Sec. VIII). We conjecture that these pictures are valid for the long-time behavior of entanglement quite generally. The FIG. 10. Infinite chain with regions A, B, C marked. B is of setup relevant to us in the non-noisy case is a quench, in length l while A and C are semi-infinite. The mutual information which the initial state is a ground state of one Hamiltonian, between A and C is nonzero so long as l< 2v t: correlations and a different Hamiltonian is used for the evolution. exist over distances up to 2v t, not v t. E E 031016-13 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) We have S ¼ S ¼ v t for all times, since the appropriate A C E minimal cuts are vertical. If l> 2v t, S is given by two E B vertical cuts, so I vanishes. When l< 2v t, S is instead AC E B dominated by a horizontal cut, so that I ¼ 2v t − l. AC E The entanglement tsunami is sometimes taken to mean that at time t,a “boundary layer” of width v t inside a given region is entangled with the exterior. If this region were maximally entangled with the exterior, this would repro- duce the correct value of the entanglement across a cut (S ¼ v t). However, this picture is not correct: the result for the mutual information shows that correlations exist over FIG. 13. Bottom: Semi-infinite chain with regions A, B (length distances up to 2v t, not v t. So although v is a speed, it E E E l and l , respectively), and C adjacent to the boundary. Top: The A B should not be thought of as the speed at which the boundary mutual information between B and C for this geometry, for the of the entangled region moves. two regimes indicated. Although the rule for calculating the entanglement is almost trivial, the consequences are not always intuitively [For a simpler example of exponentially small values of obvious. First consider the case where the regions A and C the mutual information, consider hI i at infinite times in a above are finite rather than infinite (and embedded in an AC finite system. If the system contains L qubits and A ∪ C infinite chain); see Fig. 11. When the length d of the regions contains N qubits, the mutual information is exponentially A and C exceeds their separation l, the time dependence of small whenever N< L=2,and given by Eq. (44) as the mutual information is as shown in Fig. 11. [The sequence −1 −ðL−2NÞ hI of minimal cut configurations required for calculating S i ∼ ð2 ln 2Þ 2 .] [The sequence of cuts for S AC B in this case is simply (a), (c)]. in this case is (a), (b), (c) shown in Fig. 12.] By contrast, when Finally, consider the effect of a boundary. Take a semi- the separation l exceeds the length d, the mutual information infinite chain with regions A, B, C adjacent to the boundary is always zero (or more precisely, exponentially small). as in Fig. 13 (C is semi-infinite). Consider the mutual information between B and C, I ¼ S þ S − S .We BC B C A must distinguish the case l <l =2 from the case A B l >l =2. (In the former case, the first event is that the A B minimal cut at the boundary of A goes from being vertical to being horizontal; in the latter, the first event is that the two vertical cuts at the boundary of B are replaced by a horizontal one.) The resulting expressions for I are BC plotted in Fig. 13. VI. NUMERICAL EVIDENCE FOR KPZ GROWTH We now give numerical evidence that noisy entangle- FIG. 11. Bottom: Infinite chain with finite regions A and C each ment growth in 1D is in the KPZ universality class. We of length d, separated by distance l. Top: The mutual information study the time evolution of spin-1=2 chains with open between A and C in the case d> l. In the opposite regime the boundary conditions, taking the initial state to be a pro- mutual information vanishes. duct state with all spins pointing in the same direction (either the positive y or z direction) and keeping track of the entanglement entropy across each bond during the evolu- tion. The discrete time evolution is a circuit of one- and two-site unitaries. Figure 14 shows the structure of a single time step: two layers of two-site unitaries are applied, one layer on odd and one on even bonds, together with single- site unitaries. Each unitary is chosen independently and randomly (from a certain set specified below). We use the symbol R to denote a generic single-site unitary and U to denote a two-site unitary. We consider three kinds of dynamics, distinguished by FIG. 12. Sequence of minimal cut configurations (red lines) the choice of unitaries. To begin with, we study “Clifford determining the entropy of region B in Fig. 11. (a) gives way to (b) when 2v t ¼ l and (b) gives way to (c) when 2v t þ l ¼ 2d. evolution” in which the unitaries are restricted to the set of E E 031016-14 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) In all our simulations, each two-site unitary U in the circuit is chosen with equal probability from three pos- sibilities: the two types of CNOT gate [Eq. (31)] and the identity matrix: ðLÞ ðRÞ U ∈ f1; CNOT ; CNOT g: ð46Þ In this section, we discuss the simplest Clifford dynamics, which includes only these gates, and no one-site unitaries (R ¼ 1). When the initial state is polarized in the y– direction, this set of gates is sufficient to give nontrivial entanglement evolution, with universal properties that turn FIG. 14. Schematic structure of a layer in the quantum circuits out to be the same as those for more generic gate sets. We used for simulations. also study the “full” Clifford dynamics in which all the Clifford generators are used, choosing the single-site so-called Clifford gates (Sec. IV). Clifford evolution can be unitaries randomly from the three options simulated efficiently (in polynomial time) using the stabi- lizer representation we discuss in Sec. IV. This allows us to R ∈ f1;R ;R g: ð47Þ H P access very long times and to pin down KPZ exponents accurately. Next, we study more general dynamics for Results for this case are similar and are given in which polynomial-time classical simulation is impossible, Appendix. D. giving evidence that KPZ behavior holds more generally. To begin with, Fig. 15 shows the evolution of the The two types of non-Clifford dynamics we study here are bipartite von Neumann entropy SðxÞ (in units of log 2) referred to as the phase evolution and the universal for a single realization of the noise (i.e., a particular random evolution: we give details below. For these dynamics we circuit) in a system of L ¼ 459 sites. The curves show use a matrix product representation of the state imple- successively later times. Note that the entropy saturates mented via ITensor [86]. more rapidly closer to the boundary, because the maximum The fingerprints of KPZ behavior that we search for are entanglement across a bond is proportional to its distance the two independent critical exponents β and α (Sec. II). We from the boundary. At very late times Sðx; tÞ saturates to a extract β both from the fluctuations in the von Neumann pyramidlike profile representing close-to-maximal entan- entropy and from the corrections to the mean value glement. Our interest is in the stochastic growth prior to [Eqs. (6) and (7)], and we extract α from the spatial saturation, which we show is KPZ–like. All observables in correlations in the entanglement at distances shorter than the correlation length ξðtÞ [Eq. (9)]. For Clifford circuits, we also touch on the entanglement probability distribution. A. Clifford evolution Clifford circuits, or “stabilizer circuits,” are a special class of quantum circuits that play an important role in quantum information theory. As shown by Gottesman and Knill, they can be simulated efficiently, even when the entanglement entropy grows rapidly, by representing the quantum state in terms of stabilizers [87]: see Sec. IV. The time evolution operator for a Clifford circuit belongs to the Clifford group, a subgroup of the unitary group on the full Hilbert space. This group may be generated by a small set of local Clifford gates: the two-site controlled NOT gates [Eq. (31)] and the single-site Hadamard and phase gates R and R [Eqs. (29) and (30)]. For circuits H P FIG. 15. The von Neumann entropy Sðx; tÞ for a system of built from these gates, time evolving the state on L spins up length L ¼ 459, as a function of x, for several successive times to a time t takes a computational time of order Lt, and (t ¼ 340, 690, 1024, 1365, 1707, 2048, and 4096), in the Clifford measuring the entanglement across a given bond in the final evolution. This shows that the state evolves from a product state state takes a time of order L at most. This is in sharp to a near-maximally-entangled one. Prior to saturation the contrast to the exponential scaling that is inevitable for entanglement displays KPZ-like stochastic growth. Sðx; tÞ is in more general circuits. units of log 2. 031016-15 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) FIG. 16. The von Neumann entropy Sðx; tÞ in units of log 2, far from the boundaries, in a system of length L ¼ 1025 at various times (from bottom to top t ¼ 170, 340, 512, and 682) evolved with the Clifford evolution scheme. ξ schematically shows the 1=z typical correlation length Eq. (8), which grows in time like t . the following are measured far from the boundary, in order to avoid finite-size effects associated with saturation. FIG. 17. Top: Growth of the mean entanglement with time for Figure 16 shows successive snapshots for a subregion of the Clifford evolution with only CNOT gates (in units of log 2). a larger system of L ¼ 1025 bonds (times t ¼ 170, 340, The solid red curve is a fit using Eq. (49). The exponent β is found 512, 682, from bottom to top). The maximal slope that can to be β ¼ 0.33  0.01, in agreement with the KPZ prediction appear is 1, in accord with Eq. (3). Note the gradual β ¼ 1=3. Dashed line shows asymptotic linear behavior. Bottom: roughening of the surface and the growing correlation Growth in the fluctuations in the entanglement with time. The length. dashed line shows the expected asymptotic behavior, wðtÞ ∼ t , Figure 17 shows the height and width of the growing with β ¼ 1=3. The fit includes a subleading correction: Eq. (49), surface, with β ¼ 0.32  0.02. Error bars denote the 1σ uncertainty. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðtÞ¼hS ðx; tÞi;wðtÞ¼ ⟪S ðx; tÞ⟫; ð48Þ vN vN β ¼ 0.33  0.01; β ¼ 0.32  0.02: ð50Þ h w for very long times. These quantities are averaged over at Both estimates of β are in excellent agreement with the least 10 realizations. In each realization only the entan- KPZ value β ¼ 1=3. The solid lines in Fig. 17 show the fits glement across the center bond is used (therefore all data (the fit parameters are in Table I). The dashed lines show points are uncorrelated) and the system size is L ¼ at, the slopes corresponding to the expected asymptotic power 1=3 where a is chosen to avoid finite-size effects (see Sec VI B). laws, hðtÞ ∼ t and wðtÞ ∼ t . We obtain estimates β and β of the exponent β by fitting The analysis in Sec. IV implies that v is a well-defined h w the data to the expected forms [cf. Eqs. (6) and (7)]: velocity, and v t is a sharply defined length scale character- izing the range of entanglement in the state. We may β β η h w hðtÞ¼ v t þ Bt ;wðtÞ¼ Ct þ Dt : ð49Þ confirm this by measuring this length scale directly; see the section below. Here, η (with η < β ) captures subleading corrections. Note the small value of the subleading exponent η obtained We find from the fit. This implies that finite time corrections are TABLE I. Summary of all fitting parameters to Eq. (49) used in this section. The errors set the estimated 2σ uncertainty. Evolution type v β B β C η D E h w Clifford 0.1006  0.0001 0.33  0.01 0.66  0.04 0.32  0.02 0.28  0.05 0.08  0.10.16  0.04 Phase 0.133  0.03   0.54  0.04  0.223  0.004   0.168  0.03 Universal 0.202  0.001   0.09  0.005  0.14  0.003   0.36  0.01 031016-16 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) providing further support for KPZ universality in this system. B. Numerics on speeds and scaling forms We argue in Sec. V that in addition to determining the entanglement growth rate, v can also be viewed as a speed. This is the speed of the fictitious particles in Sec. IV. Operationally, the simplest manifestation of this speed is in the saturation behavior of the entanglement. The analytical arguments imply that to leading order (at large t and l) the entanglement across a cut at position l (l ≤ L=2) has the simple scaling form given above in Eq. (24): FIG. 18. The logarithmic derivative of the width dw=d log t u for u< 1 versus time for the Clifford evolution. The universal behavior S ¼ v tfðl=v tÞ;fðuÞ¼ ð51Þ A E E 1=3 1 for u ≥ 1. with exponent t is observed at shorter time scales compared with Fig. 17. This gives simply S ¼ v t for t < l=v , and S ¼ l for A E E A t > l=v . This means that there is no influence of the reduced if we plot the numerical derivative dw=d log t rather boundary at times t < l=v . (See also the numerical results 1=3 than w itself (both quantities scale as t at long times). This in Refs. [39,40], indicating sharp saturation in circuits 1=3 is done in Fig. 18. The data fit well to the t law even at short where interactions between any pair of spins are allowed.) times. This will be useful for the more general dynamics In Fig. 20, we test this result numerically for the Clifford where long times are not available. evolution. We set l ¼ L=2 and plot To complete the check of the two independent KPZ SðL=2;tÞ L=2 exponents, Fig. 19 shows the spatial correlator GðrÞ vs ð52Þ defined in Eq. (9), as a function of separation r, for three v t v t E E successive times. For small r, the correlation grows like a α 9 10 as a function of L, for several values of the time (t ¼ 2 , 2 , r with α ≃ 1=2, in agreement with the KPZ prediction for 11 12 this exponent. For distances r ≫ ξðtÞ, the correlator satu- 2 ,and 2 ). Here, v ¼ 0.1 is taken from the fits to Fig. 17. rates to a value proportional to wðtÞ. The figure gives an According to Eq. (51), this plot should converge for large t to idea of the size of the correlation length ξðtÞ for these times. aplotof fðuÞ against u. The results are in excellent Finally, recent advances in KPZ theory have yielded an agreement with the scaling form, confirming, for the case analytical expression for the full probability distribution of of Clifford circuits, that v is a meaningful velocity. the KPZ height field [56–66]. In Appendix F we show that It is also interesting to compare the entanglement velocity this analytical result compares well with Clifford numerics, v with the butterfly velocity v .Weobtain v from the E B B average spatial extent W of a growing Pauli string (see Sec. IVA) under the unitary Clifford dynamics at time t,as v ¼ W=2t. Remarkably, we find that v ¼ v =2 within B E B 2 1=2 FIG. 19. Correlation function GðrÞ¼h½SðrÞ − Sð0Þ i at time t ¼ 512, 1024, and 2048 for the Clifford evolution, showing FIG. 20. The entropy across the center of the chain (in units of excellent agreement with the KPZ prediction GðrÞ ∼ r with log 2) divided by v t versus L=2v t for various fixed values of t. E E χ ¼ 1=2 in the regime r ≪ ξðtÞ. This plot converges nicely to the scaling form in Eq. (51). 031016-17 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) “ ” FIG. 21. The average size W of a growing Pauli string as a function of time for two protocols, CNOT evolution (subscript “CNOT”) and the full Clifford evolution (subscript “Cliff”). The correspondence with the dashed lines, showing the average entanglement entropy multiplied by four, is consistent with v ¼ v =2. (Taken from a system of size L ¼ 1024.) E B FIG. 22. Top: Growth of the mean entanglement as a function numerical precision for both the CNOT-only Clifford dynam- of time for the universal and phase gate set fitted to Eq. (49), with ics and the “full” Clifford dynamics defined above. This is β set to 1=3. The dashed line shows the expected asymptotic shown in Fig. 21, wherewe plot W starting versus time for the behavior for comparison. Error bars indicate 1 standard deviation two protocols. The initial Pauli strings we consider in this (1σ) uncertainty. simulation are single-site Y operators. [The CNOT dynamics is not ergodic on the space of Pauli strings (unlike the full Universal evolution.—This set of gates, unlike the Clifford dynamics). Nevertheless, any operator grows in size others, is “universal” in the quantum information sense at the same rate v .] We compare W with 4 times the average (any unitary acting on the full Hilbert space of the spin entanglement entropy, 4SðtÞ. The two curves lie on top of chain can be approximated, arbitrarily closely, by a product each other, consistent with v =v ¼ 1=2. E B of gates from this set). The single-site gates include the We also find the same ratio for v =v in the exactly E B eightfold rotations mentioned above, together with the solvable large-q model (Sec. II A 2). However, it is possible Hadamard gate R [. (29)]. R is applied with probability to construct non-fine-tuned random circuits, involving Haar- H H 1=2 and the rotations with probability 1=16 each. random unitaries at finite q, in which the ratio is less than 1=2 Figure 22 shows the height and width hðtÞ and wðtÞ for the [73], so this value is not universal. A natural question is two protocols (averaged over 380 realizations for the phase whether it is generic for random Clifford circuits. evolution and 200 realizations for the universal evolution, and over bonds x with 20 <x< 480). The figure shows fits C. Universal and phase evolution to the forms in Eq. (49) with β and β fixed to the KPZ value h w The phase and universal dynamics take us outside the and η fixed to zero (fit parameters are in Table I). The fits with Clifford realm, and cannot be simulated efficiently on a Eq. (49) are consistent with the data. It is not possible to classical computer (in polynomial time). We give evidence extract precise estimates for β from the slope of the log-log that the correspondence with KPZ continues to hold in this plot of wðtÞ, although for the phase evolution the slope at late more generic situation. However, our results are not as times is in reasonable agreement with the expected KPZ conclusive as in the Clifford evolution as we do not have value, shown by the dashed gray trend line. access to such long times. For an alternative attack on β, we plot the numerical The simulations are performed on spin-1=2 chains of derivative dwðtÞ=d ln t. Recall that in the Clifford case the length L ¼ 500 bonds (501 spins) using the ITensor slope of this quantity (when plotted against time on a log- package [86]. The two types of dynamics are defined as log plot) has smaller finite-size corrections than the slope follows. [The two-site unitaries are always chosen from the for wðtÞ itself. The corresponding plot is shown in Fig. 23, set in Eq. (46); the initial state is taken polarized in the y for times up to t ¼ 25 (averaging over more than 5000 1=3 direction.] realizations). The dashed gray lines are the t trend lines. Phase evolution.—Each single-site unitary is chosen Results for both types of dynamics are in good agreement randomly and uniformly from the set of eightfold rotations with the expected slope β ¼ 1=3. about the z axis in spin space: R ¼ exp ðπinσ =8Þ, with Next, we examine the spatial correlator Eq. (9) in Fig. 24. n ∈ 1; …; 8. For both types of dynamics, the behavior for r ≪ ξðtÞ 031016-18 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) VII. FREE FERMIONS ARE NONGENERIC The growth of entanglement in systems of free particles is highly nongeneric. In the presence of noise, the entan- glement of a system of free particles on the lattice grows pffiffi only as S ∼ t, in contrast to the behavior S ∼ t of generic systems. The case of spatially homogeneous noise has been discussed recently [88]. The basic point is the same when the noise varies in space: the fact that the single-particle wave functions spread diffusively in the presence of noise implies that the entanglement cannot be larger than pffiffi Oð tÞ [88]. FIG. 23. The logarithmic derivative of the width dw=d log t As a concrete example, consider a short-range hopping versus time for the phase and universal evolution protocols. For Hamiltonian for free fermions, 1=3 comparison, we plot the universal behavior with exponent t in gray (dashed line). (The derivative is calculated using three data HðtÞ¼ H ðtÞc c ; ð53Þ ij j points. Errors are estimated from maximal and minimal slopes ij obtained within 1 standard deviation from the averaged data points.) with noisy matrix elements H ðtÞ. For simplicity, take the ij initial state to consist of particles localized at sites i ∈ S for some set S; for example, we could take S to consist of all the even-numbered sites: jΨð0i¼ c j0i: ð54Þ i∈S Under the evolution, each creation operator evolves into a superposition of creation operators, † † ðiÞ c → ψ ðj; tÞc ; ð55Þ i j ðiÞ where ψ ðj; tÞ is the solution of the time-dependent Schrödinger equation for a particle initially localized at ðjÞ i. In the absence of noise, ψ spreads ballistically, but in the presence of noise. it spreads only diffusively. The fact pffiffi that each creation operator is spread out over only Oð tÞ sites after a time t immediately implies that the mean pffiffi entanglement is at most of order t. (See also Ref. [88].) Note, however, that this argument does not tell us how large the fluctuations are. (Random unitary evolution of a single wave packet is discussed in Ref. [89].However,wemust 2 1=2 FIG. 24. Correlation function GðrÞ¼h½SðrÞ − Sð0Þ i at consider the full many-body wave function, since the three values of the time for the phase (top) and universal (bottom) formalism of Ref. [90] for the free-fermion density matrix gate sets, showing good agreement with the KPZ exponent value shows that the initially occupied orbitals do not simply α ¼ 1=2. contribute additively to the entanglement.) pffiffi We have confirmed numerically that hSi ∝ t for a agrees well with the KPZ exponent value α ¼ 1=2 at the noisy 1D hopping model, using the formalism of Ref. [90] largest available time. to construct the reduced density matrix. This is much The very long times accessible in the Clifford simulation slower than the linear-in-time growth of generic interacting allow us to establish KPZ exponents with high accuracy pffiffi models. The t scaling should apply for free fermions in there. For the more generic dynamical rules we cannot any number of dimensions. In 1D it also applies to certain reach the same level of precision, but, nevertheless, the noisy spin models via the Jordan-Wigner transformation: KPZ exponent values are compatible with the data. for example, the transverse field XY model: Next, we briefly discuss a fine-tuned situation in which entanglement dynamics are not KPZ-like: namely, when y y x x HðtÞ¼ fJ ðtÞ½σ σ þ σ σ þ h ðtÞσ g: ð56Þ the system is made up of free particles. Then, in Sec. VIII, i i i iþ1 i iþ1 i we move to higher dimensions. 031016-19 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) However, any generic perturbation to the spin chain spoils the free-fermion correspondence. We then expect the generic KPZ behavior to reassert itself. VIII. HIGHER DIMENSIONS We discuss several ways of thinking about entanglement growth in 1D. One of these, the directed polymer picture, generalizes naturally to higher dimensions: the polymer is simply replaced by a d-dimensional membrane embedded in (d þ 1)-dimensional space-time. As in 1D, we think of this membrane as a coarse-grained version of a minimal FIG. 26. Minimal membrane for a disk-shaped region in d ¼ 2. cut bisecting a unitary circuit. The membrane is subject to pinning by “disorder” in space-time arising from the dynamical noise. See Fig. 25 for the two-dimensional case. A. Universal fluctuations of SðtÞ in noisy systems We can explore two kinds of questions using this picture. Consider the entanglement SðtÞ for a region A whose First, we can examine universal properties that are specific to boundary ∂A has length or area j∂Aj. In the d ¼ 2 case, the noisy scenario: as in 1D, fluctuations are governed by shown in Fig. 25, j∂Aj is the length of the spatial boundary. universal exponents. Second, we can calculate leading-order Neglecting fluctuations, the “world volume” of the minimal properties of SðtÞ that do not involve fluctuations and that membrane scales as j∂Aj × t. This gives the leading scaling are, therefore, likely to be valid even in the absence of noise, of the membrane’s energy and, hence, of the entanglement. i.e., for dynamics with a time-independent Hamiltonian. As in 1D, subleading terms encode universal data. We now In higher dimensions the behavior of SðtÞ has nontrivial consider these terms. dependence on the geometry of the region for which we The pinning of a membrane or domain wall by disorder calculate the entanglement. We suggest the “minimal mem- is well studied [47,91–94] (a brief summary is in brane in space-time” as a simple and general heuristic for Appendix E). Translating standard results into the language such calculations. Below, we discuss the case of a spherical of the entanglement in a d-dimensional noisy quantum region (Sec. VIII B) and contrast our results with an alter- system, we find that in both d ¼ 1 and d ¼ 2 there is a native simple conjecture. For other toy models for entangle- unique dynamical phase with nontrivial critical exponents. ment spreading, see Refs. [10,23]. The same is true for continuum systems [more precisely, for Denoting the region for which we wish to calculate the systems with continuous (statistical) spatial translational entropy by A, and its boundary by ∂A, the membrane lives symmetry] in d ¼ 3. However, if a lattice is present, two in a space-time slice of temporal thickness t, and terminates stable phases (and thus a dynamical phase transition) are possible in d ¼ 3; one with nontrivial exponents and one at ∂A on the upper boundary of this time slice; see Fig. 25. with trivial ones. In the trivial phase, the membrane is For simple shapes and for times shorter than the saturation “smooth” and is pinned by the lattice. In the nontrivial time, the membrane also has a boundary on the lower slice, phases, the membrane is instead pinned by disorder in a as shown in Figs. 25 and 26. In this section, we focus on “rough” configuration. We discuss the nontrivial phases entanglement growth prior to saturation. (which are the only ones possible in d< 3 and for continuum systems in d ¼ 3). Generally, fluctuations have a weaker effect in higher dimensions than in 1D. For simplicity, take a quantum system which is infinite in one direction and of size L in the other d − 1 directions, and consider the entanglement for a cut perpendicular to the infinite direction. Since A and its complement are both infinite, SðtÞ will grow indefinitely for this geometry. However, there are two regimes, t ≲ L and t ≫ L (here, we drop a dimensionful prefactor). For times t ≲ L (see Appendix E for details): d−1 θþ1−d hSðtÞi ¼ L ðv t þ Bt þÞ; ð57Þ 2 1=2 ðd−1Þ=2 θ−ðd−1Þ=2 ⟪SðtÞ ⟫ ∝ L t ; ð58Þ where the exponent θ is defined below. This reproduces the FIG. 25. Minimal membrane picture for the entanglement of two regions in d ¼ 2. 1D result with θ ¼ β. Note that when d> 1, fluctuations 031016-20 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) are suppressed with respect to the mean by a factor of However, this simplifies in the limit of interest. We first 1=2 send t; R → ∞ with t=R fixed. In this limit r_ðt Þ remains j∂Aj : distant regions of the boundary give rise to finite, but the curvature terms become negligible (see, for essentially independent fluctuations which add up incoher- example, the explicit solution below), so we can write ently. In the opposite regime t ≫ L, the temporal dimen- E ¼ Eðr_Þ. Now we make the second approximation that sion of the membrane is much larger than its spatial dimensions, so there is a crossover to the 1D directed t=R is small, meaning that we can keep only the Oðr_ Þ polymer problem. However, the exponent of the higher- correction. dimensional problem appears in the universal L depend- The boundary condition at the top of the space-time slice is rðtÞ¼ R. We consider times prior to saturation, so the ence of the growth rate: membrane also has a free boundary on t ¼ 0. In the d−1 θ−1 SðtÞ¼ðvL þ wL Þt þ  : ð59Þ relevant limit, its energy is 1=3 (The higher corrections will include the t term associated 0 0 0 2 E ¼ 2πv dt rðt Þ½1 þ ar_ðt Þ þ : ð62Þ with the 1D universality class.) Numerically, the exponent is θ ¼ 0.84ð3Þ in d ¼ 2 and θ ¼ 1.45ð4Þ in d ¼ 3 [94]. The Minimal energy requires the boundary condition r_ð0Þ¼ 0. subleading exponent in Eq. (57) is negative for d> 1,so When t=R is small, we may expand in 1=R. This gives this correction may be hard to observe numerically. 0 2 02 rðt Þ ≃ R − ðt − t Þ=ð4aRÞ, as illustrated in Fig. 26. The corresponding entropy is B. Minimal membrane picture for dynamics without noise SðtÞ¼ 2πv Rt 1 − þ  : ð63Þ In higher dimensions we can ask how SðtÞ depends on 2 12aR the geometry of region A when this geometry is nontrivial. This calculation generalizes trivially to higher dimensions, Interestingly, the membrane picture makes predictions where the correction is of the same order. Corrections due about this that do not involve the noise-induced fluctua- to fluctuations come in with negative powers of t, and are tions, and that are likely also to be valid for Hamiltonian negligible in the limit we are discussing. dynamics without noise (with the replacement S → S=s eq Note that the first correction in the brackets in Eq. (63) is we discuss in Sec. V). of order ðt=RÞ , and not of order t=R. This result differs As an instructive special case, take A to be a disk-shaped from what one might naively have expected if one guessed region of radius R in d ¼ 2. (A ball in higher dimensions is that at time t an annulus of width v~ × t inside the disk is precisely analogous.) We assume continuous rotational entangled with the outside, where v~ is a tsunami velocity. symmetry, at least on average. At short times, the leading This picture gives an entropy proportional to the area of the scaling of the entanglement is SðtÞ ≃ 2πv Rt, since the annulus, world-sheet area of the membrane is approximately 2πR × t. However, there are corrections to this arising vt ~ 2 2 SðtÞ∼πR − πðR − vt ~ Þ ¼ 2πRvt ~ 1 − ; ð64Þ from the curvature of ∂A. We consider the limit of large t and large R with a fixed ratio t=R. In this regime, the effects of fluctuations may be leading to a negative correction of order t=R. The difference neglected, and instead the energetics of the membrane are between Eqs. (63) and (64) also indicates that a picture in terms of independently spreading operators is misleading, determined by deterministic elastic effects. We write the in agreement with what we find in 1D. energy of the membrane as It is interesting to note that in the regime where t=R is of E ¼ d SE; ð60Þ order 1, the full r_ dependence of Eðr_Þ plays a role. This suggests that an infinite number of nonuniversal parameters where d S is the membrane’s area element and E is its enter the expression for SðtÞ in this regime, and that there is “energy” density. We get SðtÞ by minimizing E with no general, universal scaling form for the entanglement of a appropriate boundary conditions. sphere in d> 1. However, we do expect saturation to Next, we Taylor expand E in terms of local properties remain discontinuous, as in 1D [Eq. (24)], occurring via a of the membrane. For a flat “vertical” membrane (i.e., with transition between an optimal membrane configuration that normal perpendicular to the t axis), E ¼ v . In general, reaches the bottom of the space-time slice and one (with however, E depends on the angle φ by which the surface E ¼ πR ) that does not. locally deviates from verticality, as well as, for example, the IX. OUTLOOK local curvatures κ and κ in the spatial and temporal s t directions. Using rotational symmetry to parametrize the Quantum quenches generate complex, highly entangled membrane by the radius rðt Þ, states whose dynamics cannot usually be tracked explicitly. For this reason, analytical approaches to quenches have 2 2 2 2 4 d SE ¼v rdθdtð1þar_ þbκ þcκ þcr_ þÞ: ð61Þ E t s typically relied on additional structure in the quantum 031016-21 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) dynamics: for example, integrability, or absence of inter- As a more speculative question in the domain of actions, or conformal invariance. This paper instead studies noisy dynamics, we may ask whether there exist time- dynamics that are as unstructured as possible. We propose independent Hamiltonians that show KPZ entanglement that noisy dynamics are a useful toy model for quantum fluctuations, despite the absence of explicit noise, in some quenches in generic (nonintegrable, non-conformally- dynamical regimes. We emphasize that this seems unlikely invariant) systems. on asymptotically long time scales for a generic system Many of our results are, of course, specific to noisy (since local reduced density matrices and observables dynamics: in particular, the emergence of KPZ behavior at will eventually thermalize), but it may hold on intermediate long wavelengths in 1D, and the detailed pictures for time scales in certain systems in which some degrees of entanglement growth afforded by the “KPZ triumvirate.” freedom act effectively as chaotic classical variables and But we suggest that some of our heuristic pictures apply to provide effective noise. non-noisy entanglement growth as well (with the replace- In the text, we discuss only initial states with area-law ment S → S=s mentioned above). We propose a general entanglement. A natural extension is to initial states with, for eq directed polymer picture or minimal membrane picture for example, submaximal volume-law entanglement. The natu- the scaling of the entanglement and mutual information ral expectation, say, in 1D, is that the directed polymer picture (Secs. V, VIII B) and we use the operator spreading picture extends to this case if we glue the unitary circuit to a tensor to clarify the meaning of the “entanglement velocity” and network representation of the initial state. Then the entropy its distinction from the operator spreading velocity Sðx; tÞ would include a fluctuating part with KPZ exponents (Sec. V). “Thermalization is slower than operator spread- together with a contribution from the initial state. Another ing” in generic 1D systems (i.e., in general ,v is smaller natural direction to explore is the role of conserved quantities. than v ): by contrast, this is not true in 1+1D CFTs [2],or Turning to higher dimensions, it would also be useful to in certain toy models [23]. It would be interesting to make test the higher-dimensional membrane pictures of Sec. VIII, more detailed comparisons with holographic models [8]. perhaps exploiting Clifford circuits to reduce the numerical Many interesting questions remain. First, within the difficulty of higher-dimensional dynamics. realm of noisy systems, an analytical treatment for the There are also many further questions regarding deter- regime with weak noise would be desirable, i.e., for ministic systems for which the tools we introduce here may dynamics of the form give insight. For example, the forthcoming Ref. [77] will discuss entanglement growth in disordered or inhomo- HðtÞ¼ H þ λH ðtÞ; ð65Þ 0 1 geneous spin chains from the point of view of surface where H is a time-independent many-body Hamiltonian, growth, while Ref. [73] will give results for the spreading of H ðtÞ represents noise, and λ is small. Our conjecture is that quantum operators. KPZ exponents apply for any nonzero value of λ [unless HðtÞ is fine-tuned]—i.e., that there is no universal distinction ACKNOWLEDGMENTS between continuous time dynamics and quantum circuits. (Note that there is no distinction between these two cases at We thank M. Kardar for useful discussions. A. N. and J. R. the level of conservation laws: once noise is added, energy is acknowledge the support of fellowships from the Gordon not conserved even in the continuous time case.) However, and Betty Moore Foundation under the EPiQS initiative our derivations and numerics correspond, roughly speaking, (Grant No. GBMF4303). A. N. also acknowledges support to the large λ regime. Perhaps the opposite regime could be from EPSRC Grant No. EP/N028678/1. S. V. is supported addressed using a more explicit renormalization group (RG) partly by the KITP Graduate Fellows Program and the U.S. treatment, although it is not obvious how to set this up. DOE Office of Basic Energy Sciences, Division of Materials Such a RG treatment might also shed light on the nature of Sciences and Engineering under Award No. DE-SC0010526. the entanglement spectrum or, equivalently, the dependence J. H. is supported by the Pappalardo Fellowship in Physics of S ðtÞ on the index n. While we believe that all the Rényi at MIT. entropies execute KPZ growth in the presence of noise, we have not pinned down the n dependence of the various APPENDIX A: GROWTH OF HARTLEY constants. The solvable models suggest that the leading- ENTROPY S IN 1D order behavior may be independent of n at large times. What is the appropriate scaling form for the spectrum? Limited Consider a one-dimensional quantum spin chain of local time scales prevent us from addressing this numerically Hilbert space dimension q, prepared initially in a product (except for Clifford circuits, where all the S are trivially state, and apply a sequence of random unitaries that couple two neighboring spins. The location of the local unitary at a equal). (The entanglement spectrum is one window on the structure of the quantum states generated by the random given time step is arbitrary. In the following we fix the dynamics. We can also ask in what ways these states differ location of the unitary, but take it to be Haar random. We prove that in this situation the Hartley entropy S from ground states of random Hamiltonians, when the amount of entanglement is similar.) generically (i.e., with probability 1) obeys 031016-22 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) S ðx; t þ 1Þ¼ min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1; ðA1Þ reduces to seven different cases excluding invalid ones 0 0 0 and those related by reflection, one easily checks that there if a unitary is applied at the bond x. The logarithm is of is always a choice among I, W, E that makes Eq. (A1) true. base q. Let us treat three exemplary cases here. If s and s are L R This formula can be interpreted in matrix-product-state entangled at time t, then Sðx − 1;tÞ¼ Sðx þ 1;tÞ and language. If d is the minimal value of the local bond Sðx; tÞ¼ Sðx − 1;tÞþ 1, so one chooses the identity I. dimension required for an exact matrix-product-state rep- If s is entangled with a spin on the left of s and s L L R resentation of the state, then S ðxÞ¼ log d . A heuristic is entangled with a spin on the right of s , then 0 x R parameter-counting argument for the local bond dimension, Sðx − 1;tÞ¼ Sðx þ 1;tÞ¼ 1 þ Sðx; tÞ. One chooses the given in Appendix A3, suggests Eq. (A1). swap W to obtain Sðx; t þ 1Þ¼ Sðx; tÞþ 2.If s and s L R However, a more rigorous proof is necessary as such are both unentangled, then one applies the entangling heuristic arguments can fail. In particular, one might unitary E to obtain Sðx;t þ 1Þ− 1 ¼ Sðx;tÞ¼ Sðx− 1;tÞ¼ naively conjecture a stronger statement: namely, that for Sðx þ 1;tÞ. any state at time t, if the unitary at bond x is Haar random, For the second part, recall that for any bipartite state then Eq. (A1) is true with probability 1. This conjecture is jψi¼ M jiijji; ðA3Þ i;j false; a counterexample is given in Appendix A2.We now i;j give a proof of Eq. (A1). the number of nonzero Schmidt coefficients is equal to the number of nonzero singular values of the matrix M, which 1. Proof of Eq. (A1) is nothing but the rank of M. For any positive integer r, the Our genericity proof consists of two parts. First, we show rank of M is smaller than r if and only if every r × r that given locations of unitaries, there exist certain unitaries submatrix has determinant zero; i.e., all r × r minors such that at each time step Eq. (A1) is true. Second, we vanish. Thus, a bipartite state jψi having Hartley entropy show the negation of Eq. (A1), (log of rank of M) strictly smaller than log r is expressed by a system of polynomial equations on the coefficients of jψi. S ðx; t þ 1Þ < min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1; ðA2Þ 0 0 0 If jψi is given by U …U U j0i, where j0i is a fixed t 2 1 happens if and only if a system of polynomial equations in product state, then the coefficients are some polynomials of the entries of the unitaries is satisfied. (The inequality “>” the entries of the unitaries U , and hence the equations that never holds, as we note in the main text.) By the first part of expresses vanishing determinants are polynomial equations the proof, the zero locus of these polynomial equations in the entries of the unitaries. does not cover the entire set of unitaries. Therefore, it is Our claim Eq. (A1) completely determines the Hartley only a submanifold of strictly smaller dimension, which entropy based on the location of unitaries, and therefore the implies it has measure zero. spatial configuration of the unitaries tells us which minors we For the first part, it is sufficient to consider only three should check. Namely, the size r of the minors we turn into types of local unitaries: the identity I, the swap gate W, the polynomial equations is given by (the exponential of) and a unitary E with the property that it turns a pair of the right-hand side of Eq. (A2). In other words, given a spatial −1=2 unentangled polarized spins, j11i, into q jiii,a i¼1 configuration of unitaries, the polynomial equations that maximally entangled state. Without loss of generality we express Eq. (A2) are determined. The polynomial equations may take the initial product state to be the polarized are over tL variables, and the actual number of equations is state j…1111…i. much larger yet finite. We do not need explicit expressions We show that using these three types of unitaries at the for these polynomials, only the fact of their existence. These given locations, one can construct a state whose entangle- polynomials might a priori read 0 ¼ 0; i.e., they could be ment entropy is given by Eq. (A1). Since Eq. (A1) defines the trivially satisfied. In that case, the solution to the polynomial entropy inductively, we only have to show it inductively, too. equation would be the entire set of unitaries, and Eq. (A1) At t ¼ 0, all the spins are unentangled, so we can simply could never be satisfied. However, we just showed in the first choose E for every designated location. Clearly, Eq. (A1) is part that this cannot happen because there exists a choice of satisfied. At later times, if we do not apply E except on unitaries for which Eq. (A1) is satisfied. This implies that the an unentangled pair of spins, then a spin can be either polynomial equations are nontrivial and define a measure unentangled or maximally entangled with a single other zero subset of the entire set of unitaries. This completes the spin. Therefore, at time t> 0, the spin s that is immedi- genericity proof. ately left to the bond x can be (i) unentangled, (ii) entangled with a spin to the left of s , (iii) entangled with the spin s L R 2. Counterexample to the stronger conjecture that is immediately to the right of the bond y,or (iv) entangled with a spin to the right of s . We show above that Eq. (A1) holds when all unitaries are chosen generically and the initial state is a product state. These are exclusive possibilities, and similarly s has four options. Enumerating all 16 cases, which in fact Naively one might make the stronger conjecture: that the 031016-23 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) update rule Eq. (A1) holds whenever a generic unitary U is U and U , we find that all (5 × 5) minors vanish, implying 0 1 applied to an arbitrary—possibly fine-tuned—state jΨi. that V has rank at most 4. Therefore, for this nongeneric We construct an explicit jΨi, which is a counterexample to initial state, this stronger conjecture. Consider four degrees of freedom ABCD. The spins B S ðx; t þ 1Þ ≠ min½S ðx − 1;tÞ;S ðx þ 1;tÞ þ 1: ðA9Þ 0 0 0 and C have dimension 2 each, and A and D have dimension 3 each. (To conform with our consideration of spin chains, 3. Parameter-counting argument the subsystems A and D should be regarded as subspaces of two or more spin-1=2’s.) The most general form of a Consider a 1D state jΨi in a matrix product representa- quantum state on ABCD is tion. Labeling the states of the qubits (spins) by σ; σ ; … running from 1 to q, 2 1 X X XX T jaijbijcijdi: ðA4Þ abcd σ 0σ 0 jΨi¼ ð…A A …Þj…σσ …i: ðA10Þ a ;a a ;a x−1 x x xþ1 a;d¼0 b;c¼0 fσg fag We consider T ¼ T δ , i.e., C is in j0i, where abcd abd c;0 Since the state is not translationally invariant, we allow 0 1 0 1 the bond dimension d to vary from bond to bond 01 0 00 1 (a ¼ 1; …;d ). In an efficient representation, d is equal x x x B C B C 0 0 T ¼ @10 0A ;T ¼ @00 0A : a0d a1d to the rank of the reduced density matrix for a cut at x: 00 0 10 0 ad ad S ðxÞ d ¼ q : ðA11Þ ðA5Þ We ask how S ðxÞ changes when we apply a unitary U to (This does not give a normalized state, but we are only 0 the two spins, σ and σ , either side of bond x. This effects concerned about ranks.) the change (repeated indices are summed): The Hartley entropy for the cut A=BCD is simple to compute. As we remark in the previous section, it is the 0 0 σ 0σ τ 0τ 0 0 A A → U A A : ðA12Þ rank of the coefficient matrix. Interpreting this matrix as a a ;a a ;a σσ ;ττ a ;a a ;a x−1 x x xþ1 x−1 x x xþ1 linear map, the rank is the dimension of the image of the map from BCD to A. The image is precisely the linear span To update the matrix product representation, we must find 0 0 of columns of T and T . They have three linearly ~ ~ new matrices A and A which satisfy a0d a1d independent columns, implying that the Hartley entropy for σ 0σ A=BCD is log 3. Similarly, the rank of the coefficient τ 0τ 2 ~ ~ A A ¼ U 0 0A A : ðA13Þ a ;a a ;a σσ ;ττ a ;a a ;a x−1 x x xþ1 x−1 x x xþ1 matrix for ABC=D is the dimension of the linear span of the 0 0 rows of T and T , which reads 3. That is, the Hartley a0d a1d 0 ~ ~ In order to solve this equation for A and A , it will generally entropy for ABC=D is log 3. be necessary to increase the bond dimension at x to a new If Eq. (A1) were to be true for the generic choice of Haar 0 0 value d . Naively, the necessary value of d will generically x x random unitary on BC, then we should be able to find a be determined by equating the number of independent unitary on BC such that equations in Eq. (A13) with the number of degrees of ~ ~ freedom in A and A . (However, the previous section shows S ðAB=CDÞ¼ log 3 þ 1 ¼ log 6: ðA6Þ 0 2 2 that this expectation can fail for certain choices of A and A .) We show this cannot hold. Applying the unitary U on BC The number of equations is q d d , since this is the x−1 xþ1 the state, we obtain number of possible values for the external indices in 0 0 Eq. (A13). The numbers of degrees of freedom in A and 0 0 0 0 0 0 U T ¼ U T þ U T ; ðA7Þ b c ;bc abcd b c ;00 b c ;10 a0d a1d |fflfflffl{zfflfflffl} |fflfflffl{zfflfflffl} 0 0 0 2 b;c A are qd d and qd d , respectively. However, d of x−1 x xþ1 x x U U 0 1 these are redundant, because the state is unchanged by the 0 0 σ σ 0σ 0σ −1 ~ ~ ~ ~ transformation A → A M, A → M A , with M an where U and U are 2 × 2 matrices. The coefficient matrix 0 1 0 0 arbitrary d × d matrix. Equating the number of equations for the cut AB=CD is then x x with the number of independent degrees of freedom gives 0 0 V ¼ U ⊗ T þ U ⊗ T ; ðA8Þ 0 1 0 1 0 0 ðd − qd Þðd − qd Þ¼ 0: ðA14Þ x x−1 x xþ1 whose rank should be 6 if S ðAB=CDÞ¼ log 6. Computing all the minors of the 6 × 6 matrix V for arbitrary matrices Choosing the smallest solution, 031016-24 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) d ¼ q minfd ;d g: ðA15Þ Z jψi¼þjψi; ðC2Þ x x−1 xþ1 i This agrees with Eq. (A1) since S ðxÞ¼ log d . for all i ¼ 1; …;n. The condition that jψi is nonzero and 0 x unique is equivalent to the condition that the operator APPENDIX B: HAAR AVERAGE FOR Trρ g ðC3Þ Let ρ ðtÞ be the reduced density matrix for a cut at x, jGj g∈G obtained by tracing out the spins to the left of the cut. Each index on this matrix labels a configuration of the spins to is a projector of rank 1 [96,97]. Since jψi is in the image of the right of the cut. Let us temporarily label these spins this projector, we see 1; 2; …, and let the spin immediately to the left of the cut be denoted 0. The indices on the reduced density matrices are then jψihψj¼ g: ðC4Þ jGj g∈G σ ;σ ;σ ;… 0 1 2 σ ;σ ;… σ ;… 1 2 2 ρ ðtÞ ; ρ ðtÞ ; ρ ðtÞ : ðB1Þ μ ;μ ;μ ;… x−1 0 1 2 x μ ;μ ;… xþ1 μ ;… 1 2 2 Since this is a normalized pure density matrix, its trace is In the following we assume that repeated indices are equal to 1. But a Pauli matrix has the property that it is summed. After applying a unitary on bond x, traceless. Therefore, only the identity element on the right has nonzero trace: 0 0 σ ;σ ;σ ;… σ ;σ ;…  0 1 1 2 ρ ðt þ 1Þ ¼ U 0 0 U ρ ðtÞ : 0 0 0 0 x μ ;μ ;… τσ ;σ σ x−1 1 τμ ;μ μ μ ;μ ;μ ;… 1 2 0 1 1 0 1 0 1 1 1 2 ⊗n n 1 ¼ dimðC Þ ¼ 2 : ðC5Þ Let us average Trρ ðt þ 1Þ over the choice of unitary, for a jGj jGj fixed initial state: From this expression, it is straightforward to obtain 0 0 00 00 σ ;σ ;σ ;… μ ;μ ;μ ;… 2 2 expressions for reduced density matrices. Suppose the 2 0 1 0 1 hTrρ ðt þ 1Þ i¼ ρ ðtÞ ρ ðtÞ 0 0 00 00 x x−1 x−1 μ ;μ ;μ ;… σ ;σ ;σ ;… 2 2 0 1 0 1 n-qubit system is partitioned into two complementary × hU 0 0 U U 00 00U i: 0 0 00 00 regions A and B. Tracing out B,wehave τσ ;σ σ νμ ;μ μ 1 τμ ;μ μ 1 νσ ;σ σ 0 1 1 0 1 1 0 1 0 1 The Haar average for four elements of a UðdÞ matrix (here, ρ ¼ Tr ðgÞ: ðC6Þ A B d ¼ q , and each index on U represents a pair of spin g∈G indices) is Tr ðgÞ is nonzero if and only if the tensor component hU U 0 0U U i 0 0 a;b a ;b c;d Haar corresponding to B is identity, in which case c ;d ¼ fδ δ 0 0δ δ 0 0 þ δ 0δ 0 δ 0δ 0 g jBj a;c a ;c b;d b ;d a;c a ;c b;d b ;d Tr ðgÞ¼ 2 gj ; ðC7Þ d − 1 where gj denotes the tensor components of g correspond- − fδ δ 0 0δ 0δ 0 þ δ 0δ 0 δ δ 0 0g : ðB2Þ a;c a ;c b;d b ;d a;c a ;c b;d b ;d ing to A. The set of all gj such that Tr ðgÞ ≠ 0 can be regarded as a subgroup of G, which we denote by G . The index contractions give the result in the text: The formula for ρ now reads 2 2 −1 2 2 hTrρ ðt þ 1Þ i ¼ qðq þ 1Þ ðTrρ þ Trρ Þ: x Haar jBj x−1 xþ1 X X 2 jG j 1 ρ ¼ g ¼ g: ðC8Þ jAj ðB3Þ 2 jG j g∈G g∈G A A APPENDIX C: ENTANGLEMENT ENTROPY It is immediate that ρ is proportional to a projector since it OF STABILIZER STATES is a sum over a group. It follows that the rank of ρ is equal jAj A stabilizer state is a state of an n-qubit system defined to 2 =jG j. In particular, the (Rényi or von Neumann) by a complete set fg ; …;g g of commuting tensor entropy of ρ with base-2 logarithm is 1 n A products of Pauli matrices through equations g jψi¼þjψi: ðC1Þ Sðρ Þ¼jAj − log jG j: ðC9Þ i A 2 A The group generated by fg ; …;g g is naturally called a 1 n The subgroup G has period 2, and therefore log jG j is an stabilizer group, and denoted by G [87,95]. A trivial A A example is the all-spin-up state, defined as integer, which is equal to the number of independent 031016-25 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) stabilizers supported only on A. This expression for the and (30). respectively. (The Hadamard gate corresponds to entanglement entropy has also appeared in Refs. [82,83]. swapping the X and Z vectors, while the phase gate Now, regard the stabilizer group G as a binary vector space corresponds to adding the X vector to the Z vector.) V by ignoring the overall phase (sign) factors. Let Π be the The von Neumann entropy in units of log 2 and the truncation map retaining the components corresponding to corresponding width averaged over ∼2 × 10 realizations the region A, and similarly let Π be the truncation map for (except for the last data point, where ∼2 × 10 realizations B ¼ A. It is routine to check that V decomposes as V ⊕ are used for the average) are plotted in Fig. 27. The fit to 0 0 V ⊕ V for some subspace V ⊆V, where V and V are the the KPZ universal form Eq. (49) gives β ¼ 0.2  0.15 B A B h spans of stabilizers supported only on A and B, respectively. and β ¼ 0.3  0.04. We also obtain v ¼ 0.194  0.001, w E Both the truncation maps are injective on V . It follows that B ¼ 0.4  0.2, C ¼ 0.4  0.1, D ¼ 0.4  0.6, and η ¼ S ¼jBj − dim V ¼ dim ðΠ VÞ − jAj. This completes −0.4  0.8. These results are consistent with the KPZ A F B F A 2 2 universality and with the data presented in Fig. 17. the proof of Eq. (33). APPENDIX E: DETAILS OF STATISTICS APPENDIX D: NUMERICS FOR FULL OF MEMBRANES CLIFFORD EVOLUTION The exponents governing the membrane problem are In Sec. VI A, we present numerical results for random traditionally denoted θ and ζ, and are related by 2ζ − θ ¼ unitary evolution using only the CNOT gates Eq. (31). Here, 2 − d [47]. Consider a patch of the membrane with linear we present similar analysis using the full set of generators for dimensions scaling as l. This includes both its temporal the Clifford group, showing that the additional gates do not dimension and its internal spatial dimensions: after a modify the universal behavior. The additional single-site rescaling of time, the membrane is statistically isotropic gates are the Hadamard and phase gates defined in Eqs. (29) on large scales. The mean “energy” of this patch of d θ membrane scales as l þ const × l , with fluctuations of order l . The length scale for wandering of the membrane in the transverse direction is of order l . The numerical results we quote in the main text are in good agreement with an epsilon expansion about d ¼ 4, which gives ζ ≃ 0.208ð4 − dÞ [93] (see also Ref. [98]). The scaling forms for the entanglement we discuss in the text are easily found by regarding the membrane as made up of patches of appropriate linear size: size t for Eqs. (57) and (58) and size L for Eq. (59). Note that the geometry of the membrane, including the transverse length scale (which is Δx ∼ t for the regime t ≲ L), determines the dimensions of the space-time region around ∂A for which the final entanglement is sensitive to small changes in HðtÞ, i.e., in the history of the noise. APPENDIX F: ENTANGLEMENT PROBABILITY DISTRIBUTION As we mention in the main text, a remarkable recent advance in KPZ theory has been the derivation of the full universal probability distribution for the height of the surface at fixed position and fixed large time [56–66]; see Refs. [99–101] for reviews. In our case, this height corresponds to the entanglement S across a cut in a system undergoing noisy unitary dynamics. One may separate out FIG. 27. Top: Growth of the mean entanglement in units of log the nonuniversal growth rate v , and the nonuniversal 2 as a function of time for the random Clifford evolution (only constant D governing the scale of fluctuations, by writing CNOT gates). The red solid curve is a fit using the form Eq. (49). Dashed line shows asymptotic linear behavior. Bottom: Growth S ¼ v t þ Dt χ: ðF1Þ in the fluctuations in the entanglement with time. The exponent β E is found to be β ¼ 0.3  0.04, in agreement with the KPZ prediction β ¼ 1=3. The dashed line shows the expected asymp- The rescaled random variable χ is then expected to have totic behavior, wðtÞ ∼ t with β ¼ 1=3. a universal probability distribution PðχÞ at late times. 031016-26 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) states (Sec. IX) indicates that the lower end point of the polymer is then no longer free, and instead favors x ¼ 0.] In Fig. 28, we fit numerical data for the probability distribution of S to the expected Tracy-Widom form and, for comparison, to a Gaussian distribution. The data are for the “full” Clifford dynamics (defined in Sec. VI A) at time t ¼ 2048. Each fit involves two parameters, corresponding to the mean and the variance. The Tracy-Widom distribu- tion used is the theoretically expected one with β ¼ 1,but, in fact, the present data do not allow us to discriminate between TW and TW . The Tracy-Widom distribu- β¼1 β¼2 tion is a much better fit to the data than the Gaussian, as quantified in the caption of Fig. 28. This is further confirmation of KPZ universality in the Clifford case. A more detailed investigation of the probability distribution is beyond the scope of this paper, in view of finite-time effects at the accessible time scales. [1] P. Calabrese and J. Cardy, Evolution of Entanglement Entropy in One-Dimensional Systems, J. Stat. Mech. (2005) P04010. [2] P. Calabrese and J. Cardy, Entanglement Entropy and Conformal Field Theory, J. Phys. A 42, 504005 (2009). [3] P. Calabrese and J. Cardy, Quantum Quenches in 1 þ 1- Dimensional Conformal Field Theories, J. Stat. Mech. (2016) 064003. [4] C. T. Asplund, A. Bernamonti, F. Galli, and T. Hartman, FIG. 28. Observed probability distribution for the entanglement Entanglement Scrambling in 2D Conformal Field Theory, entropy across the central bond of chain of length L ¼ 2048 at J. High Energy Phys. 09 (2015) 110. time t ¼ 2048 under the full Clifford dynamics, fitted to two [5] V. Hubeny, M. Rangamani, and T. Takayanagi, A Covar- probability distributions. Top: Best fit to the Tracy-Widom (TW) iant Holographic Entanglement Entropy Proposal, J. High distribution with β ¼ 1. A fit to TW with β ¼ 2 is not shown, but Energy Phys. 07 (2007) 062. is indistinguishable at the scale of the figure. Bottom: Best fit to [6] J. Abajo-Arrastia, J. Aparício, and E. López, Holographic the Gaussian. Clearly, the Tracy-Widom distribution fits the data Evolution of Entanglement Entropy, J. High Energy Phys. better than Gaussian, as the latter shows systematic deviation. 11 (2010) 149. 2 −4 The 1 − R values for the fits are 2.1 × 10 for TW , [7] T. Hartman and J. Maldacena, Time Evolution of Entan- β¼1 −4 −3 glement Entropy from Black Hole Interiors, J. High Energy 2.0 × 10 for TW , and 1.6 × 10 for Gaussian. β¼2 Phys. 05 (2013) 014. [8] H. Liu and S. J. Suh, Entanglement Tsunami: Universal This probability distribution depends on the initial con- Scaling in Holographic Thermalization, Phys. Rev. Lett. dition for the surface. For a surface that is initially flat, PðχÞ 112, 011601 (2014). is the Tracy-Widom distribution with β ¼ 1. This is the [9] H. Liu and S. J. Suh, Entanglement Growth during case relevant to our setup, where Sðx; t ¼ 0Þ¼ 0. (In the Thermalization in Holographic Systems, Phys. Rev. D directed polymer interpretation, this corresponds to a setup 89, 066012 (2014). where the x–coordinate of the upper end point of the [10] H. Casini, H. Liu, and M. Mezei, Spread of Entanglement polymer is fixed but that of the lower end point is free; and Causality, J. High Energy Phys. 07 (2016) 077. [11] M. Alishahiha, A. F. Astaneh, and M. R. Mohammadi again, this is the setup relevant to our minimal cut picture.) Mozaffar, Thermalization in Backgrounds with Hyper- Other initial conditions for a growing surface can give scaling Violating Factor, Phys. Rev. D 90, 046004 (2014). different universal forms for PðχÞ—for example, the so- [12] S. Kundu and J. F. Pedraza, Spread of Entanglement for called “narrow wedge” initial condition gives the Tracy- Small Subsystems in Holographic CFTs, Phys. Rev. D 95, Widom distribution with β ¼ 2. (The latter distribution is 086008 (2017). likely to be relevant to noisy growth of entanglement [13] M. Fagotti and P. Calabrese, Evolution of Entanglement between two subsystems that are initially unentangled with Entropy Following a Quantum Quench: Analytic Results each other, but separately highly entangled.) [The gener- for the XY Chain in a Transverse Magnetic Field, Phys. alization of the directed polymer picture to entangled initial Rev. A 78, 010306 (2008). 031016-27 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) [14] I. Peschel and V. Eisler, Reduced Density Matrices and in a Quantum Many-Body System, Nature (London) 528, Entanglement Entropy in Free Lattice Models, J. Phys. A 77 (2015). [33] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. 42, 504003 (2009). Schittko, P. M. Preiss, and M. Greiner, Quantum Thermal- [15] J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. ization through Entanglement in an Isolated Many-Body Daley, Entanglement Growth in Quench Dynamics with System, Science 353, 794 (2016). Variable Range Interactions, Phys. Rev. X 3, 031015 [34] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, (2013). Measuring Entanglement Growth in Quench Dynamics of [16] P. Hauke and L. Tagliacozzo, Spread of Correlations in Bosons in an Optical Lattice, Phys. Rev. Lett. 109, 020505 Long-Range Interacting Quantum Systems, Phys. Rev. (2012). Lett. 111, 207202 (2013). [35] P. Hayden and J. Preskill, Black Holes as Mirrors: [17] A. Coser, E. Tonni, and P. Calabrese, Entanglement Quantum Information in Random Subsystems, J. High Negativity after a Global Quantum Quench, J. Stat. Mech. Energy Phys. 09 (2007) 120. (2014) P12017. [36] Y. Sekino and L. Susskind, Fast Scramblers, J. High [18] A. S. Buyskikh, M. Fagotti, J. Schachenmayer, F. Essler, Energy Phys. 10 (2008) 065. and A. J. Daley, Entanglement Growth and Correlation [37] C. Dankert, R. Cleve, J. Emerson, and E. Livine, Exact and Spreading with Variable-Range Interactions in Spin and Approximate Unitary 2-Designs and Their Application to Fermionic Tunneling Models, Phys. Rev. A 93, 053620 Fidelity Estimation, Phys. Rev. A 80, 012304 (2009). (2016). [38] D. A. Roberts, D. Stanford, and L. Susskind, Localized [19] V. Alba and P. Calabrese, Entanglement and Thermody- Shocks, J. High Energy Phys. 03 (2015) 051. namics after a Quantum Quench in Integrable Systems, [39] R. Oliveira, O. C. O. Dahlsten, and M. B. Plenio, Generic arXiv:1608.00614. Entanglement Can be Generated Efficiently, Phys. Rev. [20] A. M. Lauchli and C. Kollath, Spreading of Correlations Lett. 98, 130502 (2007). and Entanglement after a Quench in the One-Dimensional [40] O. C. O. Dahlsten, R. Oliveira, and M. B. Plenio, The Bose Hubbard Model, J. Stat. Mech. (2008) P05018. Emergence of Typical Entanglement in Two-Party Ran- [21] H. Kim and D. A. Huse, Ballistic Spreading of Entangle- dom Processes, J. Phys. A 40, 8081 (2007). ment in a Diffusive Nonintegrable System, Phys. Rev. Lett. [41] M. Žnidarič, Exact Convergence Times for Generation of 111, 127205 (2013). Random Bipartite Entanglement, Phys. Rev. A 78, 032324 [22] L. Zhang, H. Kim, and D. A. Huse, Thermalization of (2008). Entanglement, Phys. Rev. E 91, 062128 (2015). [42] F. G. S. L. Brandao, A. W. Harrow, and M. Horodecki, [23] W. W. Ho and D. A. Abanin, Entanglement Dynamics in Local Random Quantum Circuits Are Approximate Poly- Quantum Many-Body Systems, Phys. Rev. B 95, 094302 nomial-Designs, Commun. Math. Phys. 346, 397 (2016). (2017). [43] J. Gütschow, S. Uphoff, R. F. Werner, and Z. Zimborás, [24] M. Žnidarič,T. ž. Prosen, and P. Prelovšek, Many-Body Time Asymptotics and Entanglement Generation of Clif- Localization in the Heisenberg XXZ Magnet in a Random ford Quantum Cellular Automata, J. Math. Phys. 51, Field, Phys. Rev. B 77, 064426 (2008). 015203 (2010). [25] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded [44] S. Leichenauer and M. Moosa, Entanglement Tsunami in Growth of Entanglement in Models of Many-Body (1 þ 1)-Dimensions, Phys. Rev. D 92, 126004 (2015). Localization, Phys. Rev. Lett. 109, 017202 (2012). [45] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic Scaling of [26] M. Serbyn, Z. Papić, and D. A. Abanin, Universal Slow Growing Interfaces, Phys. Rev. Lett. 56, 889 (1986). Growth of Entanglement in Interacting Strongly Disor- [46] D. A. Huse, C. L. Henley, and D. S. Fisher, Respond, Phys. dered Systems, Phys. Rev. Lett. 110, 260601 (2013). Rev. Lett. 55, 2924 (1985). [27] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenom- [47] D. A. Huse and C. L. Henley, Pinning and Roughening of enology of Fully Many-Body-Localized Systems, Phys. Domain Walls in Ising Systems due to Random Impurities, Rev. B 90, 174202 (2014). Phys. Rev. Lett. 54, 2708 (1985). [28] R. Vosk and E. Altman, Many-Body Localization in One [48] D. Forster, D. R. Nelson, and M. J. Stephen, Large- Dimension as a Dynamical Renormalization Group Fixed Distance and Long-Time Properties of a Randomly Stirred Point, Phys. Rev. Lett. 110, 067204 (2013). Fluid, Phys. Rev. A 16, 732 (1977). [29] A. Chandran and C. R. Laumann, Semiclassical Limit for [49] T. Halpin-Healy, Directed Polymers versus Directed the Many-Body Localization Transition, Phys. Rev. B 92, Percolation, Phys. Rev. E 58, R4096 (1998). 024301 (2015). [50] S. Ryu and T. Takayanagi, Holographic Derivation of [30] D. J. Luitz, N. Laflorencie, and F. Alet, Extended Slow Entanglement Entropy from the Anti–de Sitter Space/ Dynamical Regime Close to the Many-Body Localization Conformal Field Theory Correspondence, Phys. Rev. Lett. Transition, Phys. Rev. B 93, 060201 (2016). 96, 181602 (2006). [31] F. Verstraete, V. Murg, and J. I. Cirac, Matrix Product [51] B. Swingle, Entanglement Renormalization and Hologra- States, Projected Entangled Pair States, and Variational phy, Phys. Rev. D 86, 065007 (2012). Renormalization Group Methods for Quantum Spin Sys- [52] F. Pastawski, B. Yoshida, D. Harlow, and J. Preskill, tems, Adv. Phys. 57, 143 (2008). Holographic Quantum Error-Correcting Codes: Toy [32] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Models for the Bulk/Boundary Correspondence, J. High Rispoli, and M. Greiner, Measuring Entanglement Entropy Energy Phys. 06 (2015) 149. 031016-28 QUANTUM ENTANGLEMENT GROWTH UNDER RANDOM … PHYS. REV. X 7, 031016 (2017) [53] P. Hayden, S. Nezami, X.-L. Qi, N. Thomas, M. Walter, [73] A. Nahum, S. Vijay, and J. Haah “Operator Spreading in and Z. Yang, Holographic Duality from Random Tensor Random Unitary Circuits” (to be published). Networks, J. High Energy Phys. 11 (2016) 009. [74] G. Evenbly and G. Vidal, Tensor Network States and [54] K. A. Takeuchi and M. Sano, Evidence for Geometry- Geometry, J. Stat. Phys. 145, 891 (2011). Dependent Universal Fluctuations of the Kardar-Parisi- [75] M. Kardar, Roughening by Impurities at Finite Temper- Zhang Interfaces in Liquid-Crystal Turbulence, J. Stat. atures, Phys. Rev. Lett. 55, 2923 (1985). [76] M Kardar, Statistical Physics of Fields (Cambridge Phys. 147, 853 (2012). University Press, Cambridge, England, 2007). [55] K. A. Takeuchi, M. Sano, T. Sasamoto, and H. Spohn, [77] A. Nahum, J. Ruhman, and D. Huse, “Dynamics of Growing Interfaces Uncover Universal Fluctuations Behind Scale Invariance, Sci. Rep. 1, 34 (2011). entanglement and transport in 1D systems with quenched [56] P. Calabrese, P. Le Doussal, and A. Rosso, Free-Energy randomness” (to be published). Distribution of the Directed Polymer at High Temperature, [78] E. H. Lieb and D. W. Robinson, in Statistical Mechanics: Selecta of Elliott H. Lieb, edited by B. Nachtergaele, Europhys. Lett. 90, 20002 (2010). J. P. Solovej, and J. Yngvason (Springer, Berlin, 2004), [57] V. Dotsenko, Bethe Ansatz Derivation of the Tracy-Widom Distribution for One-Dimensional Directed Polymers, pp. 425–431. Europhys. Lett. 90, 20003 (2010). [79] D. A. Roberts and B. Swingle, Lieb-Robinson Bound and [58] T. Sasamoto and H. Spohn, One-Dimensional Kardar- the Butterfly Effect in Quantum Field Theories, Phys. Rev. Parisi-Zhang Equation: An Exact Solution and Its Uni- Lett. 117, 091602 (2016). versality, Phys. Rev. Lett. 104, 230602 (2010). [80] D. Gottesman, The Heisenberg Representation of [59] T. Sasamoto and H. Spohn, Exact Height Distributions for Quantum Computers, arXiv:quant-ph/9807006. [81] We truncate all stabilizers to region A by throwing away all the KPZ Equation with Narrow Wedge Initial Condition, Nucl. Phys. B834, 523 (2010). the spin operators acting outsideA. In this process some of the [60] T. Sasamoto and H. Spohn, The Crossover Regime for the stabilizers become trivial, and some become redundant, i.e., Weakly Asymmetric Simple Exclusion Process, J. Stat. equal to products of the others. I is the number of stabilizers Phys. 140, 209 (2010). that are independent when truncated to A. Equivalently, I is [61] G. Amir, I. Corwin, and J. Quastel, Probability Distribu- the rank of the matrix Ψ after the rows corresponding to the tion of the Free Energy of the Continuum Directed complement of A have been deleted; this is how we calculate Random Polymer in 1 þ 1 Dimensions, Commun. Pure the entropy numerically for Sec. VI A. Appl. Math. 64, 466 (2011). [82] A. Hamma, R. Ionicioiu, and P. Zanardi, Bipartite [62] P. Calabrese and P. Le Doussal, Exact Solution for the Entanglement and Entropic Boundary Law in Lattice Spin Kardar-Parisi-Zhang Equation with Flat Initial Condi- Systems, Phys. Rev. A 71, 022315 (2005). tions, Phys. Rev. Lett. 106, 250603 (2011). [83] A. Hamma, R. Ionicioiua, and P. Zanardi, Ground State [63] S. Prolhac and H. Spohn, Height Distribution of the Entanglement and Geometric Entropy in the Kitaev Model, Kardar-Parisi-Zhang Equation with Sharp-Wedge Initial Phys. Lett. A 337, 22 (2005). Condition: Numerical Evaluations, Phys. Rev. E 84, [84] T. Halpin-Healy and Y.-C. Zhang, Kinetic Roughening 011119 (2011). Phenomena, Stochastic Growth, Directed Polymers and [64] P. Le Doussal and P. Calabrese, The KPZ Equation with All That. Aspects of Multidisciplinary Statistical Mechan- Flat Initial Condition and the Directed Polymer with One ics, Phys. Rep. 254, 215 (1995). Free End, J. Stat. Mech. (2012) P06001. [85] O. C. O. Dahlsten and M. B. Plenio, Exact Entanglement [65] T. Imamura and T. Sasamoto, Exact Solution for the Probability Distribution of Bi-Partite Randomised Stabi- Stationary Kardar-Parisi-Zhang Equation, Phys. Rev. lizer States, Quantum Inf. Comput. 6, 527 (2006). Lett. 108, 190603 (2012). [86] ITensor, http://itensor.org. [66] T. Imamura and T. Sasamoto, Stationary Correlations for [87] D. Gottesman, A Class of Quantum Error-Correcting the 1D KPZ Equation, J. Stat. Phys. 150, 908 (2013). Codes Saturating the Quantum Hamming Bound, Phys. [67] D. N. Page, Average Entropy of a Subsystem, Phys. Rev. Rev. A 54, 1862 (1996). [88] G. Roósz, R. Juhász, and F. Iglói, Nonequilibrium Dy- Lett. 71, 1291 (1993). [68] C. Nadal, S. N. Majumdar, and M. Vergassola, Statistical namics of the Ising Chain in a Fluctuating Transverse Distribution of Quantum Entanglement for a Random Field, Phys. Rev. B 93, 134305 (2016). Bipartite State, J. Stat. Phys. 142, 403 (2011). [89] L. Saul, M. Kardar, and N. Read, Directed Waves in [69] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Random Media, Phys. Rev. A 45, 8859 (1992). Matrix Product State Representations, Quantum Inf. [90] I. Peschel, Calculation of Reduced Density Matrices Comput. 7, 401 (2007). from Correlation Functions, J. Phys. A 36, L205 [70] C. Chamon and E. R. Mucciolo, Virtual Parallel Comput- (2003). ing and a Search Algorithm Using Matrix Product States, [91] T. Nattermann, Ising Domain Wall in a Random Pinning Phys. Rev. Lett. 109, 030503 (2012). Potential, J. Phys. C 18, 6661 (1985). [71] P. Meakin, P. Ramanlal, L. M. Sander, and R. C. Ball, [92] M. Kardar, Domain Walls Subject to Quenched Impurities, Ballistic Deposition on Surfaces, Phys. Rev. A 34, 5091 J. Appl. Phys. 61, 3601 (1987). (1986). [93] D. S. Fisher, Interface Fluctuations in Disordered Systems: [72] J. M. Kim and J. M. Kosterlitz, Growth in a Restricted 5 − ϵ Expansion and Failure of Dimensional Reduction, Solid-on-Solid Model, Phys. Rev. Lett. 62, 2289 (1989). Phys. Rev. Lett. 56, 1964 (1986). 031016-29 NAHUM, RUHMAN, VIJAY, and HAAH PHYS. REV. X 7, 031016 (2017) [94] A. A. Middleton, Numerical Results for the Ground-State (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Interface in a Random Medium, Phys. Rev. E 52, R3337 Dagstuhl, (1995). Germany, 2013), Vol. 22, pp. 270–284. [95] A. R. Calderbank, E. M Rains, P. W. Shor, and N. J. A. [98] T. Halpin-Healy, Disorder-Induced Roughening of Diverse Sloane, Quantum Error Correction and Orthogonal Manifolds, Phys. Rev. A 42, 711 (1990). Geometry, Phys. Rev. Lett. 78, 405 (1997). [99] T. Kriecherbauer and J. Krug, A Pedestrian’s View on [96] A. Klappenecker and M. Rotteler, Beyond Stabilizer Interacting Particle Systems, KPZ Universality and Ran- Codes II: Clifford Codes, IEEE Trans. Inf. Theory 48, dom Matrices, J. Phys. A 43, 403001 (2010). 2396 (2002). [100] I. Corwin, The Kardar-Parisi-Zhang Equation and [97] N. Linden, F. Matus, M. B. Ruskai, and A. Winter, in Universality Class, Random Matrices Theory Appl. 01, Proceedings of the 8th Conference on the Theory of 1130001 (2012). Quantum Computation, Communication and Cryptogra- [101] T. Halpin-Healy and K. A. Takeuchi, A KPZ Cocktail- phy (TQC 2013), edited by S. Severini and F. Brandao Shaken, Not Stirred…, J. Stat. Phys. 160, 794 (2015). 031016-30

Journal

Physical Review XAmerican Physical Society (APS)

Published: Jul 1, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial