Operator Spreading in Random Unitary Circuits

I. INTRODUCTION

• The paper provides exact results and coarse-grained descriptions for operator spreading in Haar-random circuits, measured through the OTOC.Random circuits serve as minimally structured models for comparing with generic quantum dynamics.Source: §I · p. 1
• The OTOC addresses operator size, butterfly velocity, spatial structure, and possible hydrodynamic descriptions in random circuits.The operator’s significant-commutator region typically grows ballistically, motivating these questions.Source: §I · p. 2

A. Hydrodynamic equation for averaged OTOC

• In 1+1D, averaged operator dynamics map to a biased diffusion process: the front travels at finite vB and broadens diffusively as t^1/2.A fictitious Markov process on operator strings yields a lattice diffusion equation and a continuum drift-diffusion description.Source: §I · p. 2, A. Hydrodynamic equation for averaged… · pp. 4, 6
• Fluctuations between circuit realizations are argued to be small, with front-position variations O(t^1/4), below the t^1/2 averaged-front broadening.The averaged OTOC is therefore parametrically close to the spatial profile in a typical realization at large times.Source: §I · p. 2, A. Hydrodynamic equation for averaged… · p. 4

I. INTRODUCTION

• In higher dimensions, the averaged OTOC maps to KPZ-class droplet growth, giving front widths t^1/3 in 2+1D and t^0.240 in 3+1D.In 2+1D without lattice anisotropy, the asymptotic OTOC has a Tracy–Widom F2 form.Source: §I · p. 2, Figure 2 · p. 2, Figure 3 · p. 3
• Lattice operator shapes are generally anisotropic because the butterfly velocity depends on front orientation, whereas continuous rotational symmetry permits spherical asymptotics.The anisotropic shape is determined by a model-specific velocity function v(n̂) and is supported by 2+1D simulations.Source: §I · p. 3

II. OPERATOR DYNAMICS IN 1+1D

• An Ising-like partition-function mapping reduces the 1+1D averaged OTOC to an exactly solvable random-walk calculation through unitarity.The same mapping gives exact average entanglement purity and clarifies that vE can be arbitrarily smaller than vB.Source: §I · p. 3

A. Hydrodynamic equation for averaged OTOC

• The hydrodynamic treatment assumes negligible interaction between left and right operator ends and equilibration of the occupation density, while the later calculation removes these approximations.A mean-field traveling-wave treatment is explicitly uncontrolled and does not capture the true short- or long-time behavior.Source: A. Hydrodynamic equation for averaged… · p. 6

• Fluctuations are parametrically smaller than the diffusive front broadening, so the averaged OTOC profile is expected to describe individual random-circuit realizations.The diffusive broadening occurs in a single realization rather than arising solely from disorder averaging.Source: B. Hydrodynamic description including… · pp. 6–7
• The stochastic hydrodynamic description retains noise in the drift because it is more relevant under renormalization than diffusion-constant noise.Microscopic noise is expected in both parameters, but the analysis restricts attention to noisy drift.Source: B. Hydrodynamic description including… · p. 6
• Without noisy drift, the OTOC profile forms a wavepacket whose width grows diffusively and whose center moves as x_cm = v_B t.Source: B. Hydrodynamic description including… · p. 6
• The wavepacket samples a spacetime volume scaling as t^3/2, yielding typical random displacement Δx_cm ∼ t^1/4.A perturbative calculation reproduces the 1/4 exponent, also found for directed waves in random media.Source: B. Hydrodynamic description including… · p. 7

• In Clifford circuits, an initial Pauli remains a single operator string and the endpoint density stays localized on one site.Source: C. Comparison with Clifford circuit… · p. 7
• Averaged endpoint density, occupation number, and OTOC obey exactly the same dynamics in uniformly random Clifford and Haar-random circuits.The equivalence follows because random Clifford operators form a unitary 2-design.Source: C. Comparison with Clifford circuit… · p. 7
• Within individual realizations, Clifford endpoint density and OTOC fluctuate strongly, whereas generic-unitary dynamics is argued to be self-averaging.Source: C. Comparison with Clifford circuit… · p. 7

