# Theory of Ergodic Quantum Processes

Ramis Movassagh<EMAIL_ADDRESS>IBM Research, MIT-IBM AI lab, Cambridge MA,
02142, USA  Jeffrey Schenker<EMAIL_ADDRESS>Department of Mathematics,
Michigan State University, East Lansing MI, 48824, USA

###### Abstract

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. Random channels appear in a wide variety of applications,
from quantum chaos hosur2016chaos  to holographic dualities in theories of
quantum gravity hayden2016holographic  to operator dynamics von2018operator ,
to random local circuits for their potential to demonstrate quantum supremacy
arute2019quantum . 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 rank-one (entanglement breaking) channel.
Using this, we derive the limiting behaviour 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 2-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ényi 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 meta-stable, 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.

####  .1 Overview and Summary of the Results

Quantum channels represent physical changes to a quantum state and 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 watrous2018theory . Any particular evolution of a quantum system
is then formulated by the application of a sequence of the suitably chosen
quantum channels. Further, the formalism of quantum channels is fundamental
for the description of the physics of quantum matter. For example, the
expectation values and correlation functions in matrix product states (MPS)
are naturally expressed in this formalism. Quantum channels also give rise to
new counter-intuitive possibilities in communication theory smith2008quantum .

The generic (average case) behavior of quantum systems has long been of
theoretical and practical interest. Random channels appear in a wide variety
of applications, from quantum chaos hosur2016chaos  to holographic dualities
in theories of quantum gravity hayden2016holographic  to operator dynamics
von2018operator . Similarly, random local circuits are currently intensely
studied for their $k-$design properties harrow2009random , potential do
demonstrate quantum supremacy movassagh2019cayley ; bouland2019complexity  ,
and understanding operator spread and entanglement growth nahum2018operator .
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
boixo2018characterizing .

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 the lab especially now that
quantum error correction schemes have not been realized.

Figure 1: Venn diagram for distributions of possible ergodic sequences of
channels. Among the many subsets, three subsets are shown and substantially
enlarged for readability.

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 skinner2019measurement . Treating this probability as the order
parameter, it is then seen that there is a critical probability below which
the entanglement entropy is extensive and above which the state obeys an area
law bao2020theory .

We take the view of dynamical systems, a rich field with seminal results such
as the Furstenberg-Kesten theroem furstenberg1960  and the multiplicative
ergodic theorem of Oseledec oseledec1968multiplicative . We consider a
sequence of channels as a trajectory of a dynamical system. By bridging
quantum information to the theory of dynamical systems we address the
behaviour 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 includes and vastly generalizes
independently and identically distributed (iid), as well as translational
invariance (Figure 1). Below we give a more formal definition.

In this paper, 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 in each channel is an
iid random Haar isometry, 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. Finally, our formalism allows us to analyze the thermodynamic
limit of ergodic MPS. We derive formulas for the expectations of a general
observables and for the bipartite entanglement entropy across an arbitrary
cut, and prove exponential decay of two-point functions. The proofs and
derivations are given in the supplementary material (SM SM ), and in
movassagh2019ergodic .

In the past ’ergodic’ quantum channels were considered, which to the best of
our knowledge, assumed special subsets of possibilities considered herein
(Figure (1)). For example, a channel was chosen at random from some ensemble
and then repeatedly applied, i.e., $B_{k}$’s were all equal
burgarth2013ergodic , or time dynamics were analyzed for a quantum system with
repeated independently chosen random interactions with an environment
bruneau2014repeated . Others instances were studied such as certain random
independent channels and their compositions (e.g., from a finite set of random
isometries) collins2010random ; collins2011random . See collins2016random  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 pseudo-randomness emerson2003pseudo
generated by quasi-periodic dynamics.

####  .2 Ergodic theorem for composition of quantum channels

A quantum channel is a completely positive and trace preserving linear map,
$\phi$, on the space of $D\times D$ matrices:

$\phi(\rho)=\sum_{i=1}^{d}B^{i}\;\rho\;B^{i\>\dagger}.$ (1)

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 $\rho$)

$\phi_{n}\circ\cdots\circ\phi_{0}(\rho)=\sum_{i_{n},\dots,i_{0}}B_{n}^{i_{n}}\dots
B_{0}^{i_{0}}\text{ }\rho\text{ }B_{0}^{i_{0}\>\dagger}\dots
B_{n}^{i_{n}\>\dagger}.$ (2)

We will consider sequences $(\phi_{0},\phi_{1},\ldots)$ of channels drawn from
an ensemble, which we denote $\Omega$. Due to their vast generality, we take
the sequence of channels to be ergodic. 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
landau2013statistical . Later the ergodic hypothesis was made more precise in
mathematical physics culminating in the works of von Neumann and Birkhoff
moore2015ergodic . In those works it was recognized that the ergodic
hypothesis of Boltzmann was a consequence of a mathematical property of the
underlying dynamics. Formally, a map $T:\Omega\rightarrow\Omega$ on a
probability space $\Omega$ is ergodic if probabilities are invariant under $T$
and, starting from a typical sequence, the dynamics $T$ generates covers
$\Omega$ with probability one. In precise mathematical terms, $T$ is ergodic
provided that $T$ is measure preserving and $\operatorname{Prob}[{\cal A}]=0$
or $1$ for any event ${\cal A}$ with $T^{-1}({\cal A})={\cal A}$
cornfeld2012ergodic . We assume that the sequence of channels is drawn from an
ensemble in such a way that the shift
$(\phi_{0},\phi_{1},\ldots)\xmapsto{T}(\phi_{1},\phi_{2},\ldots)$ is ergodic.

In general, there may be selection rules separating distinct sectors of the
Hilbert space in which the quantum state $\rho$ resides, that do not mix under
the evolution (i.e., application of channels). Further, 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: $1.$ 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. $2.$ Has no transient subspace, because if it does then the
systems eventually leaves it with a probability arbitrary close to one.

Mathematically, both of these are guaranteed by the following two assumptions:
$1.$ For some $n_{0}$:
$\operatorname{Prob}\left[\phi_{n_{0}}\circ\cdots\circ\phi_{1}\circ\phi_{0}\text{
is strictly positive }\right]>0$, and $2.$ With probability one, if for some
observable $M\geq 0$, we have $\operatorname{tr}[M\;\phi_{0}(\rho)]=0$ for all
$\rho$, then $M\equiv 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 non-negligible 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 non-negative observable $M$, if
$\phi_{0}^{*}(M)=0$ then $M=0$.

For each $m,n\in\mathbb{Z}$ with $m<n$, let the quantum channel $\Psi_{n,m}$
be $\Psi_{n,m}\equiv\phi_{n}\circ\cdots\circ\phi_{m}$. Our main result is:

######  Theorem 1.

There exists $0<\mu<1$ and an ergodic sequence of $D\times D$ density matrices
$(\ldots,Z_{m},Z_{m+1},\ldots,Z_{n},Z_{n+1},\ldots)$ such that given
$x\in\mathbb{Z}$, the following holds

$\left\|\Psi_{n,m}(\rho)-Z_{n}\right\|\ \leq\ C_{\mu,x}\ \mu^{n-m}$ (3)

for all $m\leq x$, $n\geq x$, and any density matrix
$\rho\in\mathbb{C}^{D\times D}$, where $C_{\mu,x}<\infty$ almost surely.
Furthermore, the matrices $Z_{j}$ satisfy the shift equations

$Z_{j}=\phi_{j}(Z_{j-1})\;.$ (4)

$Z_{j}$’s 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, $\rho$, 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.

A key point is that for $n-m$ sufficiently large $\Psi_{n,m}$ is exponentially
close to the entanglement breaking (rank-one) channel
$\Psi_{n,-\infty}(\rho)=Z_{n}\operatorname{tr}[\rho]=Z_{n}$. This 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 the time step $n$. By ergodicty, the
distribution of $Z_{n}$ can be determined from the frequency of occurrences in
the trajectory $Z_{n},Z_{n+1},\ldots$. A trivial example is furnished by the
i.i.d. rank-$1$ channels $\phi_{j}$:

$\phi_{j}(\rho)\ =\ \begin{cases}\rho_{0}\operatorname{tr}[\rho]&\text{with
probability $1/2$}\\\ \rho_{1}\operatorname{tr}[\rho]&\text{with probability
$1/2$}\end{cases}\quad,$

with two distinct density matrices $\rho_{0}$ and $\rho_{1}$. In this case,
each $Z_{n}$ has the distribution

$\operatorname{Prob}[Z_{n}=\rho_{0}]=\operatorname{Prob}[Z_{n}=\rho_{1}]=\frac{1}{2},$

and the individual trajectories fluctuate randomly between the two values
$\rho_{0}$ and $\rho_{1}$, each occurring with frequency $1/2$.

To fully characterize the statistical equilibrium described by the process
$\ldots,Z_{j},\ldots$, we need to solve the recursion in Eq. (4). This
solution would depend 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_{k}^{i_{k}}$ in Eq. (2) are blocks of iid Haar unitaries.

####  .3 Applications

In the translationally invariant case the channels are the same;
$\phi_{j}=\phi_{0}$ and $Z_{j}=Z_{0}$ for all $j$. Furthermore, $Z_{0}$ is the
unique eigenmatrix of $\phi_{0}$ with eigenvalue $1$, and all other
eigenvalues of $\phi_{0}$ have modulus smaller than $1-\mu$. Eq. (3) shows the
exponentially fast convergence of the density matrix to the state $Z_{0}$
under repeated application of the channel $\phi_{0}$. That is starting from
any density matrix, repeated application of the channel derives the system
exponentially fast towards the equilibrium state $Z_{0}$.

