Theory of Ergodic Quantum Processes

The generic behavior of quantum systems has long been of theoretical and practical interest. Any quantum process is represented by a sequence of quantum channels. We consider general ergodic sequences of stochastic channels with arbitrary correlations and non-negligible decoherence. Ergodicity includes and vastly generalizes random independence. We obtain a theorem which shows that the composition of such a sequence of channels converges exponentially fast to a replacement (rank-one) channel. Using this theorem, we derive the limiting behavior of translation-invariant channels and stochastically independent random channels. We then use our formalism to describe the thermodynamic limit of ergodic matrix product states. We derive formulas for the expectation value of a local observable and prove that the two-point correlations of local observables decay exponentially. We then analytically compute the entanglement spectrum across any cut, by which the bipartite entanglement entropy (i.e., R\'enyi or von Neumann) across an arbitrary cut can be computed exactly. Other physical implications of our results are that most Floquet phases of matter are metastable and that noisy random circuits in the large depth limit will be trivial as far as their quantum entanglement is concerned. To obtain these results, we bridge quantum information theory to dynamical systems and random matrix theory.


I. OVERVIEW AND SUMMARY OF THE RESULTS
How do generic quantum systems behave? The physical change of a quantum system over a unit of time is modeled by the application of a quantum channel. Quantum channels are the most general formulation of physical (quantum) processes such as various steps of a quantum computation, effects of noise and errors on the state, and measurements [1]. Any particular evolution of a quantum system can be obtained by the application of a sequence of suitably chosen quantum channels.
The generic (average) behavior of quantum systems has long been of theoretical and practical interest. Random channels appear in a wide variety of applications, from quantum chaos [2] to holographic dualities in theories of quantum gravity [3] to operator dynamics [4]. Similarly, random local circuits are currently intensely studied for their k-design properties [5], potential to demonstrate quantum supremacy [6,7], and understanding operator spread and entanglement growth [8]. Recently, motivated by the demonstration of quantum supremacy, random circuits have become central candidates to show hardness in near-term quantum computing. For example, Google recently demonstrated a 53-qubit experimental demonstration of the hardness of sampling from the output of random circuits [9].
In nature, whenever one studies a quantum system, leakage of information to the environment is at play. A key challenge is to address the effects of decoherence in natural settings or account for it in the lab especially now that quantum error correction schemes have not been realized. To model decoherence on the output of a quantum circuit, scenarios have been considered in which after every gate a measurement is performed with some probability [10]. Here, we take a different perspective in which decoherence is directly included in the channels describing the process.
In this paper, we consider a sequence of channels as a trajectory of a dynamical system. Dynamical systems theory is a rich field with seminal results such as the Furstenberg-Kesten theorem [11] and the multiplicative ergodic theorem of Oseledec [12]. By bridging quantum information to the theory of dynamical systems, we address the behavior of quantum channels in a very general setting and in the presence of decoherence. We take the underlying dynamical system and the sequence of channels to be ergodic. Ergodicity allows for arbitrarily long-range correlations among the channels, including and vastly generalizing the extreme cases of independently and identically distributed (IID), as well as (time)-translational-invariant, channels (Fig. 1). Below, we give a more formal definition.
We first present a general theorem for an ergodic sequence of quantum channels (Theorem 1), which shows that the composition of such channels converges exponentially fast to a stochastic sequence of rank-one channels. A corollary of this result is the well-known convergence in the translation-invariant case to a fixed rank-one channel. We then apply our theorem in a natural setting in which Kraus operators [B i m m in Eq. (2)] in each channel are IID random Haar isometries, where by Haar isometry we mean a section cut out of a unitary matrix that is uniformly drawn from the space of all possible unitary matrices. We analyze the asymptotics with respect to the environment or the system or both tending to infinity and prove universal Gaussianity in these limits.
In the past, "ergodic" quantum channels were considered, which, to the best of our knowledge, assumed special subsets of possibilities considered herein. For example, a channel is chosen at random from some ensemble and then repeatedly applied [13]. In other work, time dynamics are analyzed for a quantum system with repeated independently chosen random interactions with an environment [14], and convergence results are obtained for the expected behavior of a sequence of IID channels and their applications [15]. Other instances studied include certain random independent channels and their compositions (e.g., from a finite set of random isometries) [16,17]. See Ref. [18] for a review. Our work considers a general ergodic sequence and, we believe, serves as a vast generalization of the prior work. We emphasize that this work allows for correlated randomness or even pseudorandomness [19] generated by quasiperiodic dynamics.
Interestingly, the formalism of quantum channels is suited for the calculation of physically relevant quantities such as expectation values of observables and correlation functions of low-dimensional quantum systems. Such quantum systems are best described by the density matrix renormalization group (DMRG), whose natural representations are matrix product states (MPS).
Theorem 1 allows us to analyze the thermodynamic limit of ergodic MPS. We derive formulas for the expectations of a general observable and for the bipartite entanglement entropy across an arbitrary cut. The second main result of our paper is Theorem 2 on the exponential decay of twopoint functions for ergodic MPS.
Any MPS is a ground state of a nontrivial Hamiltonian. Furthermore, MPS are the natural representation of states obtained from DMRG.
The tensors in the MPS representation of the ground state of any interacting Hamiltonian (any nontrivial Hamiltonian) will be correlated. This result is so even if the underlying Hamiltonian has random independent local terms. In particular, a state obtained from a DMRG procedure, say, the ground state of a 1D spin glass, will be correlated yet fall within the class of ergodic MPS.
This work, to the best of our knowledge, is the first to shed light on the structure of correlations for a large general class of correlated MPS that one obtains using procedures such as DMRG. This class includes MPS with IID local tensors. However, as the above remarks show, IID local tensors are mostly of theoretical interest, as correlations inevitably arise in any realistic interacting system even if the underlying Hamiltonian is fully disordered.
More broadly, it is important to explore the manifold of non-translation-invariant MPS. The present work provides a well-motivated path for further study of this line of research.
Despite various efforts over about the past two decades, previous results in the same vein as this paper were obtained only for very special classes of aforementioned quantum processes, with limited applicability to physically realistic settings. We achieve our results by bridging the mathematical fields of dynamical systems and random matrix theory to quantum information science. Roughly speaking, Theorem 1 can be seen as the "quantization" of classical results on the convergence of Markov chains. The proof is based on generalizations of the Perron-Frobenius theorem and of the multiplicative ergodic theorem. The analysis of Haar channels is enabled by the use of the Weingarten calculus, with Gaussianity of the fluctuations established by verifying Wick's theorem. Full proofs and derivations are given in the appendixes and in Ref. [20].
In Sec. II A, we detail further consequences of our work for physics of matter, quantum information theory, and MPS.

II. ERGODIC THEOREM FOR COMPOSITION OF QUANTUM CHANNELS
A quantum channel is a completely positive and tracepreserving linear map ϕ on the space of D × D matrices: FIG. 1. Venn diagram for distributions of possible ergodic sequences of channels. Among the many subsets, three subsets are shown and substantially enlarged for readability.
The total change to a quantum state is obtained by the composition of suitably chosen such maps on the initial state (i.e., density matrix ρ): We consider sequences ðϕ 0 ; ϕ 1 ; …Þ of channels drawn from an ensemble which is ergodic with respect to the time shift ðϕ 0 ; ϕ 1 ; …Þ ↦ T ðϕ 1 ; ϕ 2 ; …Þ. In statistical physics, originating in the fundamental work of Boltzmann, the ergodic hypothesis is traditionally understood as the equivalence of time averages and ensemble averages [21]. Later, the ergodic hypothesis was made more precise in mathematical physics culminating in the works of von Neumann and Birkhoff [22]. In those works, it is recognized that the ergodic hypothesis of Boltzmann is a consequence of the ergodic nature of the underlying dynamics. Formally, a map T on an ensemble is ergodic if it preserves probabilities and, starting from a typical point, the dynamics T generates covers the phase space with probability one (we recall the precise mathematical definition in Appendix A). In general, there may be selection rules separating distinct sectors of the Hilbert space in which the quantum state ρ resides that do not mix under the evolution (i.e., application of channels). Furthermore, there may be a subspace that is transient for the evolution, such as the space of excited states of a system for which the long-time evolution is described by decay to a ground state. The long-time behavior of the system can be understood by restricting attention to a single selection sector orthogonal to the transient subspace. Without loss of generality, we can assume that our system (i) has no selection rules-there are no (proper) subspaces invariant under the evolution, because if there were such a subspace, we could simply study each one separately-and (ii) has no transient subspace, because if it does, then the systems eventually leaves it with a probability arbitrary close to one.
Assumptions.-Mathematically, the only assumption we need is that with probability one there exists an integer N 0 > 0 such that for all N > N 0 , the map ϕ N ∘ ϕ N−1 ∘ Á Á Á ∘ ϕ 0 is strictly positive. This is implied by the following more easily verifiable assumptions.
(1) For some n 0 , Prob½ϕ n 0 ∘ Á Á Á ∘ ϕ 1 ∘ ϕ 0 is strictly positive > 0: (2) With probability one, if, for some observable M ≥ 0, we have tr½Mϕ 0 ðρÞ ¼ 0 for all ρ, then M ≡ 0. Assumption 1 is a generalization to ergodic quantum processes of the notion of irreducibility for a finite Markov chains, for which the transition matrices can be analyzed by the Perron-Frobenius theorem.
In the absence of selection rules (or upon reduction to an appropriate selection sector), these assumptions are expected to hold for any system interacting with an environment at positive temperature (but would not hold at zero temperature, where excited states are transient). More generally, they should hold provided there is nonnegligible decoherence, leading to the possibility of a transition between any two states of the system. Note that the second assumption can alternatively be expressed as the following condition on the dual channel: For a nonnegative observable M, if ϕ Ã 0 ðMÞ ¼ 0, then M ¼ 0. For each m; n ∈ Z with m < n, let the quantum channel Ψ n;m be Ψ n;m ≡ ϕ n ∘ Á Á Á ∘ ϕ m . Our main result is the following theorem.
Theorem 1.-There exists 0 < μ < 1 and an ergodic sequence of D × D density matrices ð…; Z m ; Z mþ1 ; …; Z n ; Z nþ1 ; …Þ such that, given x ∈ Z, the following holds: for all m ≤ x, n ≥ x, and any density matrix ρ ∈ C D×D , where C x < ∞ almost surely. Furthermore, the matrices Z j satisfy the shift equations The matrices Z j are analogous to the microstates of a quantum system at equilibrium, and the recursion given by Eq. (4) is the transition from one microstate to the next dictated by the ergodic quantum process. The theorem above proves that any initial state ρ, which may be far from equilibrium, approaches equilibrium exponentially fast in the following sense. The dynamical trajectory of any initial state, after long enough evolution with respect to the application of the ergodic sequence of the quantum channels, becomes exponentially close to the trajectory of the equilibrium microstates.
The parameter μ is a characteristic of a given ergodic quantum process; it depends on the ensemble describing of the process but is nonrandom in that it does not depend on the specific realization of a sequence of channels from the ensemble. The logarithm λ ¼ − log μ is a Lyapunov exponent giving the inverse relaxation time of the process. The prefactors C x appearing in Eq. (3) are, on the other hand, random (dependent on the specific realization) and as indicated depend on the origin of time denoted by x.
A key point is that for n − m sufficiently large Ψ n;m is exponentially close to the replacement (rank-one) channel Ψ n;−∞ ðρÞ ¼ Z n tr½ρ ¼ Z n . This result is largely a consequence of non-negligible decoherence (assumption 1 above).
An important feature of this convergence is that, although the distribution of Z n is fixed in time, the matrices Z n are random objects that fluctuate with respect to the time step n.
By ergodicity, the distribution of Z n can be determined from the frequency of occurrences in the trajectory Z n ; Z nþ1 ; …. A trivial example is furnished by the IID rank-one channels ϕ j : ρ 1 tr½ρ with probability 1=2; with two distinct density matrices ρ 0 and ρ 1 . In this case, each Z n has the distribution and the individual trajectories fluctuate randomly between the two values ρ 0 and ρ 1 , each occurring with frequency 1=2.
To fully characterize the statistical equilibrium described by the process …; Z j ; …, we need to solve the recursion in Eq. (4). This solution depends on the specification of the channels. To demonstrate the theory, we consider two natural extreme cases of an ergodic sequence of channels. First, we take the translationally invariant channels in which every channel in the sequence is equal. We then consider the other extreme in which the B i k k in Eq. (2) are blocks of IID Haar unitaries.

A. Physical consequences
We have quantifiable bounds and theorems that, to the best of our knowledge, for the first time apply to correlated quantum processes and to the ground states of interacting quantum matter (see Secs. I and IV). Any realistic physical process inevitably has temporal correlations, which is so even if the underlying process is Markovian, for which two subsequent times are correlated. We emphasize that correlated channels naturally arise in the study of nontrivial systems, such as interacting quantum many-body systems (see Secs. I and IV). Therefore, an IID assumption is of purely theoretical interest. This work, however, applies to a vast class of physically realistic quantum processes.
Three physical corollaries to our theorem are as follows. (a) Nonequilibrium phases of matter that are engineered by time-periodic driven Hamiltonians, such as in Floquet systems [23][24][25], must be metastable when interactions with an environment in positive temperature are non-negligible (equivalent of the assumption above). Therefore, generically any experimentally viable Floquet phase can only be metastable. (b) These channels become trivial as far as their information content is concerned. In fact, since the channels are asymptotically replacement channels, the process cannot even convey classical information. This constraint can be seen intuitively by the fact that a unique fixed point is reached irrespective of the input initial quantum state. We emphasize that channels may be highly correlated; this possibility holds even in the time-translation-invariant case. For example, in nearterm quantum computing with noisy random circuits [9,26,27], in the limit of large depth, all initial memory of the state is lost and a unique final state results. (c) A constant gap implies an exponential decay of correlations [28,29] and in one dimension an area law for entanglement entropy [30]. Later, Brandão and Horodecki [31] proved that in one dimension exponential decay of correlations implies an area law. Below in Sec. IV, we prove a partial converse, which says that finitely correlated states with an ergodic MPS representation have an exponential decay of correlation.

III. APPLICATIONS
In the translationally invariant case, the channels are the same, ϕ j ¼ ϕ 0 for all j, and assumption 1 reduces to the requirement that ϕ n 0 0 be strictly positive for some n 0 . The resulting sequence Z j given by Theorem 2 is also translation invariant, Z j ¼ Z 0 . Indeed, in this case, one may easily show that Z 0 is the unique eigenmatrix of ϕ 0 with eigenvalue 1, and all other eigenvalues of ϕ 0 have modulus smaller than 1 − μ. Equation (3) shows the exponentially fast convergence of the density matrix to the state Z 0 under repeated application of the channel ϕ 0 . That is, starting from any density matrix, repeated application of the channel derives the system exponentially fast toward the equilibrium state Z 0 .
In general, Z j is not an eigenmatrix of ϕ j ; rather, these matrices obey Eq. (4). Since Z j−1 and Z j are drawn from the same probability distribution, Eq. (4) implies a fixed point equation for the distribution of Z j . For general ergodic channels, the precise form of the fixed point equation depends on the correlations between the channels. For IID channels, the fixed point equation takes a particularly simple form whose solution can be determined by the method of moments (see Appendix B).
To illustrate this implication in a natural example, consider the other extreme and take the sequence of channels where Q r is a pure state (rankone projection) in C r×r , which we think of as the environment, and ρ ∈ C D×D is the initial state of the system. Operators U j are a sequence of independent Haar distributed Dr × Dr dimensional unitaries. Random Haar channels are the archetypes of generic quantum channels. Properties of a single random channel of this type are considered in Ref. [32]. Here, we consider an ergodic process obtained from an IID sequence of such channels.
We show that the average of the equilibrium matrices Z j is the totally mixed state (see Appendix B) Furthermore, using the Weingarten calculus [16,33], we show that fluctuations around the average are given by where tr½W j ¼ 0 and. in the limit that D → ∞ or r → ∞ (or both), W j has a Gaussian distribution on the space of D × D matrices: where Ξ D is a D-dependent normalization (see Appendix B). As a result, for large D the eigenvalues of Z j are distributed according to Wigner's semicircle law: . This distribution allows analytical computation of spectral properties of Z j . For example, the von Neumann entropy of Z j is SðZ j Þ ≈ −tr½Z j log Z j , which is given by This expression for large D is approximately and for large D and r the limiting expression is The entropy deviates by at most an order one quantity from the maximal value log D that is attained for a maximally mixed state.

IV. MATRIX PRODUCT STATES
MPS and their generalizations [34][35][36][37] provide efficient representations of quantum states by which classical simulation of quantum many-body systems becomes viable and are the natural representation of the density matrix renormalization group [38] and its tensor network generalizations. Applications range from efficient calculation of the ground state properties of quantum matter [39] to the outputs of quantum circuits [36]. They provide the tools for proving the existence of satisfying assignments in quantum satisfiability [40] and the area law [30]. Random MPS are the candidate boundary states for recent theoretical proposals of the theory of quantum gravity [3].
For the sake of concreteness, let us introduce the (one-dimensional) MPS on 2N þ 1 qudits, denoting where i k 's are the physical indices, d is the local (physical) dimension of the Hilbert space, A i k k 's are D × D × d tensors, and D is the bond dimension [34].
So far, rigorous results on generic MPS and their generalizations mainly focus on the translational-invariant case, where all the tensors in Eq. (9) are equal [31,41]. Other works focus on statistics of random MPS [42]. Here, we consider the general case in which translational invariance is relaxed. The most meaningful extensions pertaining to states of disordered systems and outputs of random quantum circuits require that the tensors are drawn only from a distribution. For example, if a Hamiltonian has local terms that are ergodic, then one expects the MPS representation of the states also to be ergodic. A similar statement is expected for the output state of a quantum circuit if the action of the circuit on the qubits is shift invariant.
Expectations in the MPS state jψðNÞi can be related to a generalized quantum process ϕ m ðMÞ ¼ However, they may not be quantum channels; although they are completely positive, they need not be trace preserving.
A preliminary observation is that the normalization can be expressed as follows: where jαi, α ¼ 1; …; D, are the elements of the computational basis on C D . The operators jαihβj, α; β ¼ 1; …; D, form a basis for the space of D × D matrices, allowing us to write the above identity as where Tr denotes the trace of a superoperator, defined as Tr½Φ ¼ P α;β tr½jβihαjΦðjαihβjÞ. There is a similar expression for the expectation of a local observable in the MPS state jψðNÞi. A local observable O on the spins in ½m; n is a linear operator on ⊗ n j¼m C d . Following the computation leading to Eq. (10), one may easily verify that the expectation of such an observable in the state jψðNÞi is given by the expression THEORY OF ERGODIC QUANTUM PROCESSES PHYS. REV. X 11, 041001 (2021) whereÔ is the following superoperator: In Theorem A.4 in Appendix A, we state and prove a generalization of Theorem 2 to non-trace-preserving quantum operations ϕ m provided they satisfy ϕ m ðMÞ ¼ 0 for M ≥ 0 if and only if M ¼ 0 in addition to the above assumptions. The main new feature of the generalization is an ergodic analog of the left eigenvector of a channel, namely, a sequence Z 0 j satisfying In the trace-invariant case, Z 0 j is the maximally mixed state for all j. Also, the shift equation (4) is replaced by Using this result, we can directly calculate the thermodynamic limit of the expectation value of an observable with respect to an MPS [Eq. (9)]: Furthermore, the entanglement spectrum of a bipartite partition of the infinite chain across the bond j ∼ j þ 1, i.e., the spectrum of the reduced density matrix for one-half of the chain, corresponds to the eigenvalues of . For instance, the bipartite von Neumann entanglement entropy is Since Z j and Z 0 jþ1 are positive, R j is similar to a positive matrix and log R j is well defined.
Turning our attention to correlation functions, let O 1 and O 2 be local observables supported at (or near) the site x ¼ 0. Let O 1 ðxÞ and O 2 ðxÞ denote the corresponding observables translated to have support at a general site x. We have the following.
In words, the two-point correlation function between local observables decays exponentially with the distance between their supports.
To the best of our knowledge, previous rigorous results considered only translationally invariant MPS. Equations (13) and (14) and Theorem 1 generalize these to the much larger class of MPS generated by ergodic sequences of matrices. In the translationally invariant case ϕ 0 ¼ ϕ m , the sequences Z n and Z 0 m are constant; i.e., Z nþ1 ¼ Z n and Z 0 m ¼ Z 0 mþ1 for all n, m. Equation (13) reduces to the entropy SðjÞ is independent of j, and the correlation bound in Theorem 1 is independent of x.

V. CONCLUDING REMARKS
In this paper, we present rigorous results on ergodic sequences of quantum channels with non-negligible decoherence. We show that, for an arbitrary initial density matrix, the state of the system subjected to such a sequence of channels converges exponential fast to an ergodic sequence of strictly positive density matrices. The ergodic framework of our results allows for noise models that are arbitrarily correlated.
In the second part, we prove results on the expectation values and two-point correlation functions of ergodic matrix product state. In particular, we rigorously show that the latter decays exponentially with the distance between the local observables.
Following the implications outlined in Sec. II A, we envision that these results could further pave the way for quantum computation by dissipation, shed light on the stability of nonequilibrium quantum phases of matter, and better characterize the output of random circuits in connection to near-term noisy quantum computers.
From a mathematical perspective, our result about MPS can be understood as describing correlations in the class of purely generated, ergodic, finitely correlated states, which generalizes the translation-invariant states considered in Ref. [43]. It would be of considerable interest to develop this theory further by extending it to ergodic, finitely correlated states over general C Ã algebras.
For these and other applications, it would be of technical interest to obtain asymptotics for the convergence factor μ and the prefactor C in Theorems 2 and 1 with respect to the system size of channels and bond dimension of MPS, respectively. In the context of MPS, quantitative control on the convergence factor could make it possible to consider states with varying bond dimension as the system size increases.

APPENDIXES
We first lay out the basic notation and aspects of ergodic theory needed for our work. We then state the theorems. Although some of the proofs are to be found in Ref. [20], we outline the proof ideas. We then detail the calculations for independently and identically distributed random Haar channels and take the various asymptotic limits with respect to the system size and/or the environment. These prove Eqs. (5)-(10) of the main text. Lastly, the derivation of bipartite entanglement entropy [Eq. (17) in the main text] is given.

APPENDIX A: FORMAL STATEMENT OF THE ERGODIC THEOREM
In this section, we state a result (Theorem A.3) that implies Theorem 1 in the main text. A sketch of the proof of Theorem A.3 is given. The full proof of this result is in Ref. [20].
To proceed, we need to introduce some basic notation and ideas. Let M D denote the set of D × D complex matrices. We primarily use the trace norm on M D : where ffi ffi · p denotes the operator square root. We also introduce the Hilbert-Schmidt inner product hMMi ¼ tr½M † M and norm Recall that a cone in a vector space is a set that is closed under addition and multiplication by positive scalars. LetĒ denote the closed cone: The interior ofĒ is which is an open cone.
Let LðM D Þ denote the set of linear maps from M D to M D . Given ϕ ∈ LðM D Þ, the adjoint ϕ Ã is defined by tr½M † ϕðMÞ ¼ tr½½ϕ Ã ðMÞ † M: A map ϕ ∈ LðM D Þ is positive provided ϕðMÞ ∈Ē whenever M ∈Ē. If, in addition, ϕðMÞ ∈ E whenever M ∈Ē, then we say that ϕ is strictly positive. Let ðΩ; F ; ProbÞ be a probability space with T∶Ω → Ω an invertible and ergodic map. Recall that T is ergodic provided it is probability preserving and Prob½A ¼ 0 or 1 for any event A with T −1 ðAÞ ¼ A. Starting from a given completely positive map valued random variable, ϕ 0 ∶Ω → fcompletely positive mapsg, we consider the sequence of maps defined by evaluating ϕ 0 along the trajectories of T: We follow the usual convention in probability theory and suppress the independent variable ω ∈ Ω in most formulas; when it is needed-as in Eq. (A1)-we use a subscript to denote the value of a random variable at a particular ω ∈ Ω.
The second assumption, which for quantum channels is equivalent to assumption 2 of the main text, is the following.
Assumption A.2.-With probability one, ker Remark.-If ϕ 0 is a quantum channel, then tr½ϕ 0 ðMÞ ¼ tr½M for any M, so ker ϕ 0 ∩Ē ¼ 0. Thus, for quantum channels, assumption 2 is equivalent to the equality THEORY OF ERGODIC QUANTUM PROCESSES PHYS. REV. X 11, 041001 (2021) 041001-7 ker ϕ Ã 0 ∩Ē ¼ f0g (assumption 2 of the main text). On the other hand, this identity is a nontrivial assumption for quantum channels. For example, if ϕðMÞ ¼ PMP þ SMS † with P a projection onto a subspace of dimension D=2 and S a partial isometry from I − P to P, then ϕ is a channel but ϕ Ã ðI − PÞ ¼ 0.
Given n ∈ Z, let By assumption A.1, the maps Φ n preserve the cone E. In 1948, Krein and Rutman [44] derived a generalization of the classical Perron-Frobenius theorem to the context of linear maps preserving a convex cone. Evans and Høegh-Krohn [45] apply this result to obtain results for positive maps on LðM D Þ. It follows from Theorem 2.3 in Ref. [45] that if Φ n is strictly positive, then there is a unique (up to scaling) strictly positive matrix R n such that Φ n ðR n Þ ¼ λ n R n , where λ n is the spectral radius of Φ n . Recall that the spectral radius is the maximum magnitude of all eigenvalues. Similarly, there is a unique (up to scaling) strictly positive matrix L n such that Φ Ã n ðL n Þ ¼ λ n L n . We normalize R n and L n so that tr½R n ¼ tr½L n ¼ 1.
Our main technical result is that L n converges as n → ∞, while R n converges as n → −∞. This result generalizes a theorem of Hennion on the Perron-Frobenius eigenvectors of products of positive matrices [46].
Theorem A.3 [20].-The limits lim n→−∞ R n ¼ Z 0 and lim n→∞ L n ¼ Z 0 0 exist almost surely and define random matrices Z 0 and Z 0 0 that fall in E with probability one. Furthermore, if we set Z n ¼ Z 0;T n ω and Z 0 n ¼ Z 0 0;T n ω , then where · denotes the projective action of a positive map on the strictly positive D × D matrices of trace 1: ϕ n · M ≡ 1 tr½ϕ n ðMÞ ϕ n ðMÞ: Remark.-If the maps ϕ n are quantum channels, then L n ¼ ð1=DÞI D , so Z 0 n ¼ ð1=DÞI D for all n. Given m < n, let P n;m denote the rank-one operator P n;m ðMÞ ¼ tr½Z 0 m MZ n : For n − m large, the operator ϕ n ∘ Á Á Á ∘ ϕ m is well approximated by P n;m . To formulate this result precisely, we introduce the following norm for a linear map: Theorem 2 of the main text is equivalent to the following.

Sketch of the proofs
The key to proving Theorem A.3 is to show the existence of the limits lim n→−∞ R n ¼ Z 0 and lim n→∞ L n ¼ Z 0 0 . More generally, we show that exist for any sequence Y n of density matrices, and the limits are independent of the choice of the particular sequence. Indeed, once these limits are known, for any n one takes Z n ¼ Z 0;T n ω and Z 0 n ¼ Z 0 0;T n ω to obtain the full ergodic sequence. Since Φ n · R n ¼ R n , it follows that Similarly, Φ Ã n · L n ¼ L n , and we find that The key to proving the existence of the limits in Eq. (A2) is to prove that Φ n is exponentially contractive with respect to a suitable metric dðY 1 ; Y 2 Þ on the space of density matrices: The existence of the limit in Eq. (A3) is guaranteed by Kingman's subadditve ergodic theorem [47]. A pivotal estimate in this context is the following: dðϕ · X; ϕ · YÞ ≤ cðϕÞdðX; YÞ; ðA4Þ where cðϕÞ ≤ 1 is a quantity that depends on the particular completely positive map ϕ. The subadditivity needed for the application of Kingman's theorem is provided by the inequality to prove that the right-hand side of Eq. (A3) is less than zero, it is convenient to work with the specialized metric where mðY 1 ; Y 2 Þ ¼ maxfλ ∈ RjλY 2 ≤ Y 1 g. Although this metric is more convenient for the arguments, we show that it is equivalent to the trace-norm metric for strictly positive density matrices. Since the Φ n are ultimately strictly positive by assumption 1 of the main text, the estimates obtained, in fact, hold in the trace norm. For details, see Ref. [20].

APPENDIX B: INDEPENDENTLY AND IDENTICALLY DISTRIBUTED RANDOM HAAR CHANNELS
The results so far are purely existential; for a particular ergodic sequence of channels, one wishes to know the statistical properties of Z 0 . To this end, we explicitly calculate them for two canonical examples. The first is for channels that are translation invariant which was detailed in the main text. The second class corresponds to the Kraus operators being IID Haar isometries (sections of independent Haar random unitary). The calculation of the latter case is the subject of this section.
Let U j , j ∈ Z, be an IID sequence of Haar distributed L × L unitary matrices, with L ¼ rD, where one may think of r and D as the size of the environment and system, respectively. Consider the channels ϕ j ∈ LðM D Þ given by ϕ j ðMÞ ¼ tr r ½U j ðM ⊗ Q r ÞU † j , where Q r is the r × r rankone projector Q r ≡ je 1 ihe 1 j with je 1 i ¼ ð1; 0; …; 0Þ T and tr r denotes the partial trace from M L → M D . Note that the partial trace is the adjoint (with respect to the Hilbert-Schmidt inner product) of the map M ↦ M ⊗ I r . The channel ϕ j can be expressed in the Kraus form as where the matrices U j;k are D × D blocks of U j : It is easy to see that this particular sequence of channels satisfies assumptions 1 and 2. In fact, each ϕ j is itself strictly positive with probability one. Thus, Theorem 1 applies, and the resulting stationary sequence Z j , j ∈ Z, satisfies the shift equation [Eq. (4) of the main text] We study the distribution of Z j by the method of moments. As such, we are interested in computing the expectation EfZ ⊗n j g, where Z ⊗n j ¼ Z j ⊗ Á Á Á ⊗ Z j with n factors of Z j . From EfZ ⊗n j g, we can compute the expectation of any multilinear form Q n l¼1 tr½M l Z j , via the identity: We begin computing the expectation for n ¼ 1, which is the first moment EfZ j g. The key identity here is In terms of matrices, this identity can be expressed as Z for an arbitrary M ∈ C L×L . By Eq. (B1), the conditional expectation of Z j given U k for k ≠ j satisfies since tr½Z j−1 ⊗ Q r ¼ tr½Z j−1 ¼ 1 and tr r ½I L ¼ rI D . Averaging over all U k for k ≠ j, we see that which is Eq. (7) in the main text. This result essentially proves that the expectation of ϕ j is the completely depolarizing channel.
We turn now to the computation of EfZ ⊗2 j g. The key is a generalization of Eq. (B2): Z UðLÞ Uðα 1 ; β 1 ÞUðα 2 ; β 2 ÞU † ðβ 0 where Wg½ð1; 1Þ; L ¼ ð1=L 2 − 1Þ and Wg½ð2Þ; L ¼ −½1=LðL 2 − 1Þ are Weingarten functions (see Ref. [33]). In terms of tensors, Eq. (B3) can be expressed as follows: where we introduce the notation extended by linearity to all of M L ⊗ M L . Now For the factors involving E ð2Þ L , we use the identity Thus, D . It follows that Averaging over U k , we obtain where m ð2Þ ≡ Eftr½Z 2 j−1 jg ¼ Eftr½Z 2 j g; since Z j−1 and Z j are equidistributed.
To complete the calculation of EfZ ⊗2 j g, we need to compute m ð2Þ . Note that m ð2Þ ¼ Eftr½E Solving for m ð2Þ , we find that Plugging this result into Eq. (B5), we obtain The variance of Z j can now be expressed as D . To proceed, we introduce the notation W j for the rescaled deviation of Z j from the mean: Thus, RAMIS MOVASSAGH and JEFFREY SCHENKER PHYS. REV. X 11, 041001 (2021) 041001-10 Our goal is to show that W j is asymptotically Gaussian, with covariance ð1=DÞE ð2Þ D − ð1=DÞE ð1;1Þ D . In other words, the distribution of W j is given by a traceless version of the Gaussian unitary ensemble [see Eq. (7)]. To obtain this result, we verify that W j asymptotically satisfies Wick's theorem.
Turning now to the computation of the nth moment for general n, we define R ðnÞ ¼ EfW ⊗n j g. By the shift invariance property, the right-hand side is independent of j, and using Eq. (B1) we see that Because There is a generalization of Eq. (B3) to higher-order products (see Ref. [33]). In terms of tensors, this generalization can be written as where Wg denotes the Weingarten function (see Ref. [33]), S n is the permutation group on f1; …; ng, and for v k ∈ C L the operator P α L has the action Applying Eq. (B9) to the average over U j in Eq. (B8) results in the following expression for R ðnÞ : The Weingarten function Wgðα; LÞ depends only on the conjugacy class of the permutation α: Wgðτατ −1 ; LÞ ¼ Wgðα; LÞ (see Ref. [33]). Furthermore, tr½P τ L M ⊗n depends only on the conjugacy class of τ, as we now show. The conjugacy classes of S n are labeled by m-tuples ⃗ c ¼ ðc 1 ; c 2 …; c m Þ with c 1 ≥ c 2 ≥ Á Á Á ≥ c m ≥ 1 and P j c j ¼ n. A permutation in the class ðc 1 ; c 2 …; c m Þ is given by (1) a partition of f1; …; ng into m sets of sizes c 1 ; …; c m and (2) a choice of a cyclic permutation on each subset of the partition. A permutation is cyclic if it has no proper invariant subsets-there are ðj − 1Þ! cyclic permutations in S j . For any permutation τ ∈ S n in the class ðc 1 ; c 2 …; c m Þ and M ∈ M L , the trace tr½P τ L M ⊗n is given by Since tr½Z j−1 ⊗ Q r ¼ tr½Z j−1 ¼ 1, we see from Eq. (B11) that tr½P τ L ðZ j−1 ⊗Q r −ð1=LÞI L Þ ⊗n ¼0 for any permutation τ with a fixed point, in which case the conjugacy class has c m ¼ 1. In other words, the sum over τ in Eq. (B10) can be restricted to τ ∈ S 0 n , where S 0 n ¼ fτ ∈ S n jτðjÞ ≠ j; j ¼ 1; …; ng: Let CS n denote the set of conjugacy classes of S n . Given ⃗ c ∈ CS n , let With the labeling of classes described above, this choice is consistent with the definitions of E ð1;1Þ L and E ð2Þ L given above. Associated to each class ⃗ c ¼ ðc 1 ; c 2 ; …; c m Þ ∈ CS n , we define a distinguished permutation ½⃗ c which consists of the cycle of the first c 1 numbers, followed by the cycle of the next c 2 numbers, etc., up to the last group of c m numbers: That is, ½ ⃗c ¼ ð1; 2; …;c 1 Þðc 1 þ 1; …;c 1 þ c 2 ÞÁÁÁðc 1 þ ÁÁÁþ c m−1 ; …;nÞ.
Putting all of this together, we obtain the following formula: where mð⃗ cÞ denotes the number of cycles in the conjugacy class ⃗ c, i.e., mð⃗ cÞ ¼ m if ⃗ c ¼ ðc 1 ; …; c m Þ with c m ≥ 1, wð ⃗c;r;DÞ ¼ L mð ⃗cÞ ðrD 2 þ 1Þ n=2 and we observe that tr ⊗n r ½E ⃗ c L ¼ r mð⃗ cÞ E ⃗ c D . Adding and subtracting ð1=DÞI D ⊗ Q r ¼ EfZ j−1 g ⊗ Q r inside the parentheses, we obtain the following alternative expression for wð ⃗ c; r; sÞ: wð⃗ c; r; DÞ ¼ L mð⃗ cÞ Let τ ∈ S 0 n and let which is the expectation appearing in the sum on the righthand side in Eq. (B14). The conjugacy class of τ, ½τ ¼ ðc 0 1 c 0 2 …c 0 m Þ, satisfies c 0 m ≥ 2. Given commuting matrices A and B, Applying this expression to EðτÞ, we obtain the following expression for EðτÞ: where we observe that tr 1 On the right-hand side, the terms with k l ¼ 1 are excluded, as they vanish since tr½W j−1 ¼ 0.
Given m-tuples of non-negative integers ⃗ k and ⃗ c 0 , let ⃗ k ≺ ⃗ c 0 denote the relation k l ≤ c 0 l for all l. Let ½ ⃗ k denote the permutation corresponding to a k 1 -cycle of the first k 1 integers, a k 2 -cycle of the next k 2 , etc., up to the last k m integers. Zero cycles are dropped, so ½ ⃗ k is an element of Thus, for ⃗ c 0 ∈ CS 0 n and τ ∈ ⃗ c 0 , Returning to the weights for Eq. (B12), we find that wð⃗ c; r; DÞ ¼ ðrDÞ mð⃗ cÞ Note that tr½P Here, Mð ⃗ k; r; sÞ ≔ Mð ⃗ k 0 ; r; sÞ, where ⃗ k 0 consists of the nonzero elements of ⃗ k listed in decreasing order. In principle, Eq. (B16) allows us to solve for all moments Mð⃗ c; r; sÞ. Note, in particular, the triangular structure to the equations: to solve for the moments with ⃗ c ∈ CS n , we need to know only the moments of lower order. Thus, the equations can be solved recursively for increasing n.
To analyze the solution to Eq. (B16) in the limits r → ∞ and D → ∞, we need the Weingarten function asymptotics (see Ref. [33]). Given α ∈ S n , with ½α ¼ ðc 1 c 2 …c m Þ, where jαj ¼ n − m, which is the minimal length of α as a product of 2-cycles, and Note that dðα; βÞ ¼ jαβ −1 j is a metric on S n . In particular, we have the triangle inequality Furthermore, jαj is invariant on conjugacy classes; in particular, jαj ¼ jα −1 j.

Gaussianity: Limit of large r
We now consider the asymptotics of Eq. (B16) as r → ∞. This example is technically somewhat simpler than the large D limit. Let us denote by AðrÞ ∼ BðrÞ the asymptotic relation lim r→∞ AðrÞ=BðrÞ ¼ 1. By Eq. (B17), By the triangle inequality [Eq. (B18)], jσj þ jτσj ≥ jτj. Furthermore, for any τ ∈ S 0 n , we have jτj ≥ n=2. Therefore, unless j ⃗ kj ¼ 0 and jσj þ jστj ¼ jτj ¼ n=2, it follows that Wgðτσ; rDÞr mðσÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ rD 2 p ðr − 1Þ Dr n−j ⃗ kj → 0 as r → ∞. The identity jσj þ jστj ¼ n=2 holds if and only if τ is a pairing (i.e., a product of n=2 disjoint 2-cycles) and if σ is a product of disjoint 2-cycles contained in τ (denoted σ ≺ τ). In this case, we have C τσ ¼ 1 and where P n denotes the set of pairings of f1; …; ng. It follows that all moments are of the order of one as r → ∞ and given In more detail, note that Wick's theorem would imply that lim r EfW ⊗2n j g is a sum of terms obtained by pairing the matrices in the tensor product and averaging according to the variance in Eq. (B7), namely, D . Further expanding each term in the sum according to a choice of one of the two terms in R ð2Þ for each pair of factors results in an expression for EfW ⊗2n j g as a sum over permutations with only twocycles and one-cycles. In this expansion, each permutation comes with the weight # pairings of the set covered by one cycles × ð−1Þ ð1=2Þ# one-cycles 1 D #cycles : The first term counts the number of ways that a given permutation can arise in the above expansion, while the second term comes from the weights attached to E ð2Þ D and E ð1;1Þ D in the variance. As one may easily check, the result is identical with the right-hand side of Eq. (B19).

Gaussianity: Large D limit
For the large D limit, we should first examine what happens to the moments like M½ð2Þ; r; D as D → ∞. Note that M½ð2Þ; r; D ¼ tr½P ½ð2Þ D R ð2Þ ¼ D − ð1=DÞ. Thus, we expect moments to become large, and a further normalization is needed to compute the limit of Eq. (B16). Indeed, we have Wgðτσ; rDÞD mð½⃗ cσÞ D #fk l ¼0g ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ rD 2 p ðr − 1Þ Dr n−j ⃗ kj ∼ ð−1Þ jτσj C τσ r ðn−j ⃗ kjÞ=2 r jσjþjτσj D #fk l ¼0g−j½⃗ cσj−jτσj : By the triangle inequality, j½⃗ cσj þ jτσj ≥ j½⃗ cτ −1 j ≥ 0, with equality precisely if ½⃗ c ¼ τ and σ ≺ τ. Everything is maximized for ⃗ k ¼ 0 and τ ¼ ½⃗ c, suggesting that In this section, we prove the formula for the bipartite entanglement entropy for ergodic MPS given by Eq. (17) in the main text. By translation invariance in distribution, it suffices to consider j ¼ 0. To focus on the bipartite entanglement entropy across a single cut, we make use of open boundary conditions for the chain.
Remark C.1.-The thermodynamic limit of the chain with open boundary conditions is identical to that with periodic boundary conditions; this limit follows from the convergence given in Theorem 2 of the main text. where Note that in the present caseÔ −M is a D × D matrix. It follows from Eq. (C2) and Theorem 2 of the main text that lim N→∞ hψð−M; NÞjOjψð−M; NÞi Let ρ −M denote the reduced density matrix of ½−M; 0, in the limit N → ∞. By Eq. (C3), we have hj −M ; …; j 0 jρ −M ji −M ; …; i 0 i Note that Q −M is a positive matrix. Let ⃗ v 1 ; …; ⃗ v D be its orthonormal eigenvectors, with corresponding eigenvalues λ 1 ≥ Á Á Á ≥ λ D ≥ 0. Let Thus, jαi j are eigenvectors of ρ −M , and the nonzero eigenvalues of ρ −M are given by λ i = P j λ j . It follows that the bipartite entropy across the bond 0 ∼ 1 (in the limit N → ∞ with M < ∞) is given by