III. HIGHER DIMENSIONS

• The averaged OTOC in higher-dimensional Haar-random circuits maps to a classical Markov process whose droplet-growth physics is described by KPZ scaling.The mapping yields universal exponents and scaling forms for the averaged OTOC.Source: §III · p. 7, A. Higher dimensions: setup and… · p. 8, B. Classical model in 2+1D: analytical… · p. 8
• Lattice anisotropy prevents the late-time operator shape from becoming spherical, because the butterfly velocity depends on front orientation.Along a lattice axis, sufficiently large q can produce vB = 2 without roughening; this effect is tied to the chosen discrete spacetime geometry.Source: §III · p. 7, B. Classical model in 2+1D: analytical… · pp. 9–10

A. Higher dimensions: setup and mapping to classical growth

• Four successive circuit layers cycle through the square lattice’s four columnar dimer coverings, updating neighboring site occupations through a stochastic process.An initially localized operator corresponds to occupation n = 1 at the origin and n = 0 elsewhere.Source: A. Higher dimensions: setup and… · pp. 7–8, Figure 4 · p. 8
• The averaged droplet profile precisely reproduces the averaged OTOC, with the classical droplet’s binary occupation field defining the corresponding stochastic growth problem.The quantum circuit average and classical-process average are related only after averaging, while the fictitious process can also be realized by random Clifford dynamics.Source: Figure 6 · p. 8, A. Higher dimensions: setup and… · p. 8

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

• The droplet grows linearly in time with a sharp fluctuating interface, and KPZ scaling predicts interface fluctuations proportional to t^1/3.The interior rapidly equilibrates, while the average 2D shape has four-fold rather than continuous rotational symmetry.Source: B. Classical model in 2+1D: analytical… · pp. 8–9
• The 2+1D simulation finds β = 0.3305 ± 0.0269 along θ = ±π/4, consistent with the KPZ value β = 1/3.The study tracks q = 2 clusters over M = 2 × 10^3 realizations and observes linear cluster-size growth.Source: B. Classical model in 2+1D: analytical… · p. 9, Figure 7 · p. 9

C. Scaling of the OTOC in 2+1D

• In 2+1D, the OTOC front broadens as t^1/3, and the KPZ mapping supplies exact scaling functions as well as exponents.The proposed scaling forms are intended to capture universal behavior beyond random circuits.Source: C. Scaling of the OTOC in 2+1D · p. 10
• With negligible lattice anisotropy, the averaged radial OTOC has an exact functional form determined by the GUE Tracy–Widom distribution F2.The result assumes rotational symmetry and expresses the correlator through the probability distribution of the fluctuating interface radius.Source: C. Scaling of the OTOC in 2+1D · p. 10
• For a straight front from half-plane initial conditions, the scaling form instead uses the GOE Tracy–Widom distribution F1.This provides an alternative initial-condition-dependent Tracy–Widom description of operator growth.Source: C. Scaling of the OTOC in 2+1D · p. 11

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

• In 3+1D, the front broadens with KPZ exponent β ≃ 0.240, while the corresponding universal scaling function is known numerically rather than analytically.The higher-dimensional behavior otherwise parallels the 2+1D case.Source: D. Scaling of the OTOC in 3 + 1D and… · p. 11
• For dimensions 4+1 and above, KPZ theory permits rough and smooth phases depending on the nonlinearity strength, but these regimes are experimentally inaccessible.The paper leaves open whether both phases occur in suitable many-body systems.Source: D. Scaling of the OTOC in 3 + 1D and… · p. 11

E. Shape of the operator at late times

• On lattices, the late-time operator shape generally retains discrete-symmetry anisotropy rather than becoming spherical; its boundary is determined by a model-specific velocity function vB(φ).The asymptotic shape follows a geometric construction from the direction-dependent normal growth speed.Source: E. Shape of the operator at late times · p. 11
• For q > qc ≲ 2, the lattice velocity becomes nonanalytic near axis directions, producing flat facets; at q = ∞, the operator shape is a square.The deterministic-limit velocity is vq=∞(φ) = 2(|cos φ| + |sin φ|).Source: E. Shape of the operator at late times · p. 12
• At q = 2, simulations show a growing anisotropic 2D cluster, with the axis radius approaching the maximal value R(θ = 0) → 2t.The maximal-speed behavior depends on the specific discrete spacetime geometry.Source: Figure 11 · p. 12, E. Shape of the operator at late times · p. 12

