Unitary circuits of finite depth and infinite width from quantum channels
Sarang Gopalakrishnan Austen Lamacraft









              


































    #1tensor=[rectangle,draw=blue!50,fill=blue!20,thick]
        
        
        Department of Physics and Astronomy, CUNY College of Staten Island, Staten Island NY 10314 USAPhysics Program and Initiative for Theoretical Sciences, The Graduate Center, CUNY, New York NY 10016 USATCM Group, Cavendish Laboratory, University of Cambridge,
        J. J. Thomson Ave., Cambridge CB3 0HE, UK
        We introduce an approach to compute reduced density matrices for local quantum unitary circuits of finite depth and infinite width. Suppose the time-evolved state under the circuit is a matrix-product state with bond dimension D; then the reduced density matrix of a half-infinite system has the same spectrum as an appropriate D × D matrix acting on an ancilla space.
        
        We show that reduced density matrices at different spatial cuts are related by quantum channels acting on the ancilla space. This quantum channel approach allows for efficient  numerical evaluation of the entanglement spectrum and Rényi entropies and their spatial fluctuations at finite times in an infinite system.
        
        We benchmark our numerical method on random unitary circuits, where many analytic results are available, and also show how our approach analytically recovers the behaviour of the kicked Ising model at the self-dual point.
        
        We study various properties of the spectra of the reduced density matrices and their spatial fluctuations in both the random and translation-invariant cases.
        
        Unitary circuits of finite depth and infinite width from quantum channels
            Austen Lamacraft
            December 30, 2023
        =========================================================================
        
        
        
        § INTRODUCTION
        
        
        The dynamics of isolated quantum systems under generic unitary dynamics is one of the basic problems in many-body physics <cit.>; despite considerable recent work, many aspects of this problem are not fully understood. An isolated system, evolving under chaotic dynamics from an initial product state, becomes increasingly entangled over time. At sufficiently late times, any finite-size subsystem of an infinite system is well described by a thermal reduced density matrix, provided the system obeys the eigenstate thermalization hypothesis (ETH) <cit.>; this approach to a thermal local density matrix is called “thermalization.” There is considerable numerical evidence that generic many-body systems obey ETH <cit.>.
        
        
        
        
        We are concerned with the dynamics of the reduced density matrix before the system has fully thermalized. To quantify the thermalization timescale more precisely, recall that the Rényi entropies of a subsystem A are defined in terms of the reduced density matrix ρ_A of A as
        
        
            S^(n)_A(t) = 1/1-nlog[ρ_A(t)^n],
        
        where n is called the Rényi index. The Rényi entropies fully characterize the spectrum of ρ_A.
        
        The general consensus <cit.> is that – unless a system experiences many-body localization <cit.>– the entropies initially obey S^(n)_A(t)∼ v_n t, increasing linearly with time with a growth rate v_n that depends on the Rényi index n. (However, recent results suggest that for n > 1 the growth is sub-linear for generic initial states in the presence of conservation laws <cit.>.) The implication for the spectrum of the reduced density matrix is as follows. Parameterizing the eigenvalues of ρ_A in terms of an “entanglement energy” as λ_i = e^-ϵ_i <cit.>, and introducing the “density of states”ϱ(ϵ), we have
        
        
            S^(n)_A = 1/1-nlog[∫ dϵ ϱ(ϵ) e^-nϵ].
        
        The behaviour S^(n)_A(t)∼ v_n t is then consistent with a density ϱ(ϵ) having the large deviation form
        
        
            ϱ(ϵ) ∼exp[t π(ϵ/t)],
        
        for some function π(η). In the saddle point approximation we find
        
        
            v_n = S^(n)_A/t = π(η_n)-η_n/1-n
        
        where η_n is determined by n=π'(η_n). The growth rates v_n are seen to be related to the function p(η) describing the density of states of the entanglement spectrum by Legendre transformation. Note that the numerator in Eq. (<ref>) vanishes at n=1 due to the normalization of the density matrix, yielding a finite growth rate v_1 for the von Neumann entropy S^(1)_A. Rényi entropies with smaller n grow faster; the reduced density matrix for a subsystem becomes thermal when the slowest Rényi entropy, the so-called “min-entropy”S_∞≡max_i(ϵ_i) (i.e., the largest entanglement eigenvalue) has saturated.
        
        
        
        As well as the rates v_n, the time evolution is characterized by the butterfly velocity v_B at which local perturbations spread <cit.>.
        
        Although the various velocities are generically separate, they coincide in exactly solvable models (such as random circuits in the limit of large local Hilbert space dimension or Clifford gates <cit.> and the self-dual kicked Ising model <cit.>), so the entanglement spectrum evolves in a trivial way.
        
        Away from these non-generic limits, little is known analytically about the entanglement spectrum. A few Rényi entropies can be explicitly computed by mapping the circuit dynamics to random classical partition functions <cit.> but these mappings do not yield the full entanglement spectrum.
        
        The picture that emerges from numerical studies is, however, that the entanglement spectrum has nontrivial structure in generic systems, such as a bandwidth that widens linearly in time; this feature is absent in the exactly solvable limits <cit.>. However, this structure is not fully understood at present. 
        
        In the present work we develop and apply a numerical transfer-matrix approach to compute the structure of the entanglement spectrum for spatially infinite systems at early times. This approach allows us to access some aspects of entanglement for larger subsystems than were studied in previous work: for instance, we are able to compute the spatial fluctuations of entanglement in systems of 10,000 sites. Our approach works formally with infinite systems; we assume that the system of interest has been initialized in a product state and then subjected to a finite-depth quantum circuit (i.e., evolution for a finite time) consisting of on-site or nearest-neighbor quantum gates (see Fig. <ref>). We compute the spectrum of the reduced density matrix of a bipartition into two semi-infinite regions at an arbitrary point. After applying a quantum circuit of finite depth t, the reduced density matrix has rank q^t-1. We will see that the spectrum of the reduced density matrix can be interpreted as that of a q^t-1× q^t-1 matrix R acting on an ancilla space. In addition, reduced density matrices across adjacent cuts are related by quantum channels that are straightforward to construct given the circuit. These quantum channels act as transfer matrices for the entanglement spectrum.
        
        
        
        The quantum-channel perspective is helpful for a number of reasons. First, using this method one can compute the largest few eigenvalues of the entanglement spectrum for circuits that are too deep to permit direct simulation. Second, as the transfer matrix acts on formally infinite systems, spatial fluctuations of entanglement can be directly studied. Third, standard methods from quantum optics such as the stochastic unraveling of quantum channels <cit.> can be applied to simulate the dynamics of entanglement on larger scales than direct simulation permits. Finally, for translation-invariant initial-states evolving under translation-invariant circuits, one can iterate the quantum channel until it converges, and thus extract the entanglement spectrum, free of finite-size effects, without incurring the computational overhead of time-evolving a large system.
        
        These are the issues we explore in the present work. Our main results are as follows.
        
        First, we provide algorithms based on quantum channels for computing the entanglement spectrum and its low-energy tail, as well as for computing some Rényi entropies, by propagating the quantum channel in an ancilla space.
        
        For the exactly solvable model of Ref. <cit.> we analytically demonstrate that the entanglement spectrum is trivial, using the properties of the associated quantum channel.
        
        We compute the behavior of the “low-energy” (large Schmidt rank) tail of the entanglement spectrum for random unitary circuits, random circuits with a conservation law, and translation-invariant integrable circuits, acting on various initial states.
        
        
        
        We compute the distributions of the purity and of the min-entropy;
        for random circuits we find that both the second Rényi entropy and the min-entropy follow Gaussian distributions at the accessible circuit depths (the purity, therefore, follows a log-normal distribution). For circuits with conservation laws or translation-invariant circuits, the nature of this low-entanglement-energy tail is sensitive to the fluctuations in the initial state.
        
        
        We compute the spatial correlations of entanglement, and find that their correlation length grows sub-linearly in time. However, the correlation lengths are short at the accessible times and we are not able to identify a definite exponent.
        
        This paper is organized as follows. In Sec. <ref> we briefly review concepts such as unitary circuits, matrix-product states, and quantum channels, as they apply to the algorithms introduced here. In Sec. <ref> we describe how to construct quantum channels for unitary circuits, and estimate the complexity of various exact and approximate methods for extracting the entanglement spectrum (or some of its moments). In Sec. <ref> we explicitly compute the transfer matrix for the self-dual kicked Ising model, and confirm that all the Rényi entropies coincide in this model, as they are known to <cit.>. In Sec. <ref> we present results for the entanglement spectrum, the distributions of purity and min-entropy, and the growth of spatial correlations in the entanglement, in random unitary circuits and some variants of these. Finally Sec. <ref> summarizes our results and discusses future directions.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § BACKGROUND
        
        
         §.§ Quantum circuits
        
        
        In this work we consider systems that evolve under the application of discrete local unitary gates tiled in the pattern shown in Fig. <ref>. Our approach applies both to Floquet systems in which the gates are applied periodically in time, and to random circuits where each gate is drawn independently and randomly. As with the TEBD algorithm <cit.>, it can be extended to continuous time evolution under strictly local Hamiltonians by discretizing the time evolution, e.g., through a Suzuki-Trotter decomposition.
        
        
        
         §.§ Matrix Product States
        
        
        In this section we introduce some background material on matrix product states (MPS's), which will form an essential part of the following development. More detailed expositions may be found in Refs.<cit.>.
        
        We consider a quantum system described by N identical subsystems with finite Hilbert space dimension q. The (pure) quantum states of the system are therefore defined by vectors in the Hilbert space
        
        
            _N ≡ℂ^q⊗⋯⊗ℂ^q^N times.
        
        A basis of orthonormal product states has the form
        
        
            |s_1:N⟩ = |s_1⟩_q ⊗|s_2⟩_q⋯⊗|s_N⟩_q,
        
        where |s⟩_qs=1,… q is an orthonormal basis for ℂ^q, and we have introduced the sequence notation s_1:N≡ s_1,s_2,… s_N. By taking components of a vector |Ψ⟩∈
        
            Ψ_s_1:N = ⟨s_1:N|Ψ|,⟩
        
        |Ψ⟩ can be regarded as a rank-N tensor with components Ψ_s_1:N.
        
        An MPS is a tensor of the form
        
        
            Ψ_s_1:N = A^(1)_s_1A^(2)_s_2⋯ A^(N)_s_N,
        
        where A^(1)_s, … A^(N)_s ∈ℂ^D_j× D_j+1 are matrices, and the numbers D_jj=1,… N+1 are known as the bond dimensions, with D_1=D_N+1=1. Thus the product Eq. (<ref>) has the form `row vector, product of matrices, column vector', and yields a complex number for each sequence s_1:s_N.
        
        An arbitrary vector |Ψ⟩ may be approximated by an MPS with an error that decreases as D≡max_j D_j increases. We will see that the state arising from applying a unitary circuit to a product state is exactly given by an MPS with D=q^d-1, where d is the depth of the circuit.
        
        A graphical notation for MPS proves to be extremely convenient, and is discussed extensively in Ref. <cit.>. In this notation tensors – such as the matrices or vectors A^(j)– are represented as boxes, with the number of lines or edges entering a box indicating the number of indices the object bears. An edge joining two vertices indicates the (pairwise) contraction of an index. Thus the MPS in Eq. (<ref>) is denoted
        
        
            Ψ_s_1:N =
            [baseline = (X.base),every node/.style=scale=1,scale=.85]
            
            [rounded corners] (1,2) rectangle (2,1);
             (1.5,1.5) node (X) A^(1);
             (1.5,.2) node s_1;
             (2,1.5) – (3,1.5);
            [rounded corners] (3,2) rectangle (4,1);
             (3.5,1.5) node A^(2);
             (3.5,.2) node s_2;
             (4,1.5) – (4.5,1.5);
            
             (1.5,1) – (1.5,.5);  (3.5,1) – (3.5,.5);
            …[baseline = (X.base),every node/.style=scale=1,scale=.85]
             (0.5,1.5) – (1,1.5);
            [rounded corners] (1,2) rectangle (2,1);
             (1.5,1.5) node (X) A^(N);
             (1.5,.2) node s_N;
             (1.5,1) – (1.5,.5);
        
        Here, A^(2) has three lines attached because the collection of matrices A^(2)_s∈ℂ^D× D may be regarded as rank-3 tensor A^(2)∈ℂ^D× D× q. The vertical leg represents the indices s_j that live in the `physical space' while the horizontal lines represent indices in the `bond space' (which we shall also refer to as the `ancilla space').
        
        The squared norm of a vector |Ψ⟩ may be calculated by contracting all physical indices between Ψ_s_1:N and Ψ̅_s_1:N. This has the graphical representation:
        
        
            ⟨Ψ|Ψ|=⟩∑_s_1:NΨ̅_s_1:NΨ_s_1:N =
              [anchor=base, baseline,every node/.style=scale=1,scale=.75]
              [rounded corners] (1,1.25) rectangle (2,0.25);
               (1.5,0.5) node (X) A^(1);
              [rounded corners] (1,-0.25) rectangle (2,-1.25);
               (1.5,-1) node (X) A̅^(1);
               (1.5,0.25) – (1.5,-0.25);
            
               (2,0.75) – (3,0.75);
               (2,-0.75) – (3,-0.75);
            
              [rounded corners] (3,1.25) rectangle (4,0.25);
               (3.5,0.5) node (X) A^(2);
              [rounded corners] (3,-0.25) rectangle (4,-1.25);
               (3.5,-1) node (X) A̅^(2);
               (3.5,0.25) – (3.5,-0.25);
            
               (4,0.75) – (4.5,0.75);
               (4,-0.75) – (4.5,-0.75);
            
              …[anchor=base, baseline,every node/.style=scale=1,scale=.75]
               (0.5,0.75) – (1,0.75);
               (0.5,-0.75) – (1,-0.75);
            
              [rounded corners] (1,1.25) rectangle (2,0.25);
               (1.5,0.5) node (X) A^(N);
              [rounded corners] (1,-0.25) rectangle (2,-1.25);
               (1.5,-1) node (X) A̅^(N);
               (1.5,0.25) – (1.5,-0.25);
        
        
        More generally, the reduced density matrix for the leftmost n subsystems that arises from a pure state by tracing over the remaining subsystems is denoted
        
        
            ρ_s_1:n,s'_1:n≡∑_s_n+1:NΨ̅_s_1:ns_n+1:NΨ_s'_1:ns_n+1:N=
                [anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
                [rounded corners] (1,1.75) rectangle (2,0.75);
                 (1.5,1) node (X) 1;
                [rounded corners] (1,-0.75) rectangle (2,-1.75);
                 (1.5,-1.5) node (X) 1̅;
            
                 (1.5,0.75) – (1.5,0.4);
                 (1.5,0.2) node s_1;
            
                 (1.5,-0.75) – (1.5,-0.4);
                 (1.5,-0.3) node s'_1;
            
                 (2,1.25) – (2.5,1.25);
                 (2,-1.25) – (2.5,-1.25);
            
                …[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
            
                 (0.5,1.25) – (1,1.25);
                 (0.5,-1.25) – (1,-1.25);
            
                [rounded corners] (1,1.75) rectangle (2,0.75);
                 (1.5,1) node (X) n;
                [rounded corners] (1,-0.75) rectangle (2,-1.75);
                 (1.5,-1.5) node (X) n̅;
            
                 (1.5,0.75) – (1.5,0.4);
                 (1.5,0.2) node s_n;
            
                 (1.5,-0.75) – (1.5,-0.4);
                 (1.5,-0.3) node s'_n;
            
                 (2,1.25) – (3,1.25);
                 (2,-1.25) – (3,-1.25);
            
                [rounded corners] (3,1.75) rectangle (4,0.75);
                 (3.5,1) node (X) n+1;
                [rounded corners] (3,-0.75) rectangle (4,-1.75);
                 (3.5,-1.5) node (X) n+1;
            
                 (3.5,0.75) – (3.5,-0.75);
            
                 (4,1.25) – (4.5,1.25);
                 (4,-1.25) – (4.5,-1.25);
            
                ⋯[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
            
                   (0.5,1.25) – (1,1.25);
                   (0.5,-1.25) – (1,-1.25);
            
                [rounded corners] (1,1.75) rectangle (2,0.75);
                 (1.5,1) node (X) N;
                [rounded corners] (1,-0.75) rectangle (2,-1.75);
                 (1.5,-1.5) node (X) N̅;
            
                 (1.5,0.75) – (1.5,-0.75);
                ,
        
        where for simplicity we denote A^(j) by j and a'^(j) by j̅.  The above expressions may also be written
        
        
        
            ρ_s_1:n,s'_1:n =
              [anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) 1;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) 1̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_1;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_1;
            
               (2,1.25) – (2.5,1.25);
               (2,-1.25) – (2.5,-1.25);
            
              …[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
            
               (0.5,1.25) – (1,1.25);
               (0.5,-1.25) – (1,-1.25);
            
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) n;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) n̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_n;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_n;
            
            
               (3,0) circle (.5);
               (3,-0.1) node (X) R^(n);
               (3,0.5) edge[out=90,in=0] (2,1.25);
               (3,-0.5) edge[out=270,in=0] (2,-1.25);
            
              ,
        
        where the hermitian matrices R^(j)∈ℂ^D_j+1× D_j+1 are defined by
        
        
            R^(j-1) = ∑_s A^(j)_s R^(j) A^(j)†_s,      j=n+1,… N,
        
        and R^(N)=1.
        
        
        
          §.§.§ Canonical Forms
        
        
        The MPS representation of |Ψ⟩ has a redundancy sometimes referred to as `gauge freedom'. For a set X_jj=1,… N-1 of invertible matrices, the transformation
        
        
            A^(1)_s    → A^(1)_s X_1^-1,
            
              A^(j)_s    → X_j-1A^(j)_s X_j^-1,     j=2,… N-1
            
              A^(N)_s    → X_N-1A^(N)_s ,
        
        leaves Ψ_s_1:N unchanged. Further conditions may be imposed to reduce this redundancy <cit.>. We will be concerned with matrices in left canonical form, satisfying the condition
        
        
            ∑_s A^(j)†_s A^(j)_s   = _D_j+1
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=0.75,scale=.55]
              (1,-1.5) edge[out=180,in=180] (1,1.5);
             [rounded corners] (1,2) rectangle (2,1);
             [rounded corners] (1,-1) rectangle (2,-2);
              (1.5,1) – (1.5,-1);
              (1.5,1.5) node A^(j);
              (1.5,-1.5) node A̅^(j);
              (2,1.5) – (2.5,1.5);  (2,-1.5) – (2.5,-1.5);
                =
             [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=0.750,scale=.55]
              (1,-1.5) edge[out=180,in=180] (1,1.5);
             .
        
        Note that D_2=q is necessary for A^(1) to be placed in left canonical form, but this places no restriction on the state.
        
        Analogously, matrices in right canonical form satisfy
        
        
            ∑_s A^(j)_sA^(j)†_s   = _D_j
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=0.75,scale=.55]
              (2,-1.5) edge[out=0,in=0] (2,1.5);
             [rounded corners] (1,2) rectangle (2,1);
             [rounded corners] (1,-1) rectangle (2,-2);
              (1.5,1) – (1.5,-1);
              (1.5,1.5) node A^(j);
              (1.5,-1.5) node A̅^(j);
              (0.5,1.5) – (1,1.5);  (0.5,-1.5) – (1,-1.5);
                =
             [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=0.750,scale=.55]
              (1,-1.5) edge[out=0,in=0] (1,1.5);
             .
        
        (requiring D_N=q).
        
        One benefit of the canonical forms is that they may be contracted “automatically”. For example, in terms of an MPS in right canonical form the reduced density matrix ρ_s_1:n,s'_1:n in Eq. (<ref>) takes the form
        
        
            ρ_s_1:n,s'_1:n =
              [anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) 1;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) 1̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_1;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_1;
            
               (2,1.25) – (2.5,1.25);
               (2,-1.25) – (2.5,-1.25);
            
              …[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
            
               (0.5,1.25) – (1,1.25);
               (0.5,-1.25) – (1,-1.25);
            
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) n;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) n̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_n;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_n;
            
               (2,-1.25) edge[out=0,in=0] (2,1.25);
            
              .
        
        A second benefit – which is more relevant for us – is that the spectrum of the reduced density matrix ρ_s_1:n,s'_1:n coincides with the spectrum of R^(n) for an MPS in left canonical form. This may be seen by introducing the spectral representation
        
        
            R^(n) = ∑_αλ_α r_α r_α^†,
        
        in terms of the eigenvalues λ_α and eigenvectors r_α of R^(n). Upon substitution into Eq. (<ref>) this yields a spectral representation for ρ_s_1:n,s'_1:n
        
            ρ_s_1:n,s'_1:n = ∑_αλ_α[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) 1;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) 1̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_1;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_1;
            
               (2,1.25) – (2.5,1.25);
               (2,-1.25) – (2.5,-1.25);
            
              …[anchor=base, baseline,every node/.style=scale=0.9,scale=.75]
            
               (0.5,1.25) – (1,1.25);
               (0.5,-1.25) – (1,-1.25);
            
              [rounded corners] (1,1.75) rectangle (2,0.75);
               (1.5,1) node (X) n;
              [rounded corners] (1,-0.75) rectangle (2,-1.75);
               (1.5,-1.5) node (X) n̅;
            
               (1.5,0.75) – (1.5,0.4);
               (1.5,0.2) node s_n;
            
               (1.5,-0.75) – (1.5,-0.4);
               (1.5,-0.3) node s'_n;
            
            
               (2,1.25) – (2.5,1.25);
               (2,-1.25) – (2.5,-1.25);
            
               (3,1.25) circle (.5);
               (3,1.1) node (X) r_α;
            
               (3,-1.25) circle (.5);
               (3,-1.3) node (X) r̅_α;
            
              ,
        
        in terms of vectors in _n that are orthonormal by the left canonical condition <cit.>.
        
        A more direct way of seeing that the spectra of R^(n) and ρ_A coincide, for MPS's in left canonical form, is shown in Fig. <ref>. The moments of the spectrum are given by matrix products of the form ρ_A^k; one can use the left canonical property (<ref>) to eliminate the “boxes” pairwise and arrive at the result ρ_A^k = [ (R^(n))^k]. Since all moments coincide, ρ_A and R^(n) must have the same spectrum.
        
        
        
         §.§ Quantum Channels
        
        
        The above formalism has a natural interpretation in terms of quantum channels, or completely positive trace preserving (CPTP) maps <cit.>. In this interpretation the matrices R^(j) are regarded as a sequence of density matrices that represent mixed states in the bond space. The definition of R^(j) given in Eq. (<ref>) guarantees that the maps from one R^(j) to the next are completely positive (Choi's theorem), while the left canonical condition Eq.(<ref>) ensures that they are trace preserving. Therefore, contracting a physical leg (i.e., moving the entanglement cut in real space) amounts to applying a CPTP map to the ancilla. In this context the A^(j)_s are known as Kraus operators.
        
        From now on we assume without loss of generality that all D_j=D, which may be achieved by padding the matrices with zeros. Square matrices A^(j)_s∈ℂ^D× D satisfying the left canonical condition may be parameterized in terms of unitary matrices U_j∈(qD) as
        
        
            A^(j)_s,ab = ⟨s|_q⟨a|_D U_j|0⟩_q|b⟩_D.
        
        Physically, this corresponds to the amplitude for the following process: prepare the physical subsystem j in a fixed state |0⟩_q, and then act on the state |0⟩_q|b⟩_D of the subsystem j and ancilla with a unitary, arriving in state |s⟩_q|a⟩_D. While any CPTP map may be presented in this form <cit.>, we will see that this is precisely how quantum channels arise in the case of unitary circuits.
        
        
        
        § QUANTUM CHANNELS FROM UNITARY CIRCUITS
        
        
        This section is organized as follows. We first discuss how to slice up a unitary circuit into an MPS where all matrices are in canonical form. This immediately gives us a quantum channel that propagates the ancilla-space density matrix R in the spatial direction. We then discuss two controlled but approximate ways of implementing the quantum channel for larger subsystems: first, an approach for computing the entanglement spectrum based on approximating the ancilla-space density matrix R as a lower-rank object; and second, an approach for computing the purity, specifically, by a stochastic unraveling of the quantum channel (i.e., by sampling “quantum trajectories” in ancilla space <cit.>).
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ From unitary circuits to canonical-form MPS's
        
        
        We now introduce the main idea behind our approach. A planar unitary circuit that starts from a product state may be presented as an MPS by slicing it into strips in an arbitrary way.
        
        Each slice j is associated with q^2 matrices A^(j)_s_1,s_2 indexed by two physical indices. The dimension of the ancilla is D=q^d-1.
        
        
        
        
        
        
        
        For a general decomposition of a circuit, the resulting matrices A^(j)_s_1,s_2 will not be in the appropriate canonical form, so the spectra of R^(j) and ρ_A will not coincide.
        
        
        
        An MPS in canonical form is obtained for slices along the south-west to north-east diagonal
        
        
            A^(j)_s_1,s_2;a_1:3,b_1:3=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            1-3U_12-1.5U_230U_341.5U_4 (0.5,-2.75) edge[out=90,in=0] (0,-2.25);
             (1.5,-1.25) edge[out=90,in=0] (1,-0.75);
             (2.5,0.25) edge[out=90,in=0] (2,0.75);
            
             (2.5,-1.75) edge[out=-90,in=180] (3,-2.25);
             (3.5,-0.25) edge[out=-90,in=180] (4,-0.75);
             (4.5,1.25) edge[out=-90,in=180] (5,0.75);
            
            -4-30.5-4-1.51.5-1.502.502.53.51.52.54.5 (0.5,-4.2) circle (.2);
             (1.5,-4.2) circle (.2);
            
             (3.5,2.7) node (X) s_1;
             (4.5,2.7) node (X) s_2;
            
             (-0.3,-2.25) node (X) a_3;
             (0.7,-0.75) node (X) a_2;
             (1.7,0.75) node (X) a_1;
            
             (3.3,-2.25) node (X) b_3;
             (4.3,-0.75) node (X) b_2;
             (5.3,0.75) node (X) b_1;
            
            .
        
        A graphical proof is straightforward, since
        
        
            ∑_s_1,s_2(A^(j)†_s_1,s_2 A^(j)_s_1,s_2)_b'_1:3,b_1:3=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            1-3U_12-1.5U_230U_341.5U_417.5U^†_126U^†_234.5U^†_343U^†_4 (2.5,-1.75) edge[out=-90,in=180] (3,-2.25);
             (3.5,-0.25) edge[out=-90,in=180] (4,-0.75);
             (4.5,1.25) edge[out=-90,in=180] (5,0.75);
            
             (2.5,6.25) edge[out=90,in=180] (3,6.75);
             (3.5,4.75) edge[out=90,in=180] (4,5.25);
             (4.5,3.25) edge[out=90,in=180] (5,3.75);
            
            -48.50.5-48.51.5-1.562.504.53.51.534.5 (0.5,-4.2) circle (.2);
             (1.5,-4.2) circle (.2);
            
             (0.5,8.7) circle (.2);
             (1.5,8.7) circle (.2);
            
             (3.3,-2.25) node (X) b_3;
             (4.3,-0.75) node (X) b_2;
             (5.3,0.75) node (X) b_1;
            
             (3.3,6.75) node (X) b'_3;
             (4.3,5.25) node (X) b'_2;
             (5.3,3.75) node (X) b'_1;
            
            .
        
        The unitarity of the constituent blocks
        
        
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            10U_111.5U^†_1-12.50.5-12.51.5 (0.5,-1.5) node (X) b_1;
             (1.5,-1.5) node (X) b_2;
            
             (0.5,3) node (X) b'_1;
             (1.5,3) node (X) b'_2;
            
            =
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            -12.50.5-12.51.5 (0.5,-1.5) node (X) b_1;
             (1.5,-1.5) node (X) b_2;
            
             (0.5,3) node (X) b'_1;
             (1.5,3) node (X) b'_2;
            
            
            =δ_b_1,b'_1δ_b_2,b'_2
        
        guarantees that the left canonical condition
        
        
            ∑_s_1,s_2 A^(j)†_s_1,s_2 A^(j)_s_1,s_2 = _D
        
        is satisfied. Choosing the south-east to north-west diagonal gives an MPS in right canonical form.
        
        
        
         §.§ Applying the quantum channel
        
        
        We have shown that a unitary circuit may be represented as an MPS in canonical form. This observation may be used to efficiently compute the ancilla density matrices R^(j), whose spectrum coincides with that of the reduced density matrix ρ_s_1:j,s'_1:j.
        
        Recall that the R^(j) are defined by
        
        
            R^(j-1) = ∑_s_1,s_2 A^(j)_s_1,s_2 R^(j) A^(j)†_s_1,s_2,      j=n+1,… N.
        
        After cutting along the SW-NE diagonal, R^(j-1)_a_1:d-1,a'_1:d-1  has the graphical representation
        
        
            R^(j-1)_a_1:d-1,a'_1:d-1=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            1-3U_12-1.5U_230U_341.5U_417.5U^†_126U^†_234.5U^†_343U^†_4
             (2.5,-1.75) edge[out=-90,in=180] (3,-2.25);
             (3.5,-0.25) edge[out=-90,in=180] (4,-0.75);
             (4.5,1.25) edge[out=-90,in=180] (5,0.75);
            
            
             (2.5,6.25) edge[out=90,in=180] (3,6.75);
             (3.5,4.75) edge[out=90,in=180] (4,5.25);
             (4.5,3.25) edge[out=90,in=180] (5,3.75);
            
            
             (0.5,-2.75) edge[out=90,in=0] (0,-2.25);
             (1.5,-1.25) edge[out=90,in=0] (1,-0.75);
             (2.5,0.25) edge[out=90,in=0] (2,0.75);
            
            
             (0.5,7.25) edge[out=-90,in=0] (0,6.75);
             (1.5,5.75) edge[out=-90,in=0] (1,5.25);
             (2.5,4.25) edge[out=-90,in=0] (2,3.75);
            
            -4-30.5-4-1.51.57.58.50.568.51.5-1.502.54.562.504.53.51.534.5 (0.5,-4.2) circle (.2);
             (1.5,-4.2) circle (.2);
            
             (0.5,8.7) circle (.2);
             (1.5,8.7) circle (.2);
            
            
            
            
            
            
            
             (-0.3,-2.25) node (X) a_3;
             (0.7,-0.75) node (X) a_2;
             (1.7,0.75) node (X) a_1;
            
             (1.7,3.75) node (X) a'_1;
             (0.7,5.25) node (X) a'_2;
             (-0.3,6.75) node (X) a'_3;
            
            
             (7,2.25) node (X) R^(j);
            [rounded corners] (6,3) rectangle (8,1.5);
            
            
            
             (5,0.75) edge[out=0,in=-90] (6.5,1.5);
             (4,-0.75) edge[out=0,in=-90] (7,1.5);
             (3,-2.25) edge[out=0,in=-90] (7.5,1.5);
            
             (5,3.75) edge[out=0,in=90] (6.5,3);
             (4,5.25) edge[out=0,in=90] (7,3);
             (3,6.75) edge[out=0,in=90] (7.5,3);
            
            .
        
        This expression may be simplified somewhat by noting that the topmost unitary (U_4 in the above example) may be eliminated to give
        
        
            R^(j-1)_a_1:d-1,a'_1:d-1=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            1-2U_12-0.5U_231U_316.5U^†_125U^†_233.5U^†_3
             (2.5,-0.75) edge[out=-90,in=180] (3,-1.25);
             (3.5,0.75) edge[out=-90,in=180] (4,0.25);
             (4.5,2.25) edge[out=-90,in=180] (5,1.75);
            
            
             (2.5,5.25) edge[out=90,in=180] (3,5.75);
             (3.5,3.75) edge[out=90,in=180] (4,4.25);
             (4.5,2.25) edge[out=90,in=180] (5,2.75);
            
            
             (0.5,-1.75) edge[out=90,in=0] (0,-1.25);
             (1.5,-0.25) edge[out=90,in=0] (1,0.25);
             (2.5,1.25) edge[out=90,in=0] (2,1.75);
            
            
             (0.5,6.25) edge[out=-90,in=0] (0,5.75);
             (1.5,4.75) edge[out=-90,in=0] (1,4.25);
             (2.5,3.25) edge[out=-90,in=0] (2,2.75);
            
            -3-20.5-3-0.51.56.57.50.557.51.5-0.512.53.552.513.53.5 (0.5,-3.2) circle (.2);
             (1.5,-3.2) circle (.2);
            
             (0.5,7.7) circle (.2);
             (1.5,7.7) circle (.2);
            
            
            
            
            
            
            
             (-0.3,-1.25) node (X) a_3;
             (0.7,0.25) node (X) a_2;
             (1.7,1.75) node (X) a_1;
            
             (1.7,2.75) node (X) a'_1;
             (0.7,4.25) node (X) a'_2;
             (-0.3,5.75) node (X) a'_3;
            
            
             (7,2.25) node (X) R^(j);
            [rounded corners] (6,3) rectangle (8,1.5);
            
            
             (5,1.75) edge[out=0,in=-90] (6.5,1.5);
             (4,0.25) edge[out=0,in=-90] (7,1.5);
             (3,-1.25) edge[out=0,in=-90] (7.5,1.5);
            
             (5,2.75) edge[out=0,in=90] (6.5,3);
             (4,4.25) edge[out=0,in=90] (7,3);
             (3,5.75) edge[out=0,in=90] (7.5,3);
            
            .
        
        The algorithm for applying the quantum channel is therefore:
        
        
          *  Trace over the first index of R^(j)
            R^(j)_a_1:D,a'_1:D→∑_aR^(j)_a a_1:D-1,a a'_1:D-1
        
          This reduces the number of indices of R^(j) to 2(d-2), or q^2(d-2) components.
        
          
          *  Apply the unitaries U_d-1 and U̅_d-1. This increases the rank of the resulting tensor back to 2(d-1).
        
          
          *  Continue applying unitaries from the “middle out” for j=d-2,… 2.
        
          
          *  Apply U_1 and U^†_1, with the outer indices fixed.
        
        
        
        This process involves O(d) steps of matrix multiplication, where the matrices are of size O(q^d). Thus if one directly applies the channel to a density matrix, the overall complexity is O(d q^2(d-1)). This is better than the naive O(d q^3(d-1)) because each of the unitaries that make up the channel is a sparse matrix. A Python implementation of the algorithm is available at <https://github.com/AustenLamacraft/ruc>.
        
        Since our quantum channel is constructed from a diagonal cut through the unitary circuit, there will be edge effects in a rectangular circuit of finite width. Our approach is well suited to infinite width circuits: the channel is applied repeatedly to a random initial density matrix until a steady state density matrix is approached (for translationally invariant circuits) or a stationary distribution (for random circuits). This typically occurs on the scale of a number of steps roughly equal to the depth of the circuit.
        
        
        
         §.§ Low-rank approximation of R^(n)
        
        
        Away from fine-tuned points (see Sec. <ref>), finite-depth local unitary circuits give rise to entanglement spectra that are very broad; thus, the vast majority of the eigenstates of ρ_A are close to zero and do not contribute to Rényi entropies with n ≥ 1 <cit.>. This observation is implicit in the fact that different S_n have different growth rates <cit.>. This fact allows us to propagate R^(n) with negligible error using far fewer than q^d-1 basis states. We proceed as follows. We approximate
        
        
        R^(n)≈∑_k = 1^K λ_k |k⟩⟨ k|,
        
        where λ_k are the K largest eigenvalues of R and |k⟩ are the associated eigenvectors. We renormalize all the eigenvalues to preserve the trace. We now evolve each |k⟩ under each “leg” of the quantum channel. This evolution is efficient because each of the unitaries is a very sparse matrix. At the end of this process we have the expression
        
        
        R^(n+1)≈∑_k = 1^K ∑_i = 1^4 λ_k |ϕ_ik⟩⟨ϕ_ik|,
        
        where |ϕ_ik⟩ are not mutually orthogonal or normalized, but nevertheless span a 4K-dimensional space. Eq. (<ref>) is the ancilla density matrix that would result from one step of the quantum channel applied to the approximate density matrix (<ref>). It is a legitimate density matrix, since Eq. (<ref>) was. Now we can repeat this process by approximating R^(n+1) with its top K eigenvectors, renormalizing, propagating, and so on.
        
        When K is sufficiently small, the complexity of evolving the K top eigenvectors scales as O(K d q^d-1), since the unitary gates are individually sparse matrices. The diagonalization step scales as O(K^3), meanwhile. Benchmarking our results against exact diagonalization at small sizes (and against exact results for random unitary circuits at arbitrary sizes) we find that keeping 100 eigenvectors suffices to capture the quantities that are of interest here—mainly, the purity and min-entropy, and their fluctuations.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ Trajectory approach for computing the purity
        
        In this section we describe an approach based on mapping vectors in the ancilla space rather than density matrices. Formally, this method is equivalent to the trajectory approach developed to analyze master equations in quantum optics <cit.>. Applying unitaries to a vector is an O(d q^(d-1)) operation because the unitaries are sparse. While such an approach is evidently attractive, we will see that there is a trade-off in terms of the number of times the matrices A^(j)_s_1,s_2 must be applied.
        
        The ancilla density matrix can be defined in terms of averages over trajectories in the physical indices in the following way. Starting from the Kraus form
        
        
            R^(j) = ∑_s A^(j)_s R^(j-1) A^(j)†_s,      j=n+1,… N,
        
        
        In the case of unitary circuits the index s is a composite: s=(s_1,s_2), and we have changed the indexing of slices so that indices increase going right to left. We see that L updates correspond to summing over trajectories in the physical indices of length L
        
            R^(L) = ∑_s_1:s_L∈{ℤ_q}^L A^(L)_s_L⋯  A^(1)_s_1 R^(0) A^(1)†_s_1⋯ A^(L)†_s_L.
        
        If we start from a pure state R^(0)=|ψ_0⟩⟨ψ_0| we can write this as an average over trajectories with uniform distribution
        
        
            R^(L) = q^L _s_1:L∼uniform[|ψ̃_s_1:L⟩⟨ψ̃_s_1:L|],
        
        where the vectors |ψ_s_1:L⟩ are defined as
        
        
            |ψ̃_s_1:L⟩ = A^(L)_s_L⋯ A^(1)_s_1|ψ_0⟩.
        
        These vectors are unnormalized. The normalization factors
        
        
            p(s_1:s_L)≡⟨ψ̃_s_1:L|ψ̃_s_1:L|⟩
        
        are a normalized probability distribution over trajectories by virtue of the left canonical condition Eq. (<ref>). Denoting the normalized vectors as |ψ_s_1:L⟩ we can express the ancilla density matrix as
        
        
            R^(L)=_s_1:N∼ p(·) [|ψ_s_1:L⟩⟨ψ_s_1:L|].
        
        As a trajectory increases in length, the normalization factors are updated according to
        
        
            p(s_1:s_L) = p(s_1:s_L-1)⟨ψ_L-1|A^(L)†_s_LA^(L)_s_L|ψ_L-1|,⟩
        
        so that the second factor may be interpreted as a conditional probability
        
        
            p(s_L|s_1:L-1)=⟨ψ_L-1|A^(L)†_s_LA^(L)_s_L|ψ_L-1|,⟩
        
        Eq. (<ref>) expresses the ancilla density matrix in terms of vectors, but it requires an average over trajectories. The downside of this approach is that evaluating R will require roughly γ^-1 trajectories, where γ is the purity of a half-infinite region, which sets the approximate rank of the reduced density matrix. However, this approach lends itself to parallelization while the channel based approach does not.
        
        For a translationally invariant system, and assuming this random process is ergodic, we can substitute an average over the length of single long trajectory in the L→∞ limit.
        
        
            R = lim_L→∞1/L∑_l=1^L |ψ_s_1:l⟩⟨ψ_s_1:l|.
        
        In practice, long runs and multiple trajectories are used.
        
        O(D) evaluation of matrix products is not much use if we still need O(D^3) evaluation of the spectrum of R^(L) or O(D^2) evaluation of the purity.  However, we can access the purity without dealing with R^(L) directly using
        
        
            γ_L=[R^(L)2]=_s_1:N,t_1:N∼ p(·)|⟨ψ_t_1:L|ψ_s_1:L||⟩^2,
        
        which follows from Eq. (<ref>). This formula expresses the purity as the average fidelity over pairs of trajectories. In a high purity state the ancilla vectors stay close to each other as they evolve over different trajectories, whereas in a highly entangled state different trajectories explore different regions of ancilla space.
        
        Note that the expectations discussed in this section are unrelated to any random variables that may form part of the specification of the circuit. Evaluating the average purity, for example, would require an additional average of Eq. (<ref>) over these variables.
        
        
        
        § EXACTLY SOLVABLE EXAMPLE: SELF-DUAL KICKED ISING MODEL
        
        
        In order to illustrate the utility of the formalism introduced in the previous section, we now turn to an example of a unitary circuit in which the entanglement spectrum can be determined analytically. This is the kicked Ising model at the self-dual point discussed in two recent papers <cit.>. The kicked Ising model describes the evolution of a system of L spin-1/2 subsystems (qubits) for an integer time t by the unitary operator (U_KI)^t, where U_KI=K I_𝐡 is composed of the two unitaries
        
        
            I_𝐡 = e^-iH_I[𝐡],    
              K = e^-iH_K,
        
        where
        
        
            H_I[𝐡]   =∑_j=1^L[J Z_j Z_j+1 + h_j Z_j]
            
            H_K   = b∑_j=1^L X_j,
        
        and (X_j,Y_j,Z_j) are the Pauli matrices for spin j. H_I[𝐡] is the classical Ising model with arbitrary longitudinal fields h_j, while H_K describes a transverse field.
        
        In Ref. <cit.> the growth of the entanglement entropies was found exactly for some particular initial product states at the special `self-dual' values
        
        
            |J|=|b|=π/4.
        
        For a region A of size N, the authors found that when starting from an arbitrary product state in the Z_j basis the Rényi entropies (<ref>) to be exactly given by
        
        
            lim_L→∞ S^(n)_A(t) =min(2t-2,N)log 2,
        
        independent of Rényi index n. The interpretation in terms of the entanglement spectrum is striking: there are 2^min(2t-2,N) eigenvalues equal to 2^-min(2t-2,N) and the rest are zero.
        
        We now show how our quantum channel approach may be used to derive the corresponding result for the case of a semi-infinite interval
        
        
            lim_L→∞ S^(n)_A(t) =(t-1)log 2.
        
        We can present the unitary (U_KI)^t as a unitary circuit of depth t, where the layers alternate between unitaries operating between spins 2j-1 and 2j, and between spins 2j and 2j+1. A variety of decompositions are available. For reasons that will become clear, we choose the following:
        
        
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.75]
              00U_12=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.75]
            
              [rounded corners] (0,2) rectangle (3,1);
               (1.5,1.5) node (X) ;
            
              [rounded corners] (0,-2) rectangle (3,-1);
               (1.5,-1.5) node (X) ;
            
              0.502.5010.50.5-1-0.50.510.52.5-1-0.52.522.50.522.52.5-2-2.50.5-2-2.52.5 (0.5,2.8) node (X) a;
               (2.5,2.8) node (X) b;
            
               (0.5,-2.8) node (X) c;
               (2.5,-2.8) node (X) d;
            
            ,
        
        where the one qubit () and two qubit () gates have the form
        
        
            = exp[-i b X]
               = exp[-iJ Z_1 Z_2 -i (h_1 Z_1 + h_2 Z_2)/2].
        
        In the Z basis, the elements of U_12 are
        
        
            (U_12)_ab,cd =-sin 2b/2exp(-iJ [ab+cd]-iJ̃[ac+bd])
                ×exp(-ih_1[a+c]/2-ih_2[b+d]/2)
        
        where a,b,c,d∈{1,-1} and
        
        
            J̃ = -π/4-i/2logtan b.
        
        
        Eq. (<ref>) shows that one could equally well consider a circuit in which contractions are made horizontally
        
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.75]
              00U_12=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.75]
            
              [rounded corners] (0,1.5) rectangle (1,-1.5);
               (0.5,0) node (X) Ĩ;
            
              [rounded corners] (3,1.5) rectangle (4,-1.5);
               (3.5,0) node (X) Ĩ;
            
              20.75K̃2-0.75K̃11.50.752.530.7511.5-0.752.53-0.75-0.500.7544.50.75-0.50-0.7544.5-0.75 (-0.8,0.75) node (X) a;
               (4.8,0.75) node (X) b;
            
               (-0.8,-0.75) node (X) c;
               (4.8,-0.75) node (X) d;
            
            ,
        
        where Ĩ is the Ising part involving the dual coupling J̃ from Eq. (<ref>) and K̃ are the transverse field (kicking) unitaries involving the dual field
        
        
            b̃ = -π/4-i/2logtan J.
        
        
        The matrix (U_12)_ab,cd is unitary, but the matrix Ũ_12 with elements (Ũ)_ab,cd=(U_12)_ac,bd is not except at the self-dual points Eq. (<ref>). The unitarity of Ũ_12 has the graphical representation
        
        
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
              11.5U_121-1.5U_12^† (1.5,2) edge[out=90,in=90] (2.5,2);
               (2.5,2) edge[out=-90,in=90] (2.5,-2);
               (1.5,-2) edge[out=-90,in=-90] (2.5,-2);
            
              
              
              2-21.5 (0,2.5) edge[out=0,in=90] (0.5,2);
               (0,0.5) edge[out=0,in=-90] (0.5,1);
            
               (0,-2.5) edge[out=0,in=-90] (0.5,-2);
               (0,-0.5) edge[out=0,in=90] (0.5,-1);
            
              120.5-1-20.5 (-0.3,2.5) node (X) a;
               (-0.3,0.5) node (X) b;
            
               (-0.3,-2.5) node (X) a';
               (-0.3,-0.5) node (X) b';
              =
              [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
               (0,2.5) edge[out=0,in=90] (1,2);
               (0,0.5) edge[out=0,in=90] (0.5,0);
            
               (0,-2.5) edge[out=0,in=-90] (1,-2);
               (0,-0.5) edge[out=0,in=-90] (0.5,0);
            
              2-21 (-0.3,2.5) node (X) a;
               (-0.3,0.5) node (X) b;
            
               (-0.3,-2.5) node (X) a';
               (-0.3,-0.5) node (X) b';
              
              = δ_aa'δ_bb'.
        
        The unitarity of Ũ_12 has the interesting consequence that the SW-NE MPS is in left and right canonical form, so that
        
        
            ∑_s A^(j)_sA^(j)†_s = _2^t-1.
        
        To see this, we first give the graphical representation of the left hand side of Eq. (<ref>)
        
            ∑_s_1,s_2(A^(j)_s_1,s_2A^(j)^†_s_1,s_2)_a_1:3,a'_1:3=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            11.5U_1223U_2334.5U_3446U_451-1.5U_12^†2-3U_23^†3-4.5U_34^†4-6U_45^† (3.5,6.5) edge[out=90,in=90] (6.5,6.5);
             (6.5,6.5) edge[out=-90,in=90] (6.5,-6.5);
             (3.5,-6.5) edge[out=-90,in=-90] (6.5,-6.5);
            
             (4.5,6.5) edge[out=90,in=90] (5.5,6.5);
             (5.5,6.5) edge[out=-90,in=90] (5.5,-6.5);
             (4.5,-6.5) edge[out=-90,in=-90] (5.5,-6.5);
            
            6.5-6.54.56.5-6.53.55-52.53.50.51.5-3.5-0.51.520.50.5-2-0.50.5 (0.5,0.3) circle (.2);
             (1.5,0.3) circle (.2);
             (0.5,-0.3) circle (.2);
             (1.5,-0.3) circle (.2);
            
             (0,2.5) edge[out=0,in=90] (0.5,2);
             (1,4) edge[out=0,in=90] (1.5,3.5);
             (2,5.5) edge[out=0,in=90] (2.5,5);
            
             (0,-2.5) edge[out=0,in=-90] (0.5,-2);
             (1,-4) edge[out=0,in=-90] (1.5,-3.5);
             (2,-5.5) edge[out=0,in=-90] (2.5,-5);
            
             (-0.3,2.5) node (X) a_1;
             (0.7,4) node (X) a_2;
             (1.7,5.5) node (X) a_3;
            
             (-0.3,-2.5) node (X) a_1';
             (0.7,-4) node (X) a_2';
             (1.7,-5.5) node (X) a_3';
            
            =
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            11.5U_1223U_2334.5U_341-1.5U_12^†2-3U_23^†3-4.5U_34^† (3.5,5) edge[out=90,in=90] (4.5,5);
             (4.5,5) edge[out=-90,in=90] (4.5,-5);
             (3.5,-5) edge[out=-90,in=-90] (4.5,-5);
            
            
            
            5-53.55-52.53.50.51.5-3.5-0.51.520.50.5-2-0.50.5 (0.5,0.3) circle (.2);
             (1.5,0.3) circle (.2);
             (0.5,-0.3) circle (.2);
             (1.5,-0.3) circle (.2);
            
             (0,2.5) edge[out=0,in=90] (0.5,2);
             (1,4) edge[out=0,in=90] (1.5,3.5);
             (2,5.5) edge[out=0,in=90] (2.5,5);
            
             (0,-2.5) edge[out=0,in=-90] (0.5,-2);
             (1,-4) edge[out=0,in=-90] (1.5,-3.5);
             (2,-5.5) edge[out=0,in=-90] (2.5,-5);
            
             (-0.3,2.5) node (X) a_1;
             (0.7,4) node (X) a_2;
             (1.7,5.5) node (X) a_3;
            
             (-0.3,-2.5) node (X) a_1';
             (0.7,-4) node (X) a_2';
             (1.7,-5.5) node (X) a_3';
        
        
        Where we used unitarity to eliminate the top and bottom gates. Using the motif Eq. (<ref>) we can telescope telescope the circuit until
        
        
            ∑_s_1,s_2(A^(j)_s_1,s_2A^(j)^†_s_1,s_2)_a_1:3,a'_1:3=
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
            11.5U_121-1.5U_12^† (1.5,2) edge[out=90,in=90] (2.5,2);
             (2.5,2) edge[out=-90,in=90] (2.5,-2);
             (1.5,-2) edge[out=-90,in=-90] (2.5,-2);
            
            5-53.53.55-3.5320.51.5-2-0.51.520.50.5-2-0.50.5 (0.5,0.3) circle (.2);
             (1.5,0.3) circle (.2);
             (0.5,-0.3) circle (.2);
             (1.5,-0.3) circle (.2);
            
             (0,2.5) edge[out=0,in=90] (0.5,2);
            
             (1,4) edge[out=0,in=90] (3,3.5);
             (2,5.5) edge[out=0,in=90] (3.5,5);
            
             (0,-2.5) edge[out=0,in=-90] (0.5,-2);
            
             (1,-4) edge[out=0,in=-90] (3,-3.5);
             (2,-5.5) edge[out=0,in=-90] (3.5,-5);
            
             (-0.3,2.5) node (X) a_1;
             (0.7,4) node (X) a_2;
             (1.7,5.5) node (X) a_3;
            
             (-0.3,-2.5) node (X) a_1';
             (0.7,-4) node (X) a_2';
             (1.7,-5.5) node (X) a_3';
        
        Finally, we use the explicit form of Eq. (<ref>) at the self-dual point to evaluate
        
        
            [baseline=([yshift=-1ex]current bounding box.center),every node/.style=scale=1,scale=.55]
            
              11.5U_121-1.5U_12^† (1.5,2) edge[out=90,in=90] (2.5,2);
               (2.5,2) edge[out=-90,in=90] (2.5,-2);
               (1.5,-2) edge[out=-90,in=-90] (2.5,-2);
            
              20.751.5-2-0.751.520.750.5-2-0.750.5 (0.5,0.55) circle (.2);
               (1.5,0.55) circle (.2);
               (0.5,-0.55) circle (.2);
               (1.5,-0.55) circle (.2);
            
               (0,2.5) edge[out=0,in=90] (0.5,2);
               (0,-2.5) edge[out=0,in=-90] (0.5,-2);
            
               (-0.3,2.5) node (X) a_1;
               (-0.3,-2.5) node (X) a_1';
            
             (0.5,0) node (X) c;
             (1.5,0) node (X) d;
            
              =
              ∑_b(U_12)_ab,cd(U^*_12)_a'b,cd=δ_aa',
        
        which corresponds to an initial state equal to a product state in the Z_j basis. This verifies the condition Eq. (<ref>).
        
        An MPS in both left and right canonical form describes a bistochastic quantum channel: one that preserves the identity. As a result, the ancilla density matrix R^j=2^1-t_2^t-1 for all j. Evaluating the Rényi entropies yields our result Eq. (<ref>).
        
        
        
        
        
        It seems likely that a similar analysis can be performed for general Clifford circuits, which also have degenerate entanglement spectra <cit.>. However, establishing this in general requires one to carve out several special cases (such as circuits that generate no entanglement at all from a number of initial states <cit.>) and we will not pursue this here.
        
        
        
        § NUMERICAL RESULTS: RANDOM UNITARY CIRCUITS
        
        
        In this section we present numerical results for the evolution of various entanglement measures, and their spatial fluctuations, for random unitary circuits. The coarse features of the evolution of the entanglement spectrum were already discussed in Ref. <cit.>. In particular, the bandwidth of the entanglement spectrum broadens linearly in time; as noted in that work, this broadening is a natural consequence of the wide separation between the entanglement and light-cone speeds. For circuits of depth ≥ 10 this broadening implies that an appreciable fraction of the spectrum of the reduced density matrix is zero to within numerical precision. Our focus here is on the large eigenvalues of the reduced density matrix (which dominate Rényi entropies S_n, n ≥ 1). Because the vast majority of the eigenvalues are near zero, this “low-entanglement-energy” tail can be described accurately by low-rank approximations as in Sec. <ref>, allowing us to go to circuits of depth L = 14 with minor computational effort.
        
        
        
         §.§ Benchmarking the low-rank approximation
        
        
        For depths t ≤ 10 we can compare the entanglement spectra computed by low-rank approximation with the exact ones (Fig. <ref>). Although the rank of the reduced density matrix in this case is 256, we find that working with the top 20 states allows us to match the low-energy behavior of the entanglement spectrum.
        
        
        For the largest depths we have considered, exact time evolution is not feasible; however, there is an exact result for the average purity of a semi-infinte system <cit.>, viz. γ̅= (4/5)^t-1[The exponent t-1 here is one less than in <cit.> because our partition of the system lies between the two gates in the top layer of the circuit, which therefore leaves the entanglement unaffected.]. For the numerical results presented here we increase the rank of the approximation until the mean computed purity matches this exact result to within statistical error (which is about 1%). For the largest depth we have systematically considered (t = 14) we need to keep ≈ 120 states to match the mean purity. This is only about 1% of the spectrum; thus the low-rank approximation is much more efficient than direct propagation of the channel would be.
        
        
        
         §.§ Shape of the entanglement spectrum
        
        
        In this section we discuss the shape of the entanglement spectrum and its relation to the evolution of the Rényi entropies for RUCs. First, let us recall that the eigenstate thermalization hypothesis predicts that the reduced density matrix of a subsystem should take the form 𝒩exp(- β H), where H is the Hamiltonian of the subsystem and β is the inverse temperature. Therefore the entanglement Hamiltonian -logρ∝β H + log𝒩, i.e., it is a stretched and shifted version of the physical Hamiltonian. The entanglement spectrum therefore has the same shape as the physical spectrum: for a large subsystem L_A with a local Hamiltonian, it will be essentially Gaussian in the bulk, with a bandwidth that increases as √(L_A), although the extreme value statistics (corresponding to the shape of the spectrum near its ground state) are model-dependent. When β√(L_A) is large, the entanglement spectrum will have a large bandwidth, and therefore (because of the Jacobian) the reduced density matrix will have a density of eigenvalues ϱ(λ) distribution of the form ϱ(λ) ∼ 1/λ. Infinite temperature is a singular limit, as β = 0 so the entanglement spectrum is degenerate. In practice, a typical, randomly picked state deviates form infinite temperature by an amount ∼ 1/√(L_A), so the entanglement “bandwidth” is L_A-independent (up to possible logarithmic dependences that we are not concerned with here). For Floquet systems or RUCs with no conservation laws, these arguments suggest that the entanglement spectrum of a small subsystem at very late times is degenerate up to finite size effects. Once finite size effects are included we expect a Marchenko–Pastur distribution <cit.>.
        
        
        
        For the case of interest to us – large subsystems at short times – the numerical evidence <cit.> suggests that the spectrum of the reduced density matrix has the density ϱ(λ) ∼ 1/λ over many decades, at any time 1 ≪ t ≪ l_A; this is qualitatively unlike the (compact) Wishart distribution that obtains at very late times <cit.>. The bulk of the spectrum of the reduced density matrix consists of eigenvalues below machine precision whenever t  8 (Fig. <ref>). Here, we are concerned with the “low-energy” or “high-Schmidt-coefficient” edge, which governs the behavior of the Rényi entropies S^(n), n ≥ 1, which we can follow out to later times t ≈ 14 (Fig. <ref>). The evolution at early times is nontrivial, but appears to settle down into a well-defined limiting behavior for t  7: there is a threshold in the entanglement density of states, followed by a linear increase with entanglement energy that persists out to the energies we can reliably access. The coefficient of this linear growth is approximately 2^t/2 for the larger accessible t. This is exponentially slower than the growth of the total number of states, so the fraction of states in the tails thins out exponentially in time.
        
        
        
        The Rényi index dependence of S^(n) for large n follows from this behavior of the limit shape, if we further assume that the entanglement spectrum is self-averaging.
        
        A simple model for the spectral density ρ(ϵ) of the entanglement energies that is consistent with the large deviation form Eq. (<ref>) is
        
        
            ρ(ϵ) = exp[α tΘ(E-v_∞ t)],
        
        in which case
        
        
            v_n = α - nv_∞/1-n
        
        A large n result this implies v_n/v_∞≃ 1 - 1/n, which is consistent with our numerical observations [Fig. <ref>].
        
        
        
        Although the late-time entanglement spectrum is not numerically accessible, our results allow us to comment on a few possible qualitative scenarios of the entanglement spectrum. First, it is clear numerically that the probability density of states is exponentially small near the low-energy edge of the entanglement spectrum. This turns out to be necessary for the Rényi entropies to have distinct velocities. (If one considers, e.g., a box-shaped entanglement DOS, it is simple to show that all Rényi entropies with n > 0 must have the same velocity, regardless of the aspect ratio of the box. Similar results hold for Marchenko-Pastur and other possible compact shapes.) Second, one might suppose the entanglement spectrum has a Gaussian shape. Matching exact results for S_0, the normalization of ρ, and S_2 requires the Gaussian to have a linearly growing mean and variance. Numerically, the entanglement bandwidth grows linearly in time rather than as a square root, possibly because of level repulsion. Developing a theory of how the entanglement spectrum evolves is an important question for future work: at present we do not have even a phenomenological Brownian-motion model of this growth.
        
        
        
        
        
         §.§ Statistics and spatial correlations of entanglement
        
        
        
        Out to the latest times we have considered, the sample-averaged Rényi entropies have Gaussian distributions. (Thus, quantities such as the purity are log-normally distributed.) Whether these distributions become anisotropic at much larger system sizes is unclear; however, we have not seen any sign of incipient skewness out to the times we can simulate (Fig. <ref>).
        
        We now turn to fluctuations of the entanglement across spatial cuts. The prediction of Ref. <cit.>, based on a mapping to the KPZ equation, is that the entanglement fluctuations are spatially correlated, with a correlation length ξ(t) ∼ t^2/3, and that the width of the entanglement distribution scales as t^1/3 (i.e., the entanglement “roughens”). The method used here works with an infinite system at a fixed depth, and enables one to address these spatial correlations. We find that the spatial correlations of entanglement do get longer-ranged in time, as their power spectrum clearly narrows in k-space (Fig. <ref>). The Fourier transform has a characteristic width, from which we can extract a correlation length that clearly grows sub-linearly with t. However, the correlation length remains short out to the latest times we can access, so we do not have the dynamic range to extract meaningful exponents.
        
        
        
        The KPZ picture also predicts that the entanglement “roughens” with time, i.e., its standard deviation grows. This is consistent with what we see, although, again, the roughening is too weak to extract meaningful exponents.
        
        
        
        
        
        
        
        
        
        
        
         §.§ Trajectory Approach For Purity
        
        
        Finally, we demonstrate the trajectory approach described in Section <ref>. To evaluate the purity using Eq. <ref> we evolve a pair of trajectories s_1:L and t_1:L with the transition probabilities given in Eq. (<ref>). For a random unitary circuit with a purity that fluctuates with spatial position, evaluating the unaveraged purity at a point would involve averaging over many trajectories with the same set of gates. As a proof of principle we instead focus on the ensemble averaged purity, and average the fidelity |⟨ψ_t_1:l|ψ_s_1:l||⟩^2 for l=1,… L with a trajectory of L=1000. In this case we can compare with the known exact result γ̅= (4/5)^t-1 for random unitary circuits <cit.>. Since we are now evolving vectors in the ancilla space rather than density matrices we can simulate deeper circuits. Fig. <ref> shows the average purity for depths up to 18, comparing the trajectory method with the exact result, as well as with the density matrix approach for depths up to 12. Good agreement is found in all cases.
        
        
        
         §.§ Circuits with more structure
        
        
        The transfer-matrix method discussed here extends directly from random unitary circuits to any other type of circuit that can be decomposed into a “brickwork” arrangement of two-site gates. We discuss two examples here: circuits with a conservation law and translation-invariant circuits.
        
        
        
          §.§.§ Number-conserving circuits
        
        
        As an illustrative example we now turn to circuits with a single conservation law, which we choose to be the number of ↑ spins in the computational basis <cit.>. For a random circuit, all classes of product states are equivalent; however, circuits with a conservation law yield very different entanglement DOS depending on the initial state. Fig. <ref> shows the results for three classes of initial states: random product states, random bit-strings in the computational basis, and a uniform Néel state. While the Néel state behaves analogously to the random unitary circuit, at least at these depths, we see that the other two types of product states give rise to very different entanglement spectra, with substantially higher DOS at low energies. The difference can be attributed to rare states with anomalously large weight on configurations with long strings of aligned spins, which do not entangle under number-conserving dynamics <cit.>.
        
        
        
          §.§.§ Translation-invariant circuits
        
        
        Next, we consider translation-invariant circuits, in which all the gates at a given time-step are identical. Gates could be the same at different time-steps (giving a Floquet system) or random at every time-step (giving a system with perfectly spatially correlated noise). We focus here on the former case. For concreteness we focus on the integrable Trotterization of the XXZ model that was recently introduced <cit.>. This model is parameterized by two parameters (η, λ), and the dynamics consists of repeated application of the two-site gate
        
        
            U(η, λ) ≡( [                1                0                0                0;                0 sinη/sin (η + λ)  sinλ/sin(η + λ)                0;                0  sinλ/sin(η + λ) sinη/sin (η + λ)                0;                0                0                0                1 ]).
        
        When η is imaginary and λ is real, this is a Trotterized version of the Ising phase of the XXZ chain, with larger η corresponding to larger easy-axis anisotropy.
        
        One could consider the dynamics of entanglement for either random or homogeneous initial states. For random states, we expect (and find) spatial fluctuations of entanglement, which have a growing correlation length as in random unitary circuits (Fig. <ref>). At the accessible times, we are unable to extract any clear qualitative difference between the behavior of the correlation length in these integrable circuits and the random unitary case. For translation-invariant states (specifically the Néel state), the quantum channel converges to a definite steady state after a time interval on the order of the circuit depth (Fig. <ref>). Comparing the transient behavior between the XXZ circuit and a (presumably nonintegrable) circuit consisting of tiling a random two-site gate, we see that the transient behavior of the integrable case is different: the min-entropy overshoots its steady-state value in the integrable case, but not in the random case. Fig. <ref> shows the lowest 40 entanglement energies at t = 12 as a function of the anisotropy parameter η; as one would expect, increasing the anisotropy slows down the growth of entanglement. This manifests itself as the top eigenvalue in the Schmidt spectrum drifting toward zero, leading to a large gap in the Schmidt spectrum. However, the min-entropy, starting from the Néel state, grows linearly out to the circuit depths we can access, with no signs of curvature. This is consistent with the intuitive picture of ballistic entanglement growth in integrable systems <cit.>, since quasiparticles move ballistically although spin dynamics is diffusive <cit.>.
        
        An interesting quantitative difference between integrable dynamics and generic chaotic dynamics is that (for intermediate values of η) the entanglement spectrum stays “narrow” in the integrable case: Schmidt coefficients do not rapidly spread out over many decades the way they do under random unitary dynamics. This is consistent with the quasiparticle picture, which predicts that entanglement should spread with a characteristic quasiparticle velocity that does not depend on the Rényi index. How this picture extends to the case of strong anisotropy, in which the quasiparticles have a broad distribution of velocities, is an interesting open question.
        
        
        
        § CONCLUSIONS
        
        
        In concluding, let us summarise the technical achievements of the quantum channel approach:
        
        
          *  Unitary circuits are presented directly as matrix product states in canonical form.
          
          *  The spectrum of the reduced density matrix is related to that of the ancilla states exposed by the diagonal cut.
          
          *  The resulting quantum channel allows us to work in the infinite width limit at finite depth.
          
          *  Analytical results are obtained for the kicked Ising model at the self-dual point using a simple graphical calculus, simpler than the approach of Ref. <cit.>.
          
          *  The numerical evaluation of the channel may be improved by making a low rank approximation for the ancilla density matrix or by unraveling a quantum trajectory over the physical states.
        
        
        In this work, we benchmarked the quantum channel approach against exact results for the purity in random unitary circuits, and used it to compute the shape and fluctuations of the entanglement spectrum at low “entanglement energies.” With relatively little computational effort we were able to get converged results for depths t = 14 for the entanglement spectrum and t = 18 for the purity. We expect that there is room to optimize the algorithm and perform more resource-intensive computations, allowing us to access somewhat later times than we have in the present work. Our results support and extend earlier work <cit.> showing that the entanglement spectrum has a bandwidth that grows rapidly in time, with a sharp onset. Under random unitary dynamics (with or without a conservation law), most of the Schmidt coefficients that are generated are exponentially small in circuit depth, and can therefore be truncated, allowing for efficient computation of the entanglement spectrum. We were able to compute the spatial fluctuations of entanglement for large systems, and finally for translation-invariant systems we obtained converged results for the entanglement spectrum by evolving the associated quantum channel to convergence.
        
        A key distinction between the quantum channel approach and standard methods (such as time-evolving block decimation <cit.>) for evolving a matrix-product state is that the quantum channel approach constructs the entanglement spectrum without explicitly representing the physical wavefunction at time t. Our truncation and sampling schemes are also conceptually somewhat different from that in time-evolving block decimation <cit.>. Thus the quantum channel offers benefits if one wants to compute the entanglement spectrum (since it requires storing only the object of interest); however, it is not clear how one would apply our methods to compute observables in the physical (rather than ancilla) space.
        
        A number of avenues for future work present themselves. The thermodynamic shape of the entanglement spectrum of a semi-infinite circuit is still not well understood beyond its coarsest features. At a technical level, applying our approach to the entanglement of a finite region will involve considering the quantum channel at different `times', instead of the stationary (distribution of the) ancilla density matrix. Finally, our analysis of the kicked Ising model shows that the soluble self-dual point arises simply from the relation Eq. (<ref>), and suggests a criterion for searching for more models in the same class.
        
        
        
        § ACKNOWLEDGEMENTS
        
        This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. Support of EPSRC Grant No. EP/P034616/1 (AL) and NSF Grant No. DMR-1653271 (SG) is gratefully acknowledged. S.G. thanks P.-Y. Chang, X. Chen, V. Khemani, S. Parameswaran, J. Pixley, F. Pollmann, and T. Rakovszky for helpful discussions.
        
        
        41 natexlab#1#1 bibnamefont#1#1 bibfnamefont#1#1 citenamefont#1#1 url<#>1 urlprefixURL [Polkovnikov et al.(2011)Polkovnikov,
          Sengupta, Silva, and Vengalattore]polkovnikov_reviewauthorA. Polkovnikov,
          authorK. Sengupta,
          authorA. Silva, andauthorM. Vengalattore,
          journalRev. Mod. Phys.volume83,
          pages863 (year2011),
          <https://link.aps.org/doi/10.1103/RevModPhys.83.863>.
        
        [Deutsch(1991)]deutsch_ethauthorJ. M.Deutsch,
          journalPhys. Rev. Avolume43,
          pages2046 (year1991),
          <https://link.aps.org/doi/10.1103/PhysRevA.43.2046>.
        
        [Srednicki(1994)]srednicki_ethauthorM. Srednicki,
          journalPhys. Rev. Evolume50,
          pages888 (year1994),
          <https://link.aps.org/doi/10.1103/PhysRevE.50.888>.
        
        [Rigol et al.(2008)Rigol, Dunjko, and
          Olshanii]rigol2008thermalizationauthorM. Rigol,
          authorV. Dunjko, andauthorM. Olshanii,
          journalNaturevolume452,
          pages854 (year2008).
        
        [Deutsch(2018)]deutsch_reviewauthorJ. Deutsch,
          journalReports on Progress in Physics  (year2018).
        
        [Calabrese and Cardy(2005)]calabrese2005evolutionauthorP. CalabreseandauthorJ. Cardy,
          journalJournal of Statistical Mechanics: Theory and Experimentvolume2005, pagesP04010
          (year2005).
        
        [Kim and Huse(2013)]Kim:2013aaauthorH. KimandauthorD. A.Huse,
          journalPhys. Rev. Lett.volume111,
          pages127205 (year2013),
          <https://link.aps.org/doi/10.1103/PhysRevLett.111.127205>.
        
        [Nahum et al.(2017)Nahum, Ruhman, Vijay,
          and Haah]nrvhauthorA. Nahum,
          authorJ. Ruhman,
          authorS. Vijay, andauthorJ. Haah,
          journalPhys. Rev. Xvolume7,
          pages031016 (year2017),
          <https://link.aps.org/doi/10.1103/PhysRevX.7.031016>.
        
        [von Keyserlingk et al.(2018)von
          Keyserlingk, Rakovszky, Pollmann, and Sondhi]Keyserlingk2017authorC. W.von Keyserlingk,
          authorT. Rakovszky,
          authorF. Pollmann, andauthorS. L.Sondhi,
          journalPhys. Rev. Xvolume8,
          pages021013 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevX.8.021013>.
        
        [Nahum et al.(2018)Nahum, Vijay, and
          Haah]Nahum2017authorA. Nahum,
          authorS. Vijay, andauthorJ. Haah,
          journalPhys. Rev. Xvolume8,
          pages021014 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevX.8.021014>.
        
        [Žnidarič
          et al.(2008)Žnidarič, Prosen, and
          Prelovšek]Znidaric:2008aaauthorM. Žnidarič,
          authorT. Prosen, andauthorP. Prelovšek,
          journalPhysical Review Bvolume77,
          pages064426 (year2008).
        
        [Rakovszky et al.(2019)Rakovszky,
          Pollmann, and von Keyserlingk]rpv2019authorT. Rakovszky,
          authorF. Pollmann, andauthorC. von Keyserlingk,
          journalarXiv preprint arXiv:1901.10502  (year2019).
        
        [Huang(2019)]yichen2019authorY. Huang,
          journalarXiv preprint arXiv:1902.00977  (year2019).
        
        [Li and Haldane(2008)]ent_spectrumauthorH. LiandauthorF. D. M.Haldane,
          journalPhys. Rev. Lett.volume101,
          pages010504 (year2008),
          <https://link.aps.org/doi/10.1103/PhysRevLett.101.010504>.
        
        [Pollmann and Moore(2010)]pollmann2010authorF. PollmannandauthorJ. E.Moore,
          journalNew Journal of Physicsvolume12,
          pages025006 (year2010).
        
        [Ho and Abanin(2017)]ho2017authorW. W.HoandauthorD. A.Abanin,
          journalPhys. Rev. Bvolume95,
          pages094302 (year2017),
          <https://link.aps.org/doi/10.1103/PhysRevB.95.094302>.
        
        [Mezei and Stanford(2017)]mezei2017entanglementauthorM. MezeiandauthorD. Stanford,
          journalJournal of High Energy Physicsvolume2017, pages65 (year2017).
        
        [Zhou and Nahum(2018)]tianciauthorT. ZhouandauthorA. Nahum,
          journalarXiv preprint arXiv:1804.09737  (year2018).
        
        [Bertini
          et al.(2018a)Bertini, Kos, and
          Prosen]Bertini:2018aaauthorB. Bertini,
          authorP. Kos, andauthorT. Prosen,
          journalarXiv preprint arXiv:1805.00931
          (year2018a).
        
        [Bertini
          et al.(2018b)Bertini, Kos, and
          Prosen]Bertini:2018fbzauthorB. Bertini,
          authorP. Kos, andauthorT. Prosen
          (year2018b), 1812.05090.
        
        [Chang et al.(2018)Chang, Chen,
          Gopalakrishnan, and Pixley]ccgpauthorP.-Y.Chang,
          authorX. Chen,
          authorS. Gopalakrishnan,
          andauthorJ. Pixley,
          journalarXiv preprint arXiv:1811.00029  (year2018).
        
        [Carmichael(2009)]carmichael_bookauthorH. Carmichael,
          titleAn open systems approach to quantum optics: lectures
          presented at the Université Libre de Bruxelles, October 28 to November 4,
          1991, vol. volume18 (publisherSpringer Science &
          Business Media, year2009).
        
        [Vidal(2007)]vidal2007authorG. Vidal,
          journalPhys. Rev. Lett.volume98,
          pages070201 (year2007),
          <https://link.aps.org/doi/10.1103/PhysRevLett.98.070201>.
        
        [Perez-Garcia et al.(2006)Perez-Garcia,
          Verstraete, Wolf, and Cirac]Perez-Garcia:2006aaauthorD. Perez-Garcia,
          authorF. Verstraete,
          authorM. M.Wolf, andauthorJ. I.Cirac,
          journalarXiv preprint quant-ph/0608197  (year2006).
        
        [Schollwöck(2011)]Schollwock:2011aaauthorU. Schollwöck,
          journalAnnals of Physicsvolume326,
          pages96 (year2011).
        
        [Orús(2014)]Orus:2014aaauthorR. Orús,
          journalAnnals of Physicsvolume349,
          pages117 (year2014).
        
        [Lindblad(1976)]lindblad1976generatorsauthorG. Lindblad,
          journalCommunications in Mathematical Physicsvolume48, pages119 (year1976).
        
        [Gopalakrishnan and Zakirov(2018)]gz2018authorS. GopalakrishnanandauthorB. Zakirov,
          journalarXiv preprint arXiv:1802.07729  (year2018).
        
        [Gopalakrishnan(2018)]sg_og_2018authorS. Gopalakrishnan,
          journalPhys. Rev. Bvolume98,
          pages060302 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevB.98.060302>.
        
        [Note1()]Note1noteThe exponent t-1 here is one less than in <cit.> because our partition of the system lies between the two gates in
          the top layer of the circuit, which therefore leaves the entanglement
          unaffected.[Chen and Ludwig(2017)]Chen:2017aaauthorX. ChenandauthorA. W.Ludwig,
          journalarXiv preprint arXiv:1710.02686  (year2017).
        
        [Khemani et al.(2018)Khemani, Vishwanath,
          and Huse]Khemani2017authorV. Khemani,
          authorA. Vishwanath,
          andauthorD. A.Huse,
          journalPhys. Rev. Xvolume8,
          pages031057 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevX.8.031057>.
        
        [Rakovszky et al.(2018)Rakovszky,
          Pollmann, and von Keyserlingk]rpvauthorT. Rakovszky,
          authorF. Pollmann, andauthorC. W.von Keyserlingk,
          journalPhys. Rev. Xvolume8,
          pages031058 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevX.8.031058>.
        
        [Vanicat et al.(2018)Vanicat, Zadnik, and
          Prosen]Prosen_trotterization1authorM. Vanicat,
          authorL. Zadnik, andauthorT. c. v.Prosen,
          journalPhys. Rev. Lett.volume121,
          pages030606 (year2018),
          <https://link.aps.org/doi/10.1103/PhysRevLett.121.030606>.
        
        [Ljubotina
          et al.(2019a)Ljubotina, Zadnik, and
          Prosen]Prosen_trotterization2authorM. Ljubotina,
          authorL. Zadnik, andauthorT. Prosen,
          journalarXiv preprint arXiv:1901.05398
          (year2019a).
        
        [Ljubotina
          et al.(2019b)Ljubotina, Znidaric, and
          Prosen]Prosen_trotterization3authorM. Ljubotina,
          authorM. Znidaric, andauthorT. Prosen,
          journalarXiv preprint arXiv:1903.01329
          (year2019b).
        
        [Alba and Calabrese(2017)]alba_calabreseauthorV. AlbaandauthorP. Calabrese,
          journalProceedings of the National Academy of Sciencesvolume114, pages7947 (year2017).
        
        [Ljubotina et al.(2017)Ljubotina, Žnidarič, and Prosen]lzpauthorM. Ljubotina,
          authorM. Žnidarič,
          andauthorT. Prosen,
          journalNature Communicationsvolume8,
          pages16117 EP  (year2017),
          <http://dx.doi.org/10.1038/ncomms16117>.
        
        [De Nardis et al.(2018)De Nardis,
          Bernard, and Doyon]dbd2authorJ. De Nardis,
          authorD. Bernard, andauthorB. Doyon,
          journalarXiv preprint arXiv:1812.00767  (year2018).
        
        [Gopalakrishnan and Vasseur(2018)]gv_superdiffusionauthorS. GopalakrishnanandauthorR. Vasseur,
          journalarXiv preprint arXiv:1812.02701  (year2018).
        
        [Vidal(2003)]vidal2003authorG. Vidal,
          journalPhys. Rev. Lett.volume91,
          pages147902 (year2003),
          <https://link.aps.org/doi/10.1103/PhysRevLett.91.147902>.