In general, $Z_{j}$ is not an eigenmatrix of $\phi_{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 (SM SM ).

To illustrate this in a natural example, consider the other extreme and take
the sequence of channels
$\phi_{j}(\rho)=\operatorname{tr}_{r}[U_{j}\rho\otimes Q_{r}U_{j}^{\dagger}]$
where $Q_{r}$ is a pure state (rank-one projection) in ${\mathbb{C}}^{r\times
r}$ and $U_{j}$ is a sequence of independent Haar distributed $Dr\times Dr$
dimensional unitaries. Random Haar channels are the archetypes of generic
quantum channels. We show that the average of $Z_{j}$ is the totally mixed
state (SM SM )

$\mathbb{E}\left\\{Z_{j}\right\\}=D^{-1}\>I_{D}\;.$ (5)

Furthermore, using the Weingarten calculus Collins2003 ; collins2010random ,
we show that fluctuations around the average are given by

$Z_{j}=D^{-1}\>I_{D}+(1+rD^{2})^{-\frac{1}{2}}\>W_{j}\quad,$ (6)

where $\operatorname{tr}[W_{j}]=0$ and in the limit that $D\rightarrow\infty$
or $r\rightarrow\infty$ (or both), $W_{j}$ has a Gaussian distribution on the
space of $D\times D$ matrices, see SM SM . As a result, the eigenvalues of
$Z_{j}$ are distributed according to Wigner’s semi-circle law for large $D$ up
to corrections of order $1/r$. This allows analytical computation of spectral
properties of $Z_{j}$. For example, the von Neumann entropy of $Z_{j}$ has for
large $D$ and $r$ the limiting expression

$\mathcal{S}(Z_{j})\ \approx\ \log D-\frac{1}{8r}+O(r^{-2})\;.$ (7)

The entropy deviates from the maximal value $\log D$, which is attained for
maximally mixed state, by at most an order one quantity.

######  Remark.

Three physical corollaries to our theorem are:

  * •

Hastings and Koma proved that a constant gap implies an exponential decay of
correlations hastings2006spectral , and in one-dimension an area law for
entanglement entropy hastings2007area . Later Brandão and Horodecki
brandao2013area  proved that in one-dimension exponential decay of
correlations implies an area law. Below in section (.4) we prove a partial
converse, which says that finitely correlated states with an ergodic MPS
representation have an exponential decay of correlation.

  * •

Non-equilibrium phases of matter that are engineered by time-periodic driven
Hamiltonians, such as in Floquet systems li2019observation ;
lindner2011floquet ; shtanko2018stability , can only be meta-stable when
interactions with an environment in positive temperature are non-negligible
(equivalent of Assumption 1 above). Therefore, generically any experimentally
viable Floquet phase can only be meta-stable.

  * •

These channels become trivial as far as their quantum information content is
concerned. In fact, since the channels are asymptotically rank-1 and
entanglement breaking, the process cannot even convey classical information.
This 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 holds even in the time-translation invariant
case. For example, in the near-term quantum computing with noisy random
circuits boixo2018characterizing ; harrow2017quantum ; arute2019quantum , in
the limit of large depth, all initial memory of the state is lost and a unique
final state results.

####  .4 Matrix product states

MPS and their generalizations perez2006matrix ; verstraete2008matrix ;
vidal2003efficient ; vidal2008class  provide efficient representations of
quantum states by which classical simulation of quantum many-body systems
becomes viable and are the natural representation of density matrix
renormalization group white1992density  and its tensor network
generalizations. Applications range from efficient calculation of the ground
state properties of quantum matter chan2011density  to the outputs of quantum
circuits vidal2003efficient . They provide the tools for proving the existence
of satisfying assignments in qSAT movassagh2010unfrustrated  and the area law
hastings2007area . Random MPS are the candidate boundary states for recent
theoretical proposals of the theory of quantum gravity hayden2016holographic .

For the sake of concreteness let us introduce the (one-dimensional) MPS on
$2N+1$ qudits, denoting by $\vec{i}=\\{i_{-N},i_{-N+1}\dots,i_{N}\\}$:

$|\psi(N)\rangle=\sum_{\vec{i}=1}^{d}\text{Tr}[A_{-N}^{i_{-N}}A_{-N+1}^{i_{-N+1}}\cdots
A_{N}^{i_{N}}]\text{ }|i_{-N},\dots,i_{N}\rangle$ (8)

where $i_{k}$’s are the physical indices, $d$ is the local (physical)
dimension of the Hilbert space, $A_{k}^{i_{k}}$’s are $D\times D\times d$
tensors, and $D$ is the bond dimension perez2006matrix .

So far, rigorous results on generic MPS and their generalizations have mainly
focused on the translational invariant case, where all the tensors in Eq.
$\eqref{eq:MPS}$ are equal brandao2013area ; lancien2019correlation . Other
works focused on statistics of random MPS garnerone2010statistical . 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, would require that the tensors are only
drawn 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.

An observable $O$ on the spins in $[m,n]$ is a linear operator on
$\bigotimes_{j=m}^{n}\mathbb{C}^{d}$. The expectation of such an observable in
the state $|\psi(N)\rangle$ is given by the expression

$\langle\psi(N)|O|\psi(N)\rangle\ =\
\frac{\mathrm{Tr}\left[\phi_{N}\circ\cdots\circ\phi_{n+1}\circ\widehat{O}\circ\phi_{m-1}\circ\cdots\circ\phi_{-N}\right]}{\mathrm{Tr}\left[\phi_{N}\circ\cdots\circ\phi_{-N}\right]}$
(9)

where the linear maps $\phi_{m}$ are given by
$\phi_{m}(M)=\sum_{i=1}^{d}A_{m}^{i\>\dagger}\;M\;A_{m}^{i}$ and $\widehat{O}$
is the following linear operator on $\mathbb{C}^{D\times D}$

$\widehat{O}(M)=\displaystyle{\sum_{\vec{i},\vec{j}=1}^{d}}\Bigl{[}\langle
i_{m},\ldots,i_{n}|O|j_{m},\ldots,j_{n}\rangle\ A_{n}^{i_{n}\>\dagger}\cdots
A_{m}^{i_{m}\>\dagger}\;M\;A_{m}^{j_{m}}\cdots A_{n}^{j_{n}}\Bigr{]}\;.$ (10)

The linear maps $\phi_{m}$ defined above are analogous to $\phi_{i}$ in Eq.
(1), with $B^{i}=A_{m}^{i\>\dagger}$. However they may not be quantum
channels; although they are completely positive, they need not be trace
preserving. However, in Theorem I.4 of SM SM  we state and prove a
generalization of the Theorem 1 to non-trace preserving quantum operations
$\phi_{m}$ provided they satisfy $\phi_{m}^{*}(M)=0$ for $M\geq 0$ if and only
if $M=0$ in addition to the above assumptions. The main new feature of the
generalization is an ergodic analogue of the left eigenvector of a channel,
namely a sequence $Z_{j}^{\prime}$ satisfying

$Z_{j}^{\prime}=\frac{\phi_{j}^{*}(Z^{\prime}_{j+1})}{\operatorname{tr}[\phi_{j}^{*}(Z^{\prime}_{j+1})]}\
.$

In the trace invariant case, $Z_{j}^{\prime}$ is the maximally mixed state for
all $j$. Also, the shift equation Eq. (4) is replaced by
$Z_{j}=\frac{\phi_{j}(Z_{j-1})}{\operatorname{tr}[\phi_{j}(Z_{j-1})]},$
ensuring that the matrices $Z_{j}$ have trace $1$.

Using this result, we can directly calculate the thermodynamic limit of the
expectation value of an observable with respect to an MPS (Eq. (8)):

$W(O)\equiv\lim_{N\rightarrow\infty}\langle\psi(N)|O|\psi(N)\rangle=\frac{\mathrm{tr}\left[Z^{\prime}_{n+1}\>\widehat{O}(Z_{m-1})\right]}{\mathrm{tr}[Z^{\prime}_{n+1}\>\phi_{n}\circ\cdots\circ\phi_{m}(Z_{m-1})]}\;.$
(11)

Furthermore, the entanglement spectrum of a bipartite partition of the
infinite chain across the bond $j\sim j+1$, i.e., the spectrum of the reduced
density matrix for one-half of the chain, corresponds to the eigenvalues of
$R_{j}=Z_{j}Z_{j+1}^{\prime}/\operatorname{tr}[Z_{j}Z_{j+1}^{\prime}]$ (see SM
SM ). For instance, the bipartite von Neumann entanglement entropy is

$\mathcal{S}(j)\ =\ -\operatorname{tr}[R_{j}\log R_{j}]\;.$ (12)

Since $Z_{j}$ and $Z_{j+1}^{\prime}$ are positive, $R_{j}$ is similar to a
positive matrix and $\log R_{j}$ is well defined.

Turning 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:

######  Theorem 2.

movassagh2019ergodic  There exists $0<\mu<1$ such that given local observables
$O_{1}$ and $O_{2}$, the following correlation inequality holds

$\left|W\left(O_{1}(x)O_{2}(x+\ell)\right)-W(O_{1}(x))W(O_{2}(x+\ell))\right|\leq
C\ \mu^{|\ell|}$

with $C<\infty$ almost surely depending on $\mu,x,O_{1}$, and $O_{2}$.

In words, the 2-point correlation function between local observables decays
exponentially with the distance between their supports.

To the best of our knowledge, previous rigorous results only considered
translationally invariant MPS. Eq. (11,12) and Theorem (2) generalize these to
the much larger class of MPS generated by ergodic sequences of matrices. In
the translationally invariant case $\phi_{0}=\phi_{m}$, the sequences $Z_{n}$
and $Z^{\prime}_{m}$ are constant; i.e., $Z_{n+1}=Z_{n}$ and
$Z^{\prime}_{m}=Z^{\prime}_{m+1}$ for all $n,m$. And Eq. (11) reduces to

$W(O)\ =\
\frac{\mathrm{tr}\left[Z^{\prime}_{0}\>\widehat{O}(Z_{0})\right]}{\mathrm{tr}[Z^{\prime}_{0}\>\phi_{0}^{n-m}(Z_{0})]}\;,$

the entropy $\mathcal{S}(j)$ is independent of $j$, and the correlation bound
in Theorem (2) is independent of $x$.

Concluding remarks– In this letter we presented rigorous results on ergodic
sequences of quantum channels with non-negligible decoherence. The ergodic
framework of our results allows for noise models that are arbitrarily
correlated. It would be very interesting to investigate the consequences of
the results herein in the context of random circuits in connection to near-
term noisy quantum computers. Other applications would include chaos and
scrambling in connection to recent proposals for the theory of quantum gravity
and blackhole physics. Since matrix product state ansatz is intimitely
connected to finitely correlated states.

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 fannes1992finitely . 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 $\mu$ and the prefactor $C_{\mu,s}$
with respect the system size of channels and bond dimension of MPS.

###### Acknowledgements.

RM acknowledges the support of the IBM Research Frontiers Institute and
funding from the MIT-IBM Watson AI Lab under the project Machine Learning in
Hilbert Spaces. JS was supported by the National Science Foundation under
Grant No. 1500386 and Grant No. 190015.

### References

  * (1) See supplemnetary material (SM) at URL for the derivations and proofs.
  * (2) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, 2019.
  * (3) Yimu Bao, Soonwon Choi, and Ehud Altman.  Theory of the phase transition in random unitary circuits with measurements.  Physical Review B, 101(10):104301, 2020.
  * (4) Sergio Boixo, Sergei V Isakov, Vadim N Smelyanskiy, Ryan Babbush, Nan Ding, Zhang Jiang, Michael J Bremner, John M Martinis, and Hartmut Neven.  Characterizing quantum supremacy in near-term devices.  Nature Physics, 14(6):595–600, 2018.
  * (5) Adam Bouland, Bill Fefferman, Chinmay Nirkhe, and Umesh Vazirani.  On the complexity and verification of quantum random circuit sampling.  Nature Physics, 15(2):159–163, 2019.
  * (6) Fernando GSL Brandão and Michał Horodecki.  An area law for entanglement from exponential decay of correlations.  Nature Physics, 9(11):721–726, 2013.
  * (7) Laurent Bruneau, Alain Joye, and Marco Merkli.  Repeated interactions in open quantum systems.  Journal of Mathematical Physics, 55(7):075204, 2014.
  * (8) Daniel Burgarth, Giulio Chiribella, Vittorio Giovannetti, Paolo Perinotti, and Kazuya Yuasa.  Ergodic and mixing quantum channels in finite dimensions.  New Journal of Physics, 15(7):073045, 2013.
  * (9) Garnet Kin-Lic Chan and Sandeep Sharma.  The density matrix renormalization group in quantum chemistry.  Annual review of physical chemistry, 62:465–481, 2011.
  * (10) Benoît Collins.  Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability.  Int. Math. Res. Not., (17):953–982, 2003.
  * (11) Benoît Collins and Ion Nechita.  Random quantum channels i: Graphical calculus and the bell state phenomenon.  Communications in Mathematical Physics, 297(2):345–370, 2010.
  * (12) Benoît Collins and Ion Nechita.  Random quantum channels ii: Entanglement of random subspaces, rényi entropy estimates and additivity problems.  Advances in Mathematics, 226(2):1181–1201, 2011.
  * (13) Benoit Collins and Ion Nechita.  Random matrix techniques in quantum information theory.  Journal of Mathematical Physics, 57(1):015215, 2016.
  * (14) Isaac P Cornfeld, Sergej V Fomin, and Yakov Grigorévic Sinai.  Ergodic theory, volume 245 of Grundlehren der mathematischen Wissenschaften.  Springer-Verlag, 1982.
  * (15) Joseph Emerson, Yaakov S Weinstein, Marcos Saraceno, Seth Lloyd, and David G Cory.  Pseudo-random unitary operators for quantum information processing.  science, 302(5653):2098–2100, 2003.
  * (16) David E Evans and Raphael Høegh-Krohn.  Spectral properties of positive maps on c*-algebras.  Journal of the London Mathematical Society, 2(2):345–355, 1978\.
  * (17) Mark Fannes, Bruno Nachtergaele, and Reinhard F Werner.  Finitely correlated states on quantum spin chains.  Communications in mathematical physics, 144(3):443–490, 1992.
  * (18) H. Furstenberg and H. Kesten.  Products of random matrices.  Ann. Math. Statist., 31(2):457–469, 06 1960.
  * (19) Silvano Garnerone, Thiago R de Oliveira, Stephan Haas, and Paolo Zanardi.  Statistical properties of random matrix product states.  Physical Review A, 82(5):052312, 2010.
  * (20) Aram W Harrow and Richard A Low.  Random quantum circuits are approximate 2-designs.  Communications in Mathematical Physics, 291(1):257–302, 2009.
  * (21) Aram W Harrow and Ashley Montanaro.  Quantum computational supremacy.  Nature, 549(7671):203–209, 2017.
  * (22) Matthew B Hastings.  An area law for one-dimensional quantum systems.  Journal of Statistical Mechanics: Theory and Experiment, 2007(08):P08024, 2007.
  * (23) Matthew B Hastings and Tohru Koma.  Spectral gap and exponential decay of correlations.  Communications in mathematical physics, 265(3):781–804, 2006.
  * (24) Patrick Hayden, Sepehr Nezami, Xiao-Liang Qi, Nathaniel Thomas, Michael Walter, and Zhao Yang.  Holographic duality from random tensor networks.  Journal of High Energy Physics, 2016(11):9, 2016.
  * (25) H Hennion.  Limit theorems for products of positive random matrices.  The Annals of Probability, 25(4):1545–1587, 1997.
  * (26) Pavan Hosur, Xiao-Liang Qi, Daniel A Roberts, and Beni Yoshida.  Chaos in quantum channels.  Journal of High Energy Physics, 2016(2):4, 2016.
  * (27) J. F. C. Kingman.  Subadditive ergodic theory.  Ann. Probability, 1:883–909, 1973.
  * (28) M. G. Kreĭn and M. A. Rutman.  Linear operators leaving invariant a cone in a Banach space.  Amer. Math. Soc. Translation, 1950(26):128, 1950.
  * (29) Cécilia Lancien and David Pérez-García.  Correlation length in random mps and peps.  arXiv preprint arXiv:1906.11682, 2019.
  * (30) Lev D Landau and Evgeny M Lifshitz.  Statistical Physics: Volume 5, volume 5.  Elsevier, 2013.
  * (31) Jiaming Li, Andrew K Harter, Ji Liu, Leonardo de Melo, Yogesh N Joglekar, and Le Luo.  Observation of parity-time symmetry breaking transitions in a dissipative floquet system of ultracold atoms.  Nature communications, 10(1):1–7, 2019.
  * (32) Netanel H Lindner, Gil Refael, and Victor Galitski.  Floquet topological insulator in semiconductor quantum wells.  Nature Physics, 7(6):490–495, 2011.
  * (33) Calvin C Moore.  Ergodic theorem, ergodic theory, and statistical mechanics.  Proceedings of the National Academy of Sciences, 112(7):1907–1911, 2015.
  * (34) Ramis Movassagh.  Cayley path and quantum computational supremacy: A proof of average-case $\\#p-$ hardness of random circuit sampling with quantified robustness.  arXiv preprint arXiv:1909.06210, 2019.
  * (35) Ramis Movassagh, Edward Farhi, Jeffrey Goldstone, Daniel Nagaj, Tobias J Osborne, and Peter W Shor.  Unfrustrated qudit chains and their ground states.  Physical Review A, 82(1):012318, 2010.
  * (36) Ramis Movassagh and Jeffrey Schenker.  An ergodic theorem for homogeneously distributed quantum channels with applications to matrix product states.  arXiv preprint arXiv:1909.11769, 2019.
  * (37) Adam Nahum, Sagar Vijay, and Jeongwan Haah.  Operator spreading in random unitary circuits.  Physical Review X, 8(2):021014, 2018.
  * (38) V Oseledec.  A multiplicative ergodic theorem, characteristic lyapnov exponents of dynamical systems (transactions of the moscow mathematical society, 19).  American Mathematical Society, Providence, RI, 1968.
  * (39) David Perez-Garcia, Frank Verstraete, Michael M Wolf, and J Ignacio Cirac.  Matrix product state representations.  arXiv preprint quant-ph/0608197, 2006.
  * (40) Oles Shtanko and Ramis Movassagh.  Stability of periodically driven topological phases against disorder.  Physical Review Letters, 121(12):126803, 2018.
  * (41) Brian Skinner, Jonathan Ruhman, and Adam Nahum.  Measurement-induced phase transitions in the dynamics of entanglement.  Physical Review X, 9(3):031009, 2019.
  * (42) Graeme Smith and Jon Yard.  Quantum communication with zero-capacity channels.  Science, 321(5897):1812–1815, 2008.
  * (43) Frank Verstraete, Valentin Murg, and J Ignacio Cirac.  Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems.  Advances in Physics, 57(2):143–224, 2008.
  * (44) Guifré Vidal.  Efficient classical simulation of slightly entangled quantum computations.  Physical review letters, 91(14):147902, 2003.
  * (45) Guifré Vidal.  Class of quantum many-body states that can be efficiently simulated.  Physical review letters, 101(11):110501, 2008.
  * (46) CW Von Keyserlingk, Tibor Rakovszky, Frank Pollmann, and Shivaji Lal Sondhi.  Operator hydrodynamics, otocs, and entanglement growth in systems without conservation laws.  Physical Review X, 8(2):021013, 2018.
  * (47) John Watrous.  The theory of quantum information.  Cambridge University Press, 2018.
  * (48) Steven R White.  Density matrix formulation for quantum renormalization groups.  Physical review letters, 69(19):2863, 1992.

## Supplementary Material

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 [36], 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 the Eqs. (5)-(10) of the manuscript. Lastly, the
derivation of bipartite entanglement entropy (Eq. (17) in the manuscript) is
given.

###  I Formal Statement of the Ergodic Theorem

In this section we state a result (Theorem (I.3)) that implies Theorem (1) in
the manuscript. A sketch of the proof of Theorem (I.3) is given in Section
(II). The full proof of this result is in [36].

To proceed we need to introduce some basic notation and ideas. Let
${\mathbb{M}}_{D}$ denote the set of $D\times D$ complex matrices. We
primarily use the trace norm on ${\mathbb{M}}_{D}$:

$\left\|M\right\|_{1}\ =\ \operatorname{tr}[\>|M|\>]\quad.$

Recall that $|M|=\sqrt{M^{\dagger}M}$ where $\sqrt{\cdot}$ denotes the
operator square root. We also introduce the Hilbert-Schmidt inner product and
norm:

$\left\langle\widetilde{M},\
M\right\rangle=\operatorname{tr}\left[\widetilde{M}^{\dagger}M\right]\quad\text{and}\quad\left\|M\right\|_{2}^{2}=\left\langle
M,\ M\right\rangle=\operatorname{tr}\left[M^{\dagger}M\right]\quad.$

Recall that a _cone_ in a vector space is a set that is closed under addition
and multiplication by positive scalars. Let $\overline{E}$ denote the closed
cone

$\overline{E}=\left\\{\text{positive semi-definite }D\times D\text{
matrices}\right\\}\quad.$

The interior of $\overline{E}$ is

$E=\left\\{\text{positive definite }D\times D\text{ matrices}\right\\}\quad,$

which is an open cone.

Let $\mathcal{L}({\mathbb{M}}_{D})$ denote the set of all linear maps from
${\mathbb{M}}_{D}$ to ${\mathbb{M}}_{D}$. Given
$\phi\in\mathcal{L}({\mathbb{M}}_{D})$, the adjoint map $\phi^{*}$ is defined
by

$\operatorname{tr}\left[\widetilde{M}^{\dagger}\phi(M)\right]\ =\
\operatorname{tr}\left[[\phi^{*}(\widetilde{M})]^{\dagger}M\right]\quad.$

A map $\phi\in\mathcal{L}({\mathbb{M}}_{D})$ is _positive_ provided
$\phi(M)\in\overline{E}$ whenever $M\in\overline{E}$. If in addition
$\phi(M)\in E$ whenever $M\in\overline{E}$, then we say that $\phi$ is
_strictly positive_. A _completely positive_ map is one such that
$\phi\otimes{\mathbb{I}}_{N\times
N}:{\mathbb{M}}_{D}\otimes{\mathbb{M}}_{N}\rightarrow{\mathbb{M}}_{D}\otimes{\mathbb{M}}_{N}$
is positive for every $N$, where ${\mathbb{I}}_{N}$ denotes the identity map
on ${\mathbb{M}}_{N}$. Kraus’s theorem [47] states that $\phi$ is completely
positive if and only if $\phi(M)=\sum_{i=1}^{d}B^{i}\;M\;B^{i\>\dagger}$ for
some collection of matrices $B^{i}$, $i=1,\ldots,d$. A map $\phi$ is _trace
preserving_ if $\operatorname{tr}[\phi(M)]=\operatorname{tr}[M]$ for all $M$;
equivalently $\phi^{*}(I)=I$. A _quantum channel_ is a completely positive
trace preserving map.

Let $(\Omega,\mathcal{F},\operatorname{Prob})$ be a probability space with
$T:\Omega\rightarrow\Omega$ an invertible and ergodic map. Recall that $T$ is
_ergodic_ provided it is probability preserving and $\operatorname{Prob}[A]=0$
or $1$ for any event $A$ with $T^{-1}(A)=A$. Starting from a given completely
positive map valued random variable,
$\phi_{0}:\Omega\rightarrow\\{\text{completely positive maps}\\},$ we consider
the sequence of maps defined by evaluating $\phi_{0}$ along the trajectories
of $T$:

$\phi_{n;\omega}\ =\ \phi_{0;T^{n}\omega}\quad.$ (13)

We will follow the usual convention in probability theory and suppress the
independent variable $\omega\in\Omega$ in most formulas; when it is needed —
as in Eq.(13) — we will use a subscript to denote the value of a random
variable at a particular $\omega\in\Omega$. Our main focus here is to study
the action of the composition
$\phi_{n}\circ\phi_{n-1}\circ\cdots\circ\phi_{0}$ of a long sequence of these
maps. To obtain convergence, we require several assumptions. The first
assumption is Assumption 1 of the Letter, restated here for the reader’s
convenience

######  Assumption I.1.

For some $n_{0}>0$,
$\operatorname{Prob}\left[\phi_{n_{0}}\circ\cdots\circ\phi_{0}\text{ is
strictly positive }\right]\ >\ 0\ .$

######  Remark.

A map $\phi$ is strictly positive if and only if $\phi^{*}$ is strictly
positive. Indeed, if $\phi$ is strictly positive and $M\in\overline{E}$ we
have
$\operatorname{tr}[\phi^{*}(M)M^{\prime}]=\operatorname{tr}[M\phi(M^{\prime})]>0$
for any $M^{\prime}\in\overline{E}$, as $\phi(M^{\prime})>0$. Thus
$\phi^{*}(M)$ is strictly positive. Thus assumption I.1 implies that
$\phi_{0}^{*}\circ\cdots\circ\phi_{n_{0}}^{*}$ is strictly positive with
positive probability.

The second assumption, which for quantum channels is equivalent to Assumption
2 of the Letter, is the following:

######  Assumption I.2.

With probability one,
$\ker\phi_{0}\cap\overline{E}=\ker\phi_{0}^{*}\cap\overline{E}=\\{0\\}.$ That
is, if $\phi_{0}(M)=0$ or $\phi_{0}^{*}(M)=0$ with $M\in\overline{E}$, then
$M=0$.

######  Remark.

If $\phi_{0}$ is a quantum channel, then
$\operatorname{tr}[\phi_{0}(M)]=\operatorname{tr}[M]$ for any $M$, so
$\ker\phi_{0}\cap\overline{E}=0$. Thus for quantum channels, Assumption I.2 is
equivalent the equality $\ker\phi_{0}^{*}\cap\overline{E}=\\{0\\}$ (Assumption
2 of the Letter). On the other hand, this identity is a non-trivial assumption
for quantum channels. For example, if $\phi(M)=PMP+SMS^{\dagger}$ with $P$ a
projection onto a subspace of dimension $\nicefrac{{D}}{{2}}$ and $S$ a
partial isometry from $\mathbb{I}-P$ to $P$, then $\phi$ is a channel but
$\phi^{*}(\mathbb{I}-P)=0$.

Given $n\in\mathbb{Z}$, let

$\Phi_{n}\ =\ \begin{cases}\phi_{n}\circ\cdots\circ\phi_{0}&\text{if }n>0\
,\\\ \phi_{0}&\text{if }n=0\ ,\\\ \phi_{0}\circ\cdots\circ\phi_{n}&\text{if
}n<0\quad.\end{cases}$

By Assumption I.1 the maps $\Phi_{n}$ preserve the cone $E$. In 1948, Krein
and Rutman [28] derived a generalization of the classical Perron-Frobenius
theorem to the context of linear maps preserving a convex cone. Evans and
Høegh-Krohn [16] applied this result to obtain results for positive maps on
$\mathcal{L}({\mathbb{M}}_{D})$. It follows from [16, Theorem 2.3] that if
$\Phi_{n}$ is strictly positive then there is a unique (up to scaling)
strictly positive matrix $R_{n}$ such that

$\Phi_{n}(R_{n})=\lambda_{n}R_{n}\quad,$

where $\lambda_{n}$ is the spectral radius of $\Phi_{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

$\Phi_{n}^{*}(L_{n})=\lambda_{n}L_{n}\quad.$

We normalize $R_{n}$ and $L_{n}$ so that
$\operatorname{tr}[R_{n}]=\operatorname{tr}[L_{n}]=1\quad$.

Our main technical result is that $L_{n}$ converges as $n\rightarrow\infty$,
while $R_{n}$ converges as $n\rightarrow-\infty$. This generalizes a theorem
of Hennion on the Perron-Frobenius eigenvectors of products of positive
matrices [25].

######  Theorem I.3 ([36]).

The limits

$\lim_{n\rightarrow-\infty}R_{n}=Z_{0}\qquad\text{and}\qquad\lim_{n\rightarrow\infty}L_{n}=Z_{0}^{\prime}$

exist almost surely and define random matrices $Z_{0}$ and $Z_{0}^{\prime}$
that fall in $E$ with probability one. Furthermore, if we set
$Z_{n}=Z_{0;T^{n}\omega}$ and $Z^{\prime}_{n}=Z^{\prime}_{0;T^{n}\omega}$,
then

$Z_{n}=\phi_{n}\cdot Z_{n-1}\qquad\text{and}\qquad
Z_{n}^{\prime}=\phi_{n}^{*}\cdot Z_{n+1}^{\prime}\quad,$

where $\cdot$ denotes the projective action of a positive map on the strictly
positive $D\times D$ matrices of trace $1$:

$\phi_{n}\cdot M\ \equiv\
\frac{1}{\operatorname{tr}[\phi_{n}(M)]}\phi_{n}(M)\quad.$

######  Remark.

If the maps $\phi_{n}$ are quantum channels, then $L_{n}=\frac{1}{D}I_{D}$, so
$Z_{n}^{\prime}=\frac{1}{D}I_{D}$ for all $n$.

Given $m<n$, let $P_{n,m}$ denote the rank-one operator

$P_{n,m}(M)=\operatorname{tr}[Z^{\prime}_{m}M]\;Z_{n}\quad.$

For $n-m$ large, the operator $\phi_{n}\circ\cdots\circ\phi_{m}$ is well
approximated by $P_{n,m}$. To formulate this result precisely, we introduce
the following norm for a linear map $\Phi\in\mathcal{L}({\mathbb{M}}_{D})$:

$\left\|\Phi\right\|_{1}\ =\
\max\left\\{\operatorname{tr}[\>\left|\Phi(M)\right|\>]\quad\ \middle|\
\quad\operatorname{tr}[\>|M|\>]=1\right\\}.$

Theorem 2 of the Letter is equivalent to:

######  Theorem I.4 ([36]).

Given $m<n$ in $\mathbb{Z}$, let $\Psi_{n,m}=\phi_{n}\circ\cdots\phi_{m}$.
There is $0<\mu<1$ so that given $x\in{\mathbb{Z}}$, the following bound
holds:

$\left\|\frac{1}{\operatorname{tr}[\Psi_{n,m}^{*}(I_{D})]}\Psi_{n,m}-P_{n,m}\right\|_{1}\
\leq\ C_{\mu,x}\mu^{n-m}$

for all $m\leq x$ and $n\geq x$, with $C_{\mu,x}=C_{\mu,x;\omega}$ finite
almost surely.

###  II Sketch of the Proofs

The key to proving Theorem I.3 is to show the existence of the limits
$\lim_{n\rightarrow-\infty}R_{n}=Z_{0}$ and
$\lim_{n\rightarrow\infty}L_{n}=Z_{0}^{\prime}$. More generally, we show that

$\lim_{n\rightarrow-\infty}\Phi_{n}\cdot
Y_{n}\quad\text{and}\quad\lim_{n\rightarrow-\infty}\Phi_{n}^{*}\cdot Y_{n}$
(14)

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}\omega}$ and
$Z_{n}^{\prime}=Z_{0;T^{n}\omega}^{\prime}$ to obtain the full ergodic
sequence. Since $\Phi_{n}\cdot R_{n}=R_{n}$, it follows that