F. Formal viewpoint

• After coarse-graining, the cluster density ρ(C) becomes a probability distribution over closed interfaces, so KPZ growth gives its corresponding Fokker–Planck dynamics.This identifies the formal cluster description with the stochastic evolution of a growing interface.Source: F. Formal viewpoint · pp. 12–13

A. Reduction to Ising spins

• Unitarity makes this effective Ising model simpler than analogous non-unitary tensor-network models, despite initially allowing complicated and even negative weights.The resulting local weights include three-body interactions after integrating out bra variables.Source: §IV · p. 13, A. Reduction to Ising spins · p. 15, Figure 14 · p. 15

IV. EXACT CALCULATION OF OTOC IN ‘SPACETIME’ PICTURE

• The spacetime construction is proposed as a route toward higher OTOC moments, commutator powers, and direct von Neumann entropy calculations.The paper also uses the mapping for an entanglement-related quantity.Source: §IV · p. 13
• The exact finite-time OTOC expression has no approximations, while replacing its discrete function g by a Gaussian cumulative distribution recovers the large-time formulas.The Gaussian approximation is valid when t is large.Source: §IV · p. 13
• The OTOC calculation restricts relevant circuit unitaries to the rectangle formed by the future light cone of X0 and past light cone of Yx.Null coordinates u = (t + x + 1)/2 and v = (t − x + 1)/2 parameterize this intersection.Source: §IV · p. 13

A. Reduction to Ising spins

• Haar averaging converts the four-unitary expression into an Ising-spin partition function, with interactions arising within each unitary and between neighboring unitaries.The spins take values ↑ and ↓, and the q-dependence enters through bond interactions.Source: A. Reduction to Ising spins · pp. 13–15
• For the elementary light-cone intersection ℓu = ℓv = 1, the calculation gives F = −1/(q^4 − 1).The tensor contraction uses the ↑ and ↓ vectors associated with Haar averaging.Source: A. Reduction to Ising spins · p. 14, Figure 12 · p. 15

B. Partition function for two directed paths

• Fixing the observable-attached spins imposes two bulk ↓ variables and forces two directed domain walls from the bottom to the top boundary.The first local weight rule makes each wall move only down-left or down-right.Source: B. Partition function for two directed… · p. 16
• Domain-wall deformations preserve the bulk weight, so the partition function depends on the walls’ endpoint positions and their boundary interaction weights.The right endpoint is (t − v, v), while the left endpoint is (u, t − u) in null coordinates.Source: Figure 15 · p. 16, B. Partition function for two directed… · p. 16
• The nonintersection constraint is evaluated by reflection: unrestricted paths A→B and C→D are corrected by subtracting paths A→D and C→B.This converts the constrained count into a difference of products of ordinary random-walk path counts.Source: B. Partition function for two directed… · pp. 16–17
• Combining domain-wall weights and path counts reproduces the ℓu = ℓv = 1 result q^4/(q^4 − 1) and yields the general exact expression.The spacetime variables satisfy t = ℓu + ℓv − 1 and x = ℓu − ℓv.Source: B. Partition function for two directed… · p. 17

C. Bounds on Fluctuations

• Exact fluctuation calculations are unwieldy because they involve high powers of the circuit unitaries, so the analysis uses bounds in two spatial regimes.The fluctuation argument concerns randomness between circuit realizations.Source: C. Bounds on Fluctuations · p. 17
• Far outside the butterfly light cone, the variance bound is meaningful and indicates negligible operator leakage beyond the vB-defined cone.The bound applies generally but is useful when |x| ≫ vBt.Source: C. Bounds on Fluctuations · p. 17
• Deep inside the light cone, the averaged OTOC differs from one by an exponentially small amount, and fluctuation bounds likewise suppress realization dependence.The discrepancy is bounded by O(1) exp(−(vBt − |x|)^2/2σ^2).Source: C. Bounds on Fluctuations · p. 17
• The fluctuation estimate uses approximate unitary-design results and Haar-random estimates, with a modified ring geometry introduced for the analysis.The ring has 2ct spins, and qualitative conclusions for c < 1 are anticipated rather than rigorously established.Source: C. Bounds on Fluctuations · p. 17

