Diffusive Operator Spreading for Random Unitary Free Fermion Circuits
Beatriz C. Dias, Masudul Haque, Pedro Ribeiro, Paul McClarty
DDiffusive Operator Spreading for Random Unitary Free Fermion Circuits
Beatriz Dias, Masudul Haque,
2, 3
Pedro Ribeiro,
1, 4 and Paul A. McClarty CeFEMA, Instituto Superior T´ecnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal Department of Theoretical Physics, Maynooth University, Co. Kildare, Ireland Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, 01187 Dresden, Germany Beijing Computational Science Research Center, Beijing 100084, China
We study a model of free fermions on a chain with dynamics generated by random unitary gatesacting on nearest neighbor bonds and present an exact calculation of time-ordered and out-of-time-ordered correlators. We consider three distinct cases: the random circuit with spatio-temporaldisorder (i) with and (ii) without particle number conservation and (iii) the particle non-conservingcase with purely temporal disorder. In all three cases, temporal disorder causes diffusive operatorspreading and ∼ √ t entanglement growth. This is in sharp contrast to Anderson localization for thecase of static disorder and with the ballistic behavior observed in both the clean case of Hamiltonianevolution and in fully random unitary quantum circuits. Introduction − Understanding the nature of quantumdynamics in many-body interacting systems far fromequilibrium is one of the central issues in physics. Re-moved from the sanctuary of linear response theory, onemust often resort to numerical simulations but the expo-nential scaling of Hilbert space dimension typically limitsthese to small system sizes. Even so, huge progress hasbeen made over recent years driven partly by experimentson quantum simulators, numerical developments and oc-casional exact calculations.In the spirit of looking for universal features in many-body dynamics, random unitary circuits have been in-tensively studied in the last few years (see, for example,[1–7]). These models describe space and time processesthrough the application of local random unitary gates tosome underlying degrees of freedom. Since there is noHamiltonian dynamics, even energy conservation is sac-rificed in order to uncover generic features of local dy-namics. Since traditional correlators are randomized ateach discrete time step, the strength of these models is incapturing the spread of quantum information includingthe initial linear growth of bipartite entanglement froman initial product state to a fully random state with Pageentanglement. One observable that is particularly suitedto such models is the degree to which spatially separatedlocal operators commute after time evolution: C ( r, t ) ≡
12 Tr h ρ [ O ( t ) , O r ] † [ O ( t ) , O r ] i (1)where O r is an operator localized at position r [8–19].Expanding out the correlator gives both time orderedcorrelators (TOC) and out-of-time-ordered correlators(OTOC). The quantity C ( r, t ) is known to exhibit bal-listic spreading and KPZ growth at the light cone inter-face in even the simplest variant of random circuit mod-els [1, 3] and the dynamics can be mapped to a biasedrandom walk. In the presence of a conserved charge, thepicture is modified owing to the diffusion of the conservedcharges [4, 5]. Thus, one finds there is a ballistic frontthat itself spreads diffusively. As with everything else in random matrix theory, the importance attached to theOTOC in random unitary circuits is that it is conjec-tured to provide a tractable instance of universal physicsin this case tied to thermalization and the scramblingof quantum information in general non-integrable many-body interacting systems [20–24]. − . . .
06 00 . − . . . − σ ( t )2 σ ( t ) rt C ( r , t ) r − µ ) /σ ( t )NC-T C ( r , t ) · σ ( t ) / A C-STNC-ST C ( r , t ) · σ ( t ) / A t = 450100 FIG. 1. (a) Circuit scheme in the r − t plane: the random two-site unitaries u r,r +1 are applied at neighboring sites in a brickwall pattern. (b) Exact C ( r, t ) = 1 / [ C ( r, t ) − C ( r, t )] as anheight and color map for a system with 100 sites in the case ofNC-ST (and NC-T) evolution. The black curves envelop the σ ( t ) = √ t and 2 σ ( t ) = √ t regions for C , which is the lead-ing order term in C ( r, t ). (c) Rescaled TOC (left) and OTOC(right), i.e. C ( r, t ) σ ( t ) /A and C ( r, t ) σ ( t ) /A ( t ) (with A thenormalization), as a function of ( r − µ ) /σ ( t ) for three fixedtime steps t ∈ { , , } and for NC-ST (top), CS-T (mid-dle) and NC-T (bottom). Full lines correspond to exact cal-culations (shown in (a) for NC-ST) and the color filled regionto simulations which include an average of 4000 disorder re-alizations, starting from a random state with particle numberfixed to L/
2. Both collapse to the continuum limit solution,i.e. the Gaussian g ( r, t ) = 1 / ( √ π ) exp[ − ( r − µ ) / (2 σ ( t ))].For C ( r, t ), A = 2 and σ ( t ) = √ t for all cases; for C ( r, t ), σ ( t ) = √ t for all cases, A ( t ) = 1 / (2 √ πt ) for NC-ST andNC-T and A ( t ) = 2 / ( √ πt ) for C-ST. The deviation µ = 1 / C and C around r = 0. a r X i v : . [ c ond - m a t . s t r- e l ] F e b In this paper we consider driven free fermion dynam-ics within the context of random unitary circuits. Freefermion circuits first appeared as classically simulatablematchgate circuits [25], which were later shown to cor-respond to a model of free fermions in 1D [26]. Unliketheir many-body interacting analogues, there is evidencethat free fermion models with spatio-temporal noise ex-hibit diffusive dynamical features that were observed inthe growth of the von Neumann entropy [2, 27], in large-deviation statistics [28] and in the magnetization dynam-ics of the transverse field Ising model that maps to a freefermion problem. The presence of diffusive dynamics insuch cases is a nontrivial result that sits in contrast tothe Anderson localization expected for static spatial dis-order or to the ballistic behavior expected in the cleancase. The extent to which interactions change this pic-ture is not yet fully understood. Nonetheless, in generic(i.e. interacting) quantum random circuits both OTOCsand the entanglement dynamics spread ballistically.In the following, we study the effect of spatio-temporalnoise on the dynamics of free fermions in a setting wherean exact calculation of the C ( r, t ) for a random circuitis possible. In this paper, we establish through an ex-act calculation the diffusive behavior from the underlyingdynamics of C ( r, t ) for quadratic fermions. We considerthree distinct instances of free fermion evolution: a par-ticle conserving spatio-temporal random circuit (C-ST),its generalization to a non particle conserving process(NC-ST), and a spatially homogeneous case where ran-domness appears only in the time direction (NC-T). Inthe Supplementary Section, we also present numerical re-sults for a non particle conserving circuit with quenchedspatial disorder (NC-S) that Anderson localizes.The random circuit acts on a chain of L sites( r = 0 , ..., L −
1) and periodic boundary conditions asa discrete time protocol. A time step corresponds tothe action of the unitary operation U = U even U odd , ob-tained by the successive application of the half-steps U even = Q L/ − R =0 U R, R +1 and U odd = Q L/ − R =0 U R − , R .Here, U r,r +1 is a random unitary acting nontriviallyon the local 4-dimensional Hilbert space of sites r and r + 1. The evolution follows the brick wall patternshown of Fig. 1(a). For the free fermion evolution, in-stead of sampling U r,r +1 over U(4) [2], we consider asub-group that leaves invariant the algebra of single-particle fermionic creation and annihilation operators ( a † r and a r ) [25, 26], i.e. U † r,r +1 Ψ r U r,r +1 = u r,r +1 Ψ r , whereΨ r ≡ ( a r , a r +1 , a † r , a † r +1 ) T and u r,r +1 is a 4 × τ u Tr,r +1 τ = u † r,r +1 , with τ a Pauli matrix acting on Nambu space.When det u r,r +1 = 1, such unitary operator can be ob-tained by continuous time evolution under a quadratic(free) fermionic Hamiltonian during a period of time∆ t , i.e. U r,r +1 = exp( − i Ψ † r H r Ψ r ∆ t ), in which case u r,r +1 = exp( − iH r ∆ t ). Here we take u r,r +1 to be Haar- . . .
002 0 . .
01 0 . NC-ST (filled)C-ST (dashed) 0.210.002 0.5 S / S ∞ t/L L = 20 ∝ √ t FIG. 2. Growth of the entanglement entropy divided by thesaturation value
S/S ∞ as a function of t/L for random freefermions − NC-ST (filled), C-ST (dashed) and NC-T (inset) − with L = 100 and a subsystem of 50 sites showing S ( t ) ∝ √ t at the earliest times. The data results from averaging over1000 realizations up to t = 5000, starting from a randomproduct state with particle number fixed to L/
2. The curvesfor NC-ST and NC-T overlap. distributed [29]. Note that for NC-T at each time stepthe same is applied to all pairs of sites, i.e. u r,r +1 = u ,and for C-ST u r,r +1 = v r,r +1 ⊕ v ∗ r,r +1 , with v r,r +1 a 2 × Dynamics of Entanglement − For completeness, westart by analysing the dynamics of the von Neumannentropy, S = − tr(˜ ρ ln ˜ ρ ), with ˜ ρ the reduced densitymatrix of a subsystem of size L/
2, starting from aninitial product state with a well defined particle num-ber. Fig. 2 shows that for the three processes the en-tanglement grows as ∼ D S √ t (with D S a time inde-pendent constant) for small times saturating at times t sat ∼ L . Asymptotically diffusive growth of the R´enyientropy S has been discussed in the literature, coincid-ing with linear growth of the von Neumann entropy atleast in low dimensional local Hilbert spaces in modelswith conserved quantities [30, 31]. In these cases, op-erator growth has a ballistically propagating front withdiffusive growth in the rest frame of the front. The ran-dom circuit free fermion model is therefore qualitativelydifferent to these cases. For t (cid:29) t sat , the saturation value S ∞ = s L + s + O (1 /L ), coincides with the mean en-tanglement entropy of a random Gaussian state [32, 33],with s ’ .
193 for all the considered free fermion pro-cesses and s ’ .
085 for the NC-ST case, well below thePage value ( s = ln 2 / ’ . , s = − /
2) obtained byaveraging over the full Hilbert space. These results showthat the rate of increase of the entanglement is compat-ible with diffusion of quantum information. In addition,while the saturation entanglement has volume law scal-ing, free fermion dynamics cannot explore the full Hilbertspace. We now turn to the signatures of free fermion dy-namics in the spatial correlations.
Operator Spreading − The time-ordered density-density correlator becomes trivial when averaged overtemporal disorder (see Supplementary Section [29]) but C ( r, t ), introduced in Eq. (1), remains non-trivial uponaveraging. We consider an average over separable initialstates, which is equivalent to taking ρ ∝ O r = 1 / † O r Ψwhere Ψ ≡ ( a , . . . , a L , a † , . . . , a † L ) T and O r is a localsingle-particle operator. The computation of C ( r, t ) forfree fermions - in common with other correlators - canbe brought into a form where the trace need only be per-formed over 2 L × L matrices rather than over the entireHilbert space. One can show that the many-body corre-lator can be written in terms of single-body quantities as C ( r, t ) = 1 / [ C ( r, t ) − C ( r, t )] where C ( r, t ) = tr (cid:2) O ( t ) O r (cid:3) , (2) C ( r, t ) = tr [ O ( t ) O r O ( t ) O r ] , (3)and O ( t ) = U O ( t − U † . The general relation betweenthe single particle and many-body TOC and OTOC isgiven in the Supplementary Material [29]. TOC and OTOC Numerics − The TOC ( C ) andOTOC ( C ) may now be simulated efficiently by gen-erating local pseudo-random gates u r,r +1 and averag-ing the result over many realizations of the time evo-lution. As before, we average over initial states thatare product states of definite particle number. Resultsfor L = 100 and 4000 disorder realizations are shownin Fig. 1(c) for a symmetrized particle number operator O r = 1 / a † r a r − a r a † r ) = a † r a r − /
2. In contrast to theballistic spreading seen in Haar-random circuits, C ( r, t )for random free fermions diffuses. This can be seen mostclearly in Fig. 1(c) where C σ ( t ) /A and C σ ( t ) /A ( t ) areshown to collapse to a Gaussian with standard devia-tion growing respectively as √ t and √ t at multiple fixedtimes, with the horizontal axis rescaled to ( r − µ ) /σ . Ex-amining the two terms C and C separately reveals thatboth the ordinary dynamical two-time correlator C andthe true out-of-time-ordered component C are diffusive.The magnitude of C is, however, much larger than thatof C at fixed time. These results are compatible withthe observed √ t early time entanglement growth. Exact Calculation of C ( r, t ) − We now look for the ori-gin of the diffusive behavior of C ( r, t ) examining C ( r, t )and C ( r, t ) in turn and proceeding analytically by com-puting the exact averages over the random unitaries forthe NC-ST case referring to the Supplementary Sec-tion [29] for further details. In the single particle pic-ture, we may denote states of the system by | α i ≡| r α , s α i ≡ | R r α + b α , s α i , where r α = 0 , . . . , L − s α ∈ { p, h } labels the particle-hole index; R r α = 0 , . . . , L/ − r α , r α + 1) and b r α ∈ { , } such that r α = 2 R r α + b r α . Note thatthe sites are acted upon by modulo L , due to periodicboundary conditions.Starting with C ( r, t ), we use a notation in whichoperators are rendered as state vectors C ( r, t ) =tr (cid:2) O ( t ) O r (cid:3) = hh O r || O ( t ) ii , where || O r ii ≡ || O r (0) ii , || O r ( t ) ii = P αβ || αβ ii (cid:10) α (cid:12)(cid:12) O r ( t ) (cid:12)(cid:12) β (cid:11) and || αβ ii = | α i ⊗h β | T . Since O r is the symmetrized number operator, weget || O r ii = X s || r, s ; r, s ii . (4)For a single realization, O ( t + 1) = U † O ( t ) U , whichtranslates to || O ( t + 1) ii = U † ⊗ U T || O ( t ) ii . Apply-ing one layer of our circuit, i.e. U = U even U odd , and av-eraging over multiple realizations of the random circuitas summarized in the Supplementary Section one finds || O ( t + 1) ii = ˆ W || O ( t ) ii ˆ W ≡ L/ − X R =0 || φ R ii hh φ R || L/ − X R =0 || φ R − ii hh φ R − || , (5)having introduced || φ r ii ≡ X x ∈{ r,r +1 } X s || x, s ; x, s ii . (6)With the initial condition || O (0) ii = P s || , s ; 0 , s ii , thecomplete time evolution of C is determined. The sub-space spanned by || φ r ii is closed under the evolution thusrespecting the particle-hole symmetry. These vectorizedoperators have the property hh φ r k φ r ii = δ r ,r + δ r ,r ± that leads to the simplification of the recursion relation,Eq. (5), to D h φ r (cid:13)(cid:13)(cid:13) O ( t + 1) EE = 14 (cid:16) hh φ r − || + 2 hh φ r || + hh φ r +2 || (cid:17) || O ( t ) ii . (7)From this and the initial condition we obtain C (2 R, t + 1) = C (2 R + 1 , t + 1)14 (cid:0) C (2 R − , t ) + 2 C (2 R, t ) + C (2 R + 2 , t ) (cid:17) (8)for t ≥ ( δ R − , + δ R, ) for t = 0.Taking the continuum limit, lim t,L →∞ C ( r, t ) =lim a → aC ( x = ra, τ = ta ), leads to the 1D diffusionequation ∂ τ C ( x, τ ) = ∂ x C ( x, τ ) with diffusion constant D = 1 which approximates the exact discrete time evo-lution very well even for relatively small times. Exact Calculation of C ( r, t ) − The calculation of C ( r, t ) proceeds analogously to that of C ( r, t ) but ismore involved not least because the disorder average iscarried out over a product of four unitaries rather thantwo in the case of C . As before, we write the correlatorin a vectorized notation C ( r, t ) = tr [ O ( t ) O r O ( t ) O r ] = hh Q r || S || Q ( t ) ii , (9)with Q r ≡ O r ⊗ O r fixed by the choice of observable tobe || Q r ii = X s || rs, rs, rs, rs ii − || rs, r ¯ s, rs, r ¯ s ii . (10)with s = p, h and the corresponding ¯ s = h, p and || Q r ( t ) ii = X αβµν || αβµν ii h αβ | Q r ( t ) | µν i , (11) S = X αβµν || αβµν ii hh αβνµ || , (12)where || αβµν ii = | αβ i ⊗ h µν | T . The || Q ( t ) ii stateevolves as || Q ( t + 1) ii = U † ⊗ U † ⊗ U T ⊗ U T || Q ( t ) ii (13)and an average is taken over different realizations even-tually leading to the recursion relation D(cid:10)
Θ ¯Θ r,r (cid:13)(cid:13)(cid:13) Q ( t + 1) EE = ˆ W || Q ( t ) ii , (14)ˆ W ≡ tr [ M r,r Υ r,r ] , Υ r,r ≡ hh Θ ¯Θ r − ,r − || hh Θ ¯Θ r − ,r || hh Θ ¯Θ r − ,r +2 ||hh Θ ¯Θ r,r − || hh Θ ¯Θ r,r || hh Θ ¯Θ r,r +2 ||hh Θ ¯Θ r +2 ,r − || hh Θ ¯Θ r +2 ,r || hh Θ ¯Θ r +2 ,r +2 || , analogous to Eq. (7). In the above expression, M r,r are3 × M r,r for r = r, r ± , r ± r . The || Θ ¯Θ r,r ii = || Θ r,r ii − || ¯Θ r,r ii , (15)and || Θ r,r ii = 1 √ g r,r X r α ∈{ r,r +1 } r β ∈{ r ,r +1 } X s α ,s β || αββα ii , (16) || ¯Θ r,r ii = 1 √ g r,r X r α ∈{ r,r +1 } r β ∈{ r ,r +1 } X s α ,s β || α ¯ αβ ¯ β ii , (17)with g r,r = N ( N − δ r,r ). This completely determines theevolution of the single particle OTOC. With the initialcondition, Eq. (10), one finds C ( r, t = 0) = 2 δ r, , (18) C ( r, t = 1) = 118 ( δ r, − + δ r, ) , (19) and, for subsequent times, we define K r,r ( t ) ≡ √ D(cid:10)
Θ ¯Θ r,r (cid:13)(cid:13)(cid:13) Q ( t ) EE , (20)and use Eq. (14) to get C (2 R, t + 1) = C (2 R + 1 , t + 1)= tr [ M R, R K R, R ( t )] with K R, R ( t ) ≡ K R − , R − K R − , R K R − , R +2 K R, R − K R, R K R, R +2 K R +2 , R − K R +2 , R K R +2 , R +2 . (21)The evolution thus described exactly reproduces the nu-merical results described above. We may now take thecontinuum limit from Eq. (21). This highlights one im-portant distinction between the C and C cases: theevolution of C depends on the matrix K whose ele-ments are indexed by a pair of spatial coordinates. In thecontinuum limit C ( r, t ) = lim a → C ( x = ra, τ = ta ) ’ f a K x,x ( τ ), where K x,x ( τ ) can be shown to obey the 2Ddiffusion equation ∂ τ K x,x ( τ ) = (cid:0) ∂ x + ∂ x (cid:1) K x,x ( τ ) withinitial condition K x,x (0) = 2 δ ( x ) δ ( x ) and f = 1 / t →∞ C ( r, t ) /C ( r, t ) = 0, i.e.for large times the C ( r, t ) dictates the leading behaviorof the OTOC. Extensions to C-ST, NC-T and NC-S − We haveshown that both C and C spread diffusively for freefermions in 1D in the presence of spatio-temporal noise(NC-ST). We now consider exact calculations for two fur-ther cases: C-ST where the fermion particle number isconserved and each gate in the quantum circuit is chosenrandomly, and NC-T where the unitary evolution is spa-tially homogeneous but where there is temporal noise − a single gate is chosen randomly at each time step andapplied to all pairs of sites. Fig. 2 shows that the vonNeumann entropy grows like √ t for all three cases: NC-ST, C-ST, NC-T. The Supplementary Section lays out indetail exact calculations of C and C , analogous to thecalculation summarized above for NC-ST, but with thecontinuum limit of C for C-ST having a different nor-malization due to f = 2. The result is that there is diffu-sive spreading in all three cases with diffusion constantscoinciding with those found for NC-ST. In contrast, nu-merical results obtained for the temporal homogeneouscase, NC-S, presented in [29], where even and odd lay-ers of random gates are fixed and applied repeatedly intime, show that C and C remain Anderson localized,decaying exponentially around r = 0 [34]. Conclusions − There is evidence that Hamiltonianmodels of quadratic fermions in one dimension exhibitdiffusive spreading of correlations when subjected tonoise [28, 35]. Here we have found an analytically solv-able instance of this physics in the OTOC of a randomcircuit model with both number conserving and non-conserving quadratic fermionic terms. In the long timelimit, the states that result from this dynamics are ex-tended with volume law entanglement but depart sig-nificantly from random matrix eigenstates. In contrast,previously studied random unitary circuit models exhibitballistic spreading of correlations with entanglement ap-proaching the Page value asymptotically. All at once,this strongly suggests that, in Hamiltonian models of freefermions, Anderson localization is destroyed by the co-herent noise we have considered and that ballistic propa-gation, expected for spatially homogeneous systems, be-comes diffusive in the presence of noise. Both implica-tions lay bare unusual features of quadratic fermion sys-tems.BD acknowledges support by FCT through Grant No.UIDB/04540/2020. BD and PR acknowledge support byFCT through Grant No. UID/CTM/04540/2019. [1] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X , 021014(2018).[2] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys.Rev. X , 031016 (2017).[3] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, andS. L. Sondhi, Phys. Rev. X , 021013 (2018).[4] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev.X , 031057 (2018).[5] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk,Phys. Rev. X , 031058 (2018).[6] T. Rakovszky, C. W. von Keyserlingk, and F. Pollmann,Phys. Rev. B , 125139 (2019).[7] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. X , 041019 (2018).[8] A. I. Larkin and Y. N. Ovchinnikov, Soviet Journal ofExperimental and Theoretical Physics , 1200 (1969).[9] B. Swingle, Nature Physics , 988 (2018).[10] K. Hashimoto, K. Murata, and R. Yoshii, Journal ofHigh Energy Physics , 138 (2017), arXiv:1703.09435[hep-th].[11] J. Maldacena, S. H. Shenker, and D. Stanford, Journalof High Energy Physics , 1 (2016).[12] A. Kitaev, in KITP Talks in Entanglement in Strongly-Correlated Quantum Matter (2015).[13] P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, Jour-nal of High Energy Physics , 4 (2016).[14] R. Fan, P. Zhang, H. Shen, and H. Zhai, Science Bulletin , 707–711 (2017).[15] V. Khemani, D. A. Huse, and A. Nahum, Phys. Rev. B , 144304 (2018).[16] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hay-den, Phys. Rev. A , 040302 (2016).[17] M. G¨arttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall,J. J. Bollinger, and A. M. Rey, Nature Physics , 781(2017).[18] J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng,and J. Du, Phys. Rev. X , 031011 (2017).[19] X. Mi, P. Roushan, C. Quintana, S. Mandra, J. Marshall,C. Neill, F. Arute, K. Arya, J. Atalaya, R. Babbush, et al. , arXiv preprint arXiv:2101.08870 (2021).[20] J. Riddell and E. S. Sørensen, Phys. Rev. B , 024202 (2020).[21] U. Agrawal, S. Gopalakrishnan, and R. Vasseur, Phys.Rev. B , 174203 (2019).[22] C.-J. Lin and O. I. Motrunich, Phys. Rev. B , 134305(2018).[23] J. Riddell and E. S. Sørensen, Phys. Rev. B , 054205(2019).[24] L. Colmenarez and D. J. Luitz, arXiv e-prints (2020),arXiv:2005.10257 [cond-mat.str-el].[25] L. G. Valiant, SIAM Journal on Computing , 1229(2002).[26] B. M. Terhal and D. P. DiVincenzo, Physical Review A (2002), 10.1103/physreva.65.032325.[27] G. m. H. Ro´osz, R. Juh´asz, and F. Igl´oi, Phys. Rev. B , 134305 (2016).[28] D. Bernard and T. Jin, Phys. Rev. Lett. , 080601(2019).[29] See Supplemental Material for further details and addi-tional supporting data.[30] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk,Phys. Rev. Lett. , 250602 (2019).[31] M. ˇZnidariˇc, Communications Physics (2020),10.1038/s42005-020-0366-7.[32] J. M. Mag´an, Physical Review Letters (2016),10.1103/physrevlett.116.030401.[33] C. Liu, X. Chen, and L. Balents, Physical Review B (2018), 10.1103/physrevb.97.245126.[34] A. Lagendijk, B. Van Tiggelen, and D. S. Wiersma,Phys. Today , 24 (2009).[35] D. Bernard and T. Jin, arXiv preprint arXiv:2006.12222(2020). upplementary Material : Diffusive Operator Spreading for Random Unitary FreeFermion Circuits
Beatriz Dias, Masudul Haque,
2, 3
Pedro Ribeiro,
1, 4 and Paul A. McClarty CeFEMA, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal Department of Theoretical Physics, Maynooth University, Co. Kildare, Ireland Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany Beijing Computational Science Research Center, Beijing 100084, China
The following supplementary section contains complete details of the computation of the TOCand OTOC for the three types of random free fermion circuit NC-ST, C-ST and NC-T togetherwith further numerical results supporting those given in the main text.
CONTENTS
I Outline SM-2II Fermionic quadratic Hamiltonians SM-2III Setting up the random circuits SM-3IV TOC and OTOC in the single-body picture SM-3
V Average of products of Haar-distributed orthogonal and unitary matrices SM-6VI Two-point correlator SM-7VII NC-ST: averaged TOC, C ( x, t ) SM-7
S.1 Time evolution of || O ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-71 Average over disorder realizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-82 Dynamics of || O ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-8S.2 Prescription to obtain C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-9S.3 Continuum limit of C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-10 VIII NC-ST: averaged OTOC, C ( r, t ) SM-10
S.1 Time evolution of || Q ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-111 Average over disorder realizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-112 Dynamics of || Q ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-13S.2 Prescription to obtain C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-14S.3 Continuum limit of C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-15 IX NC-T: averaged TOC, C ( r, t ) SM-16X NC-T: averaged OTOC, C ( r, t ) SM-17XI C-ST: averaged TOC, C ( r, t ) SM-18XII C-ST: averaged OTOC, C ( r, t ) SM-18
S.1 Time evolution of || Q ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-181 Average over disorder realizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-192 Dynamics of || Q ( t ) ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-19S.2 Prescription to obtain C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-21S.3 Continuum limit of C ( r, t ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-21 XIII Continuum limit of C ( r, t ) SM-22XIV Numerics SM-22
S.1 NC-S: randomness in space alone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .SM-23
XV Conclusions SM-24References SM-25
M-2
I. OUTLINE
An analytical expression for the OTOC averaged over disorder realizations for evolution under generic unitarycircuits was obtained by Nahum et al. in Ref. [1]. We present analogous analytical results for three instances offree fermion (FF) evolution: a particle conserving spatio-temporal random circuit (C-ST), its generalization to aparticle non-conserving process (NC-ST) and a spatial homogeneous case where randomness appears only along thetime direction (NC-T). In the end, we obtain a diffusive picture for the spreading of both a time-ordered (TOC) andout-of-time-order correlators in the r − t plane, unlike the usual ballistic picture shown in Ref. [1].The calculation detailed in this supplementary section runs over many pages so, here, we give a summary of themain steps to guide the reader through the remainder of the section. The analytical calculation of the OTOC givenhere is made possible firstly because the evolution of quadratic fermions can be carried out within the 2 L dimensionalspace of single particle states as opposed to the full 2 L many-particle Hilbert space as outlined in Section II. Randomquadratic fermion unitary gates are applied to neighboring pairs of sites with a single time step being composedof layers of gates applied first to odd and then to even bonds (Section III). We study a correlator of commutators C ( r, t ) (Eq. (1) in the main text) acting on the many body states that is written in terms of a time-ordered correlationfunction (TOC), C ( r, t ), and an out-of-time order correlator (OTOC), C ( r, t ). For our problem these can be analysedin the local single particle basis. It turns out (Section IV) that the many-body correlator (Eq. S7) can be written interms of single particle correlators C ( r, t ) and C ( r, t ) which are similar respectively to C ( r, t ) and C ( r, t ). However,despite the similarity, the correspondence between C and C (and C and C ) is not one-to-one, with the singleparticle TOC and OTOC being a mixture of the many-body TOC and OTOC. To compute C , we take the localoperator at position r to be the symmetrized fermionic number operator.Observables for the random circuit must be disorder averaged. A summary of results on random matrix averages(needed for the average of free fermion gates) is given in Section V and further explained throughout the text. Onecan imagine building the random circuit for a given disorder realization, computing the desired correlator and thenaveraging over realizations. This computation can be done instead by studying first the action of a single layer(acting on even or odd pairs of sites) and then building up to multiple layers. Before proceeding to calculate higherorder correlators, we briefly show that the time-ordered density-density correlator becomes trivial when averaged overtemporal disorder (Section VI). Section VII presents the calculation of the time-ordered correlator C for the NC-STcase, where we begin by establishing a vectorized notation to be used throughout the text. The calculation of C involves an average over two unitaries (Section VII S.1) and results in the TOC being given by a discrete random walkin 1D (Section VII S.2). Taking the continuum limit of the discrete evolution equation for C results in a diffusionequation in 1D (Section VII S.3) and a Gaussian broadening with standard deviation σ ( t ) = √ t . In Section VIII, thecalculation of C proceeds analogously to that of C but is more involved not least because the disorder average iscarried out over a product of four unitaries rather than two in the case of C . The result is that C is the diagonal ofa two-dimensional quantity whose discrete evolution equation, detailed in Section VIII S.2, is approximated by a 2Ddiffusion equation in the continuum limit, where the OTOC is approximated by a Gaussian with broadening width σ ( t ) = √ t (Section VIII S.3). After establishing the results for the NC-ST case, we extend these to the NC-T andC-ST cases by pointing out the differences with regards to NC-ST. The different structure of the NC-T circuit playsno role and both C (Section IX) and C (Section X) are exactly equal to the ones obtained for NC-ST. The sameholds true for C in the C-ST case (Section XI). However, the conserving character of the unitaries used to build thecircuit leads C to be different from its non-conserving counterpart. The details given in Section XII do not alter thequalitative behaviour in the continuum limit, affecting only the normalization of the Gaussian obtained before.Lastly, the analytical results obtained for the three instances of free fermion evolution (NC-ST, NC-T and C-ST)are shown to agree with simulations and also with the continuum limit solutions (Section XIV). Furthermore, beforeconcluding (Section XV) we present numerical results for a non particle conserving circuit with quenched spatialdisorder (NC-S) that Anderson localizes (Section XIV S.1). II. FERMIONIC QUADRATIC HAMILTONIANS
Before proceeding to evaluate the TOC and OTOC, let us briefly introduce the system.We consider free (or quadratic) fermions on a chain with L sites. The Hamiltonian has, at most, quadratic termsin the fermionic operators, a and a † , and it can be written as H = 12 A † HA with H = (cid:18) h ∆∆ † − h T (cid:19) , (S1)where A = ( a , · · · , a L , a † , · · · , a † L ) T is the Nambu vector and h = h † and ∆ = − ∆ T . This system is invariant underM-3the exchange of creation and annihilation operators, i.e. it is particle-hole symmetric, which translates to SHS = − H T with S = (cid:18) (cid:19) . (S2)An operator H respecting (S2) is said to be particle-hole (PH) symmetric.Given the 2 L × L single-body Hamiltonian matrix, H , we may diagonalize the problem and compute arbitraryobservables in terms of the single particle states. More specifically, having the two-point correlator χ = h AA † i , themean value of any observable of the form O = A † OA can be computed using hO ( t ) i = 12 tr (cid:0) ρA † ( t ) OA ( t ) (cid:1) = −
12 tr (cid:0) Oχ ( t ) (cid:1) , (S3)with ρ = e − β H /Z and Z = tr(e − β H ) for a thermal density matrix (this is valid for our circuits with β = 0).Furthermore, we know that the fermionic operators A ( t ) obey ∂ t A ( t ) = − iHA ( t ), meaning that they evolve accordingto A ( t ) = exp( − iHt ) A . With this, and considering the initial χ = h AA † i = diag (1 − n , . . . , − n L , n , . . . , n L ), weknow how χ evolves: χ ( t ) = e − iHt χe iHt . (S4)This means that we have the tools to compute (S3). In particular, the entanglement entropy itself can be computedin a simple manner with [2] S = − tr (cid:0) ρ log ρ (cid:1) = − tr (cid:0) χ log χ (cid:1) . (S5) III. SETTING UP THE RANDOM CIRCUITS
Unlike the gates used for generic dynamics, which needed only to respect unitarity, free fermion evolution operatorsare restricted to obey the fermionic anticanonical commutation relations.In the single-particle basis, the Nambu vector A r = ( a r , a r +1 , a † r , a † r +1 ) T evolves as A r ( t + ∆ t ) = u r,r +1 A r ( t ),with u r,r +1 = exp( − iH ∆ t ). This operator respects PH symmetry, i.e. Su Tr,r +1 S = u † r,r +1 , and it can be written as u r,r +1 = V † ˜ u r,r +1 V , with ˜ u r,r +1 its representation in the Majorana basis, which respects ˜ u Tr,r +1 ˜ u r,r +1 = . Thus, the4 × u can be obtained as u = V † OV with V = 1 √ (cid:18) − i i (cid:19) (S6)and O a 4 × u r,r +1 = v r,r +1 ⊕ v ∗ r,r +1 (where ⊕ stands for direct sum), with v r,r +1 a 2 × IV. TOC AND OTOC IN THE SINGLE-BODY PICTURE
The degree of non-commutativity of two local observables W and V centred, respectively, around positions 0 and r with finite support is given by C ( r, t ) = 12 h| [ W ( t ) , V ] | i = C ( r, t ) − C ( r, t ) , (S7)where W = U † ( t ) WU ( t ) and the TOC and OTOC are, respectively, C ( r, t ) = hW ( t ) V i and C ( r, t ) = hW ( t ) VW ( t ) Vi .Quadratic many-body observables V = V (0) and W = W (0) can be written as V = 1 / A † V A and W = 1 / A † W A ,with A = ( a , · · · , a L , a † , · · · , a † L ) T the Nambu vector and a r and a † r the fermionic operators at position r . Consideringthat W and V have this form, we can express (S7) in terms of V and W , the so called single-body observables. Forthis, it is convenient to write the relation between many and single-body operators V = ∂ v (cid:12)(cid:12)(cid:12) v =0 exp (cid:16) v A † V A (cid:17) and W = ∂ w (cid:12)(cid:12)(cid:12) w =0 exp (cid:16) w A † W A (cid:17) (S8)M-4
FIG. S1: (a) Circuit scheme for the NC-ST and C-ST cases in the r − t plane: the two-site unitaries are appliedat neighbouring sites in a brickwall pattern. The two cases differ in the gates used: the NC-ST case uses two-sitefree fermion gates u r,r +1 = V † OV , where O is a 4 × u r,r +1 = diag( w, w ), with w a 2 × u is applied to all pairs of sites. (c) Circuit scheme for the NC-S case. Oneeven and one odd layers are generated as usual, these are fixed and repeated in time (i.e. the grey box is repeated intime). In (b) and (c) the gates are free fermion gates, as for NC-ST.This operator evolves as W ( t ) = U † ( t ) WU ( t ), with U ( t ) = e − i A † H t A . . . e − i A † H A . (S9)Replacing V and W by (S8) in (S7), C ( t ) and C ( t ) become C k ( r, t ) = ∂ w (cid:12)(cid:12)(cid:12) w =0 ∂ v =0 (cid:12)(cid:12)(cid:12) v =0 ∂ w (cid:12)(cid:12)(cid:12) w =0 ∂ v (cid:12)(cid:12)(cid:12) v =0 ˜ C k ( r, t ) with j = 1 , , (S10)˜ C ( r, t ) = 12 L tr h U † ( t )e w + w A † W A U ( t )e v + v A † W A i , (S11)˜ C ( r, t ) = 12 L tr h U † ( t )e w A † W A U ( t )e v A † V A U † ( t )e w A † W A U ( t )e v A † V A i (S12)These expressions can be simplified with the properties,e − A † B n A . . . e − A † B A = e − A † ˜ BA with e ˜ B = e B . . . e B n , (S13)tr h e − A † ˜ BA i = h det (cid:16) ˜ B (cid:17)i , (S14)det( A ) = exp(tr[log A ]) , (S15)with the second one being valid for ˜ B particle-hole symmetric. Applying these three expression to (S11) and (S12),we get ˜ C k ( r, t ) = 12 L tr h e − A † ˜ B k A i = 12 L h det(1 + e ˜ B k ) i = 12 L exp (cid:18)
12 tr h log(1 + e ˜ B k ) i(cid:19) , (S16)where k = 1 , ˜ B ≡ e ( v + v ) V u † t e ( w + w ) W u t , (S17)e ˜ B ≡ e v V u † t e w W u t e vV u † t e wW u t , (S18) u t ≡ e − iH t . . . e − iH . (S19)Finally, to obtain C ( r, t ) as a simple function of W ( t ) and V , we must perform the derivatives in (S10), making useM-5of (S16). After doing so, and using tr V = tr W = 0 (which greatly simplifies the result), we obtain C ( r, t ) = 12 (cid:26)
12 tr[ W ( t ) V ] + 14 tr[ W ( t )] tr[ V ] − tr[ W ( t ) V W ( t ) V ] (cid:27) , (S20) C ( r, t ) = 12 (cid:26)
12 tr[ W ( t ) V ] + 14 tr[ W ( t )] tr[ V ] + tr[ W ( t ) V W ( t ) V ] − W ( t ) V ] (cid:27) . (S21)Combining (S20) and (S21) in (S7) we obtain C ( r, t ) = 12 h C ( r, t ) − C ( r, t ) i , (S22)with the new TOC and OTOC being, respectively, C ( r, t ) = tr[ W ( t ) V ] , (S23) C ( r, t ) = tr[ W ( t ) V W ( t ) V ] . (S24)These expressions are very similar to C ( r, t ) and C ( r, t ), with the many-body observables being replaced by therespective single-body ones. Note, however, that the correspondence between the many and single-body TOC andOTOC, given in (S20) and (S21), is not one-to-one.The relations (S20) and (S21) can be simplified if we specify the observables to be the symmetrized particlenumber operator introduced below, i.e. W = O and V = O r . Since O r respects PH symmetry, one can show that1 / O ( t ) O r ] = tr[ O ( t ) O r O ( t ) O r ] and tr[ O ( t )] tr[ O r ] = 4 such that C ( r, t ) = 12 , (S25) C ( r, t ) = 12 [1 + 2 C ( r, t ) − C ( r, t )] . (S26)Knowing that the TOC and OTOC are given by (S23) and (S24) in terms of single-particle observables, we devotethe next sections to computing these for the different types of circuits considered: NC-ST, NC-T and C-ST. Webegin by carefully computing C ( x, t ) and C ( x, t ) for the NC-ST case. Then, we see that the results obtainedapply to NC-T and C-ST with the exception of C ( x, t ) in the C-ST case. Before proceeding to do this, we addresssome useful topics. Namely, we introduce the observables to be used, commenting on their evolution operators, andsome convenient notation for the single-particle states. We comment on the average of Haar-distributed orthogonalmatrices, which will be needed when computing the TOC and OTOC for the NC-ST and NC-T cases, and of unitarymatrices, needed for the C-ST case. Finally, and before proceeding to computing four-point correlators we see thatthe two-point correlator is trivially zero.
1. Symmetrized particle number observable
The ballistic / diffusive nature of the OTOC should not depend on the particular local observable used. Forconvenience, we consider the symmetrized particle number operator at position r , i.e. O r (0) = O r = 12 ( a † r a r − a r a † r ) = a † r a r − / ⇔ O r (0) = O r = diag( . . . , r , . . . , − r + L , . . . ) . (S27)Specifically, V and W are the symmetrized number operators initially localized at positions r and 0, respectively: V = O r ⇔ V = O r and W = O ⇔ W = O .
2. Single-particle states
Rewriting our observables using Nambu’s notation as O = 1 / A † OA allows us to work with the single-bodyobservable O in the single-body basis which includes a total of 2 L single-particle states. A single-particle statecorresponds to a site r being occupied or not, i.e. we have either a particle or a hole identified, respectively, by indices p and h : | r, p i = | r i or | r, h i = | r i . It is convenient for us to label these states by | α i = | r α , s α i = | R r α + b r α , s α i , (S28)M-6where r α = 0 , ..., L − s α ∈ { p, h } is the particle-hole index. In the second equality theposition index r α is decomposed as r α = 2 R r α + b r α , where R r α = 0 , ...., L/ − r α , r α + 1)and b r α ∈ { , } labels the position within the pair R r α . This decomposition will prove useful since it emphasizes thestructure of the circuit built out of two-site unitaries. Also, we use | ¯ α i to refer to | r α , ¯ s α i where ¯ p = h and ¯ h = p . Forconvenience, we will use the different notations interchangeably.
3. Single-body time evolution operators
In the single-body picture, each observable O r evolves according to O r ( t ) = U † ( t ) O r U ( t ). Each U ( t ) is composedof t gates U = U even U odd composed of one even and one odd layer. Establishing that U r,r +1 is the 2 L × L operatorwhich acts non-trivially on sites r and r + 1, an even / odd layer can be written as U even = L/ − Y R =0 U R, R +1 and U odd = L/ − Y R =0 U R − , R , (S29)where we label the position using the pair index introduced in (S28) and where the sites are implicitly acted upon bymodulo L , due to periodic boundary conditions (this will be omitted in the future). In its turn, U r,r +1 = u r,r +1 +¯ P r,r +1 ,where u r,r +1 = P s,s P x,x ∈{ r,r +1 } | x, s i h x, s | u | x , s i h x , s | is the 4 × r, r + 1), while ¯ P r,r +1 = P s P x/ ∈{ r,r } | x, s i h x, s | is the trivial projector into the complement of the pair. Taking thisinto account, we can write each layer in terms of the non-trivial u s as U even = L/ − X R =0 u R, R +1 and U odd = L/ − X R =0 u R − , R . (S30)This identity will prove to be useful. V. AVERAGE OF PRODUCTS OF HAAR-DISTRIBUTED ORTHOGONAL AND UNITARYMATRICES
The gates which compose our circuits are random and distributed according to the Haar measure. This is theunique measure invariant under group multiplication, weighting different regions of the probability space equally andthus behaving like a uniform distribution [3]. The average with respect to the Haar probability measure µ on thematrix probability state M( N ) is denoted by R M( N ) . . . dµ ( M ) and abbreviated as a line over the averaged quantity.We will address to the average of products of entries of some Haar-distributed gate u as moments or matrix integralsof u .Since non-conserving free fermion gates are obtained from real orthogonal matrices O according to u = V OV † , theirmoments can be obtained from those of Haar-distributed orthogonal matrices. The average of products of matrixelements of a single orthogonal matrix O with respect to the Haar probability measure µ on O( N ) is obtained inRefs. [4]. In particular, it is given by the Corollary 3.4. in [4], which uses the so called Weingarten function onthe permutation group defined in Theorem 3.9 and given explicitly in Sec. 6. In particular, we are interested in theproducts of two and four entries: h α | O | β i h α | O | β i = 1 N δ α α δ β β , (S31) h β | O | α i h β | O | α i h β | O | α i h β | O | α i = 1 N ( N − N + 2) (cid:26) δ β β δ β β h ( N + 1) δ α α δ α α − δ α α δ α α − δ α α δ α α i + δ β β δ β β h − δ α α δ α α + ( N + 1) δ α α δ α α − δ α α δ α α i + δ β β δ β β h − δ α α δ α α − δ α α δ α α + ( N + 1) δ α α δ α α i(cid:27) , . (S32)Equivalent expressions for free fermion gates are obtained by doing u = V OV † , being given by (S40), (S62) and (S63).These are very similar to the ones above, with the difference that indices of the type ¯ α appear.M-7For the specific particle conserving case, the average over free fermion gates u = diag( w, w ), where w are unitarygates, can be obtained from the moments of unitary matrices, which are given by the Corollary 2.4 in [4]. In particular,the second moment of a N × N unitary matrix w , h α | w | β i ∗ h α | w | β i , is equal to (S31) while the fourth momentis given by (S109). VI. TWO-POINT CORRELATOR
With the symmetrized particle number operator as our observable, it is natural to consider tr[ O ( t ) O x ] as thetwo-point correlator. Its dynamics is determined by the evolution of O ( t + 1 /
2) = U † O ( t ) U , where U alternatesbetween U even and U odd , given by (S30). This can be averaged considering h α | u | β i ∗ h α | u | β i = 1 /N δ α α δ β β ,becoming O ( t + 1 /
2) = X R X r,r ∈{ R, R +1 } X s,s | r , s i h r, s | O ( t ) | r, s i h r , s | , (S33)where we considered U = U even , with analogous results following for U = U odd . Since P s h r, s | O ( t ) | r, s i = 0 (dueto PH symmetry), it follows that tr[ O ( t ) O x ] = 0. This leads us to look for non-trivial behaviour in higher ordercorrelators. VII. NC-ST: AVERAGED TOC, C ( x, t ) We start by considering the spatio-temporal noisy free fermion circuit without particle conservation (NC-ST), drawnin Fig. S1 (a). We first compute C ( r, t ) and then C ( r, t ). The dynamics of the later will prove to be more intricateand so we employ a vectorized notation which helps clarifying its behaviour. For a matter of consistency, this is alsoapplied to C ( r, t ). Accordingly, we consider C ( r, t ) as the overlap of two vectorized operators C ( r, t ) = tr[ O ( t ) O r ] = hh O r || O ( t ) ii , (S34)with || O r ( t ) ii given by || O r ( t ) ii = L − X x,x =0 X s,s || x, s ; x , s ii h x, s | O r ( t ) | x , s i , (S35)where || x, s ; x , s ii = | x, s i ⊗ h x , s | T and the particle-hole index s sums over p and h (we usually omit these). Also,with O r the symmetrized number operator at position r given by (S27) we have || O r ii = X s || r, s ; r, s ii . (S36)Having established this vectorized formalism as our framework, we proceed to average C ( r, t ) over multiple real-izations of the random circuit and we evaluate the ensuing dynamics. S.1. Time evolution of || O ( t ) ii The dynamics of C ( r, t ) is completely contained in || O ( t ) ii . The operator O ( t ) evolves as O ( t + 1) = U † O ( t ) U which, in the vectorized notation, translates to || O ( t + 1) ii = U † ⊗ U T || O ( t ) ii . (S37)Each circuit is a succession of layers of two-site gates acting on even and odd pairs of sites. In one unit of time oneeven and one odd layers are applied, U = U odd U even , such that || O ( t + 1) ii = (cid:16) U † ⊗ U T (cid:17) even (cid:16) U † ⊗ U T (cid:17) odd || O ( t ) ii = L/ − X R,R =0 u † R, R +1 ⊗ u T R , R +1 L/ − X R,R =0 u † R − , R ⊗ u T R − , R || O ( t ) ii , (S38)where the sites are labelled using the pair indices and where we used (S30).M-8
1. Average over disorder realizations
Next, we average random realizations of the circuit to obtain the leading order behaviour of C ( r, t ). The average || O ( t ) ii reduces to the average over the two-site gates composing the circuit. Since these are randomly chosen, theyare uncorrelated and we can perform the average independently at different half-time steps. Thus, we only need tofind u † r,r +1 ⊗ u Tr ,r +1 = || α α ii h β | u r,r +1 | α i ∗ h β | u r ,r +1 | α i || β β ii . (S39)This demands that we compute the average of products of entries of some Haar-distributed free fermion gate u ,also called the moments or matrix integrals of u . A short discussion on the topic was presented in Sec. V.Summarizing, expressions for the moments of u can be obtained from equivalent expressions existing for an Haar-distributed orthogonal matrix O and only moments of an even number of entries of u are non-zero. This implies u † r,r +1 ⊗ u Tr ,r +1 = δ r,r u † r,r +1 ⊗ u Tr,r +1 . Considering that the second moment of some free fermion gate u is given by(see Sec. V) h α | u | β i ∗ h α | u | β i = 1 N δ α α δ β β , (S40)where N = rank( u ) = 4, we get u † r,r +1 ⊗ u Tr,r +1 = 1 N || αα ii hh ββ || = || φ r ii hh φ r || , (S41)i.e. the projector evolving the pair of sites ( r, r + 1) from some time t to t + 1 /
2, with || φ r ii = 12 X x ∈{ r,r +1 } X s || x, s ; x, s ii , (S42)which obeys hh φ r || φ r ii = δ r ,r + δ r ,r ± . Also, note that || φ r ii hh φ r || leads O ( t ) to be always diagonal. Pluggingthese results into the averaged equivalent of (S38) leads to || O ( t + 1) ii = L/ − X R =0 || φ R ii hh φ R || L/ − X R =0 || φ R − ii hh φ R − || || O ( t ) ii , (S43)To know || O ( t ) ii at any t we just need to start from the initial condition || O (0) ii = P s || , s ; 0 , s ii and apply theabove expression successively. Next, we do this and simplify (S43).
2. Dynamics of || O ( t ) ii Starting with || O (0) ii = P s || , s ; 0 , s ii and applying the first odd layer of the circuit we obtain || O ( t = 1 / ii = L/ − X R =0 || φ R − ii hh φ R − || || O (0) ii = || φ − ii hh φ − || O (0) ii = || φ − ii , (S44)to which we can apply a second even layer to get || O ( t = 1) ii = L/ − X R =0 || φ R ii hh φ R || || O (1 / ii = (cid:16) || φ − ii hh φ − || + || φ ii hh φ || (cid:17) || φ − ii = 12 (cid:16) || φ − ii + || φ ii (cid:17) . (S45)Also, mixing between neighbouring pairs in successive layers is mediated by || φ r ii hh φ r || φ r ± ii = 12 || φ r ii , (S46)M-9which guarantees that the subspace {|| φ r ii} is closed under time evolution. This, allied with (S45), allows thedecomposition || O ( t ) ii = (P L/ − R =0 || φ R ii hh φ R || O ( t ) ii , for t integer P L/ − R =0 || φ R − ii hh φ R − || O ( t ) ii , for t half-integer (S47)i.e. {|| φ R ii | R = 0 , . . . , L/ − } and {|| φ R − ii | R = 0 , . . . , L/ − } for t integer and half-integer, respectively,are the natural orthonormal basis for our problem. This is, they take full advantage of the structure coming fromthe average of gates, meaning that hh r, s ; r, s || O ( t ) ii = hh r + 1 , s ; r + 1 , s || O ( t ) ii , where ( r, r + 1) forms a pair, andof the PH symmetry of O ( t ), which leads to hh r, p ; r, p || O ( t ) ii = hh r, h ; r, h || O ( t ) ii . Having established the goodbasis to use, we can employ (S46) to obtain, for t ≥ / hh φ r || O ( t + 1 / ii = hh φ r || (cid:16) || φ r − ii hh φ r − || + || φ r +1 ii hh φ r +1 || (cid:17) || O ( t ) ii = 12 (cid:16) hh φ r − || + hh φ r +1 || (cid:17) || O ( t ) ii , (S48)i.e. the time-evolution of hh φ r || O ( t ) ii can be understood as an averaging process with contributions from the nearestneighbours. Applying this twice gives hh φ r || O ( t + 1) ii = 14 (cid:16) hh φ r − || + 2 hh φ r || + hh φ r +2 || (cid:17) || O ( t ) ii . (S49)This recursive expression associated to || O ( t ) ii = P L/ − R =0 || φ R ii hh φ R || O ( t ) ii is equivalent to (S43), but moretransparent than the latter. S.2. Prescription to obtain C ( r, t ) Now that we have reduced || O ( t ) ii to its core, we focus back on the TOC.For t = 0, we use the initial condition || O (0) ii = P s || , s ; 0 , s ii in (S34) to obtain C ( r,
0) = C ( r,
0) = hh O r || O (0) ii = 2 δ r, . (S50)For t ≥
1, we can decompose || O ( t ) ii using (S47) to obtain C ( r, t ) = hh O r || O ( t ) ii = hh O r || L/ − X R =0 || φ R ii hh φ R || O ( t ) ii = L/ − X R =0 ( δ r, R + δ r, R +1 ) hh φ R || O ( t ) ii , (S51)i.e. C (2 R, t ) = C (2 R + 1 , t ) = hh φ R || O ( t ) ii . The dynamics of C ( r, t ) is then determined by that of hh φ R || O ( t ) ii ,which we just saw is given by (S49), starting with (S45) as the initial condition, i.e. C (2 R, t + 1) = C (2 R + 1 , t + 1) = ( ( δ R − , + δ R, ) , for t = 0 (cid:16) C (2 R − , t ) + 2 C (2 R, t ) + C (2 R + 2 , t ) (cid:17) , for t ≥ . (S52)Having broken down the TOC given initially by (S34), we end up with a very clean picture: the circuit’s structureleads C ( r, t ) to depend only on the pair R to which r belongs and evolving it comes down to performing a weightedaverage as specified in (S52). This averaging process is represented pictorially in Fig. S2. In the panel (c) the‘brickwall’ structure of the circuit ensures that, as time evolves, C ( r, t ) spreads across the system. Although theboundary in the r − t plane between the region with zero and non-zero C ( r, t ) describes a light cone, the weightsare concentrated around r = 0 such that the process is diffusive and not ballistic. Next, we argue in favour of thisdiffusive behaviour by going to the continuum limit.M-10FIG. S2: The values in the boxes are C ( r, t ) in the r − t plane (in empty boxes C ( r, t ) = 0). Panels (a) and (b)translate Eqs. (S48) and (S49) to a schematic, representing the evolution of C ( r, t ) in one half and one full timesteps. The scheme in (b) is obtained by applying (a) twice. (c) The previous schemes can be applied to obtain C ( r, t ),starting with C ( r,
0) = 2 δ r, as the initial condition. Even layers are depicted in blue while odd layers appear in red. S.3. Continuum limit of C ( r, t ) The OTOC is described by (S52) for t ≥
1. We show that this is a discrete diffusion equation in 1D by recoveringthe continuum diffusion equation from it.We consider the scaling form aC ( x = ra, τ = ta ). For a = 1, this coincides with the discrete C ( r, t ). If we let a →
0, this approximates C ( r, t ) in the continuum limit. After making this identification in (S52), where r = 2 R , weTaylor expand it up to O ( a ), obtaining the one-dimensional continuum diffusion equation ∂ τ C ( x, τ ) = D ∂ x C ( x, τ ) , (S53)with diffusion coefficient D = 1. Using the initial condition (S50), which translates to C ( x, t = 0) = 2 δ ( x ), thesolution to (S53) approximating C ( r, t ) for large t and L is C (2 R, t ) = C (2 R + 1 , t ) ’ aC ( x = 2 Ra, τ = ta ) = A √ πσ ( t ) exp (cid:18) − (2 R ) σ ( t ) (cid:19) , (S54)i.e. a Gaussian normalized to A = 2 with standard deviation σ = √ t . Having recovered the 1D diffusion equation(S53) shows that, for large enough times, the exact (S52) diffuses in the r − t plane. VIII. NC-ST: AVERAGED OTOC, C ( r, t ) We move on to calculate the OTOC itself for the non-conserving circuit with randomness in space and time (NC-ST).Although this will prove to be more intricate than C ( r, t ), the procedure follows the same steps.M-11We start by writing C ( r, t ) as the overlap of two vectorized operators C ( r, t ) = tr [ O ( t ) O r O ( t ) O r ] = hh Q r || S || Q ( t ) ii , (S55)with Q r ≡ O r ⊗ O r and || Q r ( t ) ii = X αβµν || αβµν ii h αβ | Q r ( t ) | µν i , (S56) S = X αβµν || αβµν ii hh αβνµ || , (S57)where || αβµν ii = | αβ i ⊗ h µν | T and the operator S reorders the indices such that hh Q r || S || Q ( t ) ii is the trace presentin (S55). Also, with O r the symmetrized number operator at position r , we have || Q r ii = X s || rs, rs, rs, rs ii − || rs, r ¯ s, rs, r ¯ s ii . (S58)We observe that in the vectorized notation a parallel is established between C ( r, t ) = hh O r || O ( t ) ii and C ( r, t ) = hh Q r || S || Q ( t ) ii : || O r ii and || O ( t ) ii are to C ( r, t ) as || Q r ii and || Q ( t ) ii are to C ( r, t ). The dynamics of the latteris set by || Q ( t ) ii , which we proceed to evaluate. S.1. Time evolution of || Q ( t ) ii The dynamics of C ( r, t ) is determined by that of || Q ( t ) ii . As done for C ( r, t ), we first compute how || Q ( t ) ii evolves and then we average the result obtained.From O ( t + 1) = U † O ( t ) U it follows that Q ( t + 1) = ( U † ⊗ U † ) Q ( t )( U ⊗ U ) which, in the vectorized notation,translates to || Q ( t + 1) ii = U † ⊗ U † ⊗ U T ⊗ U T || Q ( t ) ii . (S59)In one time step two layers of the circuit are applied such that U = U odd U even . Using (S30), it follows that || Q ( t + 1) ii = ( U † ⊗ U † ⊗ U T ⊗ U T ) even ( U † ⊗ U † ⊗ U T ⊗ U T ) odd || Q ( t ) ii = L/ − X R,R ,R ,R =0 u † R, R +1 ⊗ u † R , R +1 ⊗ u T R , R +1 ⊗ u T R , R +1 L/ − X R,R ,R ,R =0 u † R − , R ⊗ u † R − , R ⊗ u T R − , R ⊗ u T R − , R || Q ( t ) ii , (S60)valid for some realization of the circuit.
1. Average over disorder realizations
Wishing to obtain || Q ( t ) ii , we move on to average (S60) over random realizations of the circuit. As before, we cantake the average of each layer independently from the others. Thus, we wish to obtain u † r,r +1 ⊗ u † r ,r +1 ⊗ u Tr ,r +1 ⊗ u Tr ,r +1 = || α α α α ii h β | u r,r +1 | α i ∗ h β | u r ,r +1 | α i ∗ h β | u r ,r +1 | α i h β | u r ,r +1 | α i hh β β β β || . (S61)To perform this average of a product of entries of u we refer again to Sec. V. Remember that only even moments of u are non-zero. The second moment is given either by (S40) or h β | u | α i h β | u | α i = h β | u | α i ∗ h β | u | α i ∗ = 1 N δ ¯ α α δ ¯ β β (S62)M-12and the fourth moment is given by h β | u | α i ∗ h β | u | α i ∗ h β | u | α i h β | u | α i = 1 N ( N − N + 2) (cid:26) δ β β δ β β h ( N + 1) δ α α δ α α − δ ¯ α α δ ¯ α α − δ α α δ α α i + δ ¯ β β δ ¯ β β h − δ α α δ α α + ( N + 1) δ ¯ α α δ ¯ α α − δ α α δ α α i + δ β β δ β β h − δ α α δ α α − δ ¯ α α δ ¯ α α + ( N + 1) δ α α δ α α i(cid:27) , (S63)with N = rank( u ) = 4. This means that, when performing the average in (S61), either two or four entries of u mustrefer to the same unitary, leaving us with four possible terms. The case where all entries refer to the same unitary(i.e. r = r = r = r ) contributes with one term while the other three come from two different unitaries appearing,for which there are three configurations corresponding to the pairings we can make out of four elements: (12)(34),(13)(24) and (14)(23). Taking (S40), (S62) and (S63) into account, the sum over (S61) becomes L/ − X r,r ,r ,r =0 u † r,r +1 ⊗ u † r ,r +1 ⊗ u Tr ,r +1 ⊗ u Tr ,r +1 = L/ − X r,r =0 v r,r , (S64)with v r,r = (1 − δ rr ) (cid:16) || Θ ii hh Θ || + || ¯Θ ii hh ¯Θ || + || ˜Θ ii hh ˜Θ || (cid:17) r,r + δ rr ( N + 2) (cid:18)h ( N + 1) || Θ ii − || ¯Θ ii − || ˜Θ ii i hh Θ || + h − || Θ ii + ( N + 1) || ¯Θ ii − || ˜Θ ii i hh ¯Θ || + h − || Θ ii − || ¯Θ ii + ( N + 1) || ˜Θ ii i hh ˜Θ || (cid:19) r,r , (S65)where we use the definitions || Θ r,r ii = 1 √ g r,r X r α ∈{ r,r +1 } r β ∈{ r ,r +1 } X s α ,s β || αββα ii , (S66) || ¯Θ r,r ii = 1 √ g r,r X r α ∈{ r,r +1 } r β ∈{ r ,r +1 } X s α ,s β || α ¯ αβ ¯ β ii , (S67) || ˜Θ r,r ii = 1 √ g r,r X r α ∈{ r,r +1 } r β ∈{ r ,r +1 } X s α ,s β || αβαβ ii , (S68)with g r,r = N ( N − δ r,r ). Taking into account that r and r have the same parity, these obey hh i rr || i qq ii = 16 g r,r δ q,r δ q ,r + 4 √ g rr g qq δ q,r ± δ q ,r ± , (S69) hh i rr || j qq ii = δ r,r (cid:16) δ q,r δ q ,r + δ q,r +1 δ q ,r +1 + δ q,r − δ q ,r − (cid:17) , (S70)for i = j ∈ { Θ , ¯Θ , ˜Θ } . Notice that, while the evolution of || O ( t ) ii in C ( r, t ) involved mixing within each pair ofsites, ( r, r + 1), through the projector || φ r ii hh φ r || , here the evolution through v r,r entails mixing between each twopairs ( r, r + 1) and ( r , r + 1). Thus, we expect the time evolution of || O ( t ) ii to be a two-dimensional process inspace.Finally, considering (S64), considering the average of (S60) becomes || Q ( t + 1) ii = L/ − X R,R =0 v R, R L/ − X R,R =0 v R − , R − || Q ( t ) ii , (S71)with v r,r given by (S65). Starting with || Q (0) ii given by (S58), we only need to apply this expression recursively toobtain || Q ( t ) ii at successive time steps. Next, we do precisely this and simplify (S71).M-13
2. Dynamics of || Q ( t ) ii Let us start by defining || Θ ¯Θ r,r ii = || Θ r,r ii − || ¯Θ r,r ii . (S72)Starting with || Q (0) ii given by (S58), applying the first odd layer of the circuit leads to || Q ( t = 1 / ii = L/ − X R,R =0 v R − , R − || Q (0) ii = v − , − || Q (0) ii = 1 √ || Θ ¯Θ − , − ii , (S73)to which we can apply a second even layer to obtain || Q ( t = 1) ii = L/ − X R,R =0 v R, R || Q (1 / ii = ( v − , − + v − , + v , − + v , ) || Q (1 / ii = 16 √ (cid:16) || Θ ¯Θ − , − ii + || Θ ¯Θ , ii (cid:17) + 16 (cid:16) || Θ ¯Θ − , ii + || Θ ¯Θ , − ii (cid:17) . (S74)Also, we can use (S69) and (S70) to see that || Θ ¯Θ r,r ii evolves under v r,r according to v r + a,r + a || Θ ¯Θ r,r ii = (cid:18) − δ r,r (cid:19) || Θ ¯Θ r + a,r + a ii , (S75) v r + a,r − a || Θ ¯Θ r,r ii = (cid:20)
14 + (cid:18) √ − (cid:19) ( δ r,r + δ r + a,r − a ) (cid:21) || Θ ¯Θ r + a,r − a ii , (S76)where a = ±
1. These relations control the mixing occurring between neighbouring pairs in successive time steps andensure that the subspace {||
Θ ¯Θ r,r ii} is closed under time evolution. This, allied with (S73), allows the decomposition || Q ( t ) ii = (P L/ − R =0 12 || Θ ¯Θ R, R ii hh Θ ¯Θ R, R || Q ( t ) ii , for t integer P L/ − R =0 12 || Θ ¯Θ R − , R − ii hh Θ ¯Θ R − , R − || Q ( t ) ii , for t half-integer (S77)where the factor 1 / hh Θ ¯Θ r,r || Θ ¯Θ q,q ii = 2 δ r,q δ r ,q (where the indices are restricted to having the sameparity). As {|| φ r ii} was the optimal basis on which to express || O ( t ) ii , so it happens that { √ || Θ ¯Θ R, R ii | R, R =0 , . . . , L/ − } and { √ || Θ ¯Θ R − , R − ii | R, R = 0 , . . . , L/ − } are the natural orthonormal basis on which toexpress || Q ( t ) ii . Let us comment on the states || αβµν ii grouped under the same || Θ ¯Θ r,r ii . Due to PH symmetry, outof the three types of states surviving the average of Haar-distributed gates − || αββα ii , || α ¯ αβ ¯ β ii and || αβαβ ii − onlythe first two survive; explicitly: hh αβαβ || Q ( t ) ii = − hh αβαβ || Q ( t ) ii ⇒ P s α hh αβαβ || Q ( t ) ii = 0. Furthermore,PH symmetry implies hh αββα || Q ( t ) ii = − hh α ¯ αβ ¯ β || Q ( t ) ii = | h α | O ( t ) | β i | such that || Θ ¯Θ r,r ii groups the states || Θ r,r ii and − || ¯Θ r,r ii which evolve equally. The structure appearing restricts the dynamics of || Q ( t ) ii such that itis simpler than one might have guessed, with only 8 L entries out of (2 L ) being non-zero.Having established the good basis to use, we can employ (S75) and (S76) to obtain, for t ≥ / hh Θ ¯Θ r,r || Q ( t + 1 / ii = hh Θ ¯Θ r,r || (cid:16) v r − ,r − + v r +1 ,r +1 + v r − ,r +1 + v r +1 ,r − (cid:17) || Q ( t ) ii = tr " m r,r (cid:18) hh Θ ¯Θ r − ,r − || hh Θ ¯Θ r − ,r +1 ||hh Θ ¯Θ r +1 ,r − || hh Θ ¯Θ r +1 ,r +1 || (cid:19) T || Q ( t ) ii , (S78)where m r,r is a coefficient matrix which assumes different values depending on r and r : m r,r =
16 12 √ √ ! , m r,r +2 = m Tr,r − = (cid:18)
14 1412 √ (cid:19) , m r,r = 14 (cid:18) (cid:19) , (S79)with the last m r,r being valid for | r − r | >
2. In the bulk, i.e. for | r − r | >
2, this evolution equation is an averageover neighbouring sites, as depicted in Fig. S3. This is similar to what is shown to happen for C ( r, t ) in Fig. S2,where the process occurring is one instead of two-dimensional.M-14FIG. S3: Scheme of the evolution of the bulk hh Θ ¯Θ r,r || Q ( t ) ii| | r − r | > (values in boxes) in one-half time step. The bulkevolution uses the third matrix m r,r in (S79) and is effectively an average of neighbouring elements hh Θ ¯Θ r,r || Q ( t ) ii whose position indices belong to the pairs ( r, r + 1) and ( r , r + 1). Depending on whether t + 1 / hh Θ ¯Θ r,r || Q ( t + 1) ii = tr M r,r hh Θ ¯Θ r − ,r − || hh Θ ¯Θ r − ,r || hh Θ ¯Θ r − ,r +2 ||hh Θ ¯Θ r,r − || hh Θ ¯Θ r,r || hh Θ ¯Θ r,r +2 ||hh Θ ¯Θ r +2 ,r − || hh Θ ¯Θ r +2 ,r || hh Θ ¯Θ r +2 ,r +2 || T || Q ( t ) ii (S80)with M r,r a coefficient matrix given by M r,r = 172 √ √ √ √ √ √ , M r,r +2 = M Tr,r − = 148 √
13 64 √ M r,r +4 = M Tr,r − = 116 √ , M r,r = 116 , (S81)where the last M r,r is valid for | r − r | >
4. This last matrix establishes the evolution of the bulk of hh Θ ¯Θ r,r || Q ( t ) ii as a weighted average resulting from applying the scheme seen in Fig. S3 twice.At last, starting with (S74) as our initial condition, || Q ( t ) ii is given by the recursive expression (S80) alongside(S77). This is equivalent to (S71) but more transparent than the latter. Now that we know the dynamics of || Q ( t ) ii ,we can obtain that of the OTOC. S.2. Prescription to obtain C ( r, t ) Now that we have reduced || Q ( t ) ii to its core, we focus back on the C ( r, t ). At t = 0, || Q (0) ii is given by (S58)such that (S55) becomes C ( r, t = 0) = C ( r, t = 0) = hh Q r || S || Q (0) ii = 2 δ r, . (S82)For t ≥
1, we can decompose || Q ( t ) ii using (S77) to obtain C ( r, t ) = hh Q r || S || Q ( t ) ii = hh Q r || S L/ − X R,R =0 || Θ ¯Θ R, R ii hh Θ ¯Θ R, R || Q ( t ) ii = 12 √ L/ − X R =0 ( δ r, R + δ r, R +1 ) hh Θ ¯Θ R, R || Q ( t ) ii . (S83)For convenience, we can define K r,r ( t ) = 12 √ hh Θ ¯Θ r,r || Q ( t ) ii , (S84)M-15whose diagonal gives the OTOC: C (2 R, t ) = C (2 R + 1 , t ) = K R, R ( t ). Thus, we can use (S74) to get the OTOCat t = 1, which is C (2 R, t = 1) = C (2 R + 1 , t = 1) = 118 (cid:16) δ R, − + δ R, (cid:17) , (S85)and at subsequent time steps the dynamics of C ( r, t ) is determined by that of hh Θ ¯Θ r,r || O ( t ) ii ⇔ K r,r ( t ) we justsaw to be given by (S80), i.e. for t ≥ C (2 R, t + 1) = C (2 R + 1 , t + 1) = diag[ K R, R ( t + 1)]= diag tr M R, R K R − , R − K R − , R K R − , R +2 K R, R − K R, R K R, R +2 K R +2 , R − K R +2 , R K R +2 , R +2 T ( t ) (S86)Note that, although the OTOC is given by the diagonal of K r,r ( t ), we must keep track of the whole matrix to learnits dynamics. Before proceeding to inspect the dynamics of K r,r ( t ) in the continuum limit, we give a pictorial viewin Fig. S4 of how K r,r ( t ) (and thus, the OTOC) evolves in half-time steps.FIG. S4: The values in the boxes correspond to K r,r ( t ) in the r − r plane, shown here for t ∈ { / , , / , } (emptyboxes have zero value). Starting with K r,r ( t ) = 2 δ r,r , as the initial condition, we evolve this by applying (S78)successively. The highlighted diagonal gives C ( r, t ). Note that the boxes group the pairs of sites ( r, r + 1) and( r , r + 1) enhancing the symmetry K r,r ( t ) = K r +1 ,r +1 ( t ) = K r +1 ,r ( t ) = K r,r +1 ( t ). Depending on whether t isinteger or half-integer (blue/red), these pairs are even or odd. S.3. Continuum limit of C ( r, t ) We will see that the dynamics of K r,r ( t ) is a diffusive process by showing that its evolution equation can beapproximated by a two-dimensional diffusion equation in the continuum limit.In the following, we obtain the continuum limit of (S78) which evolves K r,r ( t ) in half time steps. Note that thisequation (or, more specifically, the matrix m r,r ) differs depending on whether r = r , r = r ± | r − r | > L → ∞ ) the regions with r = r and r = r ± K r,r ( t ) does not diverge, they are not expected to influence the behaviour of the bulk. Thus, we start byobtaining a continuum diffusion equation for the bulk, i.e. for K r,r ( t ) in | r − r | >
2. Then, since we wish to obtainthe OTOC, which is K r,r ( t ), we see that the diagonal differs from the bulk by a constant factor.Consider the scaling form a K ( x = ra, x = r a, τ = ta ). For a = 1, this coincides with the discrete K r,r ( t ). If welet a →
0, this approximates K r,r ( t ) in the continuum limit, i.e. lim L →∞ , t →∞ K r,r ( t ) = lim a → a K ( ra, r a, ta ).After making this identification in (S80) and considering m r,r | | r − r | > to be the isotropic bulk coefficient matrixgiven in (S79), we Taylor expand K bulk ( x, x , τ ) (i.e. the quantity K ( x, x , τ ) evolving isotropically under (S80) with m r,r | | r − r | > ) up to O ( a ) such that (S78) becomes ∂ τ K bulk ( x, x , τ ) = D ( ∂ x + ∂ x ) K bulk ( x, x , τ ) , (S87)i.e. the isotropic diffusion equation in 2D, with diffusion constant D = 1.M-16We return to the question of the diagonal, K r,r ( t ). Namely, we will account for the anisotropy present in theevolution of the diagonal which will impact the way the initial condition, which is precisely in the diagonal, will bedistributed along different directions, something which was disregarded to obtain (S87). To begin with, we distinguishbetween r = r and r = r and approximate K r,r ( t ) by a K ( x, x, τ ) for r = r and by K bulk ( x, x , τ ) for r = r . Makingthis identification in (S78) gives K bulk ( x − a, x + a, τ + 1 /
2) = 1 / K bulk ( x − a, x, τ ) + K bulk ( x − a, x + 2 a, τ ) + K bulk ( x, x + 2 a, τ )] + 1 / (2 √ K ( x, x, τ ) which, after Taylor expanding K bulk up to the lowest order, O ( a ), yields K ( x, x, τ ) ≈ cK bulk ( x, x, τ ) , (S88)with c = √ /
2. Note that we approximated K ( x, x ± a, τ ) by the bulk K bulk : similarly to what is done to obtain(S88), we can Taylor expand K bulk up to O ( a ) in the evolution equation of K bulk ( x − a, x + 2 a, τ + 1 /
2) to obtain K ( x, x ± a, τ ) ∼ K bulk ( x, x ± a, τ ).Considering (S88) in the evolution equation of the diagonal (S78) leads to the weighted average K bulk ( x, x, τ +1 /
2) =tr " m r,r (cid:18) K bulk ( x − , x − , τ ) K bulk ( x − , x + 2 , τ ) K bulk ( x + 2 , x − , τ ) K bulk ( x + 2 , x + 2 , τ ) (cid:19) T with m r,r = 16 (cid:18) (cid:19) . (S89)This allows us to directly compare the diagonal with the bulk. The anisotropy in the new coefficient matrix m r,r leads to diffusion with weighting factor f = 1 / K ( x, x, τ ), we take the diagonal of K bulk ( x, x , τ ) given as the solution of (S87) and multiply it by c − whichaccounts for the different scaling of the diagonal − and by f − which accounts for the anisotropy present.Finally, we can obtain the solution of (S87) and related it to the OTOC. The initial condition (S82) translates to K ( x, x, τ = 0) ≈ cK bulk ( x, x, τ = 0) = Aδ ( x ) δ ( x ), with A = 2 and where we used (S88), such that the solution to(S87) is cK bulk ( x, x , τ ) = A πτ D exp (cid:18) − x + x D τ (cid:19) . (S90)Finally, the OTOC is approximated by (S90) times the corrective factor fC (2 R, t ) = C (2 R + 1 , t ) = K R, R ( t ) ’ f ca K bulk (2 Ra, Ra, τ ) = A ( t ) √ πσ ( t ) exp (cid:18) − (2 R ) σ ( t ) (cid:19) , (S91)i.e. in the continuum limit, the OTOC is described by a Gaussian with broadening width σ ( t ) = √ t and varyingnormalization A ( t ) = f A/ (2 √ πt ), with A = 2 and f = 1 /
2. We conclude that C ( r, t ) is the diagonal of a quantitywhich diffuses in 2D, analogously to C ( r, t ) which instead diffused in 1D. IX. NC-T: AVERAGED TOC, C ( r, t ) We now turn to the second instance of free fermion evolution, NC-T, which, like the NC-ST, case does not conserveparticle number and, unlike it, it displays randomness only in time. In this circuit, each layer U even and U odd isbuilt acting with the same randomly chosen gate u on all pairs of sites such that the system is invariant under spacetranslation by multiples of two sites, as shown drawn in Fig. S1 (b). Next, we derive the behaviour of C ( r, t ) in thiscase by pointing out the differences from the NC-ST case, which reside in the average of unitaries. We saw in Sec.VII S.1 that the average of unitaries u † r,r +1 ⊗ u Tr ,r +1 is non-zero if u r,r +1 = u r ,r +1 . While in the NC-ST case thisdemanded r = r , the spatial homogeneity of NC-T softens this condition to r = r modulo 2. As a consequence, theoperator P r || φ r ii hh φ r || which evolves || O ( t ) ii in the NC-ST is generalized to P rr || φ rr ii hh φ rr || with || φ r ii −→ || φ rr ii = 12 X b ∈{ , } X s || r + b, s ; r + b, s ii . (S92)Considering the new evolution operator, we can start from || O (0) ii = P s || , s ; 0 , s ii to obtain || Q (1 / ii = L/ − X R,R =0 || φ R, R ii hh φ R, R || || Q (0) ii = || φ − − ii = || φ − ii . (S93)M-17The interaction between neighbouring pairs in successive layers analogous to (S46) reduces to || φ r,r ii hh φ r,r || φ r + a,r + a ii = 12 || φ r,r ii , (S94)with a = ±
1. From this, we get an evolution equation analogous to (S48) hh φ r,r || O ( t + 1 / ii = hh φ r,r || (cid:16) || φ r − ,r − ii hh φ r − ,r − || + || φ r +1 ,r +1 ii hh φ r +1 ,r +1 || (cid:17) || O ( t ) ii = 12 (cid:16) hh φ r − ,r − || + hh φ r +1 ,r +1 || (cid:17) || O ( t ) ii . (S95)This, alongside the initial condition (S93), leads only terms of the type || φ r,r ii = || φ r ii to arise. Thus, the aboveexpression is effectively reduced to (S48) which, allied with the initial condition || Q (1 / ii ) = || φ − ii being the samewe had for NC-ST, leads C ( r, t ) in the NC-T case to be exactly equal to the one obtained for NC-ST. Although thespatial homogeneity of the new circuit threatens the TOC dynamics to be described by a two-dimensional process,the initial condition fixes it to the one-dimensional process already known for NC-ST. X. NC-T: AVERAGED OTOC, C ( r, t ) We continue to consider the NC-T case, focusing now on C ( r, t ). As happens with C ( r, t ), the difference from theNC-ST case lies in the average of unitaries.We saw in Sec. VIII S.1 that, for the average u † r ,r +1 ⊗ u † r ,r +1 ⊗ u Tr ,r +1 ⊗ u Tr ,r +1 to be non-zero„ either twoor four unitaries must be equal, which, in the NC-ST case, translates to r = r ∧ r = r or r = r ∧ r = r or r = r ∧ r = r or r = r = r = r . Similarly to what we just saw for C ( r, t ), these restrictions are softenedin the NC-T case: since all the gates composing one layer are the same the four unitaries can refer to different pairssuch that v r,r is generalized to v r ,r ,r ,r by transforming (S66), (S67) and (S68) according to || Θ r,r ii → || Θ r ,r ,r ,r ii = 1 √ g r ,r X b ,b ∈{ , } s ,s || r + b , s ; r + b , s ; r + b , s ; r + b , s ii , (S96) || ¯Θ r,r ii → || ¯Θ r ,r ,r ,r ii = 1 √ g r ,r X b ,b ∈{ , } s ,s || r + b , s ; r + b , ¯ s ; r + b , s ; r + b , ¯ s ii , (S97) || ˜Θ r,r ii → || ˜Θ r ,r ,r ,r ii = 1 √ g r ,r X b ,b ∈{ , } s ,s || r + b , s ; r + b , s ; r + b , s ; r + b , s ii . (S98)These obey relations equivalent to (S69) and (S70) , hh i r ,r ,r ,r || i q ,q ,q ,q ii = 16 g r ,r δ q r δ q r δ q r δ q r + 4 √ g r r g q q δ q ,r + a δ q ,r + b δ q ,r + a δ q ,r + b , (S99) hh i r ,r ,r ,r || j q ,q ,q ,q ii = δ r r (cid:16) δ q r δ q r δ q r δ q r + δ q ,r +1 δ q ,r +1 δ q ,r +1 δ q ,r +1 + δ q ,r − δ q ,r − δ q ,r − δ q ,r − (cid:17) , (S100)with a, b = ±
1, which can be used to obtain the interaction between neighbouring pairs in successive layers v r + a,r + a,r + a,r + a || Θ ¯Θ r ,r ,r ,r ii = (cid:18) − δ r,r (cid:19) || Θ ¯Θ r + a,r + a,r + a,r + a ii , (S101) v r + a,r − a,r + a,r − a || Θ ¯Θ r ,r ,r ,r ii = (cid:20)
14 + (cid:18) √ − (cid:19) ( δ r,r + δ r + a,r − a ) (cid:21) || Θ ¯Θ r + a,r − a,r + a,r − a ii , (S102)Applying this to || Q (0) ii , given by (S58), results in || Q ( t = 1 / ii = L/ − X R ,R ,R ,R =0 v R − , R − , R − , R − || Q (0) ii = v − , − , − , − || Q (0) ii = 1 √ || Θ ¯Θ − , − , − , − ii = 1 √ || Θ ¯Θ − , − ii , (S103)M-18the same initial condition we have for the NC-ST case. Since in (S101) and (S102) r and r as well as r and r are restricted to varying by the same amount, the subspace {|| Θ ¯Θ r ,r ,r ,r ii} is closed under time evolution. Thus,starting with (S103) as our initial condition, the relation || Θ ¯Θ r ,r ,r ,r ii = || Θ ¯Θ r ,r ii effectively reduces the evolutionrelations (S101) and (S102) to the ones we had for NC-ST, (S75) and (S76). It follows that C ( r, t ) for the NC-Tcase is equal to the one obtained before for the NC-ST case. Equivalently to what happens for C ( r, t ), although thespatial homogeneity of the new circuit threatens the OTOC dynamics to be described by a four-dimensional process,initial conditions fix it to the two-dimensional process already known for NC-ST. XI. C-ST: AVERAGED TOC, C ( r, t ) Lastly, we focus on the third instance of free fermion evolution, C-ST, a restriction of the NC-ST case to a particleconserving process admitting spatio-temporal noise. The difference from NC-ST lies in the new two-site unitariesconserving particle number, i.e. we have no anomalous terms a i a j or a † i a † j such that the ph and hp sectors of u arenull (by ph sector we mean those | r, s i h r, s | u | r , s i h r , s | with s = p and s = h ). The new unitaries are diagonal byblocks: u = diag( w, w ), with w × C ( r, t ) andthen C ( r, t ).Knowing that w ∗ αβ w µν = N δ αµ δ βν according to (S31) in Sec. V, the average equivalent to (S40) becomes h α | u | β i ∗ h α | u | β i = h α | w | β i ∗ h α | w | β i δ s α s β δ s α s β = 1 N δ α α δ β β δ s α s β δ s α s β , (S104)where N = rank( w ) = 2, with the difference from the NC-ST being a restriction to the non-anomalous sectors pp and hh . This results in the evolution operator (S41) becoming || φ r ii hh φ r || −→ Φ r = 12 X x ∈{ r,r +1 } x ∈{ r ,r +1 } X s || x, s ; x, s ii hh x , s ; x , s || . (S105)The initial condition || O (0) ii = P s || , s ; 0 , s ii evolves under this as || Q (1 / ii = L/ − X R =0 Φ R − || Q (0) ii = Φ − || Q (0) ii = || φ − ii , (S106)i.e. || Q (1 / ii is equal to the one we had for both NC-ST and NC-T. Also, the interaction between neighbouringpairs in successive layers reduces to Φ r ± || φ r ii = 12 || φ r ± ii , (S107)also analogous to || φ r ii hh φ r || φ r ± ii = 1 / C ( r, t ) is the same we obtainedfor the other two cases, with the restriction in the particle-hole state in (S105), when comparing to || φ r ii hh φ r || , beingcompensated by the different prefactor. XII. C-ST: AVERAGED OTOC, C ( r, t ) While the TOC and OTOC obtained previously for the NC-T and C-ST cases match the results obtained for NC-ST, the OTOC for C-ST does not. Next, we obtain precisely C ( r, t ) for the C-ST case, pointing out the differencesfrom its non-conserving counterpart. S.1. Time evolution of || Q ( t ) ii As done for the NC-ST case, we start by computing || Q ( t ) ii .M-19
1. Average over disorder realizations
Particle conserving free fermion gates vary from their non-conserving counterparts. Being built as mentioned atthe beginning of Sec. XI, their moments can be obtained directly from those of 2 × w . Thus, theirsecond moments h β | u | α i ∗ h β | u | α i are given by (S104) while its fourth moment is given by h β | u | α i ∗ h β | u | α i ∗ h β | u | α i h β | u | α i = h β | w | α i ∗ h β | w | α i ∗ h β | w | α i h β | w | α i δ s α s β δ s α s β δ s α s β δ s α s β , (S108)with the average of unitary matrices being (see Corollary 2.4 in [4]) h β | w | α i ∗ h β | w | α i ∗ h β | w | α i h β | w | α i = 1 N ( N − h δ β β δ β β ( N δ α α δ α α − δ α α δ α α ) + δ β β δ β β ( N δ α α δ α α − δ α α δ α α ) i . (S109)with N = rank( w ) = 2. With these changes, v r,r given previously by (S65) becomes v r,r → v Cr,r = 1 − δ rr N h || θ ii hh θ || ⊗ (Π + π ) + || ˜ θ ii hh ˜ θ || ⊗ (Π + ˜ π ) i r,r + δ rr N ( N − h N || θ ii hh θ || ⊗ (Π + π ) − || θ ii hh ˜ θ || ⊗ Π − || ˜ θ ii hh θ || ⊗ Π + N || ˜ θ ii hh ˜ θ || ⊗ (Π + ˜ π ) i r,r , (S110)where we used the definitions || θ r,r ii = X x ∈{ r,r +1 } X x ∈{ r ,r +1 } || x, x , x , x ii , (S111) || ˜ θ r,r ii = X x ∈{ r,r +1 } X x ∈{ r ,r +1 } || x, x , x, x ii , (S112)which respect hh i rr || i qq ii = 4 δ q,r δ q ,r + δ q,r ± δ q ,r ± , (S113) hh i rr || j qq ii = δ r,r (cid:16) δ q,r δ q ,r + δ q,r +1 δ q ,r +1 + δ q,r − δ q ,r − (cid:17) , (S114)with i = j ∈ { θ, ¯ θ, ˜ θ } . We also usedΠ = || pppp ii hh pppp || + || hhhh ii hh hhhh || , || Π ii = || pppp ii + || hhhh ii (S115) π = || phhp ii hh phhp || + || hpph ii hh hpph || , || π ii = || phhp ii + || hpph ii (S116)˜ π = || phph ii hh phph || + || hphp ii hh hphp || , || ˜ π ii = || phph ii + || hphp ii (S117)for which P P = δ P P and P || P ii = || P ii holds, with P and P ∈ { π, ¯ π, ˜ π } .
2. Dynamics of || Q ( t ) ii Defining || Θ ˜Θ r,r ii = (cid:16) || θ r,r ii + || ˜ θ r,r ii (cid:17) ⊗ || Π ii , (S118)we obtain from initial conditions || Q ( t = 1 / ii = v C − , − || Q (0) ii = 16 || Θ ˜Θ − , − ii − || ˜ θ − , − ii ⊗ || ˜ π ii , (S119)M-20to which we can apply a second even layer of gates to obtain || Q ( t = 1) ii = L/ − X R,R =0 v C R, R || Q (1 / ii = ( v C − , − + v C − , + v C , − + v C , ) || Q (1 / ii = 118 (cid:16) || Θ ˜Θ − , − ii + || Θ ˜Θ , ii (cid:17) + 124 (cid:16) || Θ ˜Θ − , ii + || Θ ˜Θ , − ii (cid:17) − (cid:16) || ˜ θ − , − ii + || ˜ θ , ii (cid:17) ⊗ || ˜ π ii − (cid:16) || ˜ θ − , ii + || ˜ θ , − ii (cid:17) ⊗ || ˜ π ii . (S120)One can check, making use of (S113) and (S114), that || Θ ˜Θ r,r ii and || ˜ θ r,r ii ⊗ || ˜ π ii evolve under v Cr,r as v Cr + a,r + a || Θ ˜Θ r,r ii = (cid:18)
14 + δ r,r (cid:19) || Θ ˜Θ r + a,r + a ii (S121) v Cr + a,r − a || Θ ˜Θ r,r ii = (cid:18) − δ r,r (cid:19) || Θ ˜Θ r + a,r − a ii (S122) v Cr ± a,r ± a || ˜ θ r,r ii ⊗ || ˜ π ii = (cid:18)
14 + δ r,r (cid:19) || ˜ θ r ± a,r ± a ii ⊗ || ˜ π ii (S123)with a = ±
1. These relations control the mixing occurring between neighbouring pairs in successive time steps andimply that both {||
Θ ˜Θ r,r ii} and {|| ˜ θ r,r ii ⊗ || ˜ π ii} are closed under time evolution, besides being orthogonal to oneanother (since hh Π || ˜ π ii = 0). Thus, we can decompose || Q ( t ) ii = L/ − X R =0 || Θ ˜Θ R, R ii hh Θ ˜Θ R, R || δ R, R ) + || ˜ θ R, R ii hh ˜ θ R, R || ⊗ || ˜ π ii hh ˜ π || ! || Q ( t ) ii , (S124)valid for t integer, becoming valid for t half-integer after transforming 2 R → R − R → R −
1. The prefactorscome from hh Θ ˜Θ r,r || Θ ˜Θ q,q ii = 8(2 + δ r,r ) δ r,q δ r ,q and hh ˜ θ r,r || ˜ θ q,q ii hh ˜ π || ˜ π ii = 8 δ r,q δ r ,q (where the indices arerestricted to having the same parity). Ahead, we will see that only states of the type || Θ ˜Θ r,r ii contribute to theOTOC and thus we can ignore {|| ˜ θ r,r ii ⊗ || ˜ π ii} . Although || Θ ˜Θ r,r ii are analogous to the states || Θ ¯Θ r,r ii appearingin the NC-ST case, the two are inherently different: while || Θ ˜Θ r,r ii includes, besides || αββα ii , states of the type || αβαβ ii which do not appear in NC-ST due to PH symmetry, it leaves out states of the kind || α ¯ αβ ¯ β ii , which do notappear when averaging over unitary gates.Having { δ r,r ) || Θ ˜Θ r,r ii} as the natural basis to use, we can employ (S121) and (S122) to obtain, for t ≥ / hh Θ ¯Θ r,r || → hh Θ ˜Θ r,r || and m r,r → m Cr,r , with m Cr,r coefficientmatrices given by m Cr,r = (cid:18)
13 1414 13 (cid:19) , m
Cr,r +2 = ( m Cr,r − ) T = (cid:18)
14 1416 14 (cid:19) , m
Cr,r = 14 (cid:18) (cid:19) , (S125)with the last m Cr,r being valid for | r − r | >
2. Applying this twice gives an expression equivalent to (S80) with hh Θ ¯Θ r,r || → hh Θ ˜Θ r,r || and M r,r → M Cr,r , with the new coefficient matrices M Cr,r being M Cr,r = 1144
16 21 921 44 219 21 16 , M Cr,r +2 = ( M Cr,r − ) T = 132
149 113
23 149 M Cr,r +4 = ( M Cr,r − ) T = 116 , M Cr,r = 116 , (S126)where the last M Cr,r is valid for | r − r | >
4. At last, starting with (S120) as our initial condition, || Q ( t ) ii isgiven by an analogue of the recursive expression (S80), with hh Θ ¯Θ r,r || → hh Θ ˜Θ r,r || and M r,r → M Cr,r , alongside(S124). Although the details surrounding both C-ST and NC-ST are different, the overall structure appearing in thedynamics of || Q ( t ) ii is similar. Namely, the dependence on the nearest and second nearest neighbours is the samealmost everywhere, i.e. M Cr,r = M r,r for | r − r | > S.2. Prescription to obtain C ( r, t ) Now that we know the dynamics of || Q ( t ) ii we can obtain that of the OTOC.For t = 0 (S82) holds. For t ≥ C ( r, t ) using (S124) to obtain C ( r, t ) = hh Q r || S || Q ( t ) ii = hh Q r || S L/ − X R,R =0 δ r,r ) || Θ ˜Θ R, R ii hh Θ ˜Θ R, R || Q ( t ) ii = 16 L/ − X R =0 ( δ r, R + δ r, R +1 ) hh Θ ˜Θ R, R || O ( t ) ii . (S127)Note that, as we already suggested, the term || ˜ θ r,r ii ⊗ || ˜ π ii coming from (S124) does not contribute to the OTOC,since the particle-hole sector of S || Q r ii is || Π + π ii and hh Π + π || ˜ π ii = 0.Defining K Cr,r ( t ) = 16 hh Θ ˜Θ r,r || Q ( t ) ii , (S128)analogous to (S84), the OTOC is then C (2 R, t ) = C (2 R + 1 , t ) = K C R, R ( t ). Finally, we can use (S120) to get theOTOC at t = 1, C (2 R, t = 1) = C (2 R + 1 , t = 1) = 29 (cid:16) δ R, − + δ R, (cid:17) , (S129)and, at subsequent time steps, the dynamics of C ( r, t ) is determined by that of hh Θ ˜Θ r,r || O ( t ) ii ⇔ K Cr,r ( t ) we justsaw to be given by (S80) with the changes hh Θ ¯Θ || → hh
Θ ˜Θ || and M → M C . This is, for t > K r,r ( t ) → K Cr,r ( t ) and M r,r → M Cr,r , given by (S118) and (S126), respectively.FIG. S5: The values in the boxes correspond to K Cr,r ( t ) in the r − r plane, shown here for t ∈ { / , , / , } (emptyboxes have zero value). Starting with K Cr,r ( t ) = 2 δ r,r , as the initial condition, we evolve this by applying (S78)successively (with hh Θ ¯Θ || → hh
Θ ˜Θ || and M → M C , given by (S118) and (S126), respectively). The highlighteddiagonal gives C ( r, t ) for the C-ST case. This picture is equivalent to Fig. S4, valid for NC-ST and NC-T. S.3. Continuum limit of C ( r, t ) Since the bulk evolution of K Cr,r ( t ) is equal to that of K r,r ( t ), i.e. m Cr,r | | r − r | > = m r,r | | r − r | > , the continuumlimit of the OTOC for the C-ST case is analogous to the one obtained for NC-ST. The differences of the diagonalwith regards to the NC-ST case translate to c → c C = 3 / , (S130) m r,r → m Cr,r = 16 (cid:18) (cid:19) , (S131) f → f C = 2 . (S132)M-22In this case, the anisotropy of m Cr,r leads to diffusion within the diagonal with weighting factor f C = 2.This results in the OTOC being given by (S91), with the changes c = √ / → c C = 3 / f = 1 / → f C = 2.Thus, in the continuum limit, the OTOC is a Gaussian with broadening width σ ( t ) = √ t , equal to the NC-ST andNC-T cases apart from the normalization which is A ( t ) → A C ( t ) = f C /f A ( t ) = f C A/ (2 √ πt ). XIII. CONTINUUM LIMIT OF C ( r, t ) Finally, we can take into account the contributions of C ( r, t ) and C ( r, t ) in the continuum limit, given respectivelyby (S54) and (S91), to obtain C (2 R, t ) = C (2 R + 1 , t ) = 12 (cid:20) A √ πt exp (cid:18) − (2 R ) t (cid:19) − f A πt exp (cid:18) − (2 R ) t (cid:19)(cid:21) , (S133)where A = 2 and f = 1 / f → f C = 2 for the C-ST case.Although C ( r, t ) is typically given by the sum of two Gaussians with different standard deviation, in the infinitetime limit it reduces to the first one, i.e. lim t →∞ C ( r, t ) = C ( r, t ) / since lim t →∞ C ( r, t ) /C ( r, t ) = 0. This occursbecause the normalization of C ( r, t ) is constant while that of C ( r, t ) is ∝ / √ t . XIV. NUMERICS
In the previous sections, we obtained analytical expressions for the TOC and OTOC for three instances of freefermion evolution: NC-ST, C-ST and NC-T. In the following, we compare these results with the respective approximateexpressions in the continuum limit and also with data obtained from simulations. We also realize a fourth instance offree fermion evolution numerically − a temporally homogeneous case where randomness appears only along the spacedirection (NC-S). While the analytical results take into account all possible evolutions of the system, i.e. all possiblecircuits, the simulations approximate this by an average over N r = 4000 disorder realizations. We look at a systemwith L = 100 using C ( r,
0) = C ( r,
0) = 2 δ r, as the initial condition. − . . . C ( r, t ) = σ ( t ) C ( r, t ) = 2 σ ( t ) r t C ( r , t ) (a) − . . .
06 NC-ST & NC-T r t C ( r , t ) (b) − . . . r t C ( r , t ) (c) FIG. S6: (a) Height and colour maps of the exact C ( r, t ) in the r − t plane for the NC-ST, NC-T and C-ST cases,given by (S52). (b) Exact C ( r, t ) for the NC-ST and NC-T cases, given by (S85) for t = 1 and (S86) for t >
1. (c)Exact C ( r, t ) for the C-ST case, given by (S129) for t = 1 and by an expression equivalent to (S86) for t >
1, withthe changes K → K C and M → M C given respectively by (S128) and (S126). These results are for a system with L = 100 and periodic boundary conditions. The results are shown for 1 ≤ t ≤ σ ( t )and 2 σ ( t ), with σ ( t ) = √ t for C and σ ( t ) = √ t for C .We begin with the analytical results, plotted in Fig. S6. We see that the profile of C ( r, t ) and C ( r, t ) at fixed timesteps is given by a Gaussian whose width broadens with time according to σ ( t ) = √ t and σ ( t ) = √ t , respectively.We see that C dominates over C . To compare these results with the data obtained from simulations we draw theprofile of both C ( r, t ) and C ( r, t ) at fixed time steps.Fig. S7 is a comparison of the analytic results and numerical simulation data for the three different cases NC-ST, NC-T and C-ST. The continuum limit in each case is tested by seeing whether C and C respect (S54) andM-23 . . − − − − − − . . . C ( r , t ) · σ ( t ) / A ( t ) ( r − µ ) /σ ( t ) C ( r , t ) · σ ( t ) / A NC-ST g ( r, t ) NC-T t = 450100 C-ST
FIG. S7: Rescaled TOC, C ( r, t ) · σ ( t ) /A , (top) and OTOC, C ( r, t ) · σ ( t ) /A ( t ), (bottom) as a function of ( r − µ ) /σ ( t )at three fixed time steps t ∈ { , , } for a system with L = 100 and for the three cases NC-ST (left), NC-T(centre) and C-ST (right). Full lines correspond to the exact calculations seen in Fig. S6 and the colored regions todata obtained from simulations averaging over N r = 4000 disorder realizations. Also, g ( r, t ) = √ π exp (cid:16) − ( r − µ ) σ ( t ) (cid:17) .For each case, the appropriate standard deviation σ ( t ) and normalization A or A ( t ) are used. For C ( r, t ), A = 2and σ ( t ) = √ t for all cases. For C ( r, t ), σ ( t ) = √ t for all cases, A ( t ) = 1 / (2 √ πt ) for NC-ST and NC-T and A ( t ) = 2 / ( √ πt ) for C-ST. The deviation µ = 1 / µ instead of 0.(S91), respectively. This translates to C ( r, t ) · σ ( t ) /A and C ( r, t ) · σ ( t ) /A ( t ) collapsing to the Gaussian g ( r, t ) = √ π exp (cid:16) − ( r − µ ) σ ( t ) (cid:17) , with A or A ( t ) and σ ( t ) given for each case in the caption of Fig. S7. We see in Fig. S7 that,for small time steps ( t = 4), the discrete data does not follow a Gaussian, as expected, while for higher time steps itdoes. Indeed, for a reasonable system size ( L = 100), the continuum limit quickly becomes a good approximation foreither the TOC or the OTOC. For time steps significantly larger than t = 100, periodic boundary conditions causethe discrete data to diverge from the continuum limit expectation.We have established that the analytical results, the data from simulations and the continuum limit expressions arein good mutual agreement. Before proceeding to analyse the numerical results obtained for the NC-S case, we providesome indication of the fluctuations present in both C ( r, t ) and C ( r, t ) by showing a single realization for each caseof FF evolution in Fig. S8. S.1. NC-S: randomness in space alone
We now present numerical results for NC-S for L = 100 and N r = 4000. The circuit for NC-S is built as shown inFig. S1 (c): we construct two random layers as usual (one even and one odd) and then apply them repeatedly, suchthat there is discrete time translation symmetry.We see in Fig. S9 that for the NC-S case both the TOC and OTOC are localized in space, decaying exponentiallyaround r = 0 (see panels (b) and (d)). This phenomenon, arising from the disordered landscape, is a discrete timeanalogue of Anderson localization. This picture is unlike the diffusive behaviour seen for quadratic fermion evolutionand the usual ballistic behaviour. Furthermore, this hints that the diffusive behaviour observed for NC-ST and NC-TM-24 −
50 0 50 −
50 0 50 −
50 0 50 −
50 0 50150100 00 . . . C t r . . . C t r . . . C t r . . . C t r . . . C t r . . . C t r . . . C t r . . . C t r FIG. S8: Single disorder realization of C ( r, t ) (top) and C ( r, t ) (bottom) as a colour map in the r − t plane forNC-ST (left), NC-T (centre left), CS-T (centre right) and NC-S (right). The black curves are r − µ = 2 σ ( t ), i.e. t = ( r − µ ) / C and t = ( r − µ ) / C . These results are for a system with L = 100 and periodic boundaryconditions. The deviation µ appears as it did in Fig. S7. − . . . rt C ( r , t ) (a) − . . . rt C ( r , t ) (b) . . . . − − − C ( r , t ) r t = 450100 (c) − − − − − − −
10 0 10 20 C ( r , t ) r t = 450100 (d) FIG. S9: (a) and (b) Height map of C ( r, t ) and C ( r, t ), respectively, in the r − t plane for 1 ≤ t ≤
50, in the NC-Scase. (c) and (d) Profile of C ( r, t ) and C ( r, t ) on a logarithmic scale, obtained by fixing t to 4, 50 and 100 in panels(a) and (b). The data results from an average over N r = 4000 simulations, for a system with L = 100 and periodicboundary conditions.(and also C-ST) is due to randomness in time alone. XV. CONCLUSIONS
We have obtained analytical results showing that the TOC and OTOC, C ( r, t ) and C ( r, t ), spread diffusivelyin the r − t plane for noisy free fermion evolution, with or without particle conservation. In total, three distinctcases were analysed − a particle non-conserving process with randomness in space and time (NC-ST), its particleconserving analogue (C-ST) and a noisy non-conserving case with spatial homogeneity (NC-T) − which differ eitherby the structure of the circuit or by the unitaries composing it (Fig. S1).For the NC-ST case, the C ( r, t ) dynamics is similar to a classical random walk in one spatial dimension: althougha single realization of C does not describe a classical random walk, the averaged C is equivalent to an average overall possible random walks respecting the geometry of the circuit. Specifically, starting from C ( r, t = 0) = 2 δ r, , eachweight C ( r, t ) is divided among its nearest neighbours C ( r − , t + 1 /
2) and C ( r + 1 , t + 1 /
2) (Section VII). Whilethe dynamics of C is a 1D diffusive process, C is given by the diagonal of a quantity, K r,r ( t ), which diffuses in 2D(Section VIII). The bulk behaviour of K r,r ( t ) can also be described by an average over all possible classical randomM-25walks, but in 2D, with each weight K r,r ( t ) splitting into its nearest neighbours K r ± ,r ± ( t + 1 / − instead of having spatial disorder, a gate is randomly chosen ateach half time step and applied to all pairs of sites (Fig. S1). In principle, this could lead the dynamics of C and C to arise from respectively 2D and 4D processes (Sections IX and X), but the initial condition reduces them to the1D and 2D processes already known for NC-ST. To obtain the particle conserving dynamics (C-ST), the free fermiongates used in the NC-ST case, built from orthogonal matrices, are replaced by gates built out of unitary matrices.Since the average of two orthogonal and two unitary gates needed to obtain C coincide, the TOC is equal to the oneobtained for the NC-ST case (Section XI). The picture is different for C − it is the only correlator which differs fromthe NC-ST results out of the cases studied. This happens because the average of four gates needed to compute C isdifferent for orthogonal and unitary matrices (Section XII). However, this does not change the fundamental diffusivenature of the process. It merely leads to a different normalization in the continuum limit.In the continuum limit, we saw in Section VII S.3 that C ( r, t ) is given by a Gaussian normalized to A = 2 whosestandard deviation broadens in time with σ ( t ) = √ t (for all cases) and C ( r, t ) is also given by a Gaussian, butwith standard deviation σ ( t ) = √ t and whose normalization decreases in time according to A ( t ) = f A/ (2 √ πt ), with f = 1 / f = 2 for C-ST (Section XII S.3). Thus, since A ( t ) ∝ / √ t , C is the leading term in the continuum limit, i.e. lim t →∞ C ( r, t ) /C ( r, t ) = 0.Finally, let us address the connection between the single and many-body correlators. We established a relationbetween the two representations in Section IV: the many-body TOC and OTOC, C ( r, t ) and C ( r, t ), are mixturesof the single-body TOC and OTOC, C ( r, t ) and C ( r, t ), and vice versa. For our particular observables, the many-body TOC is constant, C ( r, t ) = 1 / , and the non-trivial diffusive behaviour is encoded in the many-body OTOC: C ( r, t ) = [1 + 2 C ( r, t ) − C ( r, t )]. [1] A. Nahum, S. Vijay, and J. Haah, Physical Review X , 021014 (2018).[2] I. Peschel, Journal of Physics A: Mathematical and General , L205 (2003).[3] F. Mezzadri, arXiv preprint math-ph/0609050 (2006).[4] B. Collins and P. Śniady, Communications in Mathematical Physics264