$Z_{0}=\lim_{n\rightarrow-\infty}R_{n}=\lim_{n\rightarrow-\infty}\Phi_{n}\cdot
R_{n}=\phi_{0}\cdot\left(\lim_{n\rightarrow-\infty}\Phi_{n+1;T^{-1}\omega}\cdot
R_{n}\right)\ =\ \phi_{0}\cdot Z_{-1}\quad.$

Similarly, $\Phi_{n}^{*}\cdot L_{n}=L_{n}$ and we find that

$Z_{0}^{\prime}=\lim_{n\rightarrow\infty}L_{n}=\lim_{n\rightarrow\infty}\Phi_{n}^{*}\cdot
L_{n}=\phi_{0}^{*}\cdot\left(\lim_{n\rightarrow\infty}\Phi_{n-1;T^{\omega}}\cdot
L_{n}\right)\ =\ \phi_{0}^{*}\cdot Z_{1}^{\prime}\quad.$

The key to proving the existence of the limits in Eq. (14) is to prove that
$\Phi_{n}$ is exponentially contractive with respect to a suitable metric
$d(Y_{1},Y_{2})$ on the space of density matrices:

$\lim_{n\rightarrow\infty}\frac{1}{n}\log\left[d(\Phi_{n}\cdot
Y_{1},\Phi_{n}\cdot Y_{2})\right]\ =\ \log\mu\ <\ 0.$ (15)

