Locality, Quantum Fluctuations, and Scrambling

Thermalization of chaotic quantum many-body systems under unitary time evolution is related to the growth in complexity of initially simple Heisenberg operators. Operator growth is a manifestation of information scrambling and can be diagnosed by out-of-time-order correlators (OTOCs). However, the behavior of OTOCs of local operators in generic chaotic local Hamiltonians remains poorly understood, with some semiclassical and large N models exhibiting exponential growth of OTOCs and a sharp chaos wavefront and other random circuit models showing a diffusively broadened wavefront. In this paper we propose a unified physical picture for scrambling in chaotic local Hamiltonians. We construct a random time-dependent Hamiltonian model featuring a large N limit where the OTOC obeys a Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) type equation and exhibits exponential growth and a sharp wavefront. We show that quantum fluctuations manifest as noise (distinct from the randomness of the couplings in the underlying Hamiltonian) in the FKPP equation and that the noise-averaged OTOC exhibits a cross-over to a diffusively broadened wavefront. At small N we demonstrate that operator growth dynamics, averaged over the random couplings, can be efficiently simulated for all time using matrix product state techniques. To show that time-dependent randomness is not essential to our conclusions, we push our previous matrix product operator methods to very large size and show that data for a time-independent Hamiltonian model are also consistent with a diffusively-broadened wavefront.

