Classical Models of Entanglement in Monitored Random Circuits
Oles Shtanko, Yaroslav A. Kharkov, Luis Pedro García-Pintos, Alexey V. Gorshkov
CClassical Models of Entanglement in Monitored Random Circuits
Oles Shtanko, Yaroslav A. Kharkov, Luis Pedro Garc´ıa-Pintos, and Alexey V. Gorshkov
Joint Center for Quantum Information and Computer Science and Joint Quantum Institute,NIST/University of Maryland, College Park, Maryland 20742, USA
The evolution of entanglement entropy in quantum circuits composed of Haar-random gates andprojective measurements shows versatile behavior, with connections to phase transitions and com-plexity theory. We reformulate the problem in terms of a classical Markov process for the dynamicsof bipartition purities and establish a probabilistic cellular-automaton algorithm to compute entan-glement entropy in monitored random circuits on arbitrary graphs. In one dimension, we furtherrelate the evolution of the entropy to a simple classical spin model that naturally generalizes atwo-dimensional lattice percolation problem. We also establish a Markov model for the evolutionof the zeroth R´enyi entropy and demonstrate that, in one dimension and in the limit of large localdimension, it coincides with the corresponding second-R´enyi-entropy model. Finally, we extend theMarkovian description to a more general setting that incorporates continuous-time dynamics, definedby stochastic Hamiltonians and weak local measurements continuously monitoring the system.
Recent progress in creating and manipulating many-qubit devices led to their ability to arguably attain quan-tum supremacy, i.e. to exhibit dynamics that is not effi-ciently simulable on classical computers [1–3]. Existingsupremacy proposals include circuit models that utilizerandomly generated local gates to create classically irre-producible quantum multi-qubit correlations [4, 5]. Ran-dom circuits can be characterized by a rapid growth ofentanglement across the system [6–9], and even thoughhigh entanglement alone does not guarantee hardness ofclassical simulation [10, 11], the generation of entangledstates is key to achieving such hardness in many cases[12]. Different behavior can be found if a random cir-cuit is affected by measurements that project and thusdisentangle individual qudits. Under persistent randommeasurements, a competition between entanglement pro-duction and reduction mechanisms leads to a phase tran-sition [13–25] separating phases with area- and volume-law entanglement, similar to known dynamical quantumphase transitions in Hamiltonian systems [26–29]. There-fore, the evolution of entanglement in monitored ran-dom circuits exhibits versatile behavior associated witha range of fields from complexity theory to condensedmatter physics.In this work, we study the dynamics of entanglementin a wide class of random circuits under the combined ef-fect of unitary gates and measurements. Starting witha random-walk model for entanglement growth in 1Dunitary circuits, we generalize our approach to a MonteCarlo algorithm for computing the average entanglemententropy on arbitrary graphs in the presence of both uni-tary evolution and single-qudit measurements. Further-more, our approach can be mapped to a dynamical lat-tice percolation problem, which, in the limit of infinitequdit dimension, maps to a static percolation transition[15, 17, 20]. We also study a continuous-time analogue ofrandom circuits [30–33] in which the system evolves un-der a stochastic local Hamiltonian and is subject to de-phasing and continuous weak measurements [34–36]. We thus introduce the study of random quantum circuits toa rich arena of weak-measurement physics boasting bothapplied and foundational results [37–40].Mapping of monitored quantum circuits to classicalmodels has been studied recently using the replica trickfor both closed [9, 41] and monitored [23] systems. Thetheoretical description is simplified in the limit of largequdit dimension q (cid:29) q limit, the evolution of R´enyientropies in random unitary circuits can be described bya Kardar–Parisi–Zhang-type (KPZ-type) equation [43].Hydrodynamic descriptions of random unitary circuitsare also discussed in Refs. [8, 44, 45]. In the presence ofprojective measurements, mapping to spin models con-nects the entanglement phase transitions to the percola-tion criticality and conformal field theories [15, 46]. Themeasurement-induced entanglement phase transition isalso related to the encoder-decoder problem, with mea-surements playing the role of external noise corruptingthe information [20, 23]. Compared to previous ap-proaches designed for 1D systems, our method appliesto arbitrary graphs, with significant speedups over exactsimulation.We consider a generic circuit acting on q -dimensionalqudits and consisting of Haar-random gates and projec-tive single-qudit measurements. As a measure of entan-glement between a specified subsystem A and the restof the qudits, we use average R´enyi entropy S [ A ] ≡−(cid:104) log Tr A (cid:0) ρ A (cid:1) (cid:105) , where (cid:104) . . . (cid:105) denotes an average with re-spect to circuit realizations and ρ A is the reduced densitymatrix of A . In order to compute the average, we ap-ply the annealed approximation , S [ A ] (cid:39) − log P A , where P A ≡ (cid:10) Tr A (cid:0) ρ A (cid:1)(cid:11) is the (average) purity. This approx-imation is exact for Haar-random states in the limit oflarge system size [47–49]. We demonstrate in this workthat this approximation is also remarkably accurate inapplication to random circuits.To compute the R´enyi entropy, we focus on the vector P = { P G : G ∈ F } consisting of purities of all possible a r X i v : . [ c ond - m a t . d i s - nn ] A p r subsystems G of the full system F . Below we mostlyfocus on the linear discrete-time Markovian evolution P ( t + 1) = L P ( t ) , (1)where L is a time-dependent Liouvillian (see also Ref.[50]). In the following sections, we will show the relevanceof this model to random circuits. Random Circuits.—
Consider the evolution of a purestate | ψ (cid:105) under a single Haar-random gate U supportedon subset ω , | ψ (cid:48) (cid:105) = U | ψ (cid:105) . The resulting (gate-averaged)purity of an arbitrary subset G for the state | ψ (cid:48) (cid:105) is deter-mined [see Supplemental Material (SM) [51], Sec. I] bythe purities for | ψ (cid:105) via P (cid:48) G = c − P G \ ω + c + P G ∪ ω , (2)where c − = d ( d − / ( d ω − c + = d ( d − / ( d ω − d = d ω ∩ G , and d = d ω \ G [52]. For two-qudit gates with d = d = q , we obtain c − = c + = q/ ( q + 1). Coeffi-cients c + and c − are the only non-zero elements in theLiouvillian matrix L in Eq. (1) representing a single ran-dom gate. For more than one gate, the total Liouvillianis the product of individual gate Liouvillians in reversechronological order.We now consider a projective measurement performedon a single qudit Ω of an observable O = (cid:80) i o i Π i ,where Π i are projectors satisfying Π i Π j = Π i δ ij . Ameasurement with outcome o i projects the system onto | ψ (cid:48) (cid:105) = Π i | ψ (cid:105) / (cid:112) (cid:104) ψ | Π i | ψ (cid:105) . Since measurements are al-ways “sandwiched” between Haar-random gates, withoutloss of generality, we replace the purities by their averageover the basis of O , P G = (cid:104) Tr ρ G (cid:105) O . Applying the an-nealed approximation to O /O ( (cid:104) O /O (cid:105) (cid:39) (cid:104) O (cid:105) / (cid:104) O (cid:105) ),the post-measurement purity becomes P (cid:48) G (cid:39) (cid:10) Tr G (Π i | ψ (cid:105)(cid:104) ψ | Π i ) G (cid:11) O (cid:104)(cid:104) ψ | Π i | ψ (cid:105) (cid:105) O = P G ∪ Ω + P G \ Ω P Ω . (3)To express the combined non-linear evolution underEqs. (2,3) in terms of linear evolution, we define unnor-malized purities ˜ P evolving under ˜ P ( t + 1) = L ˜ P ( t ),where unitary dynamics is given by Eq. (2), whilemeasurements are treated as the transformation ˜ P (cid:48) G =˜ P G \ Ω + ˜ P G ∪ Ω . Then, after an arbitrary number of uni-tary and measurement transformations, the physical pu-rity vector is P G = ˜ P G / ˜ P F , (4)where ˜ P F is the unnormalized purity of the full system. Classical Models.—
The linearity of purity evolu-tion enables the mapping of entanglement dynamics ontoclassical problems such as multi-particle random walks.As an illustration, consider a simple brickwork 1D uni-tary circuit [Fig. 1(a) without measurements]. It fol-lows from Eq. (2) that the dynamics of bipartition puri-ties P tx parametrized by the cut position x is a modified single-particle random-walk: P t +1 x = c (cid:0) P tx − + P tx +1 (cid:1) for x = 2 n + t (mod 2), n ∈ Z , and c = q/ ( q + 1). The R´enyientropies S ( t, x ) = − log P tx in the continuous limit thensatisfy a noise-averaged KPZ equation ∂ t S = µ (cid:0) ∂ x S − ( ∂ x S ) (cid:1) + β, (5)where β = log(( q + 1) / q ) and µ = β/ (log q ) (see SMSec. III [51]). The long-time asymptotic solution of thishydrodynamic equation, subject to vanishing boundaryconditions at the endpoints of the 1D qudit chain, is S ( t → ∞ , x ) → log (cid:20) cosh ( αN/ α ( x − N/ (cid:21) , (6)where α = (cid:112) β/µ = log q and N is the number of qu-dits. The maximum value of the entropy at the centerof the chain, S ( t → ∞ , x = N/ (cid:39) αN/ − log 2 + O (exp( − αN )), exhibits volume-law scaling. The preci-sion of the hydrodynamic approximation for q = 2 isillustrated in Fig. 2(a).In more general settings, the number of partitions in-volved in the dynamics of purities P A (or ˜ P A ) typicallygrows exponentially in time. To describe such processes,we consider a sequence of n circuit layers, each describedby a sparse Liouvillian L k , sorted in reverse chronologi-cal order. Let us define a collection of neighboring sets N k [ G ] = { G (cid:48) : L kG,G (cid:48) (cid:54) = 0 } and multiple trajectories G = { G → G · · · → G n } such that G k +1 ∈ N k [ G k ],and G = A . Assuming the system is initialized in aproduct state, i.e. P G (0) = 1 for all G , the (unnormal-ized) purity after n circuit layers is˜ P A = (cid:88) G P G , P G = n − (cid:89) k =0 L kG k ,G k +1 . (7)The sum in Eq. (7) can be represented by a gener-alized random walk. Following a Markov chain MonteCarlo (MCMC) algorithm, instead of listing all possi-ble trajectories, we consider their probabilistic genera-tion, choosing G k +1 from the set N k [ G k ] with probabil-ity π kG = L kG,G k /T k , where T k = (cid:80) G ∈ N k [ G k ] L kG,G k . Inthe MCMC scheme, the purity ˜ P A in Eq. (7) is approxi-mated by the sample average ˜ P A ≈ (cid:104) P MC G (cid:105) G = (cid:104) (cid:81) k T k (cid:105) G .In practice, in order to estimate ˜ P A , one needs an averageover a large number of trajectories due to the broad dis-tribution of P MC G . Nevertheless, simulations show thatthis number can be significantly smaller than the fullnumber of trajectories G .For simulation purposes, the MCMC algorithm can bemapped to a probabilistic classical cellular automaton[53]. It is instructive to start with the monitored 1Dbrickwork circuit shown in Fig. 1(a). We describe eachunitary and measurement layer by a Liouvillian L u k and L m k , respectively. The colored pattern in Fig. 1(b) il-lustrates the dynamics of the corresponding cellular au-tomaton evolving from top to bottom generating a tra-jectory, where each cell represents a qudit. Blue cells in unitary measurement (b) A B s i m u l a t i o n t i m e G G G ... (c) A B i jk Σ = G G G ...qudit positionqudit position prohibited: Z i =1 Z i =-1 J ij =-1 J ij = 1 e v o l u t i o n t i m e (a) e v o l u t i o n t i m e ( ) ( ) FIG. 1.
1D brickwork circuits . (a) Monitored 1D brickwork circuit. Rectangles represent two-qudit Haar-random unitaries,while black circles represent projective measurements. (b) A trajectory generated by a probabilistic cellular automaton for themonitored 1D random circuit in (a). The process starts from the top (final time) and flows down following local rules describedon the right for unitaries (ovals) and measurements (dots). The 20 simulation steps shown (from G to G ) correspond to the10 brickwork layers interleaved with 10 measurement layers shown in (a). The resulting trajectory can be converted to a singlecontribution P MC G to the sum in Eq. (7). (c) One configuration, corresponding precisely to the trajectory in (b), of a classicalspin model whose partition function defines the final purity via Eq. (8). The Hamiltonian is the 2D Ising model in Eq. (9) forclassical spins Z i associated with squares ( Z i = 1 are blue, Z i = − J ij : J ij = − J ij = 1(solid) otherwise. The operator P acts non-trivially only on the top row of cells, projecting them onto the target bipartition A , and on all directed 3-sets Σ prohibiting certain configurations of edges (shown on the right) irrespective of square colors. each row denote the qudits forming G k ∈ G , where G isa trajectory formed by alternating L u k and L m k .Each unitary layer Liouvillian L u k can be converted intoa rule: all distinct-color pairs of cells affected by a Haar-random two-qudit gate in this layer switch to the samecolor, chosen by a coin toss for each gate. Every suchswitch contributes a factor T k = 2 q/ ( q + 1) to the cor-responding amplitude P MC G . A measurement layer L m k isequivalent to a color flip of each cell representing a mea-sured qudit with probability π kG = 1 /
2, with no contribu-tion to P MC G , i.e. T k = 1. The method can be used as afast sampling algorithm producing purities [see Fig. 2(b)]and applies to arbitrary graphs [see Fig. 2(c) showingMCMC results for 2D]. The MCMC method is particu-larly useful for analysing measurement-induced entangle-ment phase transitions. For example, Fig. 2(b) illustratesentropy growth in a 1D brickwork circuit with measure-ments affecting each qudit with probability p at everylayer; MCMC shows excellent agreement with exact nu-merics. The algorithm can be used to reach larger systemsizes enabling an estimate of the critical p and critical ex-ponents, as illustrated in Fig. 2(d).Instead of the MCMC algorithm, we can alternativelyrelate ˜ P A to a partition function of a classical Hamil-tonian H , with trajectories G mapped to energy levels E [ G ] ∈ spec( H ),˜ P A = (cid:88) G exp( − βE [ G ]) = Tr (cid:16) P exp( − βH ) (cid:17) , (8)where β = log (cid:0) ( q + 1) / q (cid:1) is an effective inverse tem-perature, and P is a projector onto the levels E [ G ]. For the monitored 1D brickwork circuit, such a classi-cal model can be constructed as follows (see also SM Sec.IV [51]). Let us combine each pair of unitary and mea-surement layers into a single Liouvillian L um k coupling G (cid:48) k to G (cid:48) k +1 , where the full trajectory G = { G (cid:48) → G (cid:48) · · · → G (cid:48) n } includes only odd-time bipartitions, G (cid:48) k = G k − .Then, for any trajectory G , bipartitions propagate on asquare lattice shown in Fig. 1(c), representing a classicalanalogue of the Feynman checkerboard [54]. We associatea classical variable Z i = 1 to a square if it is part of G (cid:48) k [blue squares in Fig. 1(c)], while the rest have Z i = − i, j ) separatingadjacent squares i and j and situated below a measure-ment in Fig. 1(c) is assigned an additional degree of free-dom J ij = 1 (solid line) or J ij = − J ij = 1. In this setting, the classicalHamiltonian and projector in Eq. (8) are H = 12 N − (cid:88) (cid:104) i,j (cid:105) J ij Z i Z j , P = P A (cid:89) i,j,k ∈ Σ (cid:16) −
14 ( Z k − J ik Z i )( Z k − J jk Z j ) (cid:17) , (9)where (cid:104) i, j (cid:105) is the sum over all N edges, P A = (cid:81) i ∈ A (1+ Z i ) (cid:81) j ∈ F \ A (1 − Z j ) projects the top squaresonto the bipartition configuration A , and Σ are all pos-sible oriented 3-cell groups, as shown in Fig. 1(c). Thismodel applies to periodic boundary conditions (for openboundary conditions, see SM Sec. IV [51]). The Hamil-tonian H is simply equal to the total number of solid rededges plus the total number of dashed black edges. Connection to Hartley entropy.—
The dynamicsof the Hartley entropy, defined as the zeroth R´enyi en-tropy S [ A ] = log( R A ), where R G = rank ( ρ G ), is typi-cally qualitatively different from R´enyi entropies of order n ≥ S has deep con-nections to the 2nd R´enyi entropy S . Under a genericunitary transformation supported on ω , G ∩ ω (cid:54) = ∅ and ω \ G (cid:54) = ∅ , the rank follows non-liear Markovian dynamics(see SM Sec. II [51]) R (cid:48) G = R G \ ω R G ∪ ω R G ∆ ω min( d G ∩ ω R G ∩ ω , d ω \ G R ω \ G ) , (10)where G ∆ ω ≡ ( G ∪ ω ) \ ( G ∩ ω ). The effect of a mea-surement on qudit Ω is R (cid:48) G = R G \ Ω R G ∪ Ω R G ∆Ω min( R G ∩ Ω , R Ω \ G ) . (11)In 1D systems, these expressions become linear: R (cid:48) G =min (cid:0) d G ∩ ω R G \ ω , d ω \ G R G ∪ ω (cid:1) for unitary gates [7] and R (cid:48) G = min ( R G \ Ω , R G ∪ Ω ) for measurements. It followsfrom these expressions that finding S for a 1D sys-tem initialized in a product state can be reduced toa minimization of P G in Eq. (7), P A = min P G . For1D brickwork circuits, this is equivalent to path min-imization on a percolated lattice, as previously sug-gested in Ref. [17]. Notably, S can also be written as S [ A ] = log( q ) min G E [ G ], where E [ G ] is the spectrumof the Hamiltonian in Eq. (9) for configurations satisfy-ing P [ G ] = 1. As a consequence, in the limit q → ∞ ,entropies S and S coincide, given that, in the limit β → ∞ , the partition function in Eq. (8) reduces to S [ A ] (cid:39) β min G E [ G ] with β ∼ log q .In the hydrodynamic approximation, Hartley entropyfor 1D unitary brickwork circuits evolves according to thecontinuum limit of Eq. (10) (see SM Sec. III [51]): ∂ t S = 12 ∂ x S − | ∂ x S | + α, (12)where α = log q . For product initial states, the time-dependent solutions of Eqs. (12) and (5) coincide in thelimit q → ∞ at scales x (cid:29)
1. Furthermore, the station-ary solution of Eq. (12) for an initial product state andopen boundary conditions [ S ( x = 0) = S ( x = N ) = 0]is S ( t → ∞ , x ) = α ( N/ − | x − N/ | ), exhibiting thesame volume-law entanglement as S . Continuous Evolution.—
To complete the analysisof the evolution of purity, we propose a continuous-timeversion of monitored random circuits. We consider a localstochastic Hamiltonian H ( t ) = (cid:80) i I F \ ω i ⊗ h i ( t ), whereindividual local terms h i ( t ) are stochastic Gaussian-unitary-ensemble matrices supported on respective sets ω i . The corresponding matrix elements h i,νµ of h i sat-isfy (cid:104) h ∗ i,µν ( t ) h j,µ (cid:48) ν (cid:48) ( t (cid:48) ) (cid:105) h = α i ( t ) d − ω i δ ij δ µµ (cid:48) δ νν (cid:48) δ ( t − t (cid:48) )for some positive functions α i ( t ). As a model for monitor-ing, we consider continuous weak measurements [55–58] E n t r o p y S p =0 p =0.01 p =0.05 p =0.1 T i m e ( l a y e r s ) E n t r o p y S (a) (b) (d)(c) M ea s . r a t e p E n t r o p y S P o s i t i o n x E n t r o p y S
1D 1D2D
N=50N=40N=30N=20 x y AB xy FIG. 2.
Classical simulation ( q = 2 ) . (a) Comparison be-tween the solution of Eq. (5) taken at times t = 2 n , n ∈ Z ,(solid curves) and exact numerical simulation for a 1D circuit(dots) and even cut positions x , averaged over 10 realiza-tions; system size is N = 20. (b) Entropy growth in a 1Dbrickwork circuit with measurements. The figure comparesexact numerics averaged over 10 gate-measurement realiza-tions (dots) to a MCMC simulation using 10 samples foreach of the 10 sampled measurement configurations (dashedcurves) for system size N = 14. (c) Long-time ( t = 64) valueof the entropy for a monitored circuit acting on a 2D systemof 16 ×
16 qudits with the brickwork circuit arrangement intro-duced in Ref. [42] and with measurement probability p = 0 . samples for each of the 150 sampled measurementconfigurations. (d) Scaling of the entropy for a monitored1D system as a function of the measurement rate p , as calcu-lated using MCMC with up to 10 samples. The inset showsthat all four curves collapse to a single curve when plottedas a function of ( p − p c ) N /ν for critical point p c (cid:39) . ν (cid:39) . performed for arbitrary single-qudit physical observables O j (acting on site Ω j ) that yield a combined non-linearequation ddt | ψ (cid:105) = − iH ( t ) | ψ (cid:105) (13)+ (cid:88) j (cid:20) − κ j [ δO j ( ψ )] + ξ j ( t ) (cid:112) κ j δO j ( ψ ) (cid:21) | ψ (cid:105) , where δO j ( ψ ) = O j − (cid:104) ψ | O j | ψ (cid:105) , ξ i ( t ) are independentstochastic variables satisfying (cid:104) ξ i ( t ) ξ j ( t (cid:48) ) (cid:105) ξ = δ ij δ ( t − t (cid:48) ),and κ j characterize the strength of the coupling to themeasurement apparatus [34–36]. Similar to projec-tive measurements, continuous measurements also reduceentanglement in the system, leading to a competitionwith unitary evolution [13]. We note that what followsholds even if we replace our time-independent O j with O j ( t ) = V j ( t ) O j V † j ( t ), where V j ( t ) is an arbitrary time-dependent unitary, which changes the measurement ba-sis.We are looking for a closed set of equations for thepurity vector P , which evolves according to ∂ P /∂t = L u t ( P ) + L m t ( P ). Here, the unitary Liouvillian obeys (seeSM Sec. V [51])[ L u t ( P )] G = − (cid:88) i α i (cid:16) P G + P G ∆ ω i d ω − P G \ ω i d ω i ∩ G − P G ∪ ω i d ω i \ G (cid:17) . (14)As expected, for a single time-independent term h i , thesteady-state of the Liouvillian in Eq. (14) yields thediscrete-time evolution in Eq. (2) [51].The measurement Liouvillian obeys[ L m t ( P )] G = (cid:88) j :Ω j ∈ F \ G ζ j C G, Ω j + (cid:88) j :Ω j ∈ G ζ j C F \ G, Ω j , (15)where ζ j = 8 κ j Tr Ω j O j / ( q − C S, Ω = (cid:68)(cid:13)(cid:13) ρ S ∪ Ω − ρ S ⊗ ρ Ω (cid:13)(cid:13) (cid:69) encodes the stateof the system, and (cid:107) · (cid:107) is the 2-norm. In the mean-field approximation, C S, Ω (cid:39) P S ∪ Ω − P S P Ω , leading tothe desired closed set of equations for purities. Impor-tantly, in this approximation, when only one qudit ismeasured, the steady-state of the Liouvillian in Eq. (15)yields the discrete-time annealed-approximation evolu-tion in Eq. (3) [51]. Outlook.—
In this work, we demonstrated that lineardynamics of purities arises in monitored random circuitsand yields, under the annealed approximation, an algo-rithm for computing entanglement based on a mapping todynamical percolated lattices. Many classical spin mod-els with local interactions have efficient classical solutionsrequiring only polynomial resources. Can we establishthis property for dynamical percolated lattices exploredin this work? If so, this work opens a pathway to estab-lishing the approximate classical simulability of entan-glement dynamics in monitored random circuits. Also,having established in this work a close correspondencebetween zeroth and second R´enyi entropies, we conjec-ture such correspondence for other R´enyi entropies.
Note.—
During the preparation of this manuscript, arelated preprint appeared [59] studying purities to ana-lyze entanglement phase transitions in 1D circuits usingmappings to a quantum spin model.
Acknowledgements . We thank Michael Gullans,Abhinav Deshpande, Pradeep Niroula, Zhicheng Yang,and Soonwon Choi for fruitful discussions. This workwas supported by the NSF PFCQC program, DoEASCR FAR-QC (award No. DE-SC0020312), DoE BESMaterials and Chemical Sciences Research for Quan-tum Information Science program (award No. DE- SC0019449), DoE ASCR Quantum Testbed Pathfinderprogram (award No. DE-SC0019040), AFOSR, AROMURI, ARL CDQI, and NSF PFC at JQI. [1] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin,R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao,D. A. Buell, et al. , Nature , 505 (2019).[2] A. W. Harrow and A. Montanaro, Nature , 203(2017).[3] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush,N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, andH. Neven, Nat. Phys. , 595 (2018).[4] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani,arXiv:1803.04402 (2018).[5] R. Movassagh, arXiv:1909.06210 (2019).[6] O. C. O. Dahlsten, R. Oliveira, and M. B. Plenio, J.Phys. A – Math. Theor. , 8081 (2007).[7] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys.Rev. X , 031016 (2017).[8] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, andS. L. Sondhi, Phys. Rev. X , 021013 (2018).[9] T. Zhou and A. Nahum, Phys. Rev. B , 174205 (2019).[10] B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A ,032325 (2002).[11] S. Aaronson and D. Gottesman, Phys. Rev. A , 052328(2004).[12] G. Vidal, Phys. Rev. Lett. , 040502 (2004).[13] M. Szyniszewski, A. Romito, and H. Schomerus, Phys.Rev. B , 064204 (2019).[14] Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B ,134306 (2019).[15] A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrish-nan, D. A. Huse, and J. H. Pixley, Phys. Rev. B ,060301 (2020).[16] Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B ,205136 (2018).[17] B. Skinner, J. Ruhman, and A. Nahum, Phys. Rev. X , 031009 (2019).[18] Y. Bao, S. Choi, and E. Altman, Phys. Rev. B ,104301 (2020).[19] X. Cao, A. Tilloy, and A. D. Luca, SciPost Phys. , 24(2019).[20] M. J. Gullans and D. A. Huse, arXiv:1905.05195 (2019).[21] M. J. Gullans and D. A. Huse, arXiv:1910.00020 (2019).[22] Q. Tang and W. Zhu, Phys. Rev. Research , 013022(2020).[23] S. Choi, Y. Bao, X.-L. Qi, and E. Altman,arXiv:1903.05124 (2019).[24] L. Zhang, J. A. Reyes, S. Kourtis, C. Chamon, E. R.Mucciolo, and A. E. Ruckenstein, arXiv:2001.11428(2020).[25] J. Lopez-Piqueres, B. Ware, and R. Vasseur,arXiv:2003.01138 (2020).[26] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,Adv. Phys. , 239 (2016).[27] R. Vasseur, A. C. Potter, Y.-Z. You, and A. W. W.Ludwig, Phys. Rev. B , 134203 (2019).[28] R. Nandkishore and D. A. Huse, Annu. Rev. Conden.Ma. P. , 15 (2015).[29] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. , 021001 (2019).[30] T. Zhou and X. Chen, Phys. Rev. E , 052212 (2019).[31] S. Xu and B. Swingle, Phys. Rev. X , 031048 (2019).[32] X. Chen and T. Zhou, Phys. Rev. B , 064305 (2019).[33] T. Zhou, S. Xu, X. Chen, A. Guo, and B. Swingle,arXiv:1909.08646 (2019).[34] H. M. Wiseman and G. J. Milburn, Quantum Measure-ment and Control (Cambridge University Press, 2009).[35] K. Jacobs,
Quantum Measurement Theory and its Appli-cations (Cambridge University Press, 2014).[36] K. Jacobs and D. A. Steck, Contemp. Phys. , 279–303(2006).[37] S. Weber, A. Chantasri, J. Dressel, A. N. Jordan,K. Murch, and I. Siddiqi, Nature , 570 (2014).[38] K. Murch, S. Weber, C. Macklin, and I. Siddiqi, Nature , 211 (2013).[39] M. H. Devoret and R. J. Schoelkopf, Science , 1169(2013).[40] P. Campagne-Ibarcq, P. Six, L. Bretheau, A. Sarlette,M. Mirrahimi, P. Rouchon, and B. Huard, Phys. Rev. X , 011002 (2016).[41] J. Napp, R. L. La Placa, A. M. Dalzell, F. G. Brandao,and A. W. Harrow, arXiv:2001.00021 (2019).[42] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X , 021014(2018).[43] M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. , 889 (1986).[44] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev.X , 031057 (2018).[45] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk,Phys. Rev. X , 031058 (2018).[46] C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. W. Lud-wig, Phys. Rev. B , 104302 (2020).[47] C. Nadal, S. N. Majumdar, and M. Vergassola, Phys.Rev. Lett. , 110501 (2010).[48] C. Nadal, S. N. Majumdar, and M. Vergassola, J. Stat.Phys. , 403 (2011).[49] S. N. Majumdar and G. Schehr, J. Stat. Mech.-Theory.E. , P01012 (2014).[50] W.-T. Kuo, A. Akhtar, D. P. Arovas, and Y.-Z. You,arXiv:1910.11351 (2019).[51] See Supplemental Material at X for details omitted in themain text.[52] We refer to qudit sets using standard notation for therelative complement ( G \ G (cid:48) ), union ( G ∪ G (cid:48) ), intersection( G ∩ G (cid:48) ), as well as Hilbert space dimension d G ≡ q | G | .[53] R. Fern´andez, P.-Y. Louis, and F. R. Nardi, “Overview:Pca models and issues,” in Probabilistic Cellular Au-tomata: Theory, Applications and Future Perspectives ,edited by P.-Y. Louis and F. R. Nardi (Springer Interna-tional Publishing, Cham, 2018) pp. 1–30.[54] R. Feynman, A. Hibbs, and D. Styer,
Quantum Mechan-ics and Path Integrals , Dover Books on Physics (DoverPublications, 2010).[55] C. M. Caves, Phys. Rev. D , 1643 (1986).[56] L. Di´osi, Phys. Lett. A , 419 (1988).[57] A. N. Korotkov, Phys. Rev. B , 115403 (2001).[58] H. M. Wiseman, Quantum Semicl. Opt. , 205 (1996).[59] R. Fan, S. Vijay, A. Vishwanath, and Y.-Z. You,arXiv:2002.12385 (2020). Supplemental Material for“Classical Models of Entanglement in Monitored Random Circuits”
Oles Shtanko, Yaroslav A. Kharkov, Luis Pedro Garc´ıa-Pintos, and Alexey V. GorshkovIn this Supplemental Material, we present details omitted in the main text. In particular, in Sec. I, we deriveEqs. (2,3) in the main text and analyze the accuracy of the annealed approximation. In Sec. II, we derive Eqs. (10,11)in the main text. In Sec. III, we derive Eqs. (5,6,12) in the main text. In Sec. IV, we derive Eqs. (8,9) in the maintext. Finally, in Sec. V, we derive Eqs. (14,15) in the main text, analyze the accuracy of the mean-field approximationused to simplify Eq. (15), and show that Eqs. (14,15) can be used to derive Eqs. (2,3).
SECTION I: EVOLUTION OF R´ENYI ENTROPIES
In this section, we derive the evolution Equations (2,3) for averaged subsystem purities. We also analyze theaccuracy of the annealed approximation used in the main text to connect R´enyi entropies and averaged purities.
Unitary gates . Consider a single unitary gate U supported on subset ω . We are interested in finding theeffect of U on averaged purities. To perform this calculation, we divide the system into four subsets A , B , C ,and D , as shown in Fig. S1(a), and introduce the corresponding matrix representation for the density operator ρ = (cid:80) ind ρ αβγδ,α (cid:48) β (cid:48) γ (cid:48) δ (cid:48) | αβγδ (cid:105)(cid:104) α (cid:48) β (cid:48) γ (cid:48) δ (cid:48) | and for the unitary gate U = (cid:80) γδ,γ (cid:48) δ (cid:48) U γδ γ (cid:48) δ (cid:48) | γδ (cid:105)(cid:104) γ (cid:48) δ (cid:48) | . Here labels { α, β, γ, δ } refer to the basis states of { A, B, C, D } , respectively. Then P (cid:48) G = (cid:68) Tr G (cid:0) Tr F \ G ( U ρU † ) (cid:1) (cid:69) circ ,U = (cid:88) α k ,β k ,γ k ,δ k (cid:68) ρ α β γ δ ,α β γ δ ρ α β γ δ ; α β γ δ (cid:69) circ (cid:68) U γ δ ,γ δ U ∗ γ δ ,γ δ U γ δ ,γ δ U ∗ γ δ ,γ δ (cid:69) U , (S.1)where (cid:104) . . . (cid:105) U is the Haar average over U , while the average over the circuit (cid:104) . . . (cid:105) circ includes the average over circuitevolution preceding the application of U . For random circuits, the (cid:104) . . . (cid:105) U and (cid:104) . . . (cid:105) circ averaging can be separatedbecause gate unitaries are sampled independently.Correlation functions of Haar-random unitary matrix elements satisfy (cid:104) U a,b U ∗ c,d U a (cid:48) ,b (cid:48) U ∗ c (cid:48) ,d (cid:48) (cid:105) U = 1 d ω − (cid:16) δ ac δ a (cid:48) c (cid:48) δ bd δ b (cid:48) d (cid:48) + δ ac (cid:48) δ a (cid:48) c δ bd (cid:48) δ b (cid:48) d − d ω ( δ ac δ a (cid:48) c (cid:48) δ bd (cid:48) δ b (cid:48) d + δ ac (cid:48) δ a (cid:48) c δ bd δ b (cid:48) d (cid:48) ) (cid:17) . (S.2)Denoting d ≡ d G ∩ ω and d ≡ d ω \ G , satisfying d ω = d d , we obtain P (cid:48) G = 1 d ω − (cid:16) d d P G \ ω + d d P G ∩ ω − d ω (cid:16) d d P G \ ω + d d P G ∩ ω (cid:17)(cid:17) = 1 d ω − (cid:16) d ( d − P G \ ω + d ( d − P G ∪ ω (cid:17) , (S.3)which is Eq. (2) in the main text. This expression works for both pure and mixed states of the full system. Measurements . Consider a projective measurement of qudit Ω. We express this measurement using a projectorΠ = U Π U † , where Π = Π , Tr Π = 1, and U is a Haar-random unitary supported on Ω. We assume that the fullsystem is in a pure state, ρ = | ψ (cid:105)(cid:104) ψ | . Using the annealed approximation to take the expectation value of a ratio, weobtain P (cid:48) G ≈ (cid:10) Tr G (Π | ψ (cid:105)(cid:104) ψ | Π) G (cid:11) Π (cid:104)(cid:104) ψ | Π | ψ (cid:105) (cid:105) Π = ˜ P (cid:48) G ˜ P (cid:48) F . (S.4)Let us assume first that Ω ∈ G . Then˜ P (cid:48) G ≡ (cid:68) Tr G ( ρ G Π) (cid:69) ρ, Π = (cid:68) Tr (cid:16) ρU Π U † ρU Π U † (cid:17)(cid:69) ρ,U == 1 d − (cid:68)(cid:16) Tr G ρ G (Tr Ω Π ) + Tr G \ Ω (Tr Ω ρ G ) Tr Ω Π − d Ω (cid:16) Tr G \ Ω (Tr Ω ρ ) Tr Ω Π + (Tr G ρ ) (Tr Ω Π ) (cid:17)(cid:17)(cid:69) ρ = 1 d Ω + 1 (cid:16) P G + P G \ Ω (cid:17) , (S.5)2 A B D (a) (b) N AB N (c) G ω or Ω C F α β γ 𝛿 A B C D N BD N BC N CD N AD N AC N N N β 𝛿 γα B C D A
FIG. S1.
System subdivision and matrix-product-state representation . (a) Subdivision of the full system F into foursubsets A = F \ ( G ∪ ω ), B = G \ ω , C = ω ∩ G , and D = ω \ G . (b) Matrix product structure in 1D systems correspondingto Eq. (S.10). Each square represents a matrix, while N i are local bond dimensions. (c) Most general MPS structure, where N GG (cid:48) represents the local bond dimension between sets G and G (cid:48) . where we used the correlation function in Eq. (S.2).The case when Ω / ∈ G can be approached in a similar fashion. Making use of the symmetry ˜ P (cid:48) F \ G = ˜ P (cid:48) G for purestates, it is convenient to consider the purity ˜ P (cid:48) F \ G of the complement of G instead of ˜ P (cid:48) G itself. Since Ω ∈ F \ G , wecan then use Eq. (S.5) to obtain˜ P (cid:48) G = ˜ P (cid:48) F \ G = 1 d Ω + 1 (cid:16) P F \ G + P ( F \ G ) \ Ω (cid:17) = 1 d Ω + 1 (cid:16) P G + P G ∪ Ω (cid:17) . (S.6)Both cases can be summarized as ˜ P (cid:48) G = 1 d Ω + 1 (cid:16) P G \ Ω + P G ∪ Ω (cid:17) . (S.7)The denominator of Eq. (S.4) can also be derived from Eq. (S.5) by taking G = F :˜ P (cid:48) F = 1 d Ω + 1 (cid:16) P F + P F \ Ω (cid:17) = 1 d Ω + 1 (cid:16) P Ω (cid:17) . (S.8)Inserting Eqs. (S.7) and (S.8) into Eq. (S.4), we obtain Eq. (3) in the main text.As mentioned in the main text, the approximation S [ A ] = − (cid:10) log(Tr A ρ A ) (cid:11) (cid:39) − log( (cid:10) Tr A ρ A (cid:11) ) = − log P A (S.9)gives a good estimate of the average entropy S ( A ). The accuracy of this approximation is illustrated in Fig. S2(b),using as a model a 1D system in a random matrix product state, | ψ (cid:105) = (cid:80) s ...s k Tr ( A s . . . A s k k ) | s . . . s k (cid:105) , where A s i i are N b × N b random matrices sampled from the Gaussian unitary ensemble with local bond dimension N b [ x axis inFig. S2(b)]. Fig. S2(b) suggests that the annealed approximation remains accurate even for states with small bonddimension, producing an error as low as (cid:15) ∼ − for N b ∼ SECTION II: EVOLUTION OF HARTLEY ENTROPIES
In this section, we derive the evolution Equations (10,11) for subsystem rank (whose logarithm is the Hartleyentropy).To study the evolution of the rank, it is convenient to use a tensor network representation of the system state. Wedivide the system into four subsets shown in Fig. S1(a) and associate with each subset a corresponding identicallynamed tensor A α , B β , C γ , or D δ , where indices α , β , γ , and δ run over corresponding basis states of each subsetof qudits. The structure of the tensors will depend on the dimensionality of the system and other properties of thecircuit structure.Let us start with a simple calculation for a 1D setting. For 1D systems, the tensor network is a matrix productstate (MPS) representation of the system that can be written as | ψ (cid:105) = (cid:88) n ...n ,αβγδ A αn n B βn n C γn n D δn n | αβγδ (cid:105) , (S.10)3 A Ω (a) (b) (c) ( N b )10 R e l a t i v e e rr o r , R e l a t i v e e rr o r , Annealed approx. Mean fi eld approx.Local dimension index, log ( N b ) FIG. S2.
Numerical analysis of approximations . (a) The system consists of 9 qubits forming a 1D chain with periodicboundary conditions. The system is divided into a 4-qubit subsystem of interest A , the measured qubit Ω, and the rest ofthe system. (b) Relative error of the annealed approximation in Eq. (S.9) for the R´enyi entropy. The error is quantified as (cid:15) = (cid:12)(cid:12) S [ A ] − S (cid:48) [ A ] (cid:12)(cid:12) /S [ A ], where S [ A ] = −(cid:104) log(Tr A ρ A ) (cid:105) and S (cid:48) [ A ] = − log( (cid:104) Tr A ρ A (cid:105) ). The curve shows the error as a functionof of the logarithm log ( N b ) of the local bond dimension N b . (c) Relative error of the mean-field approximation for the correlationdistance measure in Eq. (S.93). The error is quantified as (cid:15) = | C A, Ω − C A, Ω | /C A, Ω , where C A, Ω = (cid:68)(cid:13)(cid:13) ρ A ∪ Ω − ρ A ⊗ ρ Ω (cid:13)(cid:13) (cid:69) and C A, Ω = P A ∪ Ω − P A P Ω . where indices n i = 1 , . . . , N i run over the N i basis states of the bond connecting two of the subsystems, as shownin Fig. S1(b). Now consider a transformation V supported on subsystem ω composed of C and D . The MPSrepresentation of the resulting state can be written as V | ψ (cid:105) = (cid:88) n ,αβγδ A αn n B βn n ˜ C γn n (cid:48) ˜ D δn (cid:48) n | αβγδ (cid:105) , (S.11)where n (cid:48) = 1 , . . . , N (cid:48) is expressed through the new bond dimension N (cid:48) between tensors ˜ C and , and the new MPSmatrices ˜ C and ˜ D are connected to the initial ones by the relation (cid:88) n (cid:48) ˜ C γn n (cid:48) ˜ D δn (cid:48) n = (cid:88) n ,γ (cid:48) δ (cid:48) V γδ,γ (cid:48) δ (cid:48) C γ (cid:48) n n D δ (cid:48) n n . (S.12)If V is a unitary operator, this relation contains n s = N N d G ∩ ω d ω \ G equations. The matrix ˜ C has d G ∩ ω N N (cid:48) degrees of freedom (i.e indices), while the matrix ˜ D has d ω \ G N N (cid:48) degrees of freedom. Importantly, N (cid:48) of theseindices are redundant due to the existence of a gauge transformation ˜ C (cid:48) = ˜ CM and ˜ D (cid:48) = M − ˜ D for any invertible N (cid:48) × N (cid:48) matrix M . As a result, for the number of equations to coincide with the number of independent variables,the new rank N (cid:48) must satisfy d G ∩ ω d ω \ G N N − d G ∩ ω N N (cid:48) − d ω \ G N N (cid:48) + N (cid:48) = 0 . (S.13)This expression can be rewritten as ( N (cid:48) − d G ∩ ω N )( N (cid:48) − d ω \ G N ) = 0 , an equation with two positive roots. One ofthe solutions can be used to compute the rank of the subset G , R (cid:48) G = rank( B β ˜ C γ ) = N N (cid:48) . (S.14)Taking into account that R G \ ω = rank( B β ) = N N and R G ∪ ω = rank( B β C γ D δ ) = N N , we obtain R (cid:48) G = min (cid:0) d G ∩ ω R G \ ω , d ω \ G R G ∪ ω (cid:1) , (S.15)which is Eq. (10) in the main text. In the case when d G ∩ ω = d ω \ G = q the formula simplifies to R (cid:48) G = q min (cid:0) R G \ ω , R G ∪ ω (cid:1) , (S.16)leading to the following expression for the Hartley entropy: S [ G ] = min (cid:0) S [ G \ ω ] , S [ G ∪ ω ] (cid:1) + log q. (S.17)Consider now a measurement that projects the system onto state | ψ (cid:48) (cid:105) = V | ψ (cid:105) / (cid:112) (cid:104) ψ | V | ψ (cid:105) , where V is supported onsubset Ω. We consider the same subdivision of the system up to the replacement ω → Ω. Without loss of generality,4the transformation V can be chosen to be a projector on a subset of Ω of the form V = I AB ⊗ | (cid:105)(cid:104) | CD . Therefore,˜ B β = 0 and ˜ C γ = 0 for any β, γ (cid:54) = 0. The zeroth components of these tensors are connected to the initial MPSrepresentation via (cid:88) n (cid:48) ˜ C n n (cid:48) ˜ D n (cid:48) n = const × (cid:88) n C n n D n n , (S.18)where the constant arises due to normalization.Using analysis similar to the case of unitary V , we find that Eq. (S.18) contains N N equations for N N (cid:48) + N N (cid:48) − N (cid:48) variables. Hence, the new rank satisfies R (cid:48) G = min (cid:0) R G \ Ω , R G ∪ Ω (cid:1) , (S.19)which is Eq. (11) in the main text.These results can be generalized to the more general (beyond 1D) tensor structure shown in Fig. S1(c). The MPSthen contains four-index tensors capturing the most general structure of the state: | ψ (cid:105) = (cid:88) n ,αβγδ A αn AB ,n AC ,n AD B βn AB ,n BC ,n BD C γn AC ,n BC ,n CD D δn AD ,n BD ,n CD | αβγδ (cid:105) , (S.20)where indices n GG (cid:48) = 1 , . . . , N GG (cid:48) run through the N GG (cid:48) basis state of the bond between sets G and G (cid:48) .Similarly to the 1D case, we consider a transformation V that affects subsystem ω and changes tensors C and D .The new tensors are related to the original ones via (cid:88) n (cid:48) CD ˜ C γn AC ,n BC ,n (cid:48) CD ˜ D δn AD ,n BD ,r (cid:48) CD = (cid:88) n CD ,β (cid:48) γ (cid:48) V γδ,γ (cid:48) δ (cid:48) C γ (cid:48) n AC ,n BC ,n CD D δ (cid:48) n AD ,n BD ,n CD . (S.21)Analysis similar to the 1D case applies here. Equation (S.21 contains n s = d G ∩ ω d ω \ G N AC N BC N AD N BD equationsfor the d G ∩ ω N AC N BC N (cid:48) CD variables in tensor ˜ C and the d ω \ G N AD N BD N (cid:48) CD variables in tensor ˜ D . Excluding N (cid:48) CD redundant variables, we find that N (cid:48) CD satisfies N (cid:48) CD − d G ∩ ω N AC N BC N (cid:48) CD − d ω \ G N AD N BD N (cid:48) CD + d G ∩ ω d ω \ G N AC N BC N AD N BD = 0 , (S.22)which can be rewritten as ( N (cid:48) CD − d G ∩ ω N AC N BC )( N (cid:48) CD − d ω \ G N AD N BD ) = 0 , (S.23)which has two positive roots. Taking into account that R (cid:48) G = rank( B β C γ ) = N AB N BD N AC N (cid:48) CD (S.24)and R G \ ω = rank( B β ) = N AB N BC N BD , R G ∪ ω = rank( B β C γ D δ ) = rank( A α ) = N AB N AC N AD , (S.25)and taking the smaller root in Eq. (S.23), we obtain the following equation for the rank: R (cid:48) G = min( d G ∩ ω N AC R G \ ω , d ω \ G N BD R G ∪ ω ) . (S.26)To express quantities such as N AC and N BD , one needs expressions for the rank corresponding to other subsets ofthe system: R G ∩ ω = rank( C γ ) = N AC N BC N CD ,R ω \ G = rank( D δ ) = N AD N BD N CD ,R G ∆ ω = rank( A α C γ ) = rank( B β D δ ) = N AB N AD N BC N CD . (S.27)It follows from these expressions that N BD = R G \ ω R ω \ G R G ∆ ω , N AC = R ω \ G R G ∩ ω R G ∆ ω . (S.28)5Substituting these equalities into Eq. (S.26), we obtain the final expression connecting the rank of subspace G afterthe application the unitary V to ranks before the application: R (cid:48) G = R G \ ω R G ∪ ω R G ∆ ω min( d G ∩ ω R G ∩ ω , d ω \ G R ω \ G ) . (S.29)A similar expression can be obtained for the effect of a measurement: R (cid:48) G = R G \ Ω R G ∪ Ω R G ∆Ω min( R G ∩ Ω , R Ω \ G ) , (S.30)where, as before, Ω is the measured set.In 1D systems, the evolution described by Eqs. (S.15) and (S.19) can be mapped to the minimum cut problem ona square lattice. To show this, we use the trajectory representation we defined for purities, but now formulated forranks R G . We formally write the evolution in the form R (cid:48) G = min G ( L G,G (cid:48) R G (cid:48) ) , (S.31)where the Liouvillian L has the same structure as for the evolution of physical purities P given by Eq. (2) and thenumerator of Eq. (3), with c − = d G ∩ ω and c + = d ω \ G .As for purities, we consider all possible trajectories G = { G → G · · · → G n } such that G k +1 ∈ N k [ G k ], and G = A . Then, if the system is initialized in a product state, R G (0) = 1, the evolution can be presented in the form R A = min G R G , R G = n (cid:89) k =1 L kG k ,G k +1 . (S.32)The proof of this expression can be carried out by induction. First, the expression in Eq. (S.31) is equivalent toEq. (S.32) under the condition that initially R G = 1 for any G . Next, we assume that Eq. (S.32) holds for a circuit ofdepth n − ≥ {L , . . . L n } and consider different n -step trajectories G [ G (cid:48) ] starting fromthe set G (cid:48) . Then R n +1 A = min G (cid:48) (cid:16) L G,G (cid:48) n (cid:89) G k ∈G [ G (cid:48) ] L kG k ,G k +1 (cid:17) = min G R G . (S.33)In the case of a 1D brickwork configuration, the rank can be expressed using the same statistical model as in Eq. (9), R A = min G exp( − αE [ G ]) , (S.34)where α = log q and E [ G ] are the energies of all allowed configurations corresponding to G [see Eqs. (S.55)-(S.57)].The minimum configuration energy for the statistical model is then given by the minimum total length of domainwalls excluding all dashed segments on a percolated square lattice (see Sec. IV for more details). In the case of abipartition, this reduces to a minimum path configuration connecting the cut at the top of the lattice to its otherboundaries. SECTION III: HYDRODYNAMIC EQUATIONS FOR ENTROPIES IN 1D
In this section, we derive the 1D KPZ equation for R´enyi entropies [Eq. (5) in the main text] and its analyticalsolution [Eq. (6)], as well as the continuous evolution equation for Hartley entropies [Eq. (12)].Let us start from the formal solution of Eq. (5). Substituting the definition of the entropy S ( t, x ) = − log P c ( t, x ), weobtain a linear partial differential equation for the continuous approximation of the average purity P c ( t, x ) (subscript c refers to the continuous approximation), ∂ t P c = µ∂ x P c − βP c , P c ( t, x = 0) = P c ( t, x = N ) = 1 . (S.35)For simplicity, we take a product initial state with P c ( t = 0 , x ) = 1. Substituting P c = 1 + δP c we obtain ∂ t δP c = µ∂ xx δP c − βδP c − β, (S.36)6with zero initial and boundary conditions δP c ( t,
0) = δP c ( t, N ) = δP c (0 , x ) = 0. Equation (S.36) can be solved by theFourier method. Expanding the solution in a Fourier series, δP c = ∞ (cid:88) n =1 A n ( t ) sin ( k n x ) , k n = nπ/N, (S.37)and expanding the constant term β on the RHS of (S.36) in a Fourier series, β = ∞ (cid:88) n =1 β n sin k n x, where β n = 2 βN (cid:90) N sin ( k n x ) dx = (cid:40) , if n = 2 m β/πn, if n = 2 m + 1 , (S.38)we obtain a set of ordinary differential equations for coefficients A n ( t ),˙ A n ( t ) = − ( β + µ ) A n ( t ) − βπn δ n,odd . (S.39)Here δ n,m is the Kronecker delta function. The solution to (S.39) reads A n ( t ) = (cid:40) A n (0) e − ( β + µk n ) t , if n = even A n (0) e − ( β + µk n ) t + βk n ( β + µk n ) (cid:16) e − ( β + µk n ) t − (cid:17) , if n = odd . (S.40)Finally, using the initial conditions A n (0) = 0, we obtain the solutions P c ( t, x ) = 1 + 4 βN ∞ (cid:88) m =0 (cid:16) e − ( β + µk m +1 ) t − (cid:17) sin ( k m +1 x ) k m +1 ( β + µk m +1 ) , (S.41) S ( t, x ) = − log( P c ( t, x )) . (S.42)As t → ∞ , the exponent in (S.41) vanishes. The corresponding stationary purity (i.e. at t → ∞ ) obeys the followingequation: µ∂ xx P − βP = 0 . (S.43)The solution of Eq. (S.43) with the boundary conditions in Eq. (S.35) reads P c ( t → ∞ , x ) = cosh ( α ( x − N/ αN/ , (S.44)where α = (cid:112) β/µ . Equation (S.44) also follows directly from Eq. (S.41). The asymptotic value of the entropy at x = N/ S ( t → ∞ , x = N/
2) = log [cosh ( αN/ . (S.45)In the thermodynamic limit ( N → ∞ ), we obtain volume-law growth of the entropy S ( t → ∞ , x = N/
2) = αN/ − log 2 + O (exp( − αN )) . We now derive parameters β and µ in the continuous Equation (5) from the underlying discrete evolution. To findparameter β , we consider a homogeneous solution of the discrete model, P x ( t ) = P x +1 ( t ). In this case, after a singletime step, one has P x ( t + 1) = 2 cP x ( t ). Similarly, for the entropy evolution, S ( t + 1 , x ) = S ( t, x ) + β . Taking S ( t, x ) ≈ − log P c ( t, x ) and comparing these two expressions, we find β = − log(2 c ) = log (cid:16) q + 12 q (cid:17) . (S.46)To find the diffusion coefficient µ , we consider a stationary solution P x ( t ) = P x ( t + 1). Let us use the following ansatzfor the solution to the discrete equation: P x +1 ( t ) = γP x ( t ) . Then the condition of time-invariance leads us to P x ( t + 1) = P x ( t ) = c (cid:16) γ + γ (cid:17) P x ( t ) . (S.47)7This condition is satisfied if γ = 12 c (cid:16) ± (cid:112) − c (cid:17) ∈ { q, q − } , (S.48)which corresponds to positive and negative slopes for the R´enyi entropy of the first and second halves of the chain,respectively.Assuming continuity of the entropy, we find that ∂S ∂x (cid:39) S ( t, x + 1) − S ( t, x ) = (cid:115) βµ = log q. (S.49)This finally leads to the expression µ = log (cid:2) ( q + 1) / q (cid:3) (log q ) . (S.50)For q = 2, this gives µ ≈ . µ = 1 / q → ∞ , the diffusion coefficient vanishes as µ ∼ / log q . Derivation of Eq. (12):
Using the expression for the evolution of the Hartley entropy under unitary gates, S ( t + 1 , x ) = min [ S ( t, x − , S ( t, x + 1)] + log q , and using identities min( a, b ) = a + min( b − a,
0) and min( a,
0) = a −| a | , we find that ∂ t S (cid:39) S ( t + 1 , x ) − S ( t, x ) (cid:39) [ S ( t, x + 1) − S ( t, x − − | S ( t, x + 1) − S ( t, x − | ]2+ S ( t, x − − S ( t, x ) + log q (cid:39) −| ∂ x S | + 12 ∂ x S + log q, (S.51)which results in Eq. (12). In deriving (S.51), we used approximations [ S ( t, x + 1) − S ( t, x − / (cid:39) ∂ x S and S ( t, x + 1) − S ( t, x ) (cid:39) ∂ x S + ∂ x S . Solution of the hydrodynamic Equation (12) for S in 1D: In the case of a product initial state, the initialcondition for the Hartley entropy in 1D is S ( t = 0 , x ) = 0. It is easy to see by direct substitution that the solutionof the hydrodynamic equation for the Hartley entropy (see Eq. (12) and Table I) reads: S ( t, x )log q = x, if 0 ≤ x ≤ t,t, if t ≤ x ≤ N − t,t − x, if N − t ≤ x ≤ N. (S.52)Moreover, within each segment defined in Eq. (S.52), the second derivative ∂ x S vanishes, except for special pointswhere the second derivative is discontinuous. Therefore, in the special case when the initial state is a product state,the hydrodynamic equation for the Hartley entropy reduces to ∂ t (cid:101) S = 1 − | ∂ x (cid:101) S | , (S.53)where (cid:101) S = S / log q .On the other hand, we can write the hydrodynamic KPZ equation for the R´enyi entropy (cid:101) S = S / log q in thescaling limit q → ∞ : ∂ t (cid:101) S (cid:39) − ( ∂ x (cid:101) S ) + 1 , (S.54)where we used the asymptotic values for the coefficients β = log(( q + 1) / q ) → log q and µ = β/ (log q ) → / log q at q → ∞ . Note that the piecewise solutions with slopes 0, ± q → ∞ , R´enyi and Hartley entropies coincide, i.e. S ( t, x ) → S ( t, x ), provided thatthe initial state is a product state, i.e. S (0 , x ) = S (0 , x ) = 0.In Fig. S3, we compare the dynamics of second R´enyi and Hartley entropies in 1D by solving discrete (exact) andhydrodynamic (approximate) evolution equations listed in Table I. In Fig. S3(a,b), one can see that the hydrodynamicevolution (dashed lines) agrees well with the exact solution (solid lines). At asymptotically large times, the stationarysolutions for s = S /N and s = S /N have the same shape. In Fig. S3(c), we plot the values of the entropy inthe middle of the 1D chain, x = N/
2, for increasing values of the qudit dimension q . The plot shows that, in theasymptotic limit q → ∞ , the second R´enyi entropy approaches the Hartley entropy, s ( t, x ) → s ( t, x ).8 (a) (b) (c) FIG. S3. (a,b) Evolution of second R´enyi and Hartley entropies in 1D under a brickwork Haar-random unitary circuit: com-parison of discrete-time evolution and the solutions of the corresponding continuous-time hydrodynamic equations (equationsare shown in Table I). Here we plot the specific entropy (entropy per qudit) s , = S , /N for q = 2. The continuous-timeequations were solved by a finite-difference scheme with a time-step dt = 10 − . We chose an initial state where the two halvesof the 1D system are disentangled from each other but each is in a maximally entangled state. (c) Second R´enyi and Hartleyentropies for a cut in the middle of the 1D qudit chain ( x = N/
2) as a function of time t for q = 2 , , ,
64, calculated usingdiscrete time-evolution equations. The Hartley entropy has a universal scaling law s ∼ log q for arbitrary q . The 2-R´enyientropy approaches the Hartley entropy in the limit of large qudit dimension log q → ∞ : s ( t, x ) → s ( t, x ). At long times,both entanglement entropies saturate to a volume-law scaling in the thermodynamic limit: S ( t → ∞ , x = N/ → N log q and S ( t → ∞ , x = N/ → N log q − log 2.TABLE I. Discrete-time and continuous-time evolution equations for the second R´enyi entropy and for the Hartley entropy.R´enyi entropy S Hartley entropy S Discrete-time evolution P t +1 x = qq + 1 ( P tx − + P tx +1 ) ,S t +12 ,x = − log P t +1 x S t +10 ,x = min ( S t ,x − , S t ,x +1 ) + log q Hydrodynamic evolution ∂ t S = µ (cid:0) ∂ x S − ( ∂ x S ) (cid:1) + β ∂ t S = 12 ∂ x S − | ∂ x S | + log q SECTION IV: MAPPING TO A CLASSICAL MODEL IN 1D
In this section, we derive the classical Hamiltonian H and the projector P in Eqs. (8,9) in the main text.We consider combined circuit layers each consisting of a single unitary layer and a single measurement layer in reversechronological order. This combined circuit can be mapped to a collection of trajectories G = { G (cid:48) → G (cid:48) → · · · → G (cid:48) n } as described in the main text, where the Liouvillian L um k for a single combined layer connects set G (cid:48) k to G (cid:48) k +1 , andeach trajectory contributes to the purity according to Eq. (7). Our goal is now to map each trajectory to a classicalstate on a lattice with dynamical variables attached to the center of each lattice cell and to edges associated withmeasurements. At the end of this section, we cover open boundary conditions.As a toy illustrative example, let us focus on a primitive circuit containing only two qudits, labeled Ω and Ω ,and a single combined unitary-measurement layer applied to them [see top of Fig. S4(a)]. This primitive systemcorresponds to a simple lattice with only three vertices and two edges shown inside a blue rectangle at the bottom ofFig. S4(a). Therefore, our proposed classical model has three lattice cell variables ( Z , Z , Z ) and two edge variables( J and J ) [see labels in Fig. S4(b)]. To study the mapping of trajectory { G (cid:48) → G (cid:48) } to a classical state on a lattice,we recall the trajectory defined in Fig. 1(b) of the main text, where each unitary and measurement layer is treatedseparately. Our single-step trajectory { G (cid:48) → G (cid:48) } thus corresponds to a two-step trajectory { G → G → G } in thelanguage of Fig. 1(b). For both trajectories, the starting configuration is G = G (cid:48) = A . Notably, as we will see below,the final state for any trajectory is either G (cid:48) = G = F (full system) or ∅ , regardless of measurement locations.Let us start from the case of a unitary circuit without measurements. Such a circuit can be described by six possibletrajectories including all possible initial states A , as shown below in the first column of Eq. (S.55). These trajectoriescan be mapped to the classical spin configurations shown in the second column of Eq. (S.55). The contribution P G of9 G =G G G =G Z Z Z J J (a) (b) (c) (d) prohibited: prohibited: prohibited: FIG. S4.
Correspondence between trajectories G and configurations of a dynamical square lattice. (a) Top: two-qudit trajectory G = { G → G → G } = {∅ → Ω → F } illustrated by blue and white cells reflecting, respectively, whetherthe qudit belongs to G k or not. The monitored circuit generating such a trajectory is shown on the right (the black dot is asingle-qudit measurement, while the green rectangle is a two-qudit gate). Bottom: the blue rectangle marks the portion of thesquare lattice corresponding to the { G → G → G } trajectory shown at the top. Here the colors of the top left and top rightsquare cells reflect the colors of qudits in the initial configuration G (cid:48) = G , while the color of the bottom square reflects the colorof the final configuration G (cid:48) = G . The intermediate configuration G is encoded in the edge below the measurement (solidedge if the cell keeps its color in G and dashed edge if the color changes). (b-d) Illustration of the correspondence betweentrajectories (on the background) and lattice configurations (b) in the absence of measurements, (c) with a measurement of onequdit, and (d) with measurements of both qudits. The examples of prohibited lattice configurations are shown at the bottom;these configurations are projected out using operator P . each trajectory G to the purity in Eq. (7) is given in the third column in the form of the configuration energy E [ G ]:Trajectory G Classical configuration E [ G ] = − log( P G ) /β {∅ → ∅ → ∅} , { Z = − , Z = − , Z = − , J = 1 , J = 1 } { Ω → Ω → ∅} , { Z = +1 , Z = − , Z = − , J = 1 , J = 1 } { Ω → Ω → F } , { Z = +1 , Z = − , Z = +1 , J = 1 , J = 1 } { Ω → Ω → ∅} , { Z = − , Z = +1 , Z = − , J = 1 , J = 1 } { Ω → Ω → F } , { Z = − , Z = +1 , Z = +1 , J = 1 , J = 1 } { F → F → F } , { Z = +1 , Z = +1 , Z = +1 , J = 1 , J = 1 } G k as blue and the rest of the qudits as white.For unitary circuits, the edges are not treated as dynamical variables; instead, the edges are all taken to be numbersequal to unity, i.e. J ij = +1. Also, no trajectory is mapped to the configurations { Z = +1 , Z = +1 , Z = − , J =1 , J = 1 } and { Z = − , Z = − , Z = +1 , J = 1 , J = 1 } ; therefore, these configurations must be projectedout using the operator P . We show the first of these prohibited configurations at the bottom of Fig. S4(b). Eachallowed configuration contributes with energy E [ G ] = 1 equal to the total length of domain walls (i.e. the number ofred edges).Now let us assume that a measurement is applied to one of the qudits Ω k , where k = 1 or k = 2, while the otherqudit Ω k remains unmeasured (here we define 1 = 2 and 2 = 1). The measurement adds more accessible trajectories tothe ones already listed in Eq. (S.55). Here are these additional trajectories, the corresponding classical configurations,10and the corresponding energies:Trajectory G Classical configuration E [ G ] = − log( P G ) /β {∅ → Ω k → ∅} , { Z k = − , Z k = − , Z = − , J k = − , J k = 1 } {∅ → Ω k → F } , { Z k = − , Z k = − , Z = +1 , J k = − , J k = 1 } { Ω k → ∅ → ∅} , { Z k = +1 , Z k = − , Z = − , J k = − , J k = 1 } { Ω k → F → F } , { Z k = − , Z k = +1 , Z = +1 , J k = − , J k = 1 } { F → Ω k → F } , { Z k = +1 , Z k = +1 , Z = +1 , J k = − , J k = 1 } { F → Ω k → ∅} , { Z k = +1 , Z k = +1 , Z = − , J k = − , J k = 1 } k = 1 (and k = 2). As for the unitarycase, J k, = 1; however, since we measured Ω k , J k, is now a dynamical variable. Since no trajectories correspond toconfigurations { Z = − , Z = +1 , Z = − , J k = − , J k = 1 } and { Z = +1 , Z = − , Z = +1 , J k = − , J k =1 } , these configurations must be projected out by P . The first of these two prohibited configurations is shown at thebottom of Fig. S4(c).Finally, if both qudits are measured, the following trajectories become accessible, in addition to those already listedin Eq. (S.55) and Eq. (S.56):Trajectory G Classical configurations E [ G ] = − log( P G ) /β {∅ → F → F } , { Z = − , Z = − , Z = +1 , J = − , J = − } { Ω → Ω → ∅} , { Z = +1 , Z = − , Z = − , J = − , J = − } { Ω → Ω → F } , { Z = +1 , Z = − , Z = +1 , J = − , J = − } { Ω → Ω → ∅} , { Z = − , Z = +1 , Z = − , J = − , J = − } { Ω → Ω → F } , { Z = − , Z = +1 , Z = +1 , J = − , J = − } { F → ∅ → ∅} , { Z = +1 , Z = +1 , Z = − , J = − , J = − } { Z = +1 , Z = +1 , Z = +1 , J k = − , J k = − } and { Z = − , Z = − , Z = − , J k = − , J k = − } , theseconfigurations must be projected out by P . Both prohibited configurations correspond to the “lambda”-configurationof dashed lines shown at the bottom of Fig. S4(d).Using Eqs. (S.55), (S.56), and (S.57), it is straightforward to verify that the energy of each allowed configuration isdescribed by the Ising-type interaction, E [ G ] = 1 − (cid:16) J Z Z + J Z Z (cid:17) , (S.58)while all prohibited configurations satisfy 14 ( Z − J Z )( Z − J Z ) = 1 . (S.59)Therefore, the purity of set A is given by˜ P A = Tr (cid:16) P exp( − βH ) (cid:17) , H = 1 − (cid:88) (cid:104) i,j (cid:105) J ij Z i Z j , (S.60)where the projector P excludes prohibited configurations as well as sets the initial conditions for trajectory G : P = P A (cid:16) −
14 ( Z − J Z )( Z − J Z ) (cid:17) , (S.61)where P A = (cid:81) i ∈ A (1 − Z i ) (cid:81) j ∈ F \ A (1 + Z j ), where products are taken over the Z -variables at the top of the lattice(i.e. Z and Z ).Having considered the case of two qudits and one layer, we can generalize the result to an arbitrary number of quditsand layers. This generalization is possible because: (a) the process is Markovian, so contributions from different layersdo not interfere; (b) gates in the same layer do not overlap, so their contributions are also independent. Therefore, the11 prohibited: Σ b FIG. S5.
Boundary conditions . Open boundaries put additional constraints on the allowed spin-model configurations. Edgepairs of squares (green boxes) cannot exhibit configurations shown on the right, i.e. solid red and dashed black edges areprohibited in the indicated locations independently of square colors. full multi-qudit multi-layer configuration on a square lattice is obtained by glueing together 2-qudit 1-layer elementsdescribed above, while the effective energy of this full configuration is just the sum of the energies of these elements: H = N − (cid:88) (cid:104) i,j (cid:105) J ij Z i Z j , P = P A (cid:89) ( i,j,k ) ∈ Σ (cid:16) −
14 ( Z k − J ik Z i )( Z k − J jk Z j ) (cid:17) , (S.62)which is Eq. (8) in the main text.Another approach to dealing with prohibited configurations (without projecting them out) is to consider a modifiedIsing Hamiltonian with a penalty term: ˜ P A = lim Γ →∞ Tr P A exp( − βH Γ ) , (S.63)where H Γ = 12 N − (cid:88) (cid:104) i>j (cid:105) J ij Z i Z j + Γ (cid:88) ( i,j,k ) ∈ Σ ( Z k − J ik Z i )( Z k − J jk Z j ) . (S.64)In the limit of infinite penalty coefficient Γ, the two approaches are equivalent.One can compute the energies of domain-wall configurations without explicitly using the Hamiltonian. For unitaryevolution, all Ising couplings are J ij = 1; therefore, the energy of any allowed Ising configuration E [ G ] is equal to thetotal length of domain walls. In the presence of measurements, E [ G ] is equal to the total length of solid segmentsof domain walls (i.e. solid red edges) plus the total number of dashed edges outside of the domain walls (i.e. dashedblack edges). The projector P prohibits configurations shown in Fig. 1(c). Because dashed edges resemble cuts on alattice, we can draw an analogy of the model to trajectories on a dynamical percolated lattice. Boundary conditions . Boundaries result in the presence of unpaired qudits/cells at the edges that are notaffected by the unitaries at even (or odd, depending on where exactly the boundary is) layers, as can be seen fromFig. 1(a). These qudits exhibit additional prohibited trajectories at even (or odd) times. Specifically, black dashedand solid red lines are prohibited in the locations shown in Fig. S5. Such trajectories can be projected out by taking
P → P (cid:81) ij ∈ Σ b (1 + J ij Z i Z j ) for all boundary square pairs Σ b highlighted as green boxes in Fig. S5. SECTION V: CONTINUOUS EVOLUTION
In this section, we show how continuous evolution of the system under a stochastic Hamiltonian and local continuousmeasurements gives rise to Eqs. (14) and (15) in the main text. We also analyze the accuracy of the mean-fieldapproximation used to simplify Eq. (15). Finally, we show that the continuous evolution in Eqs. (14,15) can be usedto derive the discrete evolution in Eqs. (2,3)First, we prove that the eigenbasis of the monitored operators can be randomized while preserving purities. Thisrandomization will allow us to derive a differential equation for the purities. More specifically, we consider ran-dom time-dependent unitaries U † τ,i ( t ) supported on a single qudit i and drawn from the circular unitary ensem-ble. In order to have a time-continuous description, we assume that the matrices are weakly correlated in time, (cid:104) [ U τ,i ( t )] µν [ U ∗ τ,j ( t (cid:48) )] µ (cid:48) ν (cid:48) (cid:105) = q − δ µµ (cid:48) δ νν (cid:48) δ ij f [( t − t (cid:48) ) /τ ], where f [ x ] is a smooth correlation function satisfying f [0] = 1and lim | x |→∞ f [ x ] = 0. The correlation f ( x ) will be later removed by taking the limit τ →
0. In this limit, the12unitaries represent a stochastic process where, at each time, the matrix U τ =0 ,i ( t ) is drawn independently from thecircular unitary ensemble.Consider the unitary U τ ( t ) = (cid:81) i U τ,i ( t ) and the corresponding rotating-frame wavefunction | ψ (cid:48) ( t ) (cid:105) = U † τ ( t ) | ψ ( t ) (cid:105) . (S.65)Since U τ,i ( t ) are defined as single-qudit unitary operators, the purity P G for any G remains invariant under thistransformation: P G ( | ψ (cid:48) (cid:105)(cid:104) ψ (cid:48) | ) = P G ( | ψ (cid:105)(cid:104) ψ | ) . (S.66)The new wavefunction | ψ (cid:48) ( t ) (cid:105) satisfies ddt | ψ (cid:48) (cid:105) = − i (cid:16) H (cid:48) ( t ) + H U ( t ) (cid:17) | ψ (cid:48) (cid:105) + − (cid:88) j κ j (cid:16) O j ( t ) − (cid:104) ψ (cid:48) | O j ( t ) | ψ (cid:48) (cid:105) (cid:17) + (cid:88) j ξ j ( t ) (cid:112) κ j (cid:16) O j ( t ) − (cid:104) ψ (cid:48) | O j ( t ) | ψ (cid:48) (cid:105) (cid:17) | ψ (cid:48) (cid:105) , (S.67)where H U ( t ) ≡ − i (cid:88) i U † τ,i ( t ) ∂ t U τ,i ( t ) , (S.68) H (cid:48) ( t ) ≡ U † τ ( t ) H ( t ) U τ ( t ) = (cid:88) i I F \ ω i ⊗ h (cid:48) i ( t ) , (S.69) O j ( t ) ≡ U † τ ( t ) O j U τ ( t ) = U † τ,j ( t ) O j U τ,j ( t ) . (S.70)The modified stochastic Hamiltonians h (cid:48) i ( t ) ≡ U † τ,i ( t ) h i ( t ) U τ,i ( t ) have the same distribution as h i ( t ), due to theinvariance of the Gaussian unitary ensemble under unitary rotations. The new Hamiltonian term H U ( t ) is a sum ofsingle-qudit operators. Therefore, the rotation in Eq. (S.65) effectively randomizes the measurement bases: O j → U † τ,j ( t ) O j U τ,j ( t ). The penalty for this is the addition of single-qudit terms to the Hamiltonian, but with no effecton the distribution of the original stochastic Hamiltonian terms. Importantly, we show below that the additionalsingle-qudit terms do not affect (averaged) purities P G . Unitary dynamics . Let us first focus on the effect of the combined Hamiltonian H c ( t ) = H U ( t ) + (cid:88) i H (cid:48) i ( t ) , (S.71)where each local stochastic term can be represented as H (cid:48) i ( t ) = I F \ ω i ⊗ h (cid:48) i ( t ), where h (cid:48) i ( t ) are independent randommatrices supported on sets ω i as introduced in the main text. The unitary evolution operator for infinitesimal time dt is U t,dt = T exp (cid:16) − i (cid:90) t + dtt dt (cid:48) H c ( t (cid:48) ) (cid:17) . (S.72)Therefore, the evolution of the density matrix satisfies ρ ( t + dt ) = U t,dt ρU † t,dt = ρ ( t ) − i (cid:90) t + dtt [ H U ( t (cid:48) ) , ρ ] dt (cid:48) − i (cid:88) i (cid:90) t + dtt [ H (cid:48) i ( t (cid:48) ) , ρ ] dt (cid:48) − (cid:88) ij (cid:90) t + dtt dt (cid:48) (cid:90) t + dtt dt (cid:48)(cid:48) [ H (cid:48) i ( t (cid:48) ) , [ H (cid:48) j ( t (cid:48)(cid:48) ) , ρ ]] + O ( dt ) . (S.73)The evolution of purity is connected to the density matrix ρ = | ψ (cid:48) (cid:105)(cid:104) ψ (cid:48) | as follows: ddt P G = ddt (cid:104) Tr G ρ G (cid:105) = lim dt → (cid:68) Tr G ρ G ( t + dt ) − ρ G ( t ) dt (cid:69) circ ,H (cid:48) i ( t ) . (S.74)The averaging in Eq. (S.74) includes both the circuit average (cid:104) . . . (cid:105) circ over different evolution histories realized by thestochastic Hamiltonian H (cid:48) i prior to time t and the average (cid:104) . . . (cid:105) H (cid:48) i ( t ) over the stochastic Hamiltonian H (cid:48) i at the currenttime t . At this point, we are not yet averaging over U τ,i .13To calculate the quantity on the RHS of Eq. (S.74), we use the reduction (cid:104) H (cid:48) i ( t ) OH (cid:48) j ( t (cid:48) ) (cid:105) H i = 12 α i δ ij δ ( t − t (cid:48) ) (cid:104) W i OW i (cid:105) W , (S.75)where O is any operator and W i is a Wigner random matrix from a Gaussian unitary ensemble supported on ω i andsatisfying (cid:104) W ai,bj W ∗ a (cid:48) i (cid:48) ,b (cid:48) j (cid:48) (cid:105) W = 1 d ω δ aa (cid:48) δ bb (cid:48) δ ii (cid:48) δ jj (cid:48) . (S.76)Since (cid:104) H (cid:48) i ( t ) (cid:105) H i = 0, Eq. (S.74) takes the form ddt P G = − i Tr G (cid:16) { Tr F \ G [ H U , ρ ] , ρ G } (cid:17) − α i (cid:68)(cid:16) Tr G (cid:104) (Tr F \ G [ W i , ρ ]) (cid:105) +Tr G (cid:16) Tr F \ G [ W i , [ W i , ρ ] ρ G (cid:17)(cid:17)(cid:69) circ ,W = − α i (cid:68)(cid:16) Tr G (cid:16) Tr F \ G ( W i ρ )Tr F \ G ( W i ρ ) (cid:17) +Tr G (cid:16) Tr F \ G ( ρW i )Tr F \ G ( ρW i ) (cid:17) − G (cid:16) Tr F \ G ( ρW i )Tr F \ G ( W i ρ ) (cid:17) + Tr G (cid:16) Tr F \ G ( W i ρ ) ρ G (cid:17) +Tr G (cid:16) Tr F \ G ( ρW i ) ρ G (cid:17) − G (cid:16) Tr F \ G ( W i ρW i ) ρ G (cid:17)(cid:17)(cid:69) circ ,W . (S.77)The first term in (S.77) that contains H U vanishes since it acts on single qudits.Averages in Eq. (S.77) can be calculated using the following identities: (cid:68) Tr G (cid:0) Tr F \ G ( W ρ ) (cid:1) (cid:11) W = (cid:68) Tr G (cid:0) Tr F \ G ( ρW ) (cid:1) (cid:11) W = 1 d ω Tr G ∆ ω (cid:16) ρ G ∆ ω (cid:17) , (cid:68) Tr G (cid:16) Tr F \ G ( ρW )Tr F \ G ( W ρ ) (cid:17)(cid:69) W = 1 d ω \ G Tr G ∪ ω ( ρ G ∪ ω ) , (cid:68) Tr G (cid:16) Tr F \ G ( W ρ ) ρ G (cid:17)(cid:69) W = (cid:68) Tr G (cid:16) Tr F \ G ( ρW ) ρ G (cid:17)(cid:69) W = Tr G ( ρ G ) , (cid:68) Tr G (cid:16) Tr F \ G ( W ρW ) ρ G (cid:17)(cid:69) W = 1 d G ∩ ω Tr G \ ω ( ρ G \ ω ) . (S.78)Combining equations (S.77) and (S.78), we obtain ddt P G = − (cid:88) i α i (cid:16) P G + P G ∆ ω i d ω i − P G \ ω i d ω i ∩ G − P G ∪ ω i d ω i \ G (cid:17) , (S.79)which is Eq. (14) in the main text. Continuous monitoring . Let us consider now the dynamics of the system under continuous measurements ofobservables O j ( t ) applied to qudits Ω j in the randomized basis, O j ( t ) = U τ,j ( t ) O j U † τ,j ( t ): ddt | ψ (cid:48) (cid:105) = − (cid:88) j κ j (cid:16) O j ( t ) − (cid:104) ψ (cid:48) | O j ( t ) | ψ (cid:48) (cid:105) (cid:17) + (cid:88) j ξ j ( t ) (cid:112) κ j (cid:16) O j ( t ) − (cid:104) ψ (cid:48) | O j ( t ) | ψ (cid:48) (cid:105) (cid:17) | ψ (cid:48) (cid:105) . (S.80)Here coefficients κ j characterize the strength of the coupling to the measurement apparatus, and ξ j ( t ) are independentunbiased Gaussian random variables satisfying (cid:104) ξ i ( t ) ξ j ( t (cid:48) ) (cid:105) ξ = δ ij δ ( t − t (cid:48) ). The average (cid:104) . . . (cid:105) ξ over the random variables ξ amounts to averaging over the measurement back-action in the realization of the monitoring. We find that, uponaveraging over ξ j , the purities evolve as ddt P G = 2 (cid:88) j κ j (cid:68) Tr G [ ρ G , O j ( t )] (cid:69) circ +2 (cid:88) j κ j (cid:28) Tr G (cid:104) Tr F \ G (cid:0) { O j ( t ) , ρ } − ρO j ( t )) ρ (cid:1)(cid:105) (cid:29) circ , (S.81)where, as defined above, the circuit avareage (cid:104) . . . (cid:105) circ is over the stochastic Hamiltonian H (cid:48) i prior to time t .Equation (S.81) can now be rewritten using the following identities (to ease notation, we suppress time dependenceof O j ): Tr G ([ ρ G , O j ] ) = − (cid:16) Tr G ( ρ G O j ) − Tr G ( ρ G O j ) (cid:17) δ Ω j ∈ G , Tr G (cid:0) Tr F \ G { ρ, O j } (cid:1) = 2 (cid:16) Tr G ( ρ G O j ) + Tr G ( ρ G O j ) (cid:17) δ Ω j ∈ G + 4 Tr G (cid:0) Tr F \ G ( ρO j ) (cid:1) δ Ω j ∈ F \ G , Tr G (cid:0) Tr F \ G ( { ρ, O j } ) ρ G (cid:1) = 2 Tr G (cid:0) ρ G O j (cid:1) δ Ω j ∈ G + 2 Tr L (cid:0) Tr F \ G ( ρO j ) ρ G (cid:1) δ Ω j ∈ F \ G . (S.82)14where δ Ω ∈ G = 1 if Ω ∈ G and δ Ω ∈ G = 0 otherwise.Since the choice of single-qudit rotations U τ does not affect the purities [see Eq. (S.66)], we can additionally averagethe purities over U τ without loss of generality. It is also convenient to split the Liouvillian operator into two parts as ddt P G = (cid:88) j :Ω j ∈ G κ j M ( ρ, O j , G ) + (cid:88) j :Ω j ∈ F \ G κ j M ( ρ, O j , G ) , (S.83)where M and M correspond, respectively, to measurements inside and outside the set G : M ( ρ, O j , G ) = 8 (cid:10) Tr G ( ρ G O j ) − ρO j Tr G (cid:0) ρ G O j (cid:1) + (Tr ρO j ) Tr G ρ G (cid:11) circ ,U τ , (S.84) M ( ρ, O j , G ) = 8 (cid:68) Tr G (cid:0) Tr F \ G ( ρO j ) (cid:1) − ρO j Tr G (cid:0) Tr F \ G ( ρO j ) ρ G (cid:1) + (Tr ρO j ) Tr G ρ G (cid:69) circ ,U τ . (S.85)In the limit of vanishing correlation time, κ j τ →