The existence of the limit in Eq. (15) is guaranteed by Kingman’s subadditve
ergodic theorem [27]. A pivotal estimate in this context is the following

$d(\phi\cdot X,\phi\cdot Y)\ \leq\ c(\phi)\;d(X,Y)\quad,$ (16)

where $c(\phi)\leq 1$ is a quantity that depends on the particular completely
positive map $\phi$. The subadditivity needed for the application of Kingman’s
theorem is provided by the inequality

$c(\phi\circ\psi)\leq c(\phi)\;c(\psi)\quad,$ (17)

which follows if we take $c(\phi)$ to be the smallest number such that Eq.
(16) holds. To obtain Eqs. (16,17), and to prove that the right hand side of
Eq. (15) is less than zero, it is convenient to work with the specialized
metric

$d(Y_{1},Y_{2})=\frac{1-m(Y_{1},Y_{2})m(Y_{2},Y_{1})}{1+m(Y_{1},Y_{2})m(Y_{2},Y_{2})}\quad,$

where $m(Y_{1},Y_{2})=\max\left\\{\lambda\in{\mathbb{R}}\ \middle|\ \lambda
Y_{2}\leq Y_{1}\right\\}$. 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 $\Phi_{n}$ are ultimately strictly
positive by Assumption 1 of the Letter, the estimates obtained in fact hold in
the trace-norm. For details, see [36].