Information scrambling describes a process whereby information about the initial condition of a unitarily evolving system spreads over the entire system, becoming inaccessible to any local measurement [1][2][3][4]. Because it describes an effective loss of memory, scrambling is relevant for understanding quantum thermalization (e,g., [5][6][7][8]), i.e., the emergence of irreversibly from unitary time evolution, and is also tied to the black hole information problem. Scrambling is also closely related to the dynamics of initially simple Heisenberg operators, with the growth in size and complexity of these operators probing the spreading of quantum information [9][10][11][12][13][14][15][16].
Given two local operators W 0 and V r at positions 0 and r, the out-of-time ordered correlator (OTOC), provides one way to quantify scrambling by probing how arXiv:1805.05376v1 [cond-mat.str-el] 14 May 2018 the Heisenberg operator W 0 (t) grows with time. Although scrambling can also be usefully characterized in entropic terms, OTOCs are more directly measurable, with early experiments having already been carried out in a variety of platforms [17][18][19][20][21][22][23][24][25][26][27]. A closely related quantity is the squared commutator between V and W , defined as, The physical picture is that under Heisenberg dynamics, the operator W 0 expands and eventually fails to commute with V r , as manifested by the growth of C(r, t) from zero. For chaotic local Hamiltonians, W 0 (t) is expected to expand ballistically, with speed called the butterfly velocity, so that the OTOC exhibits a casual light-conelike structure in space-time. The squared commutator remains small outside the light-cone and grows rapidly the boundary of the light-cone is crossed. Inside the lightcone, C(r, t) saturates for chaotic systems regardless of the specific form of operator W and V .
A particularly interesting question concerns the specific growth form of C(r, t) near the wavefront of the light-cone. In some models, C(r, t) grows exponentially with time, a phenomenon proposed as a quantum analog of classical butterfly effect, the exponential divergence of initially nearby trajectories. This observation has led to an emphasis on probing the footprint of quantum chaos at an intermediate time-scale, especially with a view towards defining a notion of quantum Lyapunov exponent. A well-defined Lyapunov exponent λ L , i.e., purely exponential growth of C(r, t), plus the ballistic growth of the OTOC, implies that the wavefront is sharp, and sharp wavefronts have been identified in a broad class of holographic/large N models, including the O(N ) model [28], the diffusive metal [29] and the coupled SYK model [30][31][32]. On the other hand, although significant efforts have been made [33][34][35][36], a clear signature of purely exponential growth of the OTOC in more physical systems with finite on-site degrees of freedom is absent, and there are some counterexamples in random circuit models [9,10,[37][38][39].
To reconcile the many different scenarios, in a recent paper [11] we proposed a universal form for the early growth region of the squared commutator, assuming there is a well-defined butterfly velocity v B (a different ansatz is needed for localized systems [40][41][42][43]), furthered studied by [44]. The shape of the wavefront is controlled by a single parameter p, denoted as the broadening exponent, associated with the growth rate λ p . For large N /holographic models, p = 0 and the corresponding λ p is the Lyapunov exponent. However, an exact calculation in a Haar random brickwork circuit model gave p = 1 in one dimension, indicating a diffusive broadening of the wavefront. Saddle point analysis shows p = 1 2 for general non-interacting systems with translational invariance [11,36,44]. Large-scale matrix-product state simulations using the time-dependent variational principle [45] and matrix product operator simulations [11] also gave strong evidence of wavefront broadening for chaotic local Hamiltonian systems.
In this work, we make two contributions to understanding the early growth region behavior of the OTOC. First, to understand the intriguing differences between the large N models and the Haar random circuit models (HRC), we design and analyze a new random circuit model, denoted as the Brownian coupled cluster model (BCC). BCC, as an extension of the single cluster version [46,47], describes the dynamics of clusters of N spins connected in a one-dimensional array (or more generally, connected according to any graph), similar to coupled SYK cluster models but with the couplings random in both space and time. We show that in the large N limit, BCC is similar to other large N models and has a well defined Lyapunov exponent, but the finite N correction qualitatively changes the broadening exponent from p = 0 to p = 1 in one dimension. We find that finite N corrections are actually quite dramatic, with the broadening of the wavefront characterized by a diffusion constant which scales as 1/ log 3 N at large N . We also find that there is a finite region in space-time where the wavefront remains sharp, indicating strong finite-size effects on the broadening exponent.
With this new point of view, our second contribution is to push our numerical matrix product operator simulations of operator growth in a local Hamiltonian Ising system to include 200 spins in the wavefront and up to time 250 in units of the nearest neighbor Ising coupling. By directly analyzing the way contours of constant C deviate in space-time, we find that the broadening exponent indeed converges to p = 1 in the large space-time limit. Therefore we conclude that diffusive broadening of the wavefront is generic for one-dimensional chaotic systems.
In more detail, our analysis of the BCC proceeds by focusing on operator dynamics, suitably averaged over the random couplings in the Hamiltonian. Any operator may be expanded in a complete basis of operators, with the expansion coefficients called operator amplitudes and with the square of the amplitudes forming a probability distribution, the operator probability distribution. The starting point of the analysis is the derivation of an equation of motion for the circuit-averaged operator probability distribution of a Heisenberg operator. The effect of averaging over the couplings in the quantum Hamiltonian is to give a closed stochastic equation for the operator probability distribution; in essence, the operator amplitudes evolve via unitary time evolution for each choice of couplings, and the averaging dephases this dynamics to yield a master equation for the operator probability distribution. One point should be emphasized: The randomness of the couplings in the Hamiltonian, which we sometimes call "disorder", is physically distinct from the quantum randomness manifested in the operator probability distribution. The latter will, in a certain limit, be instantiated as a random process which we call "noise".
One of the key assertions of this paper is that the disorder average is a technical convenience while the noise average contains essential physics of quantum fluctuations.
Starting from the master equation for the operator probability distribution, the analysis proceeds from two limits. In a large N limit, a mean-field-like treatment of the operator distribution becomes exact, and the operator dynamics can be translated into a closed non-linear partial differential equation for the operator weight, a quantity closely related to the OTOC. The resulting dynamical equation is similar to the Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) equation [48,49] which occurs, for example, in studies of combustion waves, invasive species, and quantum chromodynamics, among others, and was recently introduced in the context of scrambling to describe the growth of OTOCs [28,50,51]. The key physical effects embodied by FKPP-type equations are unstable exponential growth, diffusion, and eventual saturation; together these lead to traveling wave solutions with a sharp wavefront which describe the spreading of local Heisenberg operators. The leading finite N correction results in a stochastic partial differential equation, a noisy FKPP-like equation, in which the noise is multiplicative and 1/N suppressed. Drawing from the noisy FKPP literature [52,53], we argue that the noisy FKPP-like in the BCC has a diffusively broadened wavefront after averaging over noise. These analytical arguments are verified by direct numerical integration of the large-but-finite-N BCC stochastic equation. It should be emphasized again that the noise in the noisy FKPP-like equation represents quantum fluctuations, not different instances of the microscopic couplings.
In the small N limit, a different analytical treatment shows that the OTOC exhibits the same broadening diffusive broadening as in the Haar random circuit model (HRC). Moreover, we show that by representing the operator probability distribution as a "stochastic" matrix product state, it is possible to numerically solve the master equation for the time dynamics. Thanks to the dephasing provided by the disorder average over couplings, one can show that the late time operator probability distribution has low "correlation/entanglement" when viewed as a matrix product state. We further find that the operator probability distribution never has high correlation/entanglement, so that matrix product state techniques can accurately capture the operator dynamics for all times. A modest bond dimension of χ = 32 is already sufficient to converge the dynamics for 200 sites for all time.
Finally, taking these lessons from the BCC, especially the crucial role of noise, meaning quantum fluctuations, we argue that the diffusive broadening of the operator growth wavefront is generic in one-dimension. This idea has been previously conjectured based on work with random circuit models [9,10,[37][38][39]. One piece of evidence is direct numerical simulation of the time-independent Hamiltonian dynamics of Heisenberg operators in a system of 200 spins for very long time. A new analysis of the space-time contours of constant squared commutator conclusively demonstrates diffusive broadening of the wavefront at the largest sizes. Another piece of evidence is the prevalence of noiseless FKPP-like equations describing OTOC dynamics in large N or weakly coupled models, including linearized FKPP-like equations obtained in resummed perturbation theory [28] and fully non-linear FKPP-like equations obtained from self-consistent Keldysh treatments [50]. We argue that, starting from these known results, quantum fluctuations should invariably be described by adding noise, specifically multiplicative noise of the type found in the BCC. Hence these models will also suffer similar dramatic finite N effects resulting in diffusively broadened wavefronts.
The remainder of the paper is organized as follows. Section II describes in detail the notions of operator dynamics used throughout the paper and their relations to OTOCs. Section III introduces and analyzes the Brownian coupled cluster model (BCC), both at small and large N . Section IV discusses implications of the results for generic Hamiltonian systems. We conclude with some outlook, including the effects of conserved quantities and going beyond one dimension, as well as open questions.

II. GENERAL FORMALISM OF OPERATOR DYNAMICS AND ITS RELATION TO THE OUT-OF-TIME-ORDERED CORRELATOR
Consider a generic quantum system consisting of L sites with N spin-1/2s per site. The dynamics of the system is governed by a local unitary circuit U (t). In the Pauli basis, a Heisenberg operator W (t) takes the form where S is a product of Pauli operators with length N L and the time dependence is encoded in the coefficients c(S, t). The normalization is tr(W (t) † W (t)) = 2 N L , so that S |c(S, t)| 2 = 1. Under Heisenberg time evolution in a chaotic quantum system, an initially-localized operator grows and eventually equilibrates as far as local probes are concerned. The initial configuration cannot be recovered from local data, a manifestation of scrambling. The OTOC is designed precisely to quantify this process since the deviation of the correlator from its initial value for certain position indicates that the Heisenberg operator has developed a nontrivial component at that site. To understand the relationship between the OTOC and the operator string picture, we consider the following averaged OTOC, where a runs from 1 to N and α from 0 to 3, representing the identity and the three Pauli matrices.
Using the decoupling channel identity, we obtain that where S r,a is the r, a Pauli operator in the string S and the weight is w(σ α ) = 1 − δ α0 . Therefore φ(r, t) measures the average number of nontrivial local operators within a single spin cluster located at r, growing from 0 and saturating to 3/4 in a thermalizing system. Similarly, the averaged squared commutator C(r, t) equals 8 3 φ(r, t). There are two implications. First, C(r, t) does not depend on the phase information of the coefficients c(S) (which are in fact real for Hermitian W (t)), since it is completely determined by the probability distribution |c(S)| 2 . Second, C(r, t) does not care about the specific operator configuration S, since only the number of non-trivial operators matters. Both features significantly simplify the calculations.
In general, dealing with the Schrödinger equation is overcomplicated for the purposes of studying operator dynamics due to the unimportant phase information and associated high operator entanglement entropy. However, deriving a closed set of dynamical equations for the operator probability |c(S)| 2 from the Schrödinger equation is difficult. Random circuit models are useful to overcome this difficulty by introducing disorder as a dephasing mechanism, making the probability distribution dynamics tractable both numerically and analytically, for example, in recent studies of the HRC. This paper introduces another class of random circuit models, the Brownian coupled cluster (BCC), describing the dynamics of a system of coupled spin clusters with interaction random both in space and time. BCC can be regarded as a smoother version of HRC, in which the interactions are only between pairs of spins (or more generally, few-body) even in the large N limit, and thus is naturally more tied to holographic and SYK models. While in the HRC, the dimension of the on-site Hilbert space does not qualitatively affect the operator dynamics, in the BCC, we expect a smooth crossover from the HRC result for small N to holographic and SYK physics physics in the large N limit. Fig. 1 shows a schematic of the BCC model. It is best described using discrete time steps dt, with the limit dt → 0 taken later. For small dt, the whole time evolution unitary breaks into pieces,

III. BROWNIAN COUPLED CLUSTER MODEL
where α and β are the Pauli matrix indices running over from 0 to 3 (including the identity for convenience), a and b from 1 to N label the spins in the cluster, r and r label clusters (sometimes called sites), and rr means nearest neighbors. At each time step, the models contains two sets of uncorrelated random variables J andJ with mean zero and variance 1 8(N −1) dt and 1 16N dt, respectively. With the help of the random couplings, one can derive a master equation for the averaged probability distribution h = |c(S)| 2 . To simplify the calculation, we assume that h only depends on the operator weight w r = a w(S r,a ) of each cluster instead of on the details of the operator configuration S. This approximation is valid after a short relaxation time even though W starts as a specific operator. To proceed further, introduce the operator weight probabilitỹ with D the number of operators with weight configuration w, The operator weight probability is a properly normalized probability distribution over the (N +1) L possible weight strings.
The derivation of the master equation forh is recorded in the appendix A, with the result being, The evolution equation manifestly conserves the total probability, {w}h = 1 for all time, independent of the specific form of the functions γ + and γ − . For this particular problem, these functions are (12) In the following subsections, we will analyze this master equation in the infinite-N limit, study its large-N expansion, and compare the result with small-N results. Complementing these analytical results are numerical simulations of the master equation for 200 spin clusters using tensor network methods.

A. The infinite N limit
In the infinite-N limit, the master equation Eq. (11) can be approximated by a Fokker-Planck equation, Using the Ito stochastic calculus, the Fokker-Planck equation can be mapped to a Langevin equation, with η r (t)η r (t ) = δ rr δ(t−t ). This mapping explicitly demonstrates that the noise η arises from the deterministic master equation for the operator weight probability as 1/N effect. It is important that this noise η is conceptually different from the randomness of the Brownian circuit introduced to obtain the master equation; it originates purely from the quantum fluctuation in the BCC.
Later we will show that the noise, although suppressed at large N , has a drastic effect on the operator dynamics. First, we study the infinite-N limit in which the noise is set to zero and the stochastic Langevin equation becomes deterministic. After taking the continuum limit of the Langevin equation, in which φ(r, t) is assumed to vary slowly with respect to r, we obtain a FKPP-type equation, describing a growth-diffusion-saturation process. For simplicity of presentation, we hereafter refer to Eq. (15) as an FKPP equation. There are two fixed points of the dynamics, an unstable solution φ(r, t) = 0 and a stable solution φ(r, t) = 3 4 . The stable solution describes the equilibrium state where every operator string is equally probable. An initially localized operator configuration translates to an initial condition for the FKPP-type equation which is the unstable solution everywhere away from the initial local operator.
Similar to the FKPP equation, Eq. (15) admits a traveling wave solution φ(r, t) = f (r − vt) when the velocity v is larger than v c = 18g 2 (1 + g 2 ). Ahead of the wavefront, r vt, the traveling wave decays exponentially with r. For initial operator profile that is sufficiently localized, the wavefront travels with the minimal velocity v c and approaches the traveling-wave solution at late times. A detailed analysis can be found in the appendix C. Ahead of the wavefront, the traveling wave decays as exp(6(1 + g 2 )(t − r/v c )), consistent with a sharp wavefront. Therefore, the infinite-N limit of the Brownian coupled cluster model exhibits a well-defined Lyapunov exponent. The butterfly velocity and the Lyapunov exponent are v B = 18g 2 (1 + g 2 ) λ L = 6(1 + g 2 ). (16) Within the framework described by Eq. (3), the infinite-N limit has broadening exponent p = 0. The existence of a Lyapunov exponent in the infinite-N limit is in sharp contrast with the brickwork random circuit model result, where the diffusive-spreading nature of the wavefront is independent of the dimension of the on-site Hilbert space. The reason is that the brickwork model has no notion of few-body interactions within an on-site cluster due to the use of Haar random unitary matrices in the circuit. In Fig. 6(b), we explicitly verify the sharp wavefront by numerically solving Eq. 15.

B. The large-N expansion
Having established the purely exponential growth of the squared commutator in the infinite-N limit, we now investigate the behavior away from this limit. Comparing

Perturbative region
FIG. 2. The scrambling "phase diagram" in space-time. An overall unimportant spatial offset is omitted. The bold line marks the wavefront. The space-time is roughly divided into three regions. In the diffusive region, the squared commutator exhibits a diffusively broadening wavefront, quantitatively different from the exponential growth in the large N limit. Away from the wavefront, there is a region where the squared commutator grows exponential with a modified Lyapunov exponent. Further ahead the wavefront, the chaotic region gives way to the perturbative region. As N increases, the chaotic region expands as indicated by dashed arrows and eventually dominates the wavefront behavior in the infinite N limit. with the infinite-N limit, the large-N expansion affects Eq. (15) in two significant ways. First, φ(r, t) in principle only takes discrete values 0, 1 N , 2 N .... Therefore, in Eq. 15, φ(r, t) is set to zero when it is below 1/N . This hard cutoff allows the traveling wave to propagate with velocity smaller than v c , and as such the cutoff is important for obtaining the correction to butterfly velocity. Second, the noise term in the Langevin equation (14) becomes important. The deterministic differential equation Eq. (15) is augmented by a multiplicative noise term (17) Due to its multiplicative nature, the noise term only affects the physics when φ(r, t) is non-zero, and therefore the noise does not violate the casual structure of the noiseless FKPP equation . The effect of the noise is most prominent near the forward edge of the wave, with its most important effect being to make the position of the wavefront a random variable described by a biased random walk. The resulting noise-averaged front spreads diffusively in addition to the drift v B (N )t. Following an analysis of the original noisy FKPP equation [52,53], which we review in the appendix C, we are able to obtain the large N correction to Eq. (16), This is a remarkable result, indicating that the system approaches the infinite-N limit very slowly. It also shows that at large but finite N , the broadening exponent becomes p = 1.
In each realization of the noise, φ(r, t) still grows exponentially. But as √ 2Dt grows larger than the width of the traveling wave, the exponential growth of φ(r, t) is smoothed out by the diffusive movement of the wavefront's position, leading to a diffusive broadening of the noise-averaged wavefront. To quantitatively understand the effect of the noise induced by finite N on the wavefront, we approximate the traveling wave solution in a single realization of the noise by the following phenomenological model, which accounts for the saturation behind the wavefront and the growth ahead of the wavefront.
Using this simple model, the noise-averaged squared commutator, which is proportional to the noise-averaged φ(r, t), is the convolution of φ(r, t) with a Gaussian distribution describing the diffusive motion, To simplify the notation, we introduce the dimensionless , with ξ describing the strength of the noise. The result of the convolution is where erf(x) is the error function 2 √ π x 0 e −t 2 dt, and z = u−u 0 −τ is the position in the traveling frame with some unimportant offset u 0 determined by the initial condition. The next step is to analyze the behavior of Eq. (21) in space-time. It exhibits a light cone structure with a butterfly velocity independent of ξ, since the butterfly velocity is entirely set by the cutoff approximation and does not depend on the diffusion constant explicitly. This can be seen from the fact that C(τ, τ ) asymptotically approaches 1 2 . The space-time of t − r plane can be approximately divided into three regions based on the behavior of Eq. (21), as illustrated in Fig. 2. The region near the wave-front is the diffusive region, where the last two terms of Eq. (21) roughly cancel each other and C(u, τ ) is dominated by the single error function, In the limit that , consistent with the universal form with broadening exponent p = 1. This clearly demonstrate that the wavefront spreads diffusively. The growth behavior near the wavefront is dominated by the noise, and the original Lyapunov exponents does not enter.
This should be contrasted with the chaotic region where the first two error functions are far away from saturation, but the last error function is already saturated. In this region, the squared commutator is C(u, τ ) ∼ 4 3 e ξτ −z , a pure growth form with a modified Lyapunov exponent . This size of this region scales as log N , and the value of C in this region can be arbitrarily small in the long-time limit since it is enclosed by two lines with a bigger velocity v = v B (1 + 2ξ) than the butterfly velocity. Therefore, it is difficult to extract this region from numerical data of finite N spin chains.
There is also a third region denoted the perturbative region, where z is the largest scale in the system. In this case the squared commutator is infinitesimally small and behaves as 8 3N ξτ πz 2 exp(−z 2 /4ξτ ).

C. The small N limit
The large N analysis presented in the last section cannot be naively generalized to the case of small N . In the small N limit, the master equation Eq. (11) is still valid, but the approach of approximating the master equation with the Fokker-Planck equation, Eq. (13), to derive the Langevin equation, Eq. (14), is not.
Instead, we take a rather different approach by considering the probability of the operator string ending on a specific site at a given time, similar to what was studied in the HRC model. This probability ρ(r, t) is defined as Note that the sum of ρ(r, t) is conserved. From Eq. (11), one can derive the rate equation for ρ(r, t) as, where a subindex on ρ indicates a restriction: the operator string must end with that particular configuration. For example, ρ 11 (r + 1, t) is the probability of the operator string with w = l on site r, w = 1 on site r + 1, and w = 0 for all sites beyond r + 1. Now we use the approximation of local equilibrium, which is crucially different from the large N case, to relate ρ l (r, t) and ρ l l (r, t) to ρ(r). We find that ρ l (r, t) = With this approximation, we obtain a closed equation for ρ(r, t), The conservation law of ρ(r, t) determines that ξ = 9N 8 In the continuum limit, the equation reads This leads to The average squared commutator is related to ρ(r, t) as The final result is consistent with the Haar random circuit result and has broadening exponent p = 1.

D. Finite N results from tensor network simulation
The large N analysis together the small N analysis convincingly demonstrates that away from infinite N , the behavior of the squared commutator is dominated by the error function near the wavefront, leading to a broadening exponent p = 1. The effect of N is mostly encoded in the butterfly velocity and the diffusion constant. But the N -dependence of v B and D obtained from the two limits are not consistent. The small N analysis suggests that both v B and d increase with N , while the large-N analysis indicates that v B saturates to a certain value and D decreases with N . Therefore it is interesting to study how the two quantities interpolate between the two limits, especially the diffusion constant which is expected to be non-monotonic.
For this purpose, we directly simulate the master equation 11 by representing the probability distribution as a matrix product state (MPS) with physical dimension N +1. Comparing with the usual TEBD method [54,55], simulating a stochastic process with MPS (S-MPS) [56] is quite different. In the appendix D, we discuss the difference and introduce several useful techniques, including the canonical form for S-MPS along with the truncation schemes that are helpful for preserving the 1-norm of the S-MPS instead of the 2-norm. Notably the truncation scheme developed by White et. al [57] can be directly applied to here to exactly preserve the 1-norm for all time.
Generally, S-MPS requires higher bond dimension to capture local observables accurately compared with unitary MPS with the same entanglement. This is because the 1-norm normalization appropriate to S-MPS tends to amplify errors. Nevertheless, the entanglement entropy is still a good measure for determining whether the probability distribution can be represented efficiently as an MPS. Initially, the probability distribution is a product state with operator weight 1 in the center cluster of system and operator weight 0 elsewhere. In the early growth region ahead of the lightcone, the probability is not affected by the stochastic evolution and continues to enjoy low entanglement. Inside the lightcone, the probability distribution reaches the steady state where every operator string is equally probable and also admits a simple product state representation. Therefore, the entanglement entropy only accumulates around the wavefront. In practice, we find that the entanglement entropy never exceeds one bit, allowing us to obtain the whole scrambling curve up to N = 10. The resulting scrambling curve can be fit with the error function with almost perfect quality to extract the butterfly velocity and the diffusion constant. The fitting result is shown in the appendix D.
The result of the butterfly velocity is shown in Fig. 3(a) together with that obtained from small N and large N analysis for g = 1. We find that at N = 1 (single spin in the spin cluster), v B from S-MPS agrees with small N analysis perfectly. As N increases, v B deviates from the linear growth predicted at small N and smoothly connects to result from the large-N analysis. Fig. 3(b) summarizes the N -dependence of the diffusion constant from the different methods of analysis. The result from the S-MPS indeed exhibits non-monotonic behavior. At N =1, it agrees with the small N analysis. It peaks at N = 3 and drops as N increases, approaching the result from the large-N analysis. Therefore the numerical results agree with the analyses above from both limits. In principle, for large enough N , one should be able to identify the chaotic region.

E. Comparing with Haar random circuit models
It is instructive to compare the Brownian circuit model studied here with the previous studied HRC models. By designation, in the HRC models, the Haar unitary matrices equilibrate the operator string on the two sites it connects to immediately if there is non-trivial weight. In this case, the analysis in section III C becomes exact, the entire region ahead the wave-front is governed by the error function. On the other hand, in the BCC model studied here, the operator string takes a finite time to reach equilibrium even locally, the time scale being ∼ log N . The direct consequence is that although near the wave-front, the behavior of the squared commutator is dominated by the error function, there is still a region in space-time ahead of the wavefront where C(r, t) grows exponentially, as illustrated in Fig. 2. As N increases, this chaotic region expands and finally dominates the wave-front in the infinite-N limit.

IV. IMPLICATIONS FOR LOCAL HAMILTONIAN SYSTEMS
The analysis so far has shown two things. First, that there exists a random circuit model with a parameter N such that, at infinite N , the model exhibits exponential growth of the squared commutator with p = 0. Second, that for any non-infinite N , the dynamics of the model inevitably crosses over to a diffusively broadened wavefront with p = 1. It is quite plausible that any sufficiently generic random circuit model with finite on-site Hilbert space will also exhibit a diffusively broadened wavefront with p = 1 (with this already being established for the BCC and the HRC). The key question is, what aspects of this analysis hold when the couplings are not random in time? We now argue that p = 1 is generic for chaotic quantum many-body systems in d = 1 with finite local Hilbert space dimension.
The argument has two thrusts. First, we directly numerically simulate a small N Ising spin chain with conserved energy. Combining large scale numerical simulations with a new analysis technique, our previous result is improved to show that the system asymptotically approaches p = 1. Second, based on previous work in energy conserving systems showing the existence of noiseless FKPP-type equations governing the spreading of chaos at large N and/or weak coupling model, we argue that quantum fluctuations inevitably introduce multiplicative noise into these equations. The physics of the noisy FKPP equation then naturally leads to p = 1. For the latter argument, recall that we were careful to distinguish the noise in the 1/N corrected FKPP equation, which was a manifestation of quantum fluctuations, from the space-time random couplings in the microscopic Hamiltonian.
Although we focus on energy conserving systems here, we conjecture that our analysis also applies to Floquet models where the couplings are not random in time, but energy is not conserved because the Hamiltonian is time-dependent. In the case of a conserved energy, it also makes sense to talk about non-infinite temperature states. We briefly discuss how the story might be modified in this case, with a focus on the physics of the chaos bound.
A. Wavefront broadening for small N energy conserving systems According to the analysis in section III B, for finite onsite Hilbert space dimension, each contour of the squared commutator intersects with the chaotic region for a limited amount of time and eventually merges into the diffusive region. This suggests a strong finite-size effect that will hinder extraction of the broadening exponent p from fitting the squared commutator with the universal growth form. This is consistent with what we observed in our earlier analysis, where the value of p drifts along the contour. To unambiguously analyze the broadening of wavefront, we improve our numerical result by pushing both the system size and the simulation time so that the wavefront travel through ∼ 200 sites. We measure the spatial difference δx between two contours as a function of time. For the two contours we choose, we check that an MPO with bond dimension χ = 8 and bond dimension χ = 16 give identical results. Then assuming the general form in Eq. (3) applies, the asymptotic value of the slope is related to the broadening exponent, The result of this analysis for the mixed-field Ising chaotic spin chain is shown in Fig. 4. The inset of Fig. 4(b) plots δx with t on a log-log plot. The slope gradually increases and approaches to 1/2 in the large space-time limit, indicating that the wave-front broadens diffusively, p = 1, at the largest sizes and times. The initial deviation may be due to the early-time microscopic physics where C(r, t) behaves as t x /x!. However, the fact that the deviation persists to an intermediate scale suggests that there may exist a chaotic region in the spacetime causing a strong finite-size effect. Extracting this region for local Hamiltonian systems is an interesting future research direction.
To further validate this method, the same analysis is performed for the transverse-field Ising model, which describes non-interacting fermions and has a broadening exponent p = 1/2. The result is shown in Fig. 5 where one can indeed see that the slope of curve on a log-log plot approaches to 1/3.

B. Conjectured wavefront broadening for large N energy conserving systems
Given that one generic energy conserving model exhibits p = 1, it is plausible that this is a universal behavior among local chaotic Hamiltonians in one dimenion. To bolster this conjecture and to give a physical picture for it in one limit, it is useful to return to the noiseless FKPP equation. Indeed, a number of different models have been shown to have operator growth described by a noiseless FKPP-type equation, in some cases linearized and in some cases fully non-linear, at large N or weak coupling. What we argue is that such equations should inevitably be augmented by a noise term describing quantum fluctuations which has the form considered in this work. This would imply that this broad class of large N models also have p = 1 at the largest sizes and times. Assuming this is true, p = 1 then occurs at large and small N and weak and strong coupling, so it is reasonable to conjecture that it is universal property of onedimensional chaotic systems.
The argument proceeds in three steps of increasing specificity. The background assumptions are that one has a closed dynamical equation governing φ ∝ C and some parameter N which measures the local degrees of freedom. First, quantum fluctuations are expected to add a noise term to any approximate set of closed equations governing the dynamics of the OTOCs of simple operators. This is because operator growth is not deterministic in a quantum system, since a single Heisenberg operator is a superposition of many different complex operator strings. Second, the specific form of the noise term must be multiplicative for local Hamiltonians, meaning proportional to some power of φ. This is because operator growth arises from the failed cancellation between U and U † due to the insertion of the perturbation W . Operator growth never spontaneously occurs far away from the current support of W (t) precisely because U and U † cancel in far away regions. Third, the form of the multiplicative noise term should be φ/N . This is necessary to ensure that the noise is most important when φ is small, which should be true since larger values of φ are more self-averaging. This also guarantees that the associated Fokker-Planck equation has a sensible 1/N expansion, i.e., with no unusual powers of 1/N appearing. Then assuming the dynamical equation for C includes saturation effects, we have all the components necessary for the noisy FKPP analysis to apply.
One caveat here is that the analysis is framed in terms of large N models. While FKPP type equations have also been derived in weak coupling approximations, it is less clear how to identify the precise role of the N parameter in that case. One simple intuition is that N should arise because the dynamics is effectively coarse-grained over a long length scale corresponding to the inelastic mean free path. Identifying N ∼ , with some kind of inelastic mean free path, would then predict butterfly velocity corrections and a diffusion constant going like δv B ∼ 1/ log 2 and D ∼ 1/ log 3 . It would interesting to better understand the situation at weak coupling.

C. Finite temperature
Here we make some comments about the dynamics of operator growth at non-infinite T . The key point is that the FKPP equation and its noisy counterpart make no particular reference to temperature, so it is reasonable to suppose that they could hold at non-infinite T . The noiseless FKPP equation has already been derived at finite T for a variety of models; temperature only enters in so far as the parameters in the FKPP equation are temperature dependent. Thus the results on the BCC from small to large N provide some hints on the interplay between quantum fluctuation and local scrambling.
The manifestation of this interplay is the distinction between the diffusive region and the chaotic region the the space-time structure of the OTOC. In the infinite-N limit, quantum fluctuations are completely suppressed and the local scrambling time ∼ log N is infinite. In this case, the chaotic region occupies the entire spacetime. For large but finite N , as soon as quantum fluctuations are present, the wavefront broadens diffusively while the chaotic region occurs ahead the front. In the small N limit, with the local equilibrium assumption implying that the local scrambling time is short, the diffusive region extends to the entire space-time. This is also consistent with the analytical results on the HRC. Now consider a local Hamiltonian system, say a spin 1/2 system. Assuming the mixed-field Ising behavior is generic, we showed that at infinite temperature the wavefront also broadens diffusively. But it is difficult to tell whether there exists a chaotic region in addition to the diffusive region ahead of the wavefront due to the microscopic details affecting the behavior of OTOC at early time. On the other hand, at finite temperature and assuming a separation of time-scales, there is a thermally regulated version of the OTOC whose growth rate is bounded by the temperature [58]. If we assume the same bound applies to the non-thermally regulated object 1 , we would have Given this chaos bound, the diffusive broadening form can only be valid up to finite distance ahead the front, because its growth rate diverges in the large r limit. Imposing the chaos bound, we can estimate maximum distance from the wavefront for which the diffusive behavior is valid. The growth rate of the diffusive growth from is Demanding that γ(r, t) ≤ 2πT predicts that the diffusive region can at most persist up to the space-time line given by One might also imagine that even at infinite temperature, the logarithmic derivative cannot be too large, for example that it should be bounded by the size of the microsopic couplings. In the BCC, this was indeed true at large N where the crossover from the diffusive to chaotic region was roughly where the rate of growth in the diffusive region was approaching the Lyapunov rate (which itself was of order the microscopic scale). A similar argument might also apply when the system has a Lyapunov exponent less than 2πT . In any event, at finite T there may generically be a region ahead of the wavefront, at least at large but finite N or weak coupling, where the Lyapunov exponent is still visible. One issue is that one could run into the perturbative region before any exponent can be extracted. More generally, one should carefully study the thermally regulated commutators to see the precise consequences of the chaos bound, something we will report on in forthcoming work.

V. CONCLUSION AND OUTLOOK
In this work, we studied a random time-dependent Hamiltonian model with a large N limit in which we could study in detail the way the infinite N Lyapunov exponent gave way to a diffusively broadened scrambling wavefront. Based on this model and an analysis of large scale MPO simulations in a time-independent Hamiltonian model, we conjectured that the local operator growth wavefront broadens diffusively in generic local chaotic Hamiltonians with finite local Hilbert space dimension. We also showed how a modified stochastic MPS formalism could be used to simulate the operator dynamics for all times after averaging over different Hamiltonian realizations in the random model. A unifying element was the emergence of a noiseless FKPP equation at infinite N and a corresponding noisy FKPP equation at finite N . The noise was an effect of quantum fluctuations and ensured that both large N and small N exhibited p = 1 dynamics.
It is straightforward to extend the BCC model to any dimension, or indeed, to any graph. In higher dimensions, there will still be a Lyapunov exponent and sharp wavefront at infinite N . Finite N corrections will then introduce noise into the FKPP-type equation. The analog of the cutoff on φ is an extended cutoff front where φ = 1/N . This front will then experience some random dynamics with a constant drift (the butterfly velocity) and noisy local dynamics. Although the general longdistance structure may be complex, e.g., in high dimensions the noise may be relevant or irrelevant, one expects KPZ-like dynamics in low dimensions. It will be interesting to analyze the higher dimensional case in more detail, and possibly also study the model on more general graphs.
In terms of future directions, we have several works in progress. One is to consider the effect of a conserved U (1) symmetry on the operator spreading. This has already been studied in the HRC [38,39], but we anticipate new interesting physics associated with the interplay with the large N effects. Another is a study of the entanglement dynamics in the model as a function of N . Another interesting direction is to directly obtaining the noise physics in a Hamiltonian large N model.
Here r labels the cluster, a, b are labels within a cluster, and α, β label Pauli matrices. The coupling constant g can be used to dial the relative strength of the within-cluster and between-cluster interactions. Some of the other factors are chosen for convenience, and the sum over a, b is unconstrained in the between-cluster term. The coefficient A is determined by demanding that U U † = 1 on average, leading to Given an operator W , the Heisenberg operator is W (t) = U W U † . We may expand W (t) in a complete basis of operators, where S is a product of cluster operators S r with each S r a product of Pauli operators within the cluster. The coefficients c(S) can be determined from c(S) = 1 2 N L tr (SW (t)) . (A3) We will study the average operator probabilities, To determine the equation of motion of h(S) we must compute dh(S). There are two kinds of terms, depending on whether dB or dB appear in the same trace or not. When they appear in the same trace, we find where q α,β r,a,b = ±1 depending on whether σ α r,a σ β r,b commutes or anti-commutes with S andq α,β r,a,b = ±1 depending on whether σ α r,a σ β r+1,b commutes or anti-commutes with S.
Let w r denote the total weight of S on site r, w r = 0, · · · , N . The sums above can be written in terms of the w r . The first sum is r,a<b,α,β q α,β r,a,b = 16 This is because, if S has a non-zero Pauli on spin r, a or r, b, then α,β q α,β r,a,b = 0, otherwise it is 16. Thus we need to count the number of pairs a, b which are both the identity operator. This number is (N −wr)(N −wr−1) Similarly, the second sum is r,a<b,α,βq α,β r,a,b = 16 Putting everything together gives When the dB or dB factors appear in different traces, then we get terms connecting h(S) to h(σσS). These terms are (A7) To proceed further, let us assume that h depends only on the total weight w r on site r and not on the particular operator S. The dh 1 terms already manifestly depend just on the total weight. The dh 2 terms connect probabilities for different weights. The dh 1 term can be written To compute the dh 2 terms, we can consider a particular operator of the desired weight. Consider first the onsite terms. Suppose S r = σ x r,1 · · · σ x r,wr I r,wr+1 · · · I r,N . Now consider all σ α r,a σ β r,b . If the σσ commutes with S, then the dh 2 contribution vanishes. If σσ anticommutes with S, then the q(1 − q) factor is −2. If S r is the identity on one of a or b, say b, then the anticommuting terms are σ y σ β and σ z σ β . Of these eight, two keep the weight the same and six increase the weight by one. There are w r (N − w r ) choices of a, b in this class. The contribution is thus − 2w r (N − w r )[2h(w)dt + 6h(w + e r )dt]. (A9) Here e r = (0, · · · , 0, 1 r , 0, · · · , 0) is a unit vector that adds one to weight w r . If S r is non-identity on both a and b, then the anticommuting terms are (y or z)(I or x) and (I or x)(y or z). Of these eight, four keep the weight the same and four decrease the weight by one. There are w r (w r − 1)/2 such choices of a, b. The total contribution from these terms is thus The total onsite contributions are thus For the non-onsite contributions, we choose S to be an operator with w r σ x s on cluster r and w r+1 σ x s on cluster r+1. If S is the identity on r, a, then the anti-commuting operators are σ α r,a σ y r+1,b and σ α r,a σ z r+1,b . Of these eight, two keep w r the same and two increase w r by one. There are (N − w r )w r+1 such a, b. The contribution is Similarly, if S is the identity on r + 1, b, then by the same logic we find a contribution of Finally, suppose S is not the identity on both r, a and r + 1, b. Then the anti-commuting operators are (y or z)(I or x) and (I or x)(y or z). Of these eight, four leave both weights unchanged, two decrease weight w r by one, and two decrease weight w r+1 by one. There are w r w r+1 such a, b. The contribution is (A13) The total non-onsite contribution is thus r,r N w r (2h(w) + 6h(w + e r ))− w r w r (6h(w + e r ) − 2h(w − e r )) + r ←→ r (A14) By combining Eqs. (A8), (A11), and (A14) together, we obtain the complete equation of motion for h(w) in the case where h(S) depends only on the weights w r of P at the different clusters. It is, however, convenient to take one more step and include degeneracy factors.
The number of operators with weights w r are While h(w) is the probability of a single operator with weights w r , the objecth(w), defined bỹ is the total probability of all operators with weight w r , i.e., the probability of weights w r . The equation of motion forh can be obtained from that of h. The dh 1 terms immediately translate to dh 1 terms since the weights are the same on both sides of the equation. However, the dh 2 must by modified. We replace each h(w) with ah(w)/D(w). Then the various rates are modified by ratios of D(w) to D(w±e r ). These ratios are and Thus we have and dh 2,non-onsite dt = g 2 2N rr 3(N − w r + 1)w r h (w − e r ) + N w r h + w r (w r + 1)h(w + e r ) + r ←→ r (A20) Combining every thing together, we obtain the master equation ofh presented in Eq. (11) in the main text.
Appendix B: More details on the Brownian coupled cluster in the infinite-N limit In this appendix, some additional details on the infinite N noiseless limit of the BCC model are presented. The continuum limit of Eq. (14) resembles a FKPP-type partial differential equation, Eq. (15).
To justify the continuum approximation, here we directly study the original discrete ordinary differential equation on the lattice. Starting with (B1) we look for the traveling the wave solution exp(λ(t−r/v) in the small φ limit, so that the nonlinearity can be safely ignored. The relation between the velocity v and the growth rate is The minimum velocity for positive growth rate is the butterfly velocity v B and the corresponding growth rate is the Lyapunov exponent. In Fig. 6(a), we plot v B and λ L as a function of g and compare them with the analytical result from the continuum limit. Overall, the result from the continuum approximation tracks that from the discrete model. The main difference occurs in the limit that g 1: Analyzing Eq. (B2) shows that v B ∼ − 3 2 log |g| + 1 while the continuum approximation predicts that v B ∼ 3 √ 2|g| and λ L ∼ 6(1 + g 2 ). Furthermore, we directly simulate the discrete nonlinear ordinary differential equation. The result is shown in Fig. 6(b). As time increases, the spatial difference between two contours of log C(r, t) saturates. This explicitly verifies the sharp wavefront and the exponential growth of the squared commutator.
Appendix C: The butterfly velocity and the diffusion constant from the noisy FKPP equation In this appendix, we discuss how to obtain Eq. (18) for large but finite N BCC by analyzing Eq. (15) with the hard-cutoff approximation and the noise term Eq. (17). This material is essentially a review of the analysis of Brunet et. al. [52,53]. To simplify the notation, we rescale φ, r and t as following, showing that the front is sharp. On a log-log plot, the slope of the curve decreases to zero. This is in the sharp contrast of the local Hamiltonian systems, where the slope increases to 1/2 and 1/3 for the chaotic case and non-interacting case respectively as shown in Fig. 4 and Fig. 5.
After the rescaling, Eq. (15) becomes and the noise term becomes In this section, we will mainly focus on Eq. (C2) and Eq. (C3).

The noiseless case
We first discuss the case without noise, corresponding to the infinite-N limit of the BCC model. Eq. (C2) is similar to the Fisher-KPP equation, describing a growth-saturation process occurring in a wide class of systems including population dynamics, combustion, and reaction-diffusion systems. One of the interesting features of the Fisher-KPP equation is that it admits traveling wave solutions φ(r, t) = w(r − vt) that the initial configures converge to [59]. We expect that our equation Eq. (C2) obtained from unitary dynamics also falls into the FKPP universality class, because the linearized version of Eq. (C2) is the linearized Fisher-KPP equation, and because Eq. (C2) also has saturation physics so that φ = 1 is a stable solution.
For the Fisher-KPP equation, given a initial condition φ(r, 0) that is sufficiently well localized, it asymptotically approaches the traveling wave solution, where m(t) is the position of the wave front defined by the equation φ(m(t), t) = constant.
Traveling waves.-The first question to answer is the shape of the traveling wave solution w v (z) for different velocities. It can be determined from the following equation, with the boundary conditions w(∞) → 0 and w(−∞) → 1. In the region that z 0 and w(z) 1, the shape can be understood from the linearized equation. For each velocity, there are two modes e −γz and e −z/γ , where γ + 1/γ = v. At the critical velocity v c = 2, γ = 1/γ = 1 and the two modes are ze −z and ze −z . On the level of the linear equation, for a given velocity, any combination of the two modes is valid. The effect of the non-linearity can be regarded as setting a boundary condition for the linearized equation, say at z = 0, where α(v) and β(v) are tied to each other based on the solution to the full nonlinear equation. The above boundary condition forces both decay modes to appear, with the slower decaying mode dominating the behavior of w(z) in the large z limit. In the case that v < 2, γ becomes complex , and both modes can appear as long as the combination is real. Therefore we have In the current context, φ is interpreted as the operator weight and is always positive. Therefore, the last case where w(z) oscillates is physically irrelevant, setting the minimal physical velocity to v c = 2. However the last case becomes important for the noisy case discussed below.
Approaching to the traveling waves.-Great efforts have been made to understand the relationship between the initial configuration and the final traveling wave it asymptotes to. It is found that for an initial perturbation which is sufficiently localized, φ(r, t) approaches the critical traveling wave in the long time limit. This can be understood from the linearized equation, Using Green's function, we can write down the general solution as where φ(r , 0) is given by the initial condition. Starting with φ(r, 0) = e −λ|r| , we obtain that The velocity of the wavefront is determined by choosing v to cancel the t dependence in the exponent so that φ approaches a constant in the traveling frame. As a result, This demonstrates that, an initial sufficiently-localized configuration travels with the minimal velocity v c = 2, i.e., the leading term of the wavefront position m(t) is 2t. The asymptotic form is, From this, one can also obtain the sub-leading term of m(t) by canceling the 1/ √ t prefactor, which gives rise to m(t) ∼ 2t − 1 2 log(t). Then φ(m(t) + z, t) → e −z . The linearized equation gives the right velocity. However, the sub-leading term in m(t) is not correct. The nonlinearity forces the waveform at the critical velocity to decay like ze −z , slower than the e −z form obtained above. To take into account the nonlinear effects, we note that ∂ r φ(r, t) is also a solution to the linearized equation, and we can combine ∂ r φ(r, t) and φ(r, t) to obtain a new solutionφ that minimizes the effect of the nonlinear term φ(r, t) 2 . By expanding Eq. (C13) to the next order, one finds that which indeed has the corrected asymptotic behavior as a function of z. We can again obtain the sub-leading term in m(t) by canceling the time dependence, The second term was found by Bramson [59] and turns out to be independent of the specific shape of the initial condition as long as it is sufficiently localized.
To incorporate a simple saturation mechanism into the linearized equations and enforce the asymptotic shape of the traveling wave solution, Berestycki, et al. [60] recently solved the linearized equation with the following moving boundary condition, in order to obtain the vanishing correction of m(t). They found that The same correction has been identified in other models and is therefore expected to be universal, applying to the original Fisher-KPP equation and Eq. C2 as well.

The noisy case
The above picture applies to the infinite-N limit of the BCC model. As stated in the main text, 1/N expansion away from the limit has two main effects. First, we need cutoff φ(r, t) to 0 whenever it is below 1/N ; second, we need include the noise term Eq. (C3) in the Eq. (C2).
Cutoff-velocity.-We first discuss the effect of the cutoff without considering the noise. As discussed above, due to the positivity of φ, the velocity of the traveling wave can never go below v c = 2. However, with the cutoff, the part of φ below 1/N is set to zero by hand, and a smaller velocity becomes possible. The scaling behavior of the velocity as a function N can be obtained following Ref. [52]. When v < 2, the tail of the traveling acquires a oscillating part with a long wavelength in addition to the decay, with γ I and γ R the real and imaginary part of γ respectively. Even with the cut-off, w(z) should still remain positive until it decays to 1/N . This imposes a constraint on the wave-length of the oscillation part which is determined by γ I , requiring γ I < π log N . In consequence, v > v c − π 2 log 2 N . Therefore, the velocity correction scales as 1/ log 2 N , which is consistent with the numerical result presented in the inset of Fig. 3(a).
To quantify the above picture, and especially to understand how different initial conditions approach the asymptotic traveling wave with the cut-off, Brunet et al. [52] introduced a third boundary condition φ(m(t) + L, t) = 0 with L ∼ log N to the linearized equation in addition to Eq. (C16). This new boundary condition is to account for the hard-cutoff occurring at the leading edge of the traveling wave. They also set α = 0 for simplicity, which seems unnatural but does not affect the shape of the traveling wave in the large z limit. In the following, we repeat their analysis in some detail. With this setup, it is convenient to go to the traveling frame by performing the substitution φ(r, t) = w(r − m(t), t). Then the boundary conditions are We first approximateṁ(t) as its asymptotic value v, to be determined. By setting the time derivative ofṁ to zero, the boundary conditions uniquely determines the form of the asymptotic traveling wave, with the velocity v = 2 1 − π 2 L 2 . Understanding the asymptotic form, we restoreṁ(t) in the equation and study the full dynamics of w(z, t).
To the leading order of L, w(z, t) can be written as a superposition of eigenmodes, w(z, t) = a n L π sin nπ L z exp −z + π 2 L 2 (1 − n 2 )t + vt − m(t) (C21) where a n is obtained from Fourier expanding the initial condition, and m(t) is tuned so that ∂ z w(0, t) = 0 for all time. All modes with n > 1 decay exponentially. In the long-time limit, since m(t) → vt, we obtain Therefore, in order to match the boundary condition. In other words, the relaxation to the asymptotic traveling wave causes an additional shift to the wavefront position.
Noise-induced diffusive motion.-Now let us analyze the role of the noise term on top of the hard cut-off. Based on Eq. (C3), the noise scales as φ/N and therefore is most prominent at the leading edge the traveling wave where φ ∼ 1 N . We will argue that the wavefront obeys where X is a random diffusive process with X = 0 and X 2 = 2Dt. Different from the deterministic model studied above, in the actual noisy equation, it is possible that φ(r, t) is on the order of 1/N even when z > L. Let L + δ denote the maximal distance that w(r, t) = 0, where δ is a random variable. Then in the region that L < z < L+δ, the noise term scales as e δ /N , significantly larger than w(z, t) itself which scales as 1/N . The wavefront shift caused by the noise in a unit time is Based on a phenomenological reaction-diffusion model, Brunet et al. [53] obtained the probability for a large δ to appear in a unit time as decaying with the natural decay constant in the system. Then the additional velocity correction and the diffusion constant can be calculated straightforwardly as δv ∼ ∆z(δ)p(δ)dδ ∼ 1 L 3 D ∼ (∆z(δ)) 2 p(δ)dδ ∼ 1 L 3 .

(C27)
This suggests that the system approaches its infinite-N limit extremely slowly. Local observables: To measure a local observable at site i, one can first right-sweep from site 1 to i − 1 and then left-sweep from site L to site i + 1. Then the probability distribution is in the so-called mixed canonical form, where γ is the index for the states at site i, α is for the states to the left of i and β is for the states to the right of i. By construction, ρ γ 1,1 is the reduced probability distribution for site i, αβ ρ αγβ , from which the calculation of local observables is straightforward. Thus we have realized the goal sketched in Fig. 7, where calculation of local observables is reduced to manipulating local data.
Simulating master equation and truncation.-To simulate a local master equation such as Eq. (11), we view the generator of the stochastic evolution as a non-hermitian Hamiltonian which explicitly conserves 1-norm. Then we can roughly follow the usual time evolving block decimation (TEBD) steps to update the S-MPS after a short time step, but with the sweeping procedures replaced by those described in the last section.
Similar to TEBD, the bond dimension of the S-MPS typically grows rapidly with time and truncation is always necessary. Consider the situation after updating the tensors at site i and site i + 1. The probability distribution is in the following form, where γ 1 and γ 2 represent the states on site i and i + 1 respectively, and by construction ρ γ1γ2 11 is the reduced probability distribution of site i and site i + 1. Assuming that the stochastic simulation starts with an S-MPS with bond dimension χ, the dimension of the middle matrix ρ γ1γ2 m l mr is χd, i.e., the bond dimension at the bond linking sites i and i + 1 is χd. The goal is to reduce it back to χ by breaking ρ γ1γ2 m l mr into two pieces, ρ γ1 m l ,m and ρ γ2 m,mr , where the dimension of m is χ. Then one can continue to right-sweep or left-sweep to update the next bond.
There are two comparable methods to achieve this. The most straightforward method is to perform an SVD on ρ γ1γ2 m l mr by regarding it as a matrix with dimension χd. After keeping only the leading χ singular values, the bond dimension is reduced back to χ.
The second method is one developed by White et. al in a recent paper studying density matrix dynamics. They first break ρ γ1γ2 m l mr into two pieces and perform a single right/left-sweeping step (Step 3 to Step 5) on the left/right piece, so that ρ γ1γ2 m l mr is expressed as, ρ γ1γ2 m l mr = mm ρ γ1 m l ,m Q mm ρ γ2 m ,mr , with the properties that ρ γ1 1,m>d = 0 (due to Step 4 in the last section), γ1 ρ γ1 1,m = δ m,1 , and similarly for ρ γ2 m ,1 . In general, the rank of the matrix Q is χd.
As they pointed out, given Eq. (D5) and Eq. (D6), the right lower (χ − 1)d × (χ − 1)d section of Q does not affect the reduced probability distribution of site 1 to site i, and site i + 1 to site L, precisely due to Step 4 in the sweeping procedure. This gives some freedom to manipulate that section of the matrix in order to reduce the rank of Q without affecting local observables (they were concerned with developing a truncation scheme that respects local observables). For example, if one performed an SVD on that section and kept the leading χ−2d singular values, then the resulting Q matrix would have rank χ as required. One could also follow White et. al and reduce the rank while minimizing the error occurring in the correlation functions between the left part (1 ∼ i) and the right part (i + 1 ∼ L) of the system.
Comparing the two methods, the second one is more appealing since a truncation step at a particular bond does not change the reduced probability distribution of the left and right parts of the system. Furthermore, the truncation scheme automatically preserves conserved quantities, like the 1-norm or the total amount of some conserved charge, for all time regardless of the bond dimension. However, the accuracy of the time-dependence of local quantities still depends on the bond dimension. After one sweeps through the system and performs the truncation on each bond, only local observables with support on up to two neighboring sites remain unaffected. Furthermore, as time increases in the simulation, the errors made in observables with larger support could feedback to the dynamics of the two-site observables.
In practice, for our goal of calculating the squared commutator from the master equation Eq. (11), we find that the results obtained from these two truncation schemes are similar, and the first scheme is 2 to 3 times faster since it requires fewer steps of SVD.  We apply the technique described above to simulate the master equation Eq. (11). We first check the convergence of the resulting squared commutator with the bond dimension of the S-MPS. As shown in Fig. 8, an S-MPS with bond dimension χ = 32 already produces converged results for all time scales for N = 10 spins within a single cluster and a total of L = 200 clusters. This allows us to access both early growth and late-time saturation of scrambling in the BCC. We then fit C(r, t) in the region near the wavefront with an error function 1+erf r−v B t−r0 √ 4Dt to extract the butterfly velocity v B and the diffusion constant D. The fitting quality is shown in Fig. 8(b).