V. ENTANGLEMENT GROWTH

• The averaged purity of half an infinite chain is calculated exactly for arbitrary time under the circuit protocol with an initial product state.The calculation parallels the OTOC calculation, differing only in boundary conditions.Source: §V · p. 18
• Unitarity reduces the purity partition function to bulk domain-wall weights, with a single domain wall whose left- or right-down moves can be summed explicitly.The boundary factors cancel, leaving the domain-wall contribution as the relevant partition-function weight.Source: §V · p. 18
• The purity speed vP bounds the growth rate of the second Renyi entropy, while vE for the von Neumann entropy is also bounded by vP.The latter follows from SvN ≥ S2; the bound does not fix the actual S2 growth rate.Source: §V · p. 18
• Rare anomalously small values of S2 could affect q^−S2 despite small typical fluctuations, limiting what the purity calculation establishes about entropy growth.A direct finite-q calculation of S2 or SvN remains more difficult because of the logarithm.Source: §V · p. 18
• The scaling limit produces a deterministic domain-wall geometry, supporting a minimal-surface picture for entanglement growth in higher dimensions.In this picture, −log P is the energy of a coarse-grained minimal surface determined by elastic minimization.Source: §V · pp. 18–19
• Staircase circuits can make vE/vB arbitrarily small, with vE/vB ≲ 4/q2 in the large-q regime.This construction uses large but finite interaction range and does not require fine tuning.Source: Figure 16 · p. 19

A. Nonuniversality of the ratio vE/vB

• The construction shows that vE/vB need not be close to 1/2, extending earlier models with vE/vB < 1.The paper contrasts this with a deterministic disordered chain where vE/vB can be zero, while retaining translation-invariant circuit statistics here.Source: A. Nonuniversality of the ratio vE/vB · p. 19
• Left- and right-oriented staircases are applied at random Poissonian times and locations, with large finite range satisfying ℓ/q2 ≫ 1.The construction preserves statistical translational invariance.Source: A. Nonuniversality of the ratio vE/vB · p. 19
• For staircase circuits, vE/vB ≲ 4/q2, so the ratio can be made arbitrarily small by increasing q and the staircase range.The estimate follows from vE ≤ 2rℓ and vB ≃ q2rℓ/2.Source: A. Nonuniversality of the ratio vE/vB · p. 19

VI. OUTLOOK

• The paper conjectures that its random-circuit scaling forms apply to asymptotically late-time OTOCs in generic nonintegrable systems and quantum field theories.This remains a conjecture requiring tests in other calculable settings.Source: §VI · p. 19
• Mean-field treatment would produce a local traveling-wave equation, but it is not valid in physical dimensionalities.This differs from the random-circuit result, where higher-dimensional OTOCs are not governed by a local differential equation.Source: §VI · p. 20
• The random-circuit results do not provide a nontrivial bound on the operator norm, and extensions incorporating conserved quantities remain open.The Frobenius-norm OTOC and operator norm differ by a large dimension factor.Source: §VI · p. 20

Appendix A: Evolving Distribution on Operator Strings

• Haar-averaged squared operator coefficients evolve linearly through a real symmetric transition matrix, defining a fictitious classical Markov process.This process simplifies the true unitary operator dynamics without being identical to it.Source: Appendix A · p. 21
• Because transition probabilities depend only on operator support, the dynamics reduce to binary occupation variables, while the operator endpoint evolves autonomously.The endpoint update probabilities are independent of whether the neighboring site behind it is occupied.Source: Appendix A · p. 21
• For a two-site update, the endpoint remains at x with probability p = 1/(q2 + 1) and moves to x + 1 with probability 1 − p.This establishes the single-step evolution used for the one-dimensional diffusion picture.Source: Appendix A · p. 21