###  III 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 Letter. The
second class corresponds to the Kraus operators being independently and
identically distributed (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\in\mathbb{Z}$ be an iid sequence of Haar distributed $L\times
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
$\phi_{j}\in\mathcal{L}({\mathbb{M}}_{D})$ given by

$\phi_{j}(M)\ =\ \operatorname{tr}_{r}[U_{j}\;(M\otimes
Q_{r})\;U_{j}^{\dagger}]\quad,$

where $Q_{r}$ is the $r\times r$ rank one projector
$Q_{r}\equiv\left|e_{1}\right>\left<e_{1}\right|$ with
$\left|e_{1}\right>=(1,0,\ldots,0)^{T}$, and $\operatorname{tr}_{r}$ denotes
the _partial trace_ from ${\mathbb{M}}_{L}\rightarrow{\mathbb{M}}_{D}$. Note
that the partial trace is the adjoint (with respect to the Hilbert-Schmidt
inner product) of the map $M\mapsto M\otimes I_{r}$. The channel $\phi_{j}$
can be expressed in the Kraus form as

$\phi_{j}(M)\ =\ \sum_{k=1}^{r}U_{j;k}\;M\;U_{j;k}^{\dagger}$

where the matrices $U_{j;k}$ are $D\times D$ blocks of $U_{j}$:

$U_{j}\ =\ \begin{pmatrix}U_{j;1}&\cdots&\cdots\\\ U_{j;2}&\ddots\\\
\vdots&&\ddots\\\ U_{j;r}&\cdots&\cdots\end{pmatrix}\quad.$

It is easy to see that this particular sequence of channels satisfies
Assumptions 1 and 2. In fact each $\phi_{j}$ is itself strictly positive with
probability one. Thus, Theorem 1 applies and the resulting stationary sequence
$Z_{j}$, $j\in{\mathbb{Z}}$ satisfies the shift equation ((Eq. 4) of the
Letter)

$Z_{j}\ =\ \operatorname{tr}_{r}\left[U_{j}\left(Z_{j-1}\otimes
Q_{r}\right)U_{j}^{\dagger}\right]\quad.$ (18)

We study the distribution of $Z_{j}$ by the method of moments. As such we are
interested in computing the expectation $\mathbb{E}\left\\{Z_{j}^{\otimes
n}\right\\}$, where

$Z_{j}^{\otimes n}\ =\ \underbrace{Z_{j}\otimes\cdots\otimes Z_{j}}_{n\text{
factors}}\quad.$

From $\mathbb{E}\left\\{Z_{j}^{\otimes n}\right\\}$ we can compute the
expectation of any multi-linear form
$\prod_{\ell=1}^{n}\operatorname{tr}[M_{\ell}Z_{j}]$, via the identity.

$\prod_{\ell=1}^{n}\operatorname{tr}[M_{\ell}Z_{j}]\ =\
\operatorname{tr}\left[\left(M_{1}\otimes\cdots\otimes
M_{n}\right)Z_{j}^{\otimes n}\right]\quad.$

We begin computing the expectation for $n=1$, which is the first moment
$\mathbb{E}\left\\{Z_{j}\right\\}$. The key identity here is:

$\int_{U(L)}U(\alpha,\beta)\;U^{\dagger}(\beta^{\prime},\alpha^{\prime})\;\mathrm{d}U\
=\
\frac{1}{L}\delta_{\alpha}^{\alpha^{\prime}}\delta_{\beta}^{\beta^{\prime}}\quad.$
(19)

In terms of matrices this can be expressed as

$\int_{U(L)}U\>M\>U^{\dagger}\;\mathrm{d}U\ =\
\operatorname{tr}[M]\,\frac{1}{L}\;I_{L}\quad,$

for an arbitrary $M\in{\mathbb{C}}^{L\times L}$. By Eq. (18), the conditional
expectation of $Z_{j}$ given $U_{k}$ for $k\neq j$ satisfies

$\mathbb{E}\left\\{Z_{j}\middle|U_{k},\ k\neq j\right\\}\ =\
\operatorname{tr}_{r}\left[\int_{U(L)}U\left(Z_{j-1}\otimes
Q_{r}\right)U^{\dagger}\;\mathrm{d}U\right]\ =\
\operatorname{tr}[Z_{j-1}\otimes
Q_{r}]\,\frac{1}{L}\operatorname{tr}_{r}\left[I_{L}\right]\ =\
\frac{1}{D}I_{D}\quad,$

since $\operatorname{tr}[Z_{j-1}\otimes Q_{r}]=\operatorname{tr}[Z_{j-1}]=1$
and $\operatorname{tr}_{r}[I_{L}]=rI_{D}$. Averaging over all $U_{k}$ for
$k\neq j$ we see that

$\mathbb{E}\left\\{Z_{j}\right\\}\ =\ \frac{1}{D}I_{D}\quad,$

which is Eq. (7) in the Letter.

We turn now to the computation of $\mathbb{E}\left\\{Z_{j}^{\otimes
2}\right\\}$. The key is the generalization of Eq. (19) to the following (see
[10]):

$\displaystyle\int_{U(L)}U(\alpha_{1},\beta_{1})U(\alpha_{2},\beta_{2})U^{\dagger}(\beta_{1}^{\prime},\alpha_{1}^{\prime})U^{\dagger}(\beta_{2}^{\prime},\alpha_{2}^{\prime})\;\mathrm{d}U\
=$ $\displaystyle\
\operatorname{Wg}((1,1),L)\left(\delta_{\alpha_{1}}^{\alpha_{1}^{\prime}}\delta_{\beta_{1}}^{\beta_{1}^{\prime}}\delta_{\alpha_{2}}^{\alpha_{2}^{\prime}}\delta_{\beta_{2}}^{\beta_{2}^{\prime}}+\delta_{\alpha_{1}}^{\alpha_{2}^{\prime}}\delta_{\beta_{1}}^{\beta_{2}^{\prime}}\delta_{\alpha_{2}}^{\alpha_{1}^{\prime}}\delta_{\beta_{2}}^{\beta_{1}^{\prime}}\right)$
(20) $\displaystyle\
+\operatorname{Wg}((2),L)\left(\delta_{\alpha_{1}}^{\alpha_{1}^{\prime}}\delta_{\beta_{1}}^{\beta_{2}^{\prime}}\delta_{\alpha_{2}}^{\alpha_{2}^{\prime}}\delta_{\beta_{2}}^{\beta_{1}^{\prime}}+\delta_{\alpha_{1}}^{\alpha_{2}^{\prime}}\delta_{\beta_{1}}^{\beta_{1}^{\prime}}\delta_{\alpha_{2}}^{\alpha_{1}^{\prime}}\delta_{\beta_{2}}^{\beta_{2}^{\prime}}\right)\quad,$

where the _Weingarten functions_ are

$\operatorname{Wg}((1,1),L)\ =\
\frac{1}{L^{2}-1}\quad\text{and}\quad\operatorname{Wg}((2),L)\ =\
-\frac{1}{L(L^{2}-1)}\quad.$

In terms of tensors, Eq. (20) can be expressed as follows

$\displaystyle\int_{U(L)}(U\otimes U)\,M\,(U^{\dagger}\otimes
U^{\dagger})\;\mathrm{d}U\ =$ $\displaystyle\
\frac{1}{L^{2}-1}\left(\operatorname{tr}[M]\,E_{L}^{(1,1)}+\operatorname{tr}[E_{L}^{(2)}M]\,E_{L}^{(2)}\right)$
$\displaystyle\
-\frac{1}{L(L^{2}-1)}\left(\operatorname{tr}[M]\,E_{L}^{(2)}+\operatorname{tr}[E_{L}^{(2)}M]\,E_{L}^{(1,1)}\right)\quad,$

where $M$ is an arbitrary element of
${\mathbb{M}}_{L}\otimes{\mathbb{M}}_{L}$, $E_{L}^{(1,1)}=I_{L}\otimes I_{L}$
is the identity, and $E_{L}^{(2)}$ denotes the exchange operator

$E_{L}^{(2)}(v\otimes w)\ =\ w\otimes v\quad.$

######  Remark III.1.

$E_{L}^{(2)}$ is the SWAP operator.

Applying this identity to $Z_{j}^{\otimes 2}=Z_{j}\otimes Z_{j}$ we find that

$\displaystyle\mathbb{E}\left\\{Z_{j}^{\otimes 2}\middle|U_{k},\ k\neq
j\right\\}=$ $\displaystyle\ \operatorname{tr}_{r}^{\otimes
2}\left[\int_{U(L)}U^{\otimes 2}\,(Z_{j-1}\otimes Q_{r})^{\otimes
2}\,U^{\dagger\;\otimes 2}\;\mathrm{d}U\right]$ $\displaystyle=$
$\displaystyle\ \frac{1}{L^{2}-1}\left(\operatorname{tr}\left[(Z_{j-1}\otimes
Q_{r})^{\otimes 2}\right]\operatorname{tr}_{r}^{\otimes
2}\left[E_{L}^{(1,1)}\right]+\operatorname{tr}\left[E_{L}^{(2)}(Z_{j-1}\otimes
Q_{r})^{\otimes 2}\right]\operatorname{tr}_{r}^{\otimes
2}\left[E_{L}^{(2)}\right]\right)$ $\displaystyle\
-\frac{1}{L(L^{2}-1)}\left(\operatorname{tr}\left[E_{L}^{(2)}(Z_{j-1}\otimes
Q_{r})^{\otimes 2}\right]\operatorname{tr}_{r}^{\otimes
2}\left[E_{L}^{(2)}\right]+\operatorname{tr}\left[(Z_{j-1}\otimes
Q_{r})^{\otimes 2}\right]\operatorname{tr}_{r}^{\otimes
2}\left[E_{L}^{(1,1)}\right]\right)\;,$

where we introduced the notation $\operatorname{tr}_{r}^{\otimes
2}[M_{1}\otimes
M_{2}]=\operatorname{tr}_{r}[M_{1}]\otimes\operatorname{tr}_{r}[M_{2}]$,
extended by linearity to all of ${\mathbb{M}}_{L}\otimes{\mathbb{M}}_{L}$. Now

$\operatorname{tr}\left[(Z_{j-1}\otimes Q_{r})^{\otimes 2}\right]\ =\
\left(\operatorname{tr}[Z_{j-1}\otimes Q_{r}]\right)^{2}=1\quad,$

and

$\operatorname{tr}_{r}^{\otimes 2}\left[E_{L}^{(11)}\right]\ =\
\operatorname{tr}_{r}[I_{L}]\otimes\operatorname{tr}_{r}[I_{L}]\ =\
r^{2}I_{D}\otimes I_{D}=r^{2}E_{D}^{(1,1)}\quad.$

For the factors involving $E^{(2)}_{L}$ we make use of the identity

$\operatorname{tr}[M_{1}M_{2}]\ =\ \operatorname{tr}[E_{L}^{(2)}M_{1}\otimes
M_{2}]\quad.$ (21)

Thus

$\operatorname{tr}\left[E_{L}^{(2)}(Z_{j-1}\otimes Q_{r})^{\otimes 2}\right]\
=\ \operatorname{tr}\left[(Z_{j-1}\otimes Q_{r})^{2}\right]\ =\
\operatorname{tr}[Z_{j-1}^{2}]\quad,$

and

$\operatorname{tr}_{r}^{\otimes 2}\left[E_{L}^{(2)}\right]\ =\
rE_{D}^{(2)}\quad.$

It follows that

$\mathbb{E}\left\\{Z_{j}^{\otimes 2}\middle|U_{k},\ k\neq j\right\\}\ =\
\frac{1}{L(L^{2}-1)}\Biggl{[}\left(L-\operatorname{tr}[Z_{j-1}^{2}]\right)r^{2}E_{D}^{(1,1)}+\left(L\operatorname{tr}[Z_{j-1}^{2}]-1\right)rE_{D}^{(2)}\Biggr{]}\quad.$

Treating $\operatorname{tr}[Z_{j-1}^{2}]$ as the random variable, and
averaging over $U_{k}$, we obtain

$\mathbb{E}\left\\{Z_{j}^{\otimes
2}\right\\}=\frac{1}{L(L^{2}-1)}\Biggl{[}\left(L-m_{(2)}\right)r^{2}E_{D}^{(1,1)}+\left(Lm_{(2)}-1\right)rE_{D}^{(2)}\Biggr{]}\quad.$
(22)

where

$m_{(2)}\ \equiv\ \mathbb{E}\left\\{\operatorname{tr}[Z_{j-1}^{2}]\right\\}\
=\ \mathbb{E}\left\\{\operatorname{tr}[Z_{j}^{2}]\right\\}\quad,$

since $Z_{j-1}$ and $Z_{j}$ are equidistributed.

To complete the calculation of $\mathbb{E}\left\\{Z_{j}^{\otimes 2}\right\\}$
we need to compute $m_{(2)}$. To this end, note that by Eq. (21)

$m_{(2)}=\mathbb{E}\left\\{\operatorname{tr}[E_{D}^{(2)}Z_{j}^{\otimes
2}]\right\\}\quad.$

Thus

$\displaystyle m_{(2)}\ =$ $\displaystyle\
\frac{1}{L(L^{2}-1)}\Biggl{(}\left(L-m_{(2)}\right)r^{2}\operatorname{tr}[E_{D}^{(2)}]+\left(Lm_{(2)}-1\right)r\operatorname{tr}[E_{D}^{(1,1)}]\Biggr{)}$
$\displaystyle=$ $\displaystyle\
\frac{1}{L(L^{2}-1)}\Biggl{(}\left(L-m_{(2)}\right)r^{2}D+\left(Lm_{(2)}-1\right)rD^{2}\Biggr{)}\quad,$

since $E_{D}^{(2)}E_{D}^{(1,1)}=E_{D}^{(2)}$ and
$E_{D}^{(2)}E_{D}^{(2)}=E_{D}^{(1,1)}$. Solving for $m_{(2)}$ we find that

$m_{(2)}\ =\ \frac{(r+1)D}{rD^{2}+1}\quad.$

Plugging this into Eq. (22) we find

$\mathbb{E}\left\\{Z_{j}^{\otimes
2}\right\\}=\frac{1}{rD^{2}+1}\Biggl{(}rE_{D}^{(1,1)}+\frac{1}{D}E_{D}^{(2)}\Biggr{)}\quad.$

The variance of $Z_{j}$ can now be expressed as

$\mathbb{E}\left\\{Z_{j}^{\otimes
2}\right\\}-\mathbb{E}\left\\{Z_{j}\right\\}\otimes\mathbb{E}\left\\{Z_{j}\right\\}\
=\
\frac{1}{rD^{2}+1}\left(\frac{1}{D}E_{D}^{(2)}-\frac{1}{D^{2}}E_{D}^{(1,1)}\right)\quad.$

Note that
$\operatorname{tr}[\frac{1}{D}E_{D}^{(2)}]=\operatorname{tr}[\frac{1}{D^{2}}E_{D}^{(1,1)}]=1$.
Thus the prefactor $(rD^{2}+1)^{-1}$ sets the scale for the fluctuations
around the mean $\mathbb{E}\left\\{Z_{j}\right\\}=\frac{1}{D}E_{D}^{(11)}$. To
proceed we introduce the notation $W_{j}$ for the rescaled deviation of
$Z_{j}$ from the mean:

$W_{j}\ =\ \sqrt{rD^{2}+1}\left(Z_{j}-\frac{1}{D}I_{D}\right)\quad.$ (23)

Thus

$\mathbb{E}\left\\{W_{j}\right\\}=0\quad\text{and}\quad\mathbb{E}\left\\{W_{j}^{\otimes
2}\right\\}\ =\ \frac{1}{D}E_{D}^{(2)}-\frac{1}{D^{2}}E_{D}^{(11)}\quad.$ (24)

Our goal is to show that $W_{j}$ is asymptotically Gaussian, with covariance
$\frac{1}{D}E_{D}^{(2)}-\frac{1}{D^{2}}E_{D}^{(11)}$. To obtains this result,
we will verify that $W_{j}$ asymptotically satisfies Wick’s theorem.

Turning now to the computation of the $n$-th moment for general $n$, we define

$R^{(n)}\ =\ \mathbb{E}\left\\{W_{j}^{\otimes n}\right\\}\quad.$

By the shift invariance property, the right hand side is independent of $j$
and we have

$R^{(n)}\ =\
\left(rD^{2}+1\right)^{\nicefrac{{n}}{{2}}}\,\mathbb{E}\left\\{(Z_{j}-\frac{1}{D}I_{D})^{\otimes
n}\right\\}\ =\
\left(rD^{2}+1\right)^{\nicefrac{{n}}{{2}}}\,\mathbb{E}\left\\{\left(\operatorname{tr}_{r}\left[U_{j}Z_{j-1}\otimes
Q_{r}U_{j}^{\dagger}\right]-\frac{1}{D}I_{D}\right)^{\otimes
n}\right\\}\quad.$ (25)

Because $\operatorname{tr}_{r}\left[U_{j}Z_{j-1}\otimes
Q_{r}U_{j}^{\dagger}\right]-\frac{1}{D}I_{D}=\operatorname{tr}_{r}\left[U_{j}\left(Z_{j-1}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)U_{j}^{\dagger}\right]$, we have

$\left(\operatorname{tr}_{r}\left[U_{j}Z_{j-1}\otimes
Q_{r}U_{j}^{\dagger}\right]-\frac{1}{D}I_{D}\right)^{\otimes n}\ =\
\operatorname{tr}_{r}^{\otimes n}\left[U_{j}^{\otimes n}\left(Z_{j-1}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)^{\otimes n}U_{j}^{\dagger\otimes
n}\right]\quad.$

There is a generalization of Eq. (20) to higher order products (see [10]):

$\int_{U(L)}\prod_{j=1}^{n}U(\alpha_{j},\beta_{j})\;U^{\dagger}(\beta_{j}^{\prime},\alpha_{j}^{\prime})\;\mathrm{d}U\
=\ \sum_{\sigma,\tau\in
S_{n}}\operatorname{Wg}(\tau\sigma^{-1},L)\prod_{j=1}^{n}\delta_{\alpha_{j}}^{\alpha_{\sigma(j)}^{\prime}}\delta_{\beta_{j}}^{\beta_{\tau(j)}^{\prime}}\quad,$
(26)

where $\operatorname{Wg}$ denotes the Weingarten function (see [10]) and
$S_{n}$ is the permutation group on $\\{1,\ldots,n\\}$. In terms of tensors,
this equation can be written as

$\int_{U(L)}U^{\otimes n}\,M\,U^{\dagger\otimes n}\;\mathrm{d}U\ =\
\sum_{\sigma,\tau\in
S_{n}}\operatorname{Wg}(\tau\sigma,L)\operatorname{tr}\left[P_{L}^{\tau}M\right]\
P_{L}^{\sigma}\quad,$ (27)

where the operator $P^{\alpha}_{L}$ has the action:

$P_{L}^{\alpha}(v_{1}\otimes\cdots\otimes v_{j}\otimes\cdots\otimes v_{n})\ =\
v_{\alpha(1)}\otimes\cdots\otimes v_{\alpha(j)}\otimes\cdots
v_{\alpha(n)}\quad,$

for $v_{k}\in{\mathbb{C}}^{L}$. Applying Eq. (27) to the average over $U_{j}$
in Eq. (25) results in the following expression for $R^{(n)}$

$R^{(n)}\ =\ \left(rD^{2}+1\right)^{\frac{n}{2}}\sum_{\sigma\in
S_{n}}\left(\sum_{\tau\in
S_{n}}\operatorname{Wg}(\tau\sigma,L)\,\mathbb{E}\left\\{\operatorname{tr}\left[P_{L}^{\tau}\left(Z_{j-1}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)^{\otimes
n}\right]\right\\}\right)\operatorname{tr}_{r}^{\otimes
n}[P_{L}^{\sigma}]\quad.$ (28)

The Weingarten function $\operatorname{Wg}(\alpha,L)$ depends only on the
conjugacy class of the permutation $\alpha$ (see [10]):

$\operatorname{Wg}(\tau\alpha\tau^{-1},L)\ =\
\operatorname{Wg}(\alpha,L)\quad.$

Furthermore, for any $M\in{\mathbb{M}}_{L}$,
$\operatorname{tr}\left[P_{L}^{\tau}M^{\otimes n}\right]$ depends only on the
conjugacy class of $\tau$, as we will now show. The conjugacy classes of
$S_{n}$ are labeled by $m$-tuples $\vec{c}=(c_{1},c_{2}\cdots,c_{m})$ with
$c_{1}\geq c_{2}\geq\cdots\geq c_{m}\geq 1$ and $\sum_{j}c_{j}=n$. A
permutation in the class $(c_{1},c_{2}\cdots,c_{m})$ is given by

  1. 1.

a partition of $\\{1,\ldots,n\\}$ into $m$ sets of sizes $c_{1},\ldots,c_{m}$,
and

  2. 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 $\tau\in S_{n}$
in the class $(c_{1},c_{2}\cdots,c_{m})$ and $M\in{\mathbb{M}}_{L}$, the trace
$\operatorname{tr}\left[P_{L}^{\tau}M^{\otimes n}\right]$ is given by

$\operatorname{tr}\left[P_{L}^{\tau}M^{\otimes n}\right]\ =\
\prod_{j=1}^{m}\operatorname{tr}\left[M^{c_{j}}\right]\quad.$ (29)

Since $\operatorname{tr}[Z_{j-1}\otimes Q_{r}]=\operatorname{tr}[Z_{j-1}]=1$,
we see from Eq. (29) that $\operatorname{tr}[P_{L}^{\tau}(Z_{j-1}\otimes
Q_{r}-\frac{1}{L}I_{L})^{\otimes n}]=0$ for any permutation $\tau$ with a
fixed point, in which case the conjugacy class has $c_{m}=1$. In other words,
the sum over $\tau$ in Eq. (28) can be restricted to $\tau\in S_{n}^{\prime}$,
where

$S_{n}^{\prime}\ =\ \left\\{\tau\in S_{n}\ \middle|\ \tau(j)\neq j,\
j=1,\ldots,n\right\\}\quad.$

Let $CS_{n}$ denote the set of conjugacy classes of $S_{n}$. Given $\vec{c}\in
CS_{n}$ let

$E_{L}^{\vec{c}}=\sum_{\sigma\in\vec{c}}P_{L}^{\sigma}.$

With the labeling of classes described above this is consistent with the
definitions of $E_{L}^{(1,1)}$ and $E_{L}^{(2)}$ given above. Associated to
each class $\vec{c}=(c_{1},c_{2},\ldots,c_{m})\in CS_{n}$ we define a
distinguished permutation $[\vec{c}]$ which consists of the cycle

$1\mapsto 2\mapsto\cdots\mapsto c_{1}\mapsto 1$

of the first $c_{1}$ numbers, followed by the cycle

$c_{1}+1\mapsto c_{1}+2\mapsto\cdots\mapsto c_{1}+c_{2}\mapsto c_{1}+1$

of the next $c_{2}$ numbers, etc. up to the last group of $c_{m}$ numbers,

$\sum_{j=1}^{m-1}c_{j}+1\mapsto\sum_{j=1}^{m-1}c_{j}+2\mapsto\cdots\mapsto
n\mapsto\sum_{j=1}^{m-1}c_{j}\quad.$

That is,

$[\vec{c}]=(1,2,\ldots,c_{1})(c_{1}+1,\ldots,c_{1}+c_{2})\cdots(c_{1}+\ldots+c_{m-1},\ldots,n)\quad.$

Putting all of this together we obtain the following formula

$R^{(n)}\ =\mathbb{E}\left\\{W_{j}^{\otimes n}\right\\}\ =\ \sum_{\vec{c}\in
CS_{n}}w(\vec{c};r,D)\,\frac{1}{D^{m(\vec{c})}}E_{D}^{\vec{c}}\quad,$ (30)

where $m(\vec{c})$ denotes the number of cycles in the conjugacy class
$\vec{c}$, i.e., $m(\vec{c})=m$ if $\vec{c}=(c_{1},\ldots,c_{m})$ with
$c_{m}\geq 1$,

$w(\vec{c};r,D)\ =\
L^{m(\vec{c})}\left(rD^{2}+1\right)^{\frac{n}{2}}\sum_{\tau\in
S_{n}^{\prime}}\operatorname{Wg}(\tau[\vec{c}],L)\>\mathbb{E}\left\\{\operatorname{tr}\left[P_{L}^{\tau}\left(Z_{j-1}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)^{\otimes n}\right]\right\\}\quad,$ (31)

and we have observed that $\operatorname{tr}_{r}^{\otimes
n}[E_{L}^{\vec{c}}]=r^{m(\vec{c})}E_{D}^{\vec{c}}$. Adding and subtracting
$\frac{1}{D}I_{D}\otimes Q_{r}\ =\ \mathbb{E}\left\\{Z_{j-1}\right\\}\otimes
Q_{r}$ inside the parentheses, we obtain the following alternative expression
for $w(\vec{c;r,s})$,

$w(\vec{c};r,D)\ =\ L^{m(\vec{c})}\sum_{\tau\in
S_{n}^{\prime}}\operatorname{Wg}(\tau[\vec{c}],L)\;\mathbb{E}\left\\{\operatorname{tr}\left[P_{L}^{\tau}\left(W_{j-1}\otimes
Q_{r}+\sqrt{rD^{2}+1}\left(\frac{1}{D}I_{D}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)\right)^{\otimes n}\right]\right\\}\quad.$ (32)

Let $\tau\in S_{n}^{\prime}$. Then its conjugacy class
$[\tau]=(c_{1}^{\prime}c_{2}^{\prime}\ldots c_{m}^{\prime})$ satisfies
$c_{m}^{\prime}\geq 2$. Given commuting matrices $A$ and $B$,

$\operatorname{tr}\left[P_{L}^{\tau}(A+B)^{\otimes
n}\right]=\prod_{\ell=1}^{m}\sum_{k_{\ell}=0}^{c_{\ell}^{\prime}}{c_{\ell}^{\prime}\choose
k_{\ell}}\operatorname{tr}\left[A^{k_{\ell}}B^{c_{\ell}^{\prime}-k_{\ell}}\right]\quad.$

Applying this to the expectations in the sum on the right hand side of Eq.
(32) we obtain

$\mathbb{E}\left\\{\operatorname{tr}\left[P_{L}^{\tau}\left(W_{j-1}\otimes
Q_{r}+\sqrt{rD^{2}+1}\left(\frac{1}{D}I_{D}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)\right)^{\otimes n}\right]\right\\}\\\ =\
\mathbb{E}\left\\{\prod_{\ell=1}^{m}\biggl{[}\left(1-(1-r)^{1-c_{\ell}^{\prime}}\right)D\left(\frac{(r-1)\sqrt{rD^{2}+1}}{L}\right)^{c_{\ell}^{\prime}}+\sum_{k_{\ell}=2}^{c_{\ell}^{\prime}}{c_{\ell}^{\prime}\choose
k_{\ell}}\left(\frac{(r-1)\sqrt{rD^{2}+1}}{L}\right)^{c_{\ell}^{\prime}-k_{\ell}}\operatorname{tr}[W_{j-1}^{k_{\ell}}]\biggr{]}\right\\}\quad,$

where we have observed that

  1. 1.

$\operatorname{tr}\left[\left(\frac{1}{D}I_{D}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)^{c}\right]=\left(\frac{1}{D}-\frac{1}{L}\right)^{c}D+\left(-\frac{1}{L}\right)^{c}(r-1)D\
=\ \left(\frac{r-1}{L}\right)^{c}\left(1-(1-r)^{1-c}\right)D\quad$, and

  2. 2.

$(W_{j-1}\otimes Q_{r})(\frac{1}{D}I_{D}\otimes Q_{r}-\frac{1}{L}I_{L})\ =\
\frac{r-1}{L}W_{j-1}\otimes Q_{r}\quad$.

On the right hand side the terms with $k_{\ell}=1$ are excluded because they
vanish since $\operatorname{tr}[W_{j-1}]=0$.

Given $m$-tuples of non-negative integers $\vec{k}$ and $\vec{c}^{\prime}$,
let $\vec{k}\prec\vec{c}^{\prime}$ denote the relation $k_{\ell}\leq
c^{\prime}_{\ell}$ for all $\ell.$ Let $[\vec{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 $[\vec{k}]$ is an element of $S_{|\vec{k}|}$, where
$|\vec{k}|=\sum_{\ell}k_{\ell}$. Note that

$\prod_{\ell=1}^{m}\operatorname{tr}[W_{j-1}^{k_{\ell}}]=\operatorname{tr}[P_{D}^{[\vec{k}]}W^{(|\vec{k}|)}_{j-1}]\quad.$

Thus, for $\vec{c}^{\prime}\in CS_{n}^{\prime}$ and $\tau\in\vec{c}^{\prime}$

$\mathbb{E}\left\\{\operatorname{tr}\left[P_{L}^{\tau}\left(W_{j-1}\otimes
Q_{r}+\sqrt{rD^{2}+1}\left(\frac{1}{D}I_{D}\otimes
Q_{r}-\frac{1}{L}I_{L}\right)\right)^{\otimes n}\right]\right\\}\\\ =\
\sum_{\begin{subarray}{c}\vec{k}\prec c^{\prime}\\\ k_{\ell}\neq
1\end{subarray}}{\vec{c}^{\prime}\choose\vec{k}}\left(\frac{(r-1)\sqrt{rD^{2}+1}}{rD}\right)^{n-|\vec{k}|}D^{\\#\\{k_{\ell}=0\\}}\prod_{\ell\
:\
k_{\ell}=0}\left(1-(1-r)^{1-c_{\ell}^{\prime}}\right)\operatorname{tr}[P_{D}^{[\vec{k}]}R^{(|\vec{k}|)}]\quad,$

where
${\vec{c}^{\prime}\choose\vec{k}}=\prod_{\ell=1}^{m}{c_{\ell}^{\prime}\choose
k_{\ell}}.$ Returning to the weights for Eq. (30) we find that

$w(\vec{c};r,D)\ =\ (rD)^{m(\vec{c})}\sum_{\tau\in
S_{n}^{\prime}}\Biggl{(}\operatorname{Wg}(\tau[\vec{c}],rD)\\\
\times\sum_{\begin{subarray}{c}\vec{k}\prec[\tau]\\\ k_{\ell}\neq
1\end{subarray}}{[\tau]\choose\vec{k}}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{rD}\right)^{n-|\vec{k}|}D^{\\#\\{k_{\ell}=0\\}}\prod_{\ell\
:\
k_{\ell}=0}\left(1-(1-r)^{1-c_{\ell}^{\prime}(\tau)}\right)\operatorname{tr}[P_{D}^{[\vec{k}]}R^{(|\vec{k}|)}]\Biggr{)}\quad.$
(33)

Note that $\operatorname{tr}[P_{s}^{[\vec{k}]}R^{(|\vec{k}|)}]$ depends only
on the non-zero elements of $\vec{k}$ and is unchanged if we permute the
entries of $\vec{k}$.

Given $\vec{c}\in CS_{n}^{\prime}$, let
$M(\vec{c};r,D)=\operatorname{tr}[P_{D}^{[\vec{c}]}R^{(n)}]$. Using Eq. (30)
to compute $M(\vec{c};r,D)$, we find

$M(\vec{c};r,D)\ =\ \sum_{\sigma\in
S_{n}}w([\sigma];r,D)\frac{1}{D^{m(\sigma)}}\operatorname{tr}[P_{D}^{[\vec{c}]}P_{D}^{\sigma}]\
=\ \sum_{\sigma\in
S_{n}}w([\sigma];r,D)\frac{D^{m([\vec{c}]\sigma)}}{D^{m(\sigma)}}\quad.$

Here $m(\sigma)=m([\sigma])$ denotes the number of cycles in the permutation
$\sigma$. Using Eq. (33) we obtain the following linear equation for the
moments $M(\vec{c};r,D)$:

$M(\vec{c};r,D)\ =\ \sum_{\begin{subarray}{c}\sigma\in S_{n}\\\ \tau\in
S_{n}^{\prime}\end{subarray}}\Biggl{(}Wg(\tau\sigma,rD)\;r^{m(\sigma)}D^{m([\vec{c}]\sigma)}\\\
\times\sum_{\begin{subarray}{c}\vec{k}\prec[\tau]\\\ k_{\ell}\neq
1\end{subarray}}M(\vec{k};r,D){[\tau]\choose\vec{k}}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{rD}\right)^{n-|\vec{k}|}D^{\\#\\{k_{\ell}=0\\}}\prod_{\ell\
:\ k_{\ell}=0}\left(1-(1-r)^{1-c_{\ell}^{\prime}(\tau)}\right)\Biggr{)}\quad.$
(34)

Here $M(\vec{k};r,s)$ denotes $M(\vec{k}^{\prime};r,s)$ where
$\vec{k}^{\prime}$ consists of the non-zero elements of $\vec{k}$ listed in
decreasing order. In principle, Eq. (34) allows us to solve for all moments
$M(\vec{c};r,s)$. Note in particular the triangular structure to the
equations: to solve for the moments with $\vec{c}\in CS_{n}$ we only need to
know the moments of lower order. Thus the equations can be solved recursively
for increasing $n$.

To analyze the solution to Eq. (34) in the limits $r\rightarrow\infty$ and
$D\rightarrow\infty$ , we need the Weingarten function asymptotics (see [10]).
Given $\alpha\in S_{n}$, with $[\alpha]=(c_{1}c_{2}\ldots c_{m})$,

$\operatorname{Wg}(\alpha,L)\ =\
C_{\alpha}L^{-n-|\alpha|}(-1)^{|\alpha|}+O(L^{-n-|\alpha|-2})\quad,$ (35)

where $|\alpha|=n-m$, which is the minimal length of $\alpha$ as a product of
$2$-cycles, and

$C_{\alpha}=\prod_{i=1}^{m}\frac{(2c_{i})!}{c_{i}!(c_{i}+1)!}\quad.$

Note that $d(\alpha,\beta)=|\alpha\beta^{-1}|$ is a metric on $S_{n}$. In
particular, we have the triangle inequality

$|\alpha\gamma^{-1}|\ \leq\ |\alpha\beta^{-1}|+|\beta\gamma^{-1}|\quad.$ (36)

Furthermore $|\alpha|$ is invariant on conjugacy classes; in particular,
$|\alpha|=|\alpha^{-1}|$.

####  III.1 Gaussianity: limit of large $r$

We now consider the asymptotics of Eq. (34) as $r\rightarrow\infty$. This is
technically somewhat simpler than the large $D$ limit. Let us denote by
$A(r)\sim B(r)$ the asymptotic relation
$\lim_{r\rightarrow\infty}\nicefrac{{A(r)}}{{B(r)}}=1$. By Eq. (35)

$\operatorname{Wg}(\tau\sigma,rD)\;r^{m(\sigma)}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{sD}\right)^{n-|\vec{k}|}\
\sim\
(-1)^{|\tau\sigma|}C_{\tau\sigma}D^{-n-|\tau\sigma|}\frac{r^{\frac{n-|\vec{k}|}{2}}}{r^{|\sigma|+|\tau\sigma|}}\quad.$

