Noisy Coupled Qubits: Operator Spreading and the Fredrickson-Andersen Model
NNoisy coupled qubits: operator spreading and the Fredrickson–Andersen model
Daniel A. Rowlands and Austen Lamacraft
TCM Group, Cavendish Laboratory, University of Cambridge,J. J. Thomson Ave., Cambridge CB3 0HE, UK ∗ (Dated: March 16, 2020)We study noise-averaged observables for a system of exchange-coupled quantum spins (qubits),each subject to a stochastic drive, by establishing mappings onto stochastic models in the strong-noise limit. Averaging over noise yields Lindbladian equations of motion; when these are subjected toa strong-noise perturbative treatment, classical master equations are found to emerge. The dynamicsof noise averages of operators displays diffusive behaviour or exponential relaxation, depending onwhether the drive conserves one of the spin components or not. In the latter case, the secondmoment of operators – from which the average subsystem purity and out-of-time-order correlationfunctions can be extracted – is described by the Fredrickson–Andersen model, originally introducedas a model of cooperative relaxation near the glass transition. It is known that fluctuations ofa ballistically propagating front in the model are asymptotically Gaussian in one dimension. Weextend this by conjecturing, with strong numerical evidence, that in two dimensions the long-timefluctuations are in the Kardar–Parisi–Zhang universality class, complementing a similar observationin random unitary circuits. I. INTRODUCTION
The success of microscopic models of matter hingesupon the assumption that the resulting macroscopic de-scription is relatively insensitive to the precise dispositionof the constituent particles. When applied to dynamicalphenomena – collective motion – this assumption seemsat odds with our usual understanding of generic dynam-ical systems: that they display chaos and an exponentialsensitivity to initial conditions.That we can derive the (deterministic) laws of hydro-dynamics from the motion of gas particles and the as-sumption of molecular chaos shows that this contradic-tion is not as severe as it may at first seem. By focussingon coarse-grained variables like the average local velocity,the underlying chaotic motion fades into the backgroundand serves only to give rise to the pressure, viscosity andother parameters of the effective description.Nevertheless, if one believes that the butterfly effect ismore than a figure of speech, something must have beenlost along the way. By focusing on average quantities,the growth of fluctuations from the microscopic to themacroscopic is obscured. Long a part of statistical fluiddynamics [1–3], these questions have only recently beentaken up in quantum field theory [4–11] and many-bodyphysics [12–18].Traditionally, these fields have been concerned with av-erages (cid:104)O j ( t ) (cid:105) , and response functions i (cid:104) [ O j ( t ) , O k (0)] (cid:105) ofHeisenberg picture observables O j ( t ). However, the actof taking expectations in these quantities obscures thepossibility that in a given experiment we may observe avery different response in observable O j ( t ) to a pertur-bation coupled to observable O j (0). The variance of theresponse function defines the out-of-time-order correla- ∗ [email protected] tion function (OTOC) C jk ( t ) ≡ (cid:104) [ O j ( t ) , O k (0)] † [ O j ( t ) , O k (0)] (cid:105) , (1)that has been suggested as a diagnostic of many-bodyquantum chaos [9–11, 19]. In the light of the above dis-cussion, it is convenient to think of the OTOC in termsof a supersystem consisting of two independent copies ofthe system under consideration, and extract it from theoperator O j ( t ) ⊗ O j ( t ). The duplicate system is some-times known as the thermofield double [20].In recent years, OTOCs have been calculated ina variety of models, including the Sachdev–Ye–Kitaev(SYK) model [21–23], the many-body localized phase ofone-dimensional spin models [24–28], weakly interactingfermions [29, 30], as well as chaotic single-particle sys-tems [31].In this work we will be concerned with the dynamicsof a system of coupled qubits (spin-1/2 objects), sub-ject to classical noise described by a stochastic process η t [32]. We will be concerned with the first two operatormoments: O j ≡ E η [ O j ] , O j ⊗ O j ≡ E η [ O j ⊗ O j ] . (2)where O j ∈ { X j , Y j , Z j } is one of the Pauli matricesdescribing qubit j . The motivations for this study are:1. Conventional wisdom suggests that noise is anti-thetical to quantum coherence. On the other hand,the evolution of a quantum system in the presenceof classical noise is still unitary. We will see thatthe expected loss of coherence is only true on aver-age: the first and second moments have completelydifferent behaviour.2. The limit of strong noise provides a controlled ap-proximation in which we can obtain a tractable dy-namics of the moments. a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r FIG. 1. Growth of fronts in the one dimensional FA modelwith a finite average ‘butterfly’ velocity v B .
3. In the era of real noisy intermediate-scale quantumcomputers [33], there is need to understand the dy-namics of quantum information in the presence ofstrong noise.The stochastic models we introduce may be regardedas continuous-time analogues of the random unitary cir-cuit model studied in many recent works [12–15, 34–41].Though they share some phenomenology, our analysisof the stochastic models is completely different, beingbased on Lindblad equations. Expectations over stochas-tic trajectories are taken at the first step, and the analysisof the strong-noise limit is based on conventional manybody perturbation theory. This allows any (determinis-tic) coupling between qubits to be taken into account –though we focus on a Heisenberg coupling for simplicityof presentation – and allows models with conservation ofone of the spin components to be handled on the samefooting.As is the case for the random unitary circuit mod-els, the dynamics of moments in certain cases can beidentified with the probability distribution of a classicalstochastic process [42–45] (see Table I). For the secondmoment (or OTOC) in a model without any conservedquantities, this is the Fredrickson–Andersen (FA) model,originally introduced to describe dynamics at the glasstransition [46]. The FA model is an example of a kineti-cally constrained model – see [47] for a recent review – andhas a rich phenomenology that we apply to the stochasticspin model [48]. Specifically, we will see that with appro-priate initial conditions the FA model describes ballisti-cally growing fronts that are associated with the spread-ing of the support of local operators in the Heisenbergpicture (see Fig.1). A characteristic speed for operatorspreading in many-body systems was first identified in[7], and it has since become known as the ‘butterfly ve-locity’ v B . A. Outline
The outline of this work is as follows. Section II de-scribes the models we use, and presents our derivationsof the mappings to stochastic models in the strong-noiselimit. In Section III we study the behaviour of frontsin the FA model in one and two dimensions. In Sec-tion IV we then apply the resulting phenomenology toOTOCs and the behaviour of subsystem purity. Ourconclusions are discussed in Section V. Technical detailsrelated to numerical simulations and calculation of theeffective Lindbladians may be found in the appendices.
II. MODELS AND MAPPINGA. Models
We consider a system of N spin-1/2 objects (qubits),and the closely related cases of (1) Heisenberg pictureevolution of observables O ( t ): ∂ t O ( t ) = i [ H ( η t ) , O ( t )]under a Hamiltonian H ( η t ) depending on a stochasticprocess η t , and (2) the von Neumann equation for thedensity matrix ρ ( t ): ∂ t ρ ( t ) = − i [ H ( η t ) , ρ ( t )].We will be concerned with the first two operator mo-ments: O ≡ E η [ O ] , O ⊗ O ≡ E η [ O ⊗ O ] . (3)The second (and higher) moments of an operator carry agreat deal more information about the dynamics than theaverage alone. In particular, the second moment gives us:1. The average purity γ ≡ tr (cid:2) ρ A (cid:3) , (4)where ρ A is the reduced density matrix of a sub-system A . Labelling bases of the subsystem A andits complement A c by the (multi-)indices A and A c ,we have γ = (cid:88) A , A , A c , A c ρ A A c , A A c ρ A A c , A A c , (5)which may be extracted from ρ ⊗ ρ .2. The average OTOC C ( j − k, t ) ≡
12 tr (cid:104) ρ [ O j ( t ) , O k (0)] † [ O j ( t ) , O k (0)] (cid:105) , (6)where O i,j are operators on qubit j and k . C ( x, t )may be extracted from O j ( t ) ⊗ O j ( t ) by contract-ing indices with ρ and O k (0).The generalization of our approach to higher momentsis straightforward, and we will return to this point in theconclusion.We will consider two models, one with conserved total z -component of spin (C) and one without (NC). In bothcases, the Hamiltonian is the sum of deterministic andstochastic terms H C = H + N (cid:88) j =1 η jt Z j (C) H NC = H + N (cid:88) j =1 η jt · σ j , (NC)where σ = ( X, Y, Z ) are the usual Pauli matrices and η t is assumed to be delta-correlated noise, E η (cid:104) η jt η kt (cid:48) (cid:105) = gδ ( t − t (cid:48) ) δ j,k , (7)while η t = ( η xt , η yt , η zt ) contains three independent com-ponents (assumed identical for ease of notation). H is atime-independent Hamiltonian describing local couplingbetween the spins. Our approach is not particularly sen-sitive to the precise form of H , but for simplicity we willbegin with the Heisenberg chain H = J N (cid:88) j =1 σ j · σ j +1 . (8)The generalization to other models, including those inhigher dimension, will be evident.If we formally express the white noise η jt as the Itˆodifferential of a Brownian motion B jt , (cid:90) t (cid:48) η jt dt = (cid:90) t (cid:48) dB jt , (9)then one finds that the generator dH C t of infinitesimalstochastic unitary evolution in model C is given by [32] dH C t = Hdt + √ g (cid:88) j Z j dB jt . (10)A straightforward exercise in Itˆo calculus yields the equa-tion of motion of the density matrix [49] dρ t = e − idH C t ρ t e idH C t − ρ t = − i [ H, ρ ] dt − g (cid:88) j [ Z j , [ Z j , ρ ]] dt − i √ g (cid:88) j [ Z j , ρ Ψ ] dB jt , (11)the expectation of which is thus of Lindblad form [50] ∂ t ¯ ρ = L H (¯ ρ ) (cid:122) (cid:125)(cid:124) (cid:123) − i [ H, ¯ ρ ] − D (¯ ρ ) (cid:122) (cid:125)(cid:124) (cid:123) g (cid:88) j [ Z j , [ Z j , ¯ ρ ]] . (12)The corresponding equation for model NC, obtained byan analogous procedure, is ∂ t ¯ ρ = − i [ H, ¯ ρ ] − g (cid:88) j (cid:88) a =1 [ σ aj , [ σ aj , ¯ ρ ]] . (13) The derivation for the first moment of an operator is ob-tained similarly, the only differences being ρ → O and dH t → − dH t . We also note that this calculation can al-ternatively be done, albeit less concisely, by interpretingthe stochastic process η t in the Stratonovich sense (seee.g. Appendix A.4 of [51]). B. Model C at strong noise: Symmetric ExclusionProcess for O A key simplification occurs in the limit of strong noise( g large), where the dynamics of the moments is re-stricted certain slow subspaces . In the case of ¯ O , this isthe kernel of the dissipator D ( ¯ O ). For example, in ModelC the dynamics of O is restricted to the 2 N -dimensionaldiagonal matrix elements in the basis of Z j eigenstates | z : z N (cid:105) [52] with z j = ±
1. That is, only the matrixelements (cid:104) z : z N |O| z : z N (cid:105) (14)survive the averaging as they are unaffected by the de-phasing noise. We now show that the evolution of theprobability distribution of a spin configuration z : z N isdescribed by the symmetric exclusion process (SEP) [53](see Table I). Model NC Model C O Exponential decay Symmetric Exclusion Process
O ⊗ O
Fredrickson–Andersen ‘Octahedral’ modelTABLE I. The behaviour of the first and second operatormoments in Models NC and C in the strong-noise limit.
The dynamics on the slow subspace S can be analyzedperturbatively in g − (as is done in [51] and [54]); theeffective Liouvillian to leading order takes the form L eff = −P S L H D − L H P S , where P S is the projector onto S and D − is the inverse of the restriction of D to its coimage.Explicit evaluation leads to L eff ¯ ρ = − J g (cid:88) j [ σ + j σ − j +1 + h.c. , [ σ + j σ − j +1 + h.c. , P S ¯ ρ ]] . (15)Let us regard ¯ ρ ∈ S as an element (cid:126)ρ of a vector spaceover R (sometimes referred to as ‘superspace’) with basis {| (cid:105) (cid:104) | , | (cid:105) (cid:104) |} ⊗ N . Then L eff acts as a matrix L on (cid:126)ρ ,giving the master equation ∂ t (cid:126)ρ = L (cid:126)ρ, (16)with L = J g (cid:80) i ( (cid:126)σ i · (cid:126)σ i +1 − − L (if we think of the master equa-tion as an imaginary-time Schr¨odinger equation) coin-cides with that of a Heisenberg ferromagnet. Up to a con-stant, L is thus seen to be the generator of the SEP [55].In the one-dimensional case, we have an alternative routeto this result as the model is found to be integrable bymeans of a mapping to an imaginary-U Hubbard model[56]. In the strong-noise limit, the Bethe ansatz equationsreduce to those of the spin-1/2 ferromagnetic Heisenbergmodel, from which the quadratically dispersing Liouvil-lian spectrum and consequent diffusive relaxation follow(as was established earlier in [57] by analytic evaluationof the single-particle Green’s function). C. Model NC at strong noise:Fredrickson–Andersen model for
O ⊗ O
The dynamics of the first moment in model NC is triv-ial: the slow subspace is one-dimensional (i.e. containsonly the identity) and so we find fast local relaxationrather than any hydrodynamics as in model C.The second moment of Model NC does have interest-ing dynamics in the presence of strong noise, however.Because the noise in this model randomizes all compo-nents of the spins,
O ⊗ O lives in the tensor product ofthe space spanned by the rotationally invariant single-sitefactors | j (cid:105) ≡ j ⊗ j | j (cid:105) ≡
16 [ X j ⊗ X j + Y j ⊗ Y j + Z j ⊗ Z j ] . (17)Any O ⊗ O of this form has the expansion
O ⊗ O = (cid:88) n : n N ∈{ , } N C O n : n N | n : n N (cid:105) . (18)Using the properties of the Pauli matrices it is easy toshow tr (cid:2) O (cid:3) = (cid:88) n : n N ∈{ , } N C O n : n N . (19)Since the trace of any operator product tr [ O ( t ) O ( t )] isconserved under Heisenberg evolution, we may think of C O n : n N as a probability distribution (up to overall nor-malization) and its evolution equation ∂ t C O = LC O (20)as a (classical) master equation. In Refs. [42–45], a re-lated discrete-time Markov chain was obtained for thedynamics of operator moments due to randomly chosentwo-qubit unitary transformations. This Markov chainwas the basis of the calculations of OTOCs and purity inthe random unitary circuit model in Ref. [14].What stochastic process is described by L ? We willsee that it is the Fredrickson–Andersen (FA) model [46].The FA model is defined on a lattice with sites that mayeither be in state 1 or 0, with pairs of neighbouring sites j and k undergoing the transitions1 j k Γ − (cid:42)(cid:41) − Γ j k (21) with rates Γ , . In the stationary state, sites are inde-pendent with probability p = Γ / (Γ + Γ ) to be 1. Wefind Γ = Γ / J / g for model NC, i.e., 1s are threetimes more common than 0s: C O n : n N | stationary = 14 N (cid:89) j n j . (22)Two further comments: (1) The dynamics of C A and C ρ are identical because the rates are quadratic in J . (2)Individual trajectories of the FA model have no meaning,as only the probability distribution C O appears in themoment O ⊗ O .The derivation of the effective dynamics for
O × O fol-lows the same pattern as for the first moment. Noise av-eraging the stochastic differential equation for the secondmoment of ρ in model NC leads to the Lindblad equationfor the replicated system ∂ t ρ ⊗ ρ = − i [ H , ρ ⊗ ρ ] − g (cid:88) j,a [Σ aj , [Σ aj , ρ ⊗ ρ ]] (23)where we have introduced the operators H = H ⊗ ⊗ H and Σ aj = σ aj ⊗ ⊗ σ aj .The kernel of the dissipator determines the slow sub-space S = span (cid:0) {| (cid:105) , | (cid:105)} ⊗ N (cid:1) , where the {| (cid:105) , | (cid:105)} stateswere defined in (17). The effective Liouvillian to leadingorder (see Appendix B for the derivation for model C; themodel NC result is obtained similarly) acts on elementsof S according to L eff ( · · · ) = − g P S [ H , [ H , · · · ] . (24)A matrix representation for L eff can again be found (seeAppendix C for details) if we take {| (cid:105) , | (cid:105)} ⊗ N as a vectorspace basis. This matrix is given by L = (cid:80) j L j,j +1 , where L j,j +1 = 4 J g − /
30 0 − /
30 1 1 − / , (25)can be identified as the transition rate matrix of acontinuous-time Markov process: the one-spin facilitatedFA model with rates given in Eq. (21).A similar analysis for the second moment of modelC does not appear to lead to a mapping to a classicalstochastic model, but we nevertheless make some obser-vations about the effective dynamics in Appendix B. III. PHENOMENOLOGY OF FRONTSA. Fronts in the FA model
The FA model was originally introduced to describedynamics at the glass transition, and is an example of a kinetically constrained model , see [47] for a recent review.The model has a spectral gap [58], indicating that equilib-rium fluctuations are generically exponentially decayingin time. Our main interest, however, is in the nonequi-librium dynamics of the model, in particular in initialconditions with only a few 1s, or regions devoid of 1s. Inthis case a nonzero density of 1s grows into the emptyregion with a finite front velocity, see Fig. 1.The dynamics of a front in the FA model in one di-mension was recently analyzed rigorously in Ref. [59] fora variant of the model [60] in which the transition rate isindependent of the number of neighbours (see Ref.[61] forthe related case of the East model). There it was shownthat if the rightmost 1 starts at site 0, its displacement X t after time t is asymptotically given by the normaldistribution X t − v B t √ t d −→ t →∞ N (0 , s ) , (26)for some v B and s . This chimes with the arguments givenin Refs. [13–15] for the random unitary circuit model thatthe probability distribution of X t is that of a biased ran-dom walk. B. Fronts in two dimensions
The derivation of the FA model in the strong-noiselimit holds in any dimension. Ballistic motion of the frontin kinetically constrained models in higher dimensions isdiscussed in Refs. [62–64]. It is natural to ask how thefront distribution in (26) generalizes to higher dimen-sions. For the random unitary circuit model, Ref. [14]proposed – and provided numerical evidence – that thefluctuations of the front at long times are in the univer-sality class of the Kardar–Parisi–Zhang (KPZ) equation[65, 66]. In the 1 + 1-dimensional case, relevant for thegrowth of a front in two dimensions, this equation hasthe form ∂ t h = c + ν∂ x h + λ ∂ x h ) + ζ ( x, t ) . (27)Here h ( x, t ) denotes the displacement of the front in thedirection of growth, as a function of transverse coordinate x . The first term in (27) is a contribution to the ballisticgrowth rate; the second describes diffusive motion of thesurface; the third captures a quadratic dependence of thelocal growth rate on the tilt of the surface; the last is aspatially uncorrelated white noise. The quadratic termis a relevant perturbation below two spatial dimensionsthat is responsible for novel scaling behaviour. For theone-dimensional case considered here, fluctuations of thesurface have a dynamical critical exponent z – describingthe relative scaling of spatial and temporal fluctuations as t ∼ x z – of z = 3 /
2, and growth exponent β – describingthe growth of interface fluctuations as h ∼ t β – of β =1 / FIG. 2. (Top) Upward growth of a front in the 2D FA modelinto a region of 0s (black), starting from a row of 1s. (Bottom)Growth of front variance with time. Dashed line is the powerlaw 0 . t / , consistent with the KPZ growth exponent β =1 / growth from a row of 1s, corresponding to flat initialconditions, rather than from a single 1, which leads toa rounded cluster. This option was not available to theauthors of Ref. [14], as the peculiarities of the circuitmodel mean there is no roughening in a lattice direction.In a simulation of 10 time steps, we observe nearly twodecades of scaling with the KPZ exponent β ∼ / h ( x, t ) d −→ t →∞ ct + αt / χ, (28)where χ is a random variable with known distribution,and c and α are constants. In the case of flat initialconditions, χ is drawn from the Tracy-Widom distribu-tion corresponding to the Gaussian Orthogonal Ensemble(GOE) [68–70].In Fig. 3 we show a comparison between the proba-bility distribution of the front obtained at the end ofour numerical simulation, and the best fit Tracy-WidomGOE and Gaussian distributions. The superiority of theTracy-Widom fit is evident, in particular in capturingthe skew of the distribution and the differing behaviour FIG. 3. Front probability distribution at t = 10 comparedwith best fit Tracy-Widom density F (cid:48) ( s ) (dashed) and Gaus-sian (dash-dotted), where the width and shift of the distribu-tions are fitting parameters. of the left and right tailslog p ( h ) ∝ (cid:40) −| h | h → −∞− h / h → + ∞ . (29) IV. APPLIED PHENOMENOLOGYA. OTOCs
The mapping to the FA model yields a simple expres-sion for the OTOC (6) in Model NC with infinite tem-perature initial density matrix ρ C jk ( t ) = 2 − N (cid:18) (cid:19) (cid:88) n nN ∈{ , } Nnk =1 C O j n : n N ( t ) , (30)where the normalization follows from tr (cid:2) O j ( t ) (cid:3) = 2 N .For a FA model starting with only site j in state 1, theOTOC is then (4 / k has value 1at time t .Let us further make the reasonable assumption that,after the front arrives at site k , the probability to be instate 1 quickly approaches the equilibrium value of 3 / C jk ( t ) = Pr [ k in cluster seeded by j ] , (31)or in other words, the cumulative distribution functionof the front. The resulting behaviour of the OTOC isillustrated in Fig. 4. B. Purity decay
Consider a partition of the qubits into sets A and A c ,of sizes | A | and | A c | . The average purity of a region A is FIG. 4. The averaged OTOC C jk ( t ) is identified with theprobability for sites j and k to be in the same cluster. Here,we have illustrated the Gaussian case – the functional formbeing that of an error function, the cumulative distributionfunction of a Gaussian – appropriate to one dimension. expressed as (c.f. Ref. [15, 34]) γ = tr (cid:2) ρ A (cid:3) = 2 | A c | (cid:88) nj =0 , j ∈ Anj =0 for j ∈ Ac C ρn : n N . (32)The purity is (2 | A c | times) the probability that A c containonly 0s. Consider taking as an initial condition a randompure product state, described by a density matrix ρ = 12 N (cid:89) j (cid:2) j + (cid:37) j · σ j (cid:3) (33)with unit vectors (cid:37) j . Projected into the slow subspace ofModel NC, this gives C n : n N ( t = 0) = 2 − N . Comparingwith (32), we see γ ( t = 0) = 1, as required. Note that theoverall probability of A c being empty is 1 / | A c | , but thisexponentially small factor is cancelled by the prefactorin (32).If A c is empty (i.e. contains only 0s) at time t , weexpect the fronts to be within A at earlier times. To seehow this picture leads to the decay of purity, consider thegrowth of a single front in one dimension. The positionof the front is described by (26). To find the most likelytrajectory, the probability of finding a front at a distance X inside A at time 0 must be combined with the proba-bility of the front propagating to the boundary between A and A c at time t (see Fig.5). The joint probability ofthese two events is then P ( X t = 0 , X = − X ) ∼ | A c | + X exp (cid:32) − [ X − v B t ] s t (cid:33) . (34)For large t , it suffices to find the optimum value X ∗ ofthe initial front position, giving a ‘purity front’ velocity v PF < v B v PF ≡ X ∗ t = v B − s ln 2 . (35) FIG. 5. A ‘purity front’ velocity v PF < v B arises from thejoint probability of the propagating front deviating from thevelocity v B and an initial condition with X extra 0s. The fronts move slower than the butterfly velocity v B .Note that a similar argument appears in Ref. [34], thoughwith ad hoc assumptions about the statistics of front mo-tion.Substituting the optimal value X ∗ in (34) gives theexponential decay of the purity γ ( t ) ∼ exp (cid:18) − v B t ln 2 + s t (cid:19) , (36)which enables us to define the ‘purity velocity’ v P = v B − s ln 2, such that γ ( t ) ∼ γ (0) e − v P t ln 2 .Applied to a region A of finite size, a simple generaliza-tion of the above argument implies that two fronts movetowards each other at ± v PF . However, the two fronts never touch (see Fig. 6), for when t > | A | / v P , the mostlikely initial configuration is completely empty, and thepurity saturates. Thus we have γ ( t ) ∼ (cid:40) e − v P t ln 2 t < | A | / v P12 | A | t > | A | / v P (37)These results are valid for large t and | A | , where the op-timum dominates the probability. The fact that puritydecays on a time scale linear in the subsystem size isconsistent with our local Hamiltonian, and is to be con-trasted with the fast scrambling (i.e. in a time logarith-mic in the size) possible in systems with highly nonlocalcoupling [5, 6].If | A c | < | A | , the situation is slightly different. Once t > | A c | / v P , the most probable way for an empty A c to arise is from the stationary distribution (22), assum-ing this distribution is approached exponentially quicklyfrom the initial state, giving γ ( t > | A | / v P ) = 12 | A c | . (38)The purity dynamics we have found (see Fig. 6) areconsistent with the expectation, largely based on exact FIG. 6. (Top) Schematic representation of the propagationof two fronts in 1+1 spacetime. White represents a region of0s and black dots an active region of 1s and 0s. (Bottom) − ln ¯ γ ( t ) for the case when ¯ γ ( t ) is computed for a finite sub-system. The red lines at fixed spatial positions in the upperright figure demarcate a finite subsystem of size | A | . Beforesaturation, − ln ¯ γ ( t ) grows at a rate (shown) controlled by thepurity velocity. diagonalization studies [71] and toy models [12], that bal-listic entanglement growth is a universal feature of quan-tum chaotic many-body systems. By Jensen’s inequality,the growth rate of − ln ¯ γ ( t ) that we have calculated is alower bound on the growth of the averaged second R´enyientanglement entropy. Moreover, our continuous-time re-sults supplement the analytic discrete-time calculationsin random unitary circuits [14], which have also con-firmed this phenomenology (with a purity velocity sat-isfying v P < v B ), via a mapping of the average purityonto the partition function of a directed random walk.The approach, which has been extended to obtain thehigher R´enyi entropies from a correspondence with a hi-erarchy of classical statistical mechanics models [72], hasmotivated the suggestion of a “minimal membrane” pic-ture of entanglement spreading in generic nonintegrablequantum systems [37]. V. CONCLUSIONS
We have provided a precise account of operator spread-ing for a system of interacting qubits undergoing con-tinuous time evolution, with each qubit independentlycoupled to a stochastic drive. By averaging over noise,Lindblad equations for the first and second operator mo-ments were derived and studied perturbatively in thestrong-noise limit; the central result being the identifi-cation of a mapping to the Fredrickson–Andersen modelfor the second moment dynamics in the case of noise thatdoes not conserve a spin component. Considering thephenomenology of front growth in this model then en-abled us to determine the implications for the behaviourof OTOCs and the decay of subsystem purity, which werefound to be in line with the results established in randomunitary circuit models. Although the mapping holds inarbitrary dimension, we restricted our attention to theone- and two-dimensional case: in one spatial dimension,we exploited the known exact Gaussian asymptotics ofthe front, whilst in two dimensions we conjectured, withnumerical support, that front fluctuations exhibit (1+1)-dimensional KPZ universality, thus giving us access toexact results for the front shape in terms of Tracy-Widomdistributions.Our approach generalizes naturally to higher operatormoments [73]. We expect that the identification of theslow subspaces and projection of the dynamics into thosespaces will be more involved but tractable, and will allowthe study of higher entanglement entropies and the fulldistribution of operator statistics.
Acknowledgments . DAR and AL gratefully acknowl-edge the EPSRC for financial support, under Grant Nos.EP/M506485/1 and EP/P034616/1 respectively. Wethank Oriane Blondel, Juanpe Garrahan, Robert Jackand Katarzyna Macieszczak for useful discussions. Thispaper was finished while AL was a participant in the pro-gramme ‘Quantum Paths’ at the Erwin Schr¨odinger In-ternational Institute for Mathematics and Physics (ESI).
Appendix A: Details of numerical simulation
For convenience we study the version of the FA modelin which the probability of a site to flip its state dependsonly on having neighbours, not their number [60].We use multispin coding [74–76], whereby 64 configu-rations of the model are represented as an array of (un-signed) 64-bit integers, and updated by bitwise opera-tions. This allows 64 trajectories to be simulated simul- taneously on a single core. Our simulation consisted of asingle run on each of 16 virtual cores, corresponding to1024 trajectories.We initialized a L × H lattice with L
1s in the firstrow, enforcing periodic boundary conditions along therows. The front is defined as the height of the highest 1 ineach column. As the front grows in the vertical direction,it must be periodically reset so that it remains roughlycentered. This is achieved by calculating the mean heightof the front across the 64 configurations every 10 updatesand moving the configuration downward by 10 sites whenthe mean exceeds H/ L = 10 , H = 200,and T = 10 timesteps. A height of H = 100 resulted ina breakdown of scaling behaviour at the longest times,presumably a consequence of the fluctuations of the in-terface being bounded by the finite height window.The simulation code and data analysis are writtenin Julia and are available as a Jupyter notebook, to-gether with the simulation data, at https://github.com/AustenLamacraft/FA-front . Appendix B: Second moment dynamics of Model Cin the strong-noise limit
Application of Itˆo’s lemma enables us to write downthe stochastic differential equation for the second mo-ment of the density matrix in model C d ( ρ ⊗ ρ ) = − i [ H , ρ ⊗ ρ ] dt − i √ g (cid:88) j [Σ zj , ρ ⊗ ρ ] dB jt − g (cid:88) j [Σ zj , [Σ zj , ρ ⊗ ρ ]] dt, (B1)which upon averaging leaves us with ∂ t ρ ⊗ ρ = − i [ H , ρ ⊗ ρ ] − g (cid:88) j [Σ zj , [Σ zj , ρ ⊗ ρ ]] (B2)where we have adopted analogous notation to that of(23), i.e., Σ zj = Z j ⊗ ⊗ Z j . In the strong-noise limit,the dynamics is projected onto the 6 N slow subspace S spanned by {| z (cid:105) (cid:104) z | ⊗ | z (cid:105) (cid:104) z | , | z (cid:105) (cid:104)− z | ⊗ |− z (cid:105) (cid:104) z | : z i ∈ {− , }} ⊗ N . The first nonvanishing term in per-turbation theory for the generator of the strong-noisedynamics of the second moment of the density matrix,which we also refer to as an effective Liouvillian, is givenby L eff = −P S L H D − L H P S . If we consider the actionof DL H on a single site, we have DL H ( | z (cid:105) (cid:104) z (cid:48) | ⊗ | z (cid:105) (cid:104) z (cid:48) | ) = ig (cid:32) [ H, | z (cid:105) (cid:104) z (cid:48) | ] ⊗ | z (cid:105) (cid:104) z (cid:48) |× (cid:32) ( z j − z j (cid:48) ) + 4 (cid:32) − ( z j − z j (cid:48) ) (cid:33) − z j (cid:32) − ( z j − z j (cid:48) ) (cid:33) ( z j − z j (cid:48) ) (cid:33) + | z (cid:105) (cid:104) z (cid:48) | ⊗ [ H, | z (cid:105) (cid:104) z (cid:48) | ] × (cid:32) ( z j − z j (cid:48) ) + 4 (cid:32) − ( z j − z j (cid:48) ) (cid:33) − z j (cid:32) − ( z j − z j (cid:48) ) (cid:33) ( z j − z j (cid:48) ) (cid:33) (cid:33) , (B3)from which it follows that L H P S (cid:16) ρ ⊗ ρ (cid:17) is an eigenstateof D with eigenvalue 4. The effective Liouvillian can thusbe written L eff = − g P S ad H , (B4)where we have used the adjoint action notationad H ( · · · ) := [ H , · · · ].In the strong-noise limit, the dynamics is projectedonto the 6 N -dimensional slow subspace S spannedby {| z (cid:105) (cid:104) z | ⊗ | z (cid:105) (cid:104) z | , | z (cid:105) (cid:104)− z | ⊗ |− z (cid:105) (cid:104) z | : z i ∈{− , }} ⊗ N . It is helpful to partition the single-site fac-tors of S into three types of pairs of states1 . (11 ,
11) ( − − , − − . (11 , − −
1) ( − − , . (1 − , −
11) ( − , − , (B5)where we have represented the state | z (cid:105) (cid:104) z | ⊗ | z (cid:105) (cid:104) z | by the tuple ( z z , z z ). It is helpful to visualise eachpair of states as occupying antipodal vertices of an oc-tahedron. If we consider the effective Liouvillian, whichdiffers from that of model C only by a multiplicative con-stant and the fact that the projector is into a differentslow subspace, we identify three classes of matrix element(with the possible values given in parenthesis, and theirrepresentation on the octahedron given in Fig. 7):1. Pair changing ( ± ± FIG. 7. An example of exchange (left) and pair changing(right) terms, as visualised on a square cross section of the‘octahedron’ of single-site states that span the slow subspaceof model C. The filled and empty circles represent the stateson adjacent sites.
Appendix C: Effective Liouvillian for the secondmoment of model NC