Appendix B: Velocity and diffusion constant for lattice diffusion equation

• Alternating-bond dynamics make the endpoint position equivalent to counting right and left moves, yielding its lattice probability distribution and Gaussian cumulative form.The endpoint coordinate is determined by u − v, with u + v = t.Source: Appendix B · pp. 21–22

Appendix C: Noisy diffusion equation

• Noise shifts the wavepacket center by order t1/4, smaller than its diffusive width, so density fluctuations are subleading at late times.The scaling is self-consistent because ρ − ρ0 is smaller than ρ near the peak.Source: Appendix C · p. 22

Appendix D: Random Clifford operators

• For prime q, the Clifford group is a unitary 2-design, so its second-moment average matches the Haar average.The proof uses the Pauli basis, symplectic transformations, and invariance of the Clifford and Haar averages.Source: Appendix D · pp. 22–23
• The Pauli group is generated by X and Z operators, and the Clifford group is its normalizer acting through symplectic transformations on the quotient Pauli labels.This supplies the algebraic structure used in the 2-design argument.Source: Appendix D · p. 22
• The Haar-averaged image is constrained to the span of the identity and swap operators, enabling the comparison with Clifford averaging.The identity and swap commute with U ⊗ U, which fixes the relevant invariant subspace.Source: Appendix D · p. 23
• When q is not prime, no probability distribution over the Clifford group forms a unitary 2-design.For q = 6, Pauli operators with different orders produce incompatible Clifford averages.Source: Appendix D · p. 23

Appendix E: Anomalous behaviour of the front for φ = 0

• For p > pc, the lattice-aligned front travels at the causal speed v(φ = 0) = 2 and remains unrough, with width O(1).At p = pc, the aligned front is logarithmically rough; for the lattice model, pc ≲ 2.Source: Appendix E · p. 23
• The aligned-front transition maps to directed percolation: finite light-front occupation indicates an attached, narrow front, while extinction signals detachment and KPZ roughening.The occupation dynamics form a birth-death process with active and extinct regimes.Source: Appendix E · pp. 23–24

Appendix F: Shape of a spreading droplet for weakly varying v(φ)

• For weakly varying velocity functions, the droplet-shape equations yield periodic radial profiles because the expansion reduces to total derivatives of periodic functions.The argument extends formally to general v(φ) when the solution remains smooth during the interpolation.Source: Appendix F · p. 24

Appendix G: Haar average formula

• The Haar average decomposes the two-copy representation into symmetric and antisymmetric subspaces, where Schur’s lemma fixes each block’s proportionality constant.The corresponding projectors are (I ± F)/2, with F the swap operator.Source: Appendix G · p. 24

Appendix H: Proof of Eq. (82)

• For local Haar layers, exponentially small Haar expectations remain exponentially small after sufficiently many convolutions, using the Brandão–Harrow–Horodecki result and supporting lemmas.The cited theorem does not directly cover the setting, but related lemmas suffice for the proof.Source: Appendix H · p. 25

Appendix I: A mean field approximation

• The mean-field traveling wave has a finite butterfly speed and does not broaden, unlike the random-circuit growth picture with roughening fronts.A seeded correlator grows toward one on a timescale of order λ^-1 before saturating.Source: Appendix I · p. 25
• Because the approximation is uncontrolled, any traveling-wave regime it produces is not valid for the studied systems at short or large times.A small parameter could justify mean field over finite long times in variant models, but asymptotic roughening is still expected.Source: Appendix I · p. 25
• The mean-field approximation factorizes the joint occupation distribution, replacing crucial intersite correlations with independent site probabilities.The authors state that this approximation does not capture the true behavior at short or asymptotically long times.Source: Appendix I · p. 25
• The resulting lattice equation has death, spreading, and overcounting terms and becomes a traveling-wave equation after a formal second-order lattice-spacing expansion.The continuous-time Poissonian protocol is introduced to simplify the equations without changing the basic point.Source: Appendix I · pp. 25–26