By the triangle inequality Eq. (36), $|\sigma|+|\tau\sigma|\geq|\tau|$.
Furthermore, for any $\tau\in S_{n}^{\prime}$ we have $|\tau|\geq n/2$.
Therefore unless $|\vec{k}|=0$ and $|\sigma|+|\sigma\tau|=|\tau|=n/2$, it
follows that

$\operatorname{Wg}(\tau\sigma,rD)\;r^{m(\sigma)}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{Dr}\right)^{n-|\vec{k}|}\
\rightarrow\ 0\quad\text{as }r\rightarrow\infty\quad.$

The identity $|\sigma|+|\sigma\tau|=n/2$ holds if and only if $\tau$ is a
pairing (i.e., a product of $n/2$ disjoint $2$-cycles), and if
$\sigma\prec\tau$ in the sense that $\sigma$ is a product of disjoint
$2$-cycles contained in $\tau$. In this case we have $C_{\tau\sigma}=1$ and

$\operatorname{Wg}(\tau\sigma,rD)\;r^{m(\sigma)}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{Dr}\right)^{n-|\vec{k}|}\
\rightarrow\ (-1)^{|\tau\sigma|}D^{-n-|\tau\sigma|}\quad.$

Since all terms on the right hand side of Eq. (34) tend to zero except for
when $\vec{k}=0$ and $\tau$ is a pairing, this equation reduces to

$M(\vec{c};r\rightarrow\infty,D)\ =\ \sum_{\begin{subarray}{c}\tau\in P_{n}\\\
\sigma\prec\tau\end{subarray}}(-1)^{|\tau\sigma|}D^{\frac{n}{2}-|[\vec{c}]\sigma|-|\tau\sigma|}\quad,$

where $P_{n}$ denotes the set of pairings of $\\{1,\ldots,n\\}$. It follows
that all moments are order one as $r\rightarrow\infty$ and given by the above
expression. In particular, $M(\vec{c};r\rightarrow\infty,D)$ vanishes if
$\vec{c}\in CS_{n}$ for $n$ odd.

Furthermore, we see that $w(\vec{c};r\rightarrow\infty,D)=0$ unless $n$ is
even and

$\vec{c}=(\underbrace{2,\cdots,2}_{j},\underbrace{1,\cdots,1}_{n-2j})\equiv(2^{j},1^{n-2j})$

for some $j$, for which we have

$w((2^{j},1^{n-2j});r\rightarrow\infty,D)\ =\
(-1)^{\frac{n}{2}-j}\frac{(n-2j)!}{2^{\frac{n}{2}-j}\left(\frac{n}{2}-j\right)!}\quad,$

where $\frac{(n-2j)!}{2^{\frac{n}{2}-j}\left(\frac{n}{2}-j\right)!}$ is the
number of pairings of $\\{1,\ldots,n-2j\\}$. Thus

$\lim_{r\rightarrow\infty}R^{(n)}=\begin{cases}0&\text{ if }n\text{ is
odd,}\\\
\displaystyle{\sum_{j=0}^{\nicefrac{{n}}{{2}}}}(-1)^{\frac{n}{2}-j}\frac{(n-2j)!}{2^{\frac{n}{2}-j}\left(\frac{n}{2}-j\right)!}\frac{1}{D^{n-j}}E_{D}^{(2^{j},1^{n-2j})}&\text{
if }n\text{ is even.}\end{cases}$ (37)

This result reduces to Eq. (24) for $n=1$ and $n=2$. For higher $n$ the result
it gives is the same as Wick’s theorem for
$\lim_{r}{\mathbb{E}}\\{W_{j}^{\otimes n}\\}$, showing that these matrices
have a Gaussian distribution in the large $r$ limit. For example,

$\lim_{r\rightarrow\infty}R^{(4)}\ =\
\frac{1}{D^{2}}E_{D}^{(2,2)}-\frac{1}{D^{3}}E_{D}^{(2,1,1)}+\frac{3}{D^{4}}E^{(1,1,1,1)}.$

In more detail, note that Wick’s theorem would imply that
$\lim_{r}{\mathbb{E}}\\{W_{j}^{\otimes 2n}\\}$ is a sum of terms obtained by
pairing the matrices in the tensor product and averaging according to the
variance in Eq. (24), namely
$R^{(2)}=\frac{1}{D}E^{(2)}_{D}-\frac{1}{D^{2}}E^{(1,1)}_{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
${\mathbb{E}}\\{W_{j}^{\otimes 2n}\\}$ as a sum over permutations with only
two cycles and one-cycles. In this expansion, each permutation comes with the
weight

$\\#\text{ pairings of the set covered by one cycles}\ \times\
(-1)^{\frac{1}{2}\\#\text{ one cycles}}\frac{1}{D^{\\#\text{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. (37).

####  III.2 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\rightarrow\infty$. Note that

$M((2);r,D)\ =\ \operatorname{tr}[P^{[(2)]}_{D}R^{(2)}]\ =\
D-\frac{1}{D}\quad.$

Thus we expect moments to become large, and a further normalization is needed
to compute the limit of Eq. (34). Indeed, we have

$\operatorname{Wg}(\tau\sigma,rD)\;D^{m([\vec{c}]\sigma)}D^{\\#\\{k_{\ell}=0\\}}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{Dr}\right)^{n-|\vec{k}|}\
\sim\
(-1)^{|\tau\sigma|}C_{\tau\sigma}\;\frac{r^{\frac{n-|\vec{k}|}{2}}}{r^{|\sigma|+|\tau\sigma|}}D^{\\#\\{k_{\ell}=0\\}-|[\vec{c}]\sigma|-|\tau\sigma|}\quad.$

By the triangle inequality,
$|[\vec{c}]\sigma|+|\tau\sigma|\geq|[\vec{c}]\tau^{-1}|\geq 0$, with equality
precisely if $[\vec{c}]=\tau$ and $\sigma\prec\tau$. Everything is maximized
for $\vec{k}=0$ and $\tau=[\vec{c}]$ suggesting that

$M((2^{j});r,D)\ \sim\ \text{const.}\ D^{m(\vec{c})}\quad.$

To proceed, we define

$\widetilde{M}(\vec{c};r,D)=D^{-m(\vec{c})}M(\vec{c};r,D)\quad,$

so $\widetilde{M}((2);r,D\rightarrow\infty)=1.$ Working from Eq. (34) we find
that

$\widetilde{M}(\vec{c};r,D)\ =\ \sum_{\begin{subarray}{c}\sigma\in S_{n}\\\
\tau\in
S_{n}^{\prime}\end{subarray}}\Biggl{(}Wg(\tau\sigma,rs)\;r^{m(\sigma)}D^{m([\vec{c}]\sigma)-m(\vec{c})+m(\tau)}\\\
\times\sum_{\begin{subarray}{c}\vec{k}\prec[\tau]\\\ k_{\ell}\neq
1\end{subarray}}\widetilde{M}(\vec{k};r,D){[\tau]\choose\vec{k}}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{Dr}\right)^{n-|\vec{k}|}\prod_{\ell\
:\ k_{\ell}=0}\left(1-(1-r)^{1-c_{\ell}^{\prime}(\tau)}\right)\Biggr{)}\quad.$

The large $D$ asymptotics of the coefficients are now governed by

$\operatorname{Wg}(\tau\sigma,rD)\;r^{m(\sigma)}D^{m([\vec{c}]\sigma)-m(\vec{c})+m(\tau)}\
\sim\
(-1)^{|\tau\sigma|}C_{\tau\sigma}\frac{r^{\frac{n-|\vec{k}|}{2}}}{r^{|\sigma|+|\tau\sigma|}}D^{|[\vec{c}]|-|\tau|-|[\vec{c}]\sigma|-|\tau\sigma|}\quad.$

Again by the triangle inequality,
$|\tau|+|[\vec{c}]\sigma|+|\tau\sigma|\geq|\tau|+|[\vec{c}]\tau^{-1}|\geq|[\vec{c}]|$,
so the limiting equation for $\widetilde{M}(\vec{c};r,s)$ has order one
coefficients as $D\rightarrow\infty$. Thus the resulting solution is of order
one.

Let us now look at Eq. (30) for $R^{(n)}$, with the knowledge that
$\widetilde{M}(\vec{c};r,s)$ are of order one. Rewriting Eq. (33) for the
weights $w(\vec{c};r,D)$ in terms of $\widetilde{M}$, we find that

$w(\vec{c};r,D)\ =\ \sum_{\tau\in
S_{n}^{\prime}}\Biggl{(}\operatorname{Wg}(\tau[\vec{c}],rD)(rD)^{m(\vec{c})}D^{m(\tau)}\\\
\times\sum_{\begin{subarray}{c}\vec{k}\prec[\tau]\\\ k_{\ell}\neq
1\end{subarray}}\widetilde{M}(\vec{k};r,D){[\tau]\choose\vec{k}}\left(\frac{\sqrt{1+rD^{2}}(r-1)}{rD}\right)^{n-|\vec{k}|}\prod_{\ell\
:\ k_{\ell}=0}\left(1-(1-r)^{1-c_{\ell}^{\prime}(\tau)}\right)\Biggr{)}\quad.$

By Eq. (35), the first factor of each term in the summation is asymptotically
given by

$\operatorname{Wg}(\tau[\vec{c}],rD)(rD)^{m(\vec{c})}D^{m(\tau)}\ \sim\
(-1)^{|\tau[\vec{c}]|}C_{\tau[\vec{c}]}\frac{r^{\frac{n-|\vec{k}|}{2}}}{r^{|[\vec{c}]|+|\tau\sigma|}}D^{n-|\tau|-|[\vec{c}]|-|\tau[\vec{c}]|}\quad.$

Now $|\tau|+|[\vec{c}]|+|\tau[\vec{c}]|\geq 2|\tau|\geq n$, with equality if
and only if $\tau$ is a pairing and $[\vec{c}]\prec\tau$. Thus
$w(\vec{c};r,D\rightarrow\infty)=0$ unless $n$ is even and
$\vec{c}=(2^{j},1^{n-2j})$ for some $j$; for such $\vec{c}$ we have

$w((2^{j},1^{n-2j});r,D\rightarrow\infty)=(-1)^{n-2j}p_{n-2j}\sum_{\begin{subarray}{c}\vec{k}\prec(2^{\frac{n}{2}})\\\
k_{\ell}\neq
1\end{subarray}}r^{-\frac{|\vec{k}|}{2}}\left(\frac{r-1}{r}\right)^{\\#\\{k_{\ell}=0\\}}\widetilde{M}(\vec{k};r,D\rightarrow\infty)\quad.$

Thus

$\widetilde{M}((2^{m});r,D\rightarrow\infty)\\\ \sim\
\sum_{j=0}^{m}\sum_{\sigma\in(2^{j}1^{2m-2j})}(-1)^{2m-2j}p_{2m-2j}\sum_{\begin{subarray}{c}\vec{k}\prec(2^{m})\\\
k_{\ell}\neq
1\end{subarray}}r^{-\frac{|\vec{k}|}{2}}\left(\frac{r-1}{r}\right)^{\\#\\{k_{\ell}=0\\}}\widetilde{M}(\vec{k};r,D\rightarrow\infty)\frac{1}{D^{3m-j}}\operatorname{tr}[P_{D}^{(2^{m})}P_{D}^{\sigma}]\;.$

The trace $\operatorname{tr}[P_{D}^{(2^{m})}P_{D}^{\sigma}]$ is bounded by
$D^{m+j}$, with equality precisely if $\sigma\prec[(2^{m})]$. As a result the
only terms that survive in the limit correspond to $j=m$ and
$\sigma=[(2^{m})]$. Thus,

$\widetilde{M}((2^{m});r,D\rightarrow\infty)\ =\
\sum_{\begin{subarray}{c}\vec{k}\prec(2^{m})\\\ k_{\ell}\neq
1\end{subarray}}r^{-\frac{|\vec{k}|}{2}}\left(\frac{r-1}{r}\right)^{\\#\\{k_{\ell}=0\\}}\widetilde{M}(\vec{k};r,D\rightarrow\infty)\quad.$

Introducing $q(m,r)=\widetilde{M}((2^{m});r,D\rightarrow\infty)$, this
equation reduces to

$q(m,r)\ =\ \sum_{j=0}^{m}{m\choose
j}\left(\frac{r-1}{r}\right)^{m-j}r^{-j}q(j,r)\ =\
\frac{1}{r^{m}}\sum_{j=0}^{m}{m\choose j}\left(r-1\right)^{m-j}q(j,r)\quad.$

We have $q(0,r)=1$ (by definition). It follows that $q(m,r)=1$ for all $m$ by
induction. Thus

$w((2^{j},1^{n-2j});r,D\rightarrow\infty)\ =\
(-1)^{n-2j}\frac{(n-2j)!}{2^{\frac{n}{2}-j}\left(\frac{n}{2}-j\right)!}\sum_{\begin{subarray}{c}\vec{k}\prec(2^{\frac{n}{2}})\\\
k_{\ell}\neq
1\end{subarray}}r^{-\frac{|\vec{k}|}{2}}\left(\frac{r-1}{r}\right)^{\\#\\{k_{\ell}=0\\}}\
=\
(-1)^{n-2j}\frac{(n-2j)!}{2^{\frac{n}{2}-j}\left(\frac{n}{2}-j\right)!}\quad,$

which proves that the limit is Gaussian in this case as well.

###  IV Bipartite entanglement entropy of an ergodic Matrix Product States

In this section we prove the formula for the bipartite entanglement entropy
for ergodic MPS given by Eq. (17) in the Letter. By translation invariance in
distribution it suffices to consider $j=0$. To focus on the bipartite
entanglement entropy across a single cut, we will make use of open boundary
conditions for the chain.

######  Remark IV.1.

The thermodynamic limit of the chain with open boundary conditions is
identical to that with periodic boundary conditions; this follows from the
convergence given in Theorem 2 of the Letter.

Consider an MPS defined over $[-M,N]$, with $M$, $N>0$ and open boundary
conditions. For this state, Eq. (11) of the Letter becomes:

$\left|\psi(-M,N)\right>\ =\
\sum_{\vec{i}=1}^{d}\,\vec{b}_{-M}^{i_{-M};\dagger}A_{-M+1}^{i_{-M+1}}\cdots
A_{N-1}^{i_{N-1}}\vec{b}_{N}^{i_{N}}\,\left|i_{-M},\ldots,i_{N}\right>\quad,$
(38)

where $\vec{b}_{-M}^{i}$ and $\vec{b}_{N}^{i}$, $i=1,\ldots,d$, are non-zero
vectors in ${\mathbb{C}}^{d}$ specifying the boundary conditions at each end
of the chain.

We consider the thermodynamic limit in two steps: first $N\rightarrow\infty$
and then $M\rightarrow\infty$. For an observable $O$ on the spins in $[-M,0]$,
it follows from Eq. (38) that the expectation of $O$ is given by the following
modification of Eq. (12) of the Letter:

$\left<\psi(-M,N)\right|O\left|\psi(-M,N)\right>=\
\frac{\operatorname{tr}[B_{N}\phi_{N-1}\circ\cdots\phi_{1}(\widehat{O}_{-M})]}{\operatorname{tr}[B_{N}\phi_{N-1}\cdots\circ\cdots\phi_{-M+1}(B_{-M})]}$
(39)

where

$B_{N}\ =\ \sum_{i=1}^{d}\vec{b}_{N}^{i}\vec{b}_{N}^{i;\dagger}\quad,\quad
B_{-M}\ =\ \sum_{i=1}^{d}\vec{b}_{-M}^{i}\vec{b}_{-M}^{i;\dagger}\quad,$

and

$\widehat{O}_{-M}\ =\
\sum_{\vec{i},\vec{j}=1}^{d}\left[\left<i_{-M},\ldots,i_{0}\right|O\left|j_{-M},\ldots,j_{0}\right>A_{0}^{i_{0}\dagger}\cdots
A_{-M+1}^{i_{-M+1}\dagger}\vec{b}_{-M}^{i_{-M}}\vec{b}_{-M}^{j_{-M}\dagger}A_{-M+1}^{j_{-M+1}}\cdots
A_{0}^{j_{0}}\right]\quad.$

Note that in the present case $\widehat{O}_{-M}$ is a $D\times D$ matrix.

It follows from Eq. (39) and Theorem 2 of the Letter that

$\lim_{N\rightarrow\infty}\left<\psi(-M,N)\right|O\left|\psi(-M,N)\right>=\
\frac{\operatorname{tr}[Z_{1}^{\prime}\widehat{O}_{-M}]}{\operatorname{tr}[Z_{1}^{\prime}\phi_{0}\circ\cdots\phi_{-M+1}(B_{-M})]}\quad.$
(40)

Let $\rho_{-M}$ denote the reduced density matrix of $[-M,0]$, in the limit
$N\rightarrow\infty$. By Eq. (40), we have

$\left<j_{-M},\ldots,j_{0}\right|\rho_{-M}\left|i_{-M},\ldots,i_{0}\right>\ =\
\frac{\operatorname{tr}[Z_{1}^{\prime}A_{0}^{i_{0}\dagger}\cdots
A_{-M+1}^{i_{-M+1}\dagger}\>\vec{b}_{-M}^{i_{-M}}\;\vec{b}_{-M}^{j_{-M}\dagger}\>A_{-M+1}^{j_{-M+1}}\cdots
A_{0}^{j_{0}}]}{\operatorname{tr}[Z_{1}^{\prime}\phi_{0}\circ\cdots\phi_{-M+1}(B_{-M})]}\quad.$

Let

$Q_{-M}\ \equiv\
\sqrt{Z_{1}^{\prime}}\;\\{\phi_{0}\circ\cdots\circ\phi_{-M+1}(B_{-M})\\}\sqrt{Z_{1}^{\prime}}\quad.$

Note that $Q_{-M}$ is a positive matrix. Let $\vec{v}_{1},\ldots,\vec{v}_{D}$
be its orthonormal eigenvectors, with corresponding eigenvalues
$\lambda_{1}\geq\cdots\geq\lambda_{D}\geq 0$. Let

$\left|\alpha_{j}\right>\ =\
\frac{1}{\sqrt{\lambda_{j}}}\sum_{\vec{i}=1}^{d}\vec{b}_{-M}^{i_{-M}\dagger}A_{-M+1}^{i_{-M+1}}\cdots
A_{0}^{i_{0}}\sqrt{Z_{1}^{\prime}}\vec{v}_{j}.$

Then

$\rho_{-M}\ =\
\frac{\sum_{j=1}^{D}\lambda_{j}\left|\alpha_{j}\right>\left<\alpha_{j}\right|}{\operatorname{tr}[Q_{-M}]}\
=\
\frac{\sum_{j=1}^{D}\lambda_{j}\left|\alpha_{j}\right>\left<\alpha_{j}\right|}{\sum_{j=1}^{D}\lambda_{j}},$

and

$\left<\alpha_{j}\middle|\alpha_{k}\right>\ =\
\frac{1}{\sqrt{\lambda_{j}\lambda_{k}}}\vec{v}_{j}^{\dagger}Q_{-M}\vec{v}_{k}\
=\ \delta_{j}^{k}\quad.$

Thus $\left|\alpha\right>_{j}$ are eigenvectors of $\rho_{-M}$ and the non-
zero eigenvalues of $\rho_{-M}$ are given by
$\frac{\lambda_{i}}{\sum_{j}\lambda_{j}}$.

It follows that the bipartite entropy across the bond $0\sim 1$ (in the limit
$N\rightarrow\infty$ with $M<\infty$) is given by

$S_{-M}(0)\ =\
-\operatorname{tr}[\widetilde{Q}_{-M}\log\widetilde{Q}_{-M}]\quad,$

with $\widetilde{Q}_{-M}=\frac{1}{\operatorname{tr}[Q_{-M}]}Q_{-M}$. Eq. (17)
of the Letter follows since

$\lim_{M\rightarrow\infty}\widetilde{Q}_{-M}\ =\
\frac{1}{\operatorname{tr}\left[\sqrt{Z_{1}^{\prime}}\;Z_{0}\;\sqrt{Z_{1}^{\prime}}\right]}\sqrt{Z_{1}^{\prime}}\;Z_{0}\;\sqrt{Z_{1}^{\prime}}\quad.$

by Theorem (2) of the Letter.