# Dynamics of Pseudoentanglement

Xiaozhou Feng  Matteo Ippoliti  Department of Physics, The University of Texas
at Austin, Austin, TX 78712, USA

###### Abstract

The dynamics of quantum entanglement plays a central role in explaining the
emergence of thermal equilibrium in isolated many-body systems. However,
entanglement is notoriously hard to measure, and can in fact be “forged”:
recent works have introduced a notion of pseudoentanglement describing
ensembles of many-body states that, while only weakly entangled, cannot be
efficiently distinguished from states with much higher entanglement, such as
random states in the Hilbert space. In this work we initiate the study of the
dynamical generation and propagation of pseudoentanglement. As generic quantum
dynamics tends to maximize actual entanglement, we consider constrained models
of time evolution: automaton (i.e. reversible classical) circuits that, when
fed suitable input states, provably produce the “standard models” of
pseudoentangled ensembles—uniformly random subset-phase states and subset
states—at late times, a phenomenon we name pseudothermalization. We examine
(i) how a pseudoentangled (random subset-phase) ensemble on a small subsystem
spreads to the whole system as a function of time, and (ii) how a
pseudoentangled (random subset) ensemble is generated from an initial product
state. We map the above problems onto a family of classical Markov chains on
subsets of the computational basis. The mixing times of such Markov chains are
related to the time scales at which the states produced from the dynamics
become indistinguishable from Haar-random states at the level of each
statistical moment (or number of copies). Based on a combination of rigorous
bounds and conjectures supported by numerics, we argue that each Markov
chain’s relaxation time and mixing time have different asymptotic behavior in
the limit of large system size. This is a necessary condition for a cutoff
phenomenon: an abrupt dynamical transition to equilibrium. We thus conjecture
that our random circuits give rise to asymptotically sharp
pseudothermalization transitions.

###### Contents

  1. I Introduction.
  2. II Review of key concepts
    1. II.1 Pseudoentanglement
    2. II.2 Equilibration of Markov chains and cut-off phenomenon
    3. II.3 Random walks on permutation groups
  3. III Pseudothermalization in automaton circuits
  4. IV Spreading of Pseudoentanglement
    1. IV.1 Motivation and setup
    2. IV.2 All-to-all circuit model
    3. IV.3 Local circuit model
      1. IV.3.1 First moment
      2. IV.3.2 Higher moments: Relaxation time
      3. IV.3.3 Higher moments: Mixing time
      4. IV.3.4 Higher dimension
  5. V Generation of Pseudoentanglement
  6. VI Discussion
    1. VI.1 Summary
    2. VI.2 Outlook
  7. A Proof of Proposition 1
    1. A.1 Preliminaries
    2. A.2 Proof
  8. B Properties of the $\Phi$ maps
  9. C Proof of Proposition 2
  10. D Random walk on the hypercube
    1. D.1 Single particle
    2. D.2 Many particles: multipole eigenmodes
    3. D.3 Two particles: relative coordinate eigenmodes
  11. E Properties of the local automaton gate set
  12. F Upper bound on relaxation time with local gate set
  13. G Upper bound on pseudothermalization time for subset states

##  I Introduction.

Figure 1: Main ideas of this work. (a) Quantum thermalization: chaotic quantum
dynamics on an isolated many-body systems (modeled by a Haar-random unitary
circuit, bottom) produces random states at late times, with near-maximal
entanglement entropy described by the Page curve (top). (b) In analogy with
thermalization, we consider families of random ‘automaton’ circuits (bottom,
see text) that, acting on suitable ensembles of initial states, produce
pseudoentangled state ensembles at late times. These have relatively low
entanglement entropy (top), but are computationally indistinguishable from
Page states. (c) We study the dynamical process by which pseudoentanglement is
generated and/or propagated over time, quantified by the distinguishability of
the output state ensembles from the Haar ensemble: can an observer, by making
joint measurements on $m$ copies of an output state, efficiently tell it apart
from a random state? For each value of $m$, this defines a time scale dubbed
($m$-copy) pseudothermalization time (vertical dashed lines). This time scale
may become sharp in the thermodynamic limit if the problem features a cutoff
phenomenon.

Entanglement is a fundamental concept in quantum information science and in
quantum many-body physics, being both a necessary resource for quantum
advantage in information processing Horodecki and Horodecki (2009); Chitambar
and Gour (2019) and a framework to characterize different phases of matter
beyond conventional order parameters Kitaev and Preskill (2006); Li and
Haldane (2008); Jiang  _et al._ (2012). Furthermore, the generation and
propagation of entanglement in many-body systems is a key ingredient in our
understanding of how thermal equilibrium emerges from the otherwise reversible
unitary dynamics of isolated systems Rigol  _et al._ (2008); Nandkishore and
Huse (2015); Kaufman  _et al._ (2016); Mezei and van der Schee (2020). Due to
entanglement, local subsystems may be described by mixed states even when the
global state is pure. Such local mixed states can approach universal forms at
late times in generic dynamics—namely the Gibbs ensemble $\rho_{A,{\rm
th}}\propto{\rm Tr}_{\bar{A}}(e^{-\beta H})$ Alhambra (2022) and its
generalizations Vidmar and Rigol (2016)—thus recovering predictions of
thermodynamics. In the absence of any conservation laws (e.g. in driven
systems without symmetries) the predicted late-time equilibrium state for each
local subsystem is maximally mixed: $\rho_{A,{\rm th}}\propto\mathbb{I}_{A}$.
This requires each local subsystem to become maximally entangled with the
remainder of the system. More specifically, in systems modeled by random
unitary circuits Nahum  _et al._ (2017); Bertini  _et al._ (2019); Zhou and
Nahum (2020); Brandao  _et al._ (2021); Fisher  _et al._ (2022); Potter and
Vasseur (2022); Mittal and Hunter-Jones (2023) (and even in some quasi-
periodically driven systems Pilatowsky-Cameo  _et al._ (2023, 2024)) one can
show that the global pure state of the system resembles a typical Haar-random
state in the Hilbert space, also known as a Page state, with a characteristic
scaling of entanglement known as the Page curve Page (1993), sketched in Fig.
1(a).

Despite its ubiquity and far-reaching importance, entanglement is generally
difficult to control, quantify, and measure Horodecki and Horodecki (2009),
especially beyond few-body settings Amico  _et al._ (2008). The last decade
has seen substantial progress in these directions, particularly thanks to
experimental advances in quantum simulators and theoretical developments in
randomized measurements Islam  _et al._ (2015); Kaufman  _et al._ (2016);
Elben  _et al._ (2018); Brydges  _et al._ (2019); Huang  _et al._ (2020);
Elben  _et al._ (2023); nonetheless, the cost of measuring the entanglement
entropy of an $N$-qubit system remains exponential in $N$, ultimately limiting
the reach of all these approaches.

Recent works Aaronson  _et al._ (2023); Giurgica-Tiron and Bouland (2023)
highlight the intrinsic subtlety of measuring entanglement in a dramatic way,
by introducing a notion of pseudoentanglement—a forgery of entanglement that
cannot be efficiently distinguished from the real thing. More specifically,
they construct ensembles of quantum states that, while only weakly entangled,
are indistinguishable from Page states to any observer with a polynomial
amount of resources. These ensembles are interesting from the point of view of
computation and cryptography Ji  _et al._ (2018); Lu  _et al._ (2023), and can
be efficiently constructed through the use of quantum-secure pseudorandom
permutations Aaronson  _et al._ (2023); Zhandry (2016).

Motivated by the role of entanglement in quantum thermalization, in this work
we consider the dynamical generation and propagation of pseudoentanglement in
many-body systems. Clearly, to avert the generic tendency of interacting many-
body systems to maximize entanglement, we must specialize to non-generic
models of time evolution; however, within these constraints (to be discussed
below), we introduce random circuit models that prepare pseudoentangled state
ensembles as their equilibrium steady states, much like the Page states for
generic quantum dynamics. This is sketched in Fig. 1(b). We name this
dynamical approach to pseudoentangled ensembles pseudothermalization.
Depending on one’s perspective, this may be viewed as either a spoofing of
quantum thermalization that is practically indistinguishable from the genuine
phenomenon, or as a perfectly valid instance of quantum thermalization based
on an operational definition (i.e., if a deviation from thermal equilibrium
cannot be efficiently detected, one may argue it does not count).

Specifically, we consider the random subset-phase state ensemble Aaronson  _et
al._ (2023) and the random subset state ensemble Giurgica-Tiron and Bouland
(2023), which are fixed points of ‘automaton’ quantum circuits Gopalakrishnan
(2018); Iaconis  _et al._ (2019); Iaconis (2021); Klobas  _et al._ (2021)
(i.e., permutations of the computational basis, not to be confused with
‘quantum cellular automata’). We therefore introduce families of automaton
circuits and study the convergence of certain non-equilibrium initial states
to these equilibrium fixed points, i.e., the phenomenon of
pseudothermalization. We focus on two scenarios modeled after paradigmatic
problems in ‘true’ quantum thermalization: (i) the spreading of
pseudoentanglement from a small subsystem placed in contact with a larger,
disentangled system, and (ii) the generation of pseudoentanglement from an
initial product state.

The key technical step in both cases is a mapping of the quantum
pseudothermalization problem onto the equilibration of certain classical
Markov chains over subsets of the computational basis—with the size of subsets
related to the number of state copies $m$ available to the observer. This
mapping allows us to make contact with a rich body of mathematical work on the
equilibration of Markov chains on finite spaces, including the so-called
cutoff phenomenon Diaconis and Shahshahani (1981); Diaconis (1996); Chen and
Saloff-Coste (2008); Lacoin and Leblond (2011); Chen (2021). This phenomenon
refers to the sudden equilibration that can be seen in random walks over
highly-connected graphs (for example permutations of a deck of cards induced
by a riffle shuffle Aldous and Diaconis (1986)), in contrast with e.g.
diffusion on a Euclidean lattice where equilibration is gradual and
characterized by a dynamical scaling collapse $\sim t/L^{z}$, $z=2$.

Based on this connection, we argue that pseudothermalization in our random
circuit models happens in an asymptotically sharp fashion: our output states
are easily distinguishable from Page states for a long time, until they
abruptly (i.e., over a parametrically shorter window of time) become
indistinguishable at the level of $m$ copies, Fig. 1(c). This dynamical
transition in the distinguishability of low-entanglement states from Page
states is a new fundamental phenomenon in the dynamics of isolated quantum
many-body systems, that enriches and complements our growing understanding of
thermal equilibrium at the level of higher statistical moments Roberts and
Yoshida (2017); Hunter-Jones (2019); Cotler  _et al._ (2023); Ippoliti and Ho
(2022); Pilatowsky-Cameo  _et al._ (2023); Fava  _et al._ (2023) and opens a
number of exciting directions for future work.

The rest of the paper is organized as follows. In Sec. II we provide a brief,
self-contained review of some useful background: pseudoentanglement (Sec.
II.1), Markov chains on finite sets and the cutoff phenomenon (Sec. II.2), and
certain random walks over permutation groups (Sec. II.3). We then introduce
the general framework of pseudothermalization via automaton circuit dynamics
in Sec. III, and apply it to the problem of spreading and generation of
pseudoentanglement in Sec. IV and V respectively. Finally we summarize our
results and list outstanding questions and future directions in Sec. VI.

##  II Review of key concepts

Asymptotic notation

In this work we make extensive use of the following asymptotic notation. Given
two functions $f,g:\mathbb{N}\mapsto\mathbb{R}$, we write:

  * •

$f=o(g)$ if $\lim_{N\to\infty}f(N)/g(N)=0$;

  * •

$f=O(g)$ if $f(N)<Cg(N)$ for all $N$ for some constant $C>0$;

  * •

$f=\Theta(g)$ if $f=O(g)$ and $g=O(f)$;

  * •

$f=\Omega(g)$ if $g=O(f)$;

  * •

$f=\omega(g)$ if $\lim_{N\to\infty}f(N)/g(N)=\infty$.

With a slight abuse of notation, we also write asymptotic upper and lower
bounds as e.g. $\Omega(f)\leq g\leq O(h)$ to denote that $g=O(h)$ and
$g=\Omega(f)$. We also use a tilde to denote when one of these conditions
applies up to a ${\rm polylog}(N)$ factor: e.g. we may write $f\leq
O(g(N)/\log(N))$ as $f\leq\tilde{O}(g)$.

###  II.1 Pseudoentanglement

An ensemble of quantum states on $N$ qubits is said to be pseudoentangled if
it cannot be distinguished from the unitarily invariant (Haar) random ensemble
by looking at any ${\rm poly}(N)$ number of copies Aaronson  _et al._ (2023).
More precisely, defining the $m$-th moment operator

$\rho^{(m)}_{\mathcal{E}}=\mathbb{E}_{\psi\sim\mathcal{E}}[\outerproduct{\psi}{\psi}^{\otimes
m}],$ (1)

an ensemble $\mathcal{E}$ is pseudoentangled111The definition in Ref. Aaronson
_et al._ (2023) additionally requires the ability to prepare the ensemble
efficiently through suitable pseudorandom primitives. We will not focus on
that aspect in this work. if

$\left\|\rho_{\mathcal{E}}^{(m)}-\rho_{{\rm Haar}}^{(m)}\right\|_{\rm
tr}=o\left(1/{\rm poly}(N)\right)$ (2)

whenever $m=O({\rm poly}(N))$. The trace norm $\|A\|_{\rm tr}={\rm
Tr}(\sqrt{A^{\dagger}A})$ characterizes the distinguishability between two
quantum states: if $\|\rho-\sigma\|_{\rm tr}=\epsilon$, then there is an
observable222 One can write the observable as $\mathcal{O}=\Pi_{+}-\Pi_{-}$,
where $\Pi_{\pm}$ are the projectors on positive (negative) eigenvalue
subspaces of $\rho-\sigma$.  $\mathcal{O}$, with $-I\leq\mathcal{O}\leq I$,
that gives ${\rm Tr}(\mathcal{O}\rho)-{\rm Tr}(\mathcal{O}\sigma)=\epsilon$;
then an observer can tell apart $\rho$ from $\sigma$ by measuring
$\mathcal{O}$ a number of times $\sim 1/\epsilon^{2}$. This is the optimal way
to distinguish two states Nielsen and Chuang (2010). Thus the meaning of Eq.
(2) is that, when making joint measurements on $m$ copies of states drawn from
$\mathcal{E}$, an observer would need a superpolynomial number of experiments
to tell the two ensembles apart. Since we require Eq. (2) to apply for all
$m\leq O({\rm poly}(N))$, the states from $\mathcal{E}$ necessarily appear to
be Haar-random (and thus near-maximally entangled) to any observer with
polynomial resources, even though they may in fact be much less entangled.

A simple pseudoentangled ensemble is given by random subset phase states
Aaronson  _et al._ (2023)

$\ket{\psi_{S,f}}=\frac{1}{|S|^{1/2}}\sum_{\mathbf{z}\in
S}(-1)^{f(\mathbf{z})}\ket{\mathbf{z}},$ (3)

where $f$ is a random Boolean function on $\\{0,1\\}^{N}$ chosen uniformly at
random, and $S$ is a subset of bitstrings, $S\subset\\{0,1\\}^{N}$, chosen
uniformly at random out of all subsets of some fixed cardinality $|S|=K$. We
call this space of subsets $\Sigma_{K}$. Here $K$ must scale suitably with the
system size $N$: $\omega({\rm poly}(N))\leq K\leq o(2^{N})$. The role of $K$
is to determine the actual entanglement of the states in the ensemble: we have
the average von Neumann entropy

$S_{\rm vN}(A)\sim\min\\{|A|\log 2,\log K\\}$ (4)

for all subsystems $A$ (regardless of locality), due to the fact that $K$
upper bounds the Schmidt rank of $\ket{\psi_{S,f}}$ on any bipartition of the
$N$ qubits. Thus the restriction on $K\geq\omega({\rm poly}(N))$ ensures that,
e.g., the half-cut entanglement entropy can scale faster than $\log N$
(necessary so that the states cannot be efficiently distinguished from Haar-
random ones by measuring their purity), and the restriction $K\leq o(2^{N})$
ensures that it scales more slowly than the Page value $\simeq N\log(2)/2$
Page (1993) (see Fig. 1(a,b)).

Very recently, Ref. Giurgica-Tiron and Bouland (2023) showed that the random
phase factors are not necessary to have a pseudoentangled ensemble. Random
subset states suffice:

$\ket{\psi_{S}}=\frac{1}{|S|^{1/2}}\sum_{\mathbf{z}\in S}\ket{\mathbf{z}},$
(5)

with $S$ drawn uniformly at random from $\Sigma_{K}$ (subsets of the
computational basis of cardinality $K$). This ensemble obeys

$\|\rho_{\mathcal{E}}^{(m)}-\rho_{{\rm Haar}}^{(m)}\|_{{\rm tr}}\leq
O\left(\frac{m^{2}}{K}\right)+O\left(\frac{mK}{D}\right),$ (6)

with $D=2^{N}$ the Hilbert space dimension. This is superpolynomailly small if
$\omega({\rm poly}(N))\leq K\leq o(2^{N}/{\rm poly}(N))$ and $m=O({\rm
poly}(N))$.

###  II.2 Equilibration of Markov chains and cut-off phenomenon

Figure 2: Schematic of the cutoff phenomenon for Markov chains on finite sets.
(a) A family of Markov chains on systems of size $N$ has a cutoff if it
exhibits “sudden” equilibration as $N\to\infty$, in the following sense: the
total variation distance from equilibrium remains close to maximal for a time
scale $t_{\rm mix}$ (mixing time), then decays to zero exponentially with a
time constant $t_{\rm rel}$ (relaxation time) with $\lim_{N\to\infty}t_{\rm
rel}/t_{\rm mix}=0$. (b) Even though the relaxation time itself may be growing
with $N$, upon rescaling time by $t\mapsto t/t_{\rm mix}$ the approach to
equilibrium becomes step-like as $N\to\infty$. In this sense there is a
dynamically sharp equilibration transition.

A key technical step in this work will be to map the evolution of quantum
distinguishability measures such as Eq. (2) to the equilibration of classical
Markov chains over finite sets. For this purpose it is useful to briefly
review some notations and facts about this topic, including the cut-off
phenomenon which will play an important role in our results.

A Markov chain on a finite set can be simply defined by a stochastic matrix
$\Gamma_{ij}$, with $i,j$ ranging over the set. ‘Stochastic’ means that
$\Gamma_{ij}\geq 0$, with $\sum_{i}\Gamma_{ij}=1$ for all $j$. $\Gamma_{ij}$
represents the probability of a transition from $j$ to $i$. A state of the
Markov chain is given by a probability distribution $p(i)$. In each time step,
$p$ is updated according to $p(i)\mapsto\sum_{j}\Gamma_{ij}p(j)$. We say that
the Markov chain equilibrates if $\Gamma^{t}p$ converges to some late-time
limit $\pi$ as $t\to\infty$.

The equilibration of finite-sized Markov chains is characterized by two
fundamental time scales Levin and Peres (2017). The first is the mixing time,

$t_{\rm mix}(\epsilon)=\min\\{t:\|p_{t}-\pi\|_{\rm tv}\leq\epsilon\\},$ (7)

describing the approach to the steady distribution. Here $\|\cdots\|_{\rm tv}$
is the total variation norm, $\|p\|_{\rm tv}=\frac{1}{2}\sum_{i}|p(i)|$, and
one should maximize over all possible initial states. The second important
time scale is the relaxation time, defined as the inverse spectral gap of the
transition matrix $\Gamma_{ij}$:

$t_{\rm rel}=\frac{1}{1-\lambda_{1}},$ (8)

where the spectrum of $\Gamma_{ij}$ is given by
$1=\lambda_{0}>\lambda_{1}\geq\lambda_{2}\geq\dots>-1$. The relaxation time
describes the asymptotic late-time decay towards the steady state.

In many cases of interest in physics, these two time scales have the same
behavior: e.g., for the simple random walk on a line of length $L$, both time
scales are $\propto L^{2}$. This gives a scale-invariant behavior, with the
approach to equilibrium depending on the single ratio $t/L^{2}$. However, in
highly-connected spaces (e.g. expander graphs Hoory  _et al._ (2006)) the two
time scales can be different, with $t_{\rm rel}/t_{\rm mix}\to 0$ in the
large-system limit. Such a separation of time scales, sketched in Fig. 2, is
conjectured Peres (2004); Chen and Saloff-Coste (2008) to give rise to the
cut-off phenomenon Diaconis and Shahshahani (1981); Diaconis (1996); Chen and
Saloff-Coste (2008); Lacoin and Leblond (2011); Chen (2021): an abrupt
transition to equilibrium, defined by

$\lim_{N\to\infty}\|p_{t(N)}-\pi\|_{\rm tv}=\chi(\theta),\quad t(N)=t_{\rm
mix}(N)+\theta t_{\rm rel}(N).$ (9)

Here $\chi(\theta)$ is an $N$-independent function that goes to 1 as
$\theta\to-\infty$ and to 0 as $\theta\to+\infty$. Informally, the total
variation distance from equilibrium exhibits a scaling collapse as
$N\to\infty$ on the form $\chi[(t-t_{\rm mix})/t_{\rm rel}]$, where the
crucial observation is that the mid-point of the collapse ($t_{\rm mix}$) and
its width ($t_{\rm rel}$) scale differently in the large-system limit: $t_{\rm
rel}=o(t_{\rm mix})$. In other words, if we express time in units of $t_{\rm
mix}$, the distance from equilibrium becomes step-like in the infinite-system
limit, going from almost-maximal to almost-zero over a vanishingly narrow
window of time, $1\pm t_{\rm rel}/t_{\rm mix}$. This is sketched in Fig. 2(b).

###  II.3 Random walks on permutation groups

In this work we will be especially interested in Markov chains over
permutation groups and related processes. Here we give a brief review of some
essential concepts Levin and Peres (2017).

Consider a graph $(V,E)$ with vertices $V=\\{v_{i}\\}$ and edges
$E=\\{(v_{i},v_{j})\\}$. The interchange process (IP) is a Markov chain over
the permutation group $\mathcal{S}_{|V|}$ whose unit step consists of drawing
an edge $(v_{i},v_{j})\in E$ uniformly at random and transposing its vertices
$v_{i}$, $v_{j}$. The composition of $t$ such transpositions $\tau_{i,j}$
produces a random element
$\sigma_{t}=\tau_{i_{t},j_{t}}\circ\cdots\circ\tau_{i_{1},j_{1}}$ of the
permutation group $\mathcal{S}_{|V|}$. The IP is the Markov chain
$\\{\sigma_{t}\\}_{t=0}^{\infty}$ on $\mathcal{S}_{|V|}$.

The IP induces other Markov chains of interest. The random walk (RW) with
initial state $v_{i}$ is the Markov chain
$\\{\sigma_{t}(v_{i})\\}_{t=0}^{\infty}$. This has the intuitive meaning we
expect for a random walk on the graph $(V,E)$: a particle initialized at
vertex $v_{i}$ can hop with equal probability to any of its neighboring
vertices. Additionally, the exclusion process with $m$ particles
($\textsf{EX}(m)$) is defined as the Markov chain
$\\{\sigma_{t}(A)\\}_{t=0}^{\infty}$ where $A\subset V$ is a collection of $m$
distinct vertices, and $\sigma_{t}(A)\equiv\\{\sigma_{t}(v_{i}):v_{i}\in A\\}$
is its image under permutation $\sigma_{t}$. This process describes the
evolution of $m$ indistinguishable, hard-core particles on the graph, where at
each time step the occupation number of two neighboring vertices get swapped.
When $m=1$, this clearly reduces to the random walk:
$\textsf{EX}(m)=\textsf{RW}$.

These three Markov chains are connected by a remarkable result known as
Aldous’ spectral gap conjecture333Despite the name, this is a proven
theoremCaputo  _et al._ (2010)., which states that

$t_{\rm rel}^{\textsf{IP}}=t_{\rm rel}^{\textsf{RW}}=t_{\rm
rel}^{\textsf{EX}(m)}\quad\forall\,m.$ (10)

The inequality $t_{\rm rel}^{\textsf{IP}}\geq t_{\rm rel}^{\textsf{EX}(m)}\geq
t_{\rm rel}^{\textsf{RW}}$ is a simple consequence of the fact that RW and
$\textsf{EX}(m)$ are sub-processes of IP Caputo  _et al._ (2010) (informally,
they are projections of IP onto a smaller state space), and thus it holds also
for other Markov chains on $\mathcal{S}_{|V|}$ beyond IP. However, the other
direction of the inequality is highly nontrivial and is specific to IP, where
it was conjectured in 1992 and proved in 2009 Caputo  _et al._ (2010).

Beyond IP, in this work we will consider other Markov chains over permutation
groups. In particular, we will consider a setting where the unit time step
does not just transpose two neighboring vertices, like IP, but rather performs
one of a more general set of “elementary permutations” in $\mathcal{S}_{|V|}$
sampled randomly at each time step. Such a process also gives a Markov chain
over $\mathcal{S}_{|V|}$ and also induces “generalized exclusion processes”
through the action of $\mathcal{S}_{|V|}$ on the subset spaces $\Sigma_{m}$,
in the same way as IP induces $\textsf{EX}(m)$. These “generalized exclusion
processes” differ from the standard one as they can allow multiple particles
to hop in a correlated way in each time step. Aldous’ spectral gap conjecture,
Eq. (10), does not apply to this more general case. In Sec. IV.3 we will in
fact encounter a case where the $m$-particle “generalized exclusion process”
relaxes more slowly than the single-particle walk.

##  III Pseudothermalization in automaton circuits

Figure 3: (a) Subset states and their transformation under automaton gates,
illustrated here for $N=3$ qubits and a CCX gate. (b) Spaces of subsets of the
computational basis, $\Sigma_{m}=\\{S\subset\mathbb{Z}_{2}^{N}:\,|S|=m\\}$,
shown for $N=3$ and $m=1,2,3$. Vertices represent subsets $S\in\Sigma_{m}$,
edges represent the action of an automaton gate set $\mathcal{G}$ (self-edges
are omitted for simplicity).

We aim to introduce models of quantum dynamics that achieve pseudoentanglement
as their equilibrium state. However, fully-generic quantum dynamics tends to
maximize actual entanglement. The easiest way to avoid this fate is to
constrain the dynamics to preserve the ‘subset(-phase) state’ property, Sec.
II.1. Reversible classical circuits, also known as automaton circuits in the
quantum context444Note that some works use a different definition of automaton
gates, which includes phases Iaconis (2021). Gopalakrishnan (2018); Iaconis
_et al._ (2019); Iaconis (2021); Klobas  _et al._ (2021), have this property:
an automaton circuit $U_{t}$ corresponds to a permutation
$\sigma_{t}\in\mathcal{S}_{D}$ of the computational basis,
$U_{t}\ket{\mathbf{z}}\equiv\ket{\sigma_{t}\mathbf{z}}$; thus, its action on a
subset-phase state $\ket{\psi_{S_{0},f_{0}}}$ is

$\displaystyle U_{t}\ket{\psi_{S_{0},f_{0}}}$
$\displaystyle=\frac{1}{\sqrt{K}}\sum_{\mathbf{z}\in
S_{0}}(-1)^{f_{0}(\mathbf{z})}\ket{\sigma_{t}(\mathbf{z})}=\ket{\psi_{S_{t},f_{t}}},$
(11)

where $S_{t}=\\{\sigma_{t}(\mathbf{z}):\,\mathbf{z}\in
S_{0}\\}\equiv\sigma_{t}(S_{0})$ is another subset of cardinality $K$, and
$f_{t}\equiv f_{0}\circ\sigma_{t}^{-1}$ is another phase function. In
particular, if $f_{0}$ is uniformly random then so is $f_{t}$, and if $f_{0}$
is constant then so is $f_{t}$. It follows that automaton dynamics preserves
the ‘subset(-phase) state’ property. Further, it preserves the cardinality of
the subsets, $K$. The action of automaton gates on subsets is sketched in Fig.
3(a). An immediate consequence of Eq. (11) is that, under arbitrary automaton
dynamics, the entanglement entropy about any cut cannot grow beyond $\log K$
(as $K$ upper-bounds the Schmidt rank of the state). The maximum amount of
entanglement is thus controlled by the choice of initial state, and saturation
to the Page curve can be avoided.

In this work we consider the following type of dynamics. A distribution of
initial states $\ket{\psi_{S_{0},f_{0}}}$ is specified (these may be subset
states or subset-phase states depending on the allowed values of $f_{0}$);
then, each step of the time evolution consists of sampling an automaton gate
$u_{i}$ uniformly at random from a fixed gate set $\mathcal{G}$ (we will
discuss requirements on $\mathcal{G}$ in the following) and applying it to the
state. After $t$ time steps, this results in an ensemble of time-evolved
states $\\{U_{t}\ket{\psi_{S_{0},f_{0}}}=\ket{\psi_{S_{t},f_{t}}}\\}$ defined
by all possible initial states and choices of
$U_{t}=u_{i_{t}}u_{i_{t-1}}\cdots u_{i_{1}}\in\mathcal{G}^{\circ t}$. The
relevant notion of equilibration for this state ensemble is whether, at late
times, it becomes uniformly distributed on the appropriate (subset or subset-
phase) state space. If that is the case, and subject to usual constraints on
the subset size $K$ (see Sec. II.1), then our dynamics achieves
pseudoentanglement at late times. We call this condition pseudothermalization:
while the upper limit $K$ on the Schmidt rank prevents full thermalization
(defined, for a model without conservation laws such as this one, by the
approach to a Page state), this failure of thermalization may be practically
undetectable, in a sense that is rigorously formalized through the concept of
pseudoentanglement.

A first question to address is whether such equilibrium is achieved at all for
a given automaton gate set $\mathcal{G}$. Let us focus on the subsets and
neglect the phases for now. The action of the quantum circuits
$U_{t}\ket{\psi_{S_{0}}}=\ket{\psi_{S_{t}}}$ defines a Markov chain
$\\{S_{t}\\}$ over the space of size-$K$ subsets of the computational basis,

$\Sigma_{K}\equiv\\{S\subset\\{0,1\\}^{N}:|S|=K\\}.$ (12)

The transition matrix of the Markov chain, $\Gamma_{S^{\prime},S}$, is
specified by the gate set $\mathcal{G}$:

$\Gamma_{S,S^{\prime}}=\frac{1}{|\mathcal{G}|}\sum_{u\in\mathcal{G}}\delta_{S,u(S^{\prime})}.$
(13)

Viewing $\Gamma$ as an adjacency matrix, we can thus picture the $\Sigma_{K}$
spaces as graphs; some minimal examples are shown in Fig. 3(b). To achieve
pseudothermalization, the Markov chain $\\{S_{t}\\}$ must converge at late
times to the uniform distribution on $\Sigma_{K}$, denoted here as $\pi_{K}$:

$\lim_{t\to\infty}p_{t}=\pi_{K},$ (14)

where $p_{t}(S)$ is the probability distribution over subsets after $t$ time
steps. To guarantee the uniform distribution $\pi_{K}$ as a unique steady
state, the Markov chain must be Levin and Peres (2017):

  1. (i)

Reversible: the Markov chain transition matrix is symmetric,
$\Gamma_{S^{\prime},S}=\Gamma_{S,S^{\prime}}$ for all $S,S^{\prime}$. This
property is ensured e.g. by picking a gate set $\mathcal{G}$ such that
$u^{2}=I$ for all $u\in\mathcal{G}$.

  2. (ii)

Irreducible: for any two states $S,S^{\prime}\in\Sigma_{K}$, there exists a
sequence of gates in $\mathcal{G}$ that maps $S$ to $S^{\prime}$.
Irreducibility will be discussed separately for each gate set $\mathcal{G}$ in
the following.

  3. (iii)

Aperiodic: the chain can connect any two states $S,S^{\prime}$ with paths of
any (sufficiently large) length555The prototypical counterexample is given by
walks on bipartite graphs, where all paths connecting $S$ to $S^{\prime}$ have
the same length modulo 2. In that case the populations of the two partitions
can keep oscillating indefinitely.. One way to ensure aperiodicity is to make
the chain ‘idle’, i.e. allow it to skip a step. This can be done by including
the identity gate in $\mathcal{G}$.

For a gate set $\mathcal{G}$ with all three properties above, we can conclude
that $p_{t}\to\pi_{K}$ as $t\to\infty$. This suffices to establish
pseudoentanglement at late times in two cases:

  * •

when the initial state distribution is restricted to a constant phase
function, e.g. $f_{0}=0$, we obtain the random subset ensemble Giurgica-Tiron
and Bouland (2023);

  * •

when the initial state distribution includes a uniformly random phase function
$f_{0}$, we obtain the random subset-phase ensemble Aaronson  _et al._ (2023).

Both of these ensembles are pseudoentangled, subject to certain constraints on
$K$, as we reviewed in Sec. II.1.

Having established that the steady-state ensembles in these two cases are
pseudoentangled, it remains to understand how the system approaches this
condition over the course of its dynamics, and what are the time scales
involved. As already mentioned, we will call this process
pseudothermalization, in analogy with how conventional thermalization in
quantum many-body systems can be defined from the saturation of entanglement
of local subsystems. Two key technical results will allow us to make progress
on these questions by connecting pseudothermalization with the equilibration
of certain classical Markov chains.

######  Proposition 1 (subset-phase states).

Consider an ensemble of subset-phase states
$\mathcal{E}_{p}=\\{\ket{\psi_{S,f}}:\ S\sim p,\ f\sim{\rm unif.}\\}$,
specified by a uniformly random distribution over phase functions and a
general probability distribution $p$ over $\Sigma_{K}$ (subsets). We have, for
all integers $1\leq m\leq K$,

$\left\|\rho^{(m)}_{\mathcal{E}_{p}}-\rho^{(m)}_{\rm Haar}\right\|_{\rm
tr}=2\|\Phi_{m\leftarrow K}[p]-\pi_{m}\|_{\rm
tv}+O\left(\frac{m^{2}}{K}\right)$ (15)

Here $\|\cdots\|_{\rm tr}$ and $\|\cdots\|_{\rm tv}$ are the trace norm and
total variation norm, respectively; $\rho^{(m)}$ denotes the $m$-th
statistical moment of a state ensemble, i.e.
$\rho^{(m)}_{\mathcal{E}_{p}}=\mathbb{E}_{S\sim
p}\mathbb{E}_{f}[\outerproduct{\psi_{S,f}}{\psi_{S,f}}^{\otimes m}]$ and
$\rho^{(m)}_{\rm Haar}=\int{\rm d}\psi_{\rm
Haar}\outerproduct{\psi}{\psi}^{\otimes m}$; and $\Phi_{m\leftarrow K}$ is a
stochastic map from distributions on $\Sigma_{K}$ to distributions on
$\Sigma_{m}$, defined by

$\Phi_{m\leftarrow
K}[p](S^{\prime})=\binom{K}{m}^{-1}\sum_{S\in\Sigma_{K}:\,S^{\prime}\subset
S}p(S).$ (16)

###### Proof.

See Appendix A. ∎

If we take $m={\rm poly}(N)$ and $K=\omega({\rm poly}(N))$, then the error
term on the RHS of Eq. (15) is negligible, and pseudoentanglement (defined by
the LHS being small) is tied to the non-uniformity of the distribution
$\Phi_{m\leftarrow K}[p]$ on $\Sigma_{m}$. If we view $p$ as the distribution
$\Sigma_{K}$ induced by the action of a random automaton circuit, then
$\Phi_{m\leftarrow K}[p]$ is the distribution induced by the same process on
$\Sigma_{m}$ (up to an average over initial states). We show this in Appendix
B, along with other useful properties of the stochastic maps $\Phi$. Then Eq.
(15) maps pseudoentanglement onto the mixing of a classical Markov chain on
$\Sigma_{m}$.

######  Proposition 2 (subset states).

Consider an ensemble of subset states
$\mathcal{E}^{\prime}_{p}=\\{\ket{\psi_{S}}:\ S\sim p\\}$, specified by a
general probability distribution $p$ over $\Sigma_{K}$ (subsets). We have, for
all integers $1\leq m\leq K/2$,

$\left\|\rho_{\mathcal{E}^{\prime}_{p}}^{(m)}-\rho_{\rm
Haar}^{(m)}\right\|_{\rm tr}=\left\|\mathcal{M}^{(m)}_{p}\right\|_{\rm
tr}+O\left(\frac{m}{\sqrt{K}}\right).$ (17)

Here $\mathcal{M}^{(m)}_{p}$ is a $|\Sigma_{m}|\times|\Sigma_{m}|$ matrix
defined by

$\displaystyle(\mathcal{M}^{(m)}_{p})_{S^{\prime},S^{\prime\prime}}$
$\displaystyle=\binom{K}{m}^{-1}\binom{K}{m+\delta}f^{(m+\delta)}_{p}(S^{\prime}\cup
S^{\prime\prime}),$ (18) $\displaystyle f^{(m+\delta)}_{p}$
$\displaystyle=\Phi_{m+\delta\leftarrow K}[p]-\pi_{m+\delta},$ (19)

where we use subsets $S^{\prime},S^{\prime\prime}\in\Sigma_{m}$ to index the
matrix elements, and define $\delta\equiv|S^{\prime}\setminus
S^{\prime\prime}|$.

###### Proof.

See Appendix C. ∎

We again relate pseudoentanglement (the LHS of Eq. (17) being small) to the
non-uniformity of certain probability distributions over subsets; except
unlike the previous case, Eq. (15), here we have contributions from a family
of distributions $\\{\Phi_{m+\delta\leftarrow K}[p]\\}_{\delta=0}^{m}$, rather
than from the single distribution $\Phi_{m\leftarrow K}[p]$. If we view $p$ as
arising from a Markov chain on $\Sigma_{K}$ induced by a random automaton
circuit dynamics, then the emergence of pseudoentanglement is again
related—albeit in a less straightforward way—to the mixing of the associated
Markov chains $\Phi_{m+\delta\leftarrow K}[p]$ on $\Sigma_{m+\delta}$, for all
$\delta=0,1,\dots,m$.

Having reduced the problem of pseudothermalization to that of equilibration of
certain Markov chains (namely random walks over the subset spaces
$\\{\Sigma_{m}\\}$, Fig. 3(b)), we are now in a position to study the problem
quantitatively by leveraging known results on these processes, including those
reviewed in Sec. II.2 and II.3. In the following sections we consider several
automaton gate sets $\mathcal{G}$ and initial state distributions, and analyze
the associated classical Markov chains to characterize the process of
pseudothermalization.

##  IV Spreading of Pseudoentanglement

###  IV.1 Motivation and setup

In this section, we study the problem of how pseudoentanglement spreads in
space, focusing on the setting of a smaller pseudoentangled system placed in
contact with a larger, trivial system. This problem is reminiscent of
analogous settings in the dynamics of entanglement: e.g., the growth of
entanglement after quantum quenches where thermalization may happen due to the
spreading of entangled quasiparticles in space Calabrese and Cardy (2005,
2009), or the equilibration of quantum spin chains connected to thermal baths
at their edges Prosen and Žnidarič (2009); Znidaric  _et al._ (2011);
Morningstar  _et al._ (2022); Sels (2022).

We consider a situation in which we are given a random subset-phase ensemble
$\mathcal{E}_{0}=\\{\ket{\psi_{S_{0},f_{0}}}:\ S_{0}\sim{\rm unif.},\
f_{0}\sim{\rm unif.}\\}$ on a system $A$ of $N_{A}$ qubits. We take
$|S_{0}|=K$ with $\omega({\rm poly}(N_{A}))<K<o(2^{N_{A}})$, so that the
initial ensemble on $A$ is pseudoentangled. We then introduce a system
$\bar{A}$ of $N-N_{A}$ ancillas initialized in the product state
$\ket{0}^{\otimes N-N_{A}}$, and form a larger system $A\bar{A}$ of $N$
qubits. This overall system is still in a subset-phase state, though not a
uniformly-random one; the initial state indeed is of the form
$\ket{\psi_{S_{0},f_{0}}}\otimes\ket{0}^{\otimes N-N_{A}}$. we can
straigthforwardly redefine $S_{0}$ as a subset of the computational basis of
the extended system $\\{0,1\\}^{N}$ (rather than the original system
$\\{0,1\\}^{N_{A}}$), by padding each bitstring with $N-N_{A}$ trailing zeros:
$\mathbf{z}=(z_{1},\dots z_{N_{A}})\mapsto(z_{1},\dots z_{N_{A}},0,\dots 0)$.
The phase function $f_{0}$ is never evaluated on any bitstrings outside the
original ones, so it can be trivially extended to a uniformly random binary
function on $\\{0,1\\}^{N}$. Thus after including the ancillas, we have a
subset-phase ensemble with uniformly random phases, but highly-correlated
subsets. These states are manifestly not pseudoentangled—one can efficiently
distinguish them from Haar-random states by, e.g., measuring $\langle
Z_{j}\rangle$ for any $j>N_{A}$.

We then place the original system and ancillas into contact by allowing them
to evolve together under automaton circuit dynamics, as discussed in Sec. III.
The question we aim to address is if and how pseudoentanglement is achieved on
the entire system of $N$ qubits over the course of the evolution. As the
phases are already randomized, it is sufficient to focus on the subsets.

A first apparent constraint is that the value of $K$ is conserved, even as we
add the ancillas: even assuming that the subsets indeed equilibrate to the
uniform distribution, pseudoentanglement requires $\omega({\rm
poly}(N))<K<o(2^{N})$ (see Sec. II.1). Combining these constraints with the
analogous ones from the initial state ($N\mapsto N_{A}$) yields the overall
constraint666Some examples of parameters satisfying these constraints are as
follows: $N_{A}=N/2$, $K=\exp(N^{1/2})$; $N_{A}=N^{1/2}$, $K=\exp(N^{1/4})$;
$N_{A}=(\log N)^{2}$, $K=N^{\log\log N}$. These examples are purely
illustrative; in the following we will not assume anything about $N$, $N_{A}$,
$K$ other than Eq. (20).

$\omega({\rm poly}(N))<K<o(2^{N_{A}}).$ (20)

In particular this implies that ${N_{A}}$ must grow faster than $\log(N)$.
Intuitively this means that the number of ancillas added to the original
system cannot be too large, otherwise one would “dilute” the initial
pseudoentanglement too much.

Assuming Eq. (20), it remains to analyze the equilibration of subsets toward
the uniform distribution under the chosen model of automaton dynamics. In the
rest this section, we discuss how to achieve this under two automaton gate
sets $\mathcal{G}$, one made of all-to-all gates and the other of local gates.

###  IV.2 All-to-all circuit model

Figure 4: All-to-all automaton circuit model for pseudoentanglement spreading.
(a) Sketch of the circuit, showing a subset-phase initial state
$\ket{\psi_{S_{0},f_{0}}}$ on subsystem $A$ along with $\ket{0}^{\otimes
N-N_{A}}$ on the rest of the system, and a sequence of random globally-
controlled NOT gates, $\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}$. Note that the
controls come in two variations, open and solid circles, representing control
on 0 and 1 respectively. (b) Sketch of the exact result Eq. (25) for the
pseudothermalization transitions in this model. Consecutive
pseudothermalization times $t_{\rm PTherm}^{(m)}$, $t_{\rm PTherm}^{(m+1)}$
are not sharply separated.

We start by analyzing a model with a high degree of mathematical tractability,
that comes at the cost of all-to-all interactions and exponentially-long
equilibration time scales. We emphasize that these latter aspects are not at
all necessary to pseudothermalization in general, and are in fact absent from
another model that we study next, in Sec. IV.3.

We consider a gate set made of globally-controlled NOT gates:
$\mathcal{G}=\\{\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}\\}$, where
$\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}$ flips qubit $i$ if and only if the
remaining $N-1$ qubits are in computational basis state $\ket{\mathbf{a}}$.
There are $ND/2$ gates in $\mathcal{G}$: $N$ choices for the target $i$ and
$2^{N-1}=D/2$ choices for the control string $\mathbf{a}$. The resulting
circuit is illustrated in Fig. 4(a).

Viewed as a permutation of the hypercube $\mathbb{Z}_{2}^{N}$, the gate
$\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}$ is simply a transposition of two
neighboring vertices,

$(\mathbf{a}_{1:i-1},0,\mathbf{a}_{i+1:N})\xleftrightarrow{\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}}(\mathbf{a}_{1:i-1},1,\mathbf{a}_{i+1:N})$
(21)

(we use the shorthand $\mathbf{a}_{i:j}=(a_{i},\dots a_{j})$). Clearly by a
sequence of such neighbor transpositions one can transpose two arbitrary
vertices777Given a path from $\mathbf{x}$ to $\mathbf{y}$ on the hypercube,
the sequence of neighbor transpositions along the path and back has the effect
of transposing $\mathbf{x}$ and $\mathbf{y}$ while leaving all other sites
fixed. $\mathbf{x}\leftrightarrow\mathbf{y}$; with these one can generate the
whole permutation group $\mathcal{S}_{D}$. It follows that $\mathcal{G}$
induces an irreducible Markov chain on $\Sigma_{m}$ for all values of $m$. The
Markov chains induced by $\mathcal{G}$ on the subset spaces $\\{\Sigma_{m}\\}$
are also clearly reversible ($\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}$ is its
own inverse) and aperiodic (it is always possible to draw a value of
$\mathbf{a},i$ such that $\textsf{C}^{(N-1)}_{\mathbf{a}}{X}_{i}$ acts
trivially on all bitstrings in $S$, making the walk idle). By the discussion
in Sec. III we have that the unique steady state of the Markov chain on
$\Sigma_{K}$ is the uniform distribution $\pi_{K}$. Thus at late times this
dynamics produces the random subset-phase state ensemble over the entire
system, which in the relevant regime of $(K,N)$ is pseudoentangled.

Having addressed the $t\to\infty$ limit of the dynamics, it remains to
understand how pseudoentanglement is achieved as a function of time. To this
end, we can use Proposition 1, which relates the $m$-copy distinguishability
between the time-evolved state ensemble $\mathcal{E}_{t}$ and the Haar
ensemble to the mixing of a Markov chain $\Phi_{m\leftarrow K}[p_{t}]$ on
$\Sigma_{m}$.

We note that the Markov chain on the permutation group $\mathcal{S}_{D}$
induced by this gate set $\mathcal{G}$ is precisely the interchange process
(IP) on the hypercube $\mathbb{Z}_{2}^{N}$, reviewed in Sec. II.3: each gate
picks an edge of the hypercube uniformly at random and transposes its two
vertices. It follows that the Markov chains $\Phi_{m\leftarrow K}[p_{t}]$ over
the spaces of subsets $\Sigma_{m}$ are exactly the exclusion processes
$\textsf{EX}(m)$. These are well-studied problems for which we can rely on
powerful known results.

First of all, Aldous’ spectral gap conjecture Caputo  _et al._ (2010), Eq.
(10), states that the relaxation time of $\textsf{EX}(m)$ is equal, for all
$m$, to the relaxation time of the simple random walk (RW). The simple random
walk on the hypercube can be solved straightforwardly (see Appendix D.1), and
one finds888 The exponential scaling in $N$ is due to the fact that each
elementary transposition of hypercube vertices is exponentially unlikely to
cause the random walker to hop: only $N$ gates out of $ND/2$ cause a
nontrivial hop. The conventional result for a random walker who hops once per
unit time is $t_{\rm rel}^{\textsf{RW}}=N/2$; rescaling this by the average
number of gates between hops ($D/2$) yields Eq. (22) (see also more detailed
derivation in Appendix D.1).  $t_{\rm rel}^{\textsf{RW}}=ND/4$; therefore

$t_{\rm rel}^{\textsf{EX}(m)}=t_{\rm
rel}^{\textsf{RW}}=\frac{ND}{4}\quad\forall\ m.$ (22)

Secondly, the mixing time for $\textsf{EX}(m)$ was recently solved Hermon and
Pymar (2020) and was found to scale as

$t_{\rm mix}^{\textsf{EX}(m)}\sim ND\log(Nm).$ (23)

Due to Proposition 1, this mixing time is also the time taken to achieve
$m$-copy pseudothermalization.

Based on Eq. (22) and (23), we see that

$t_{\rm mix}^{\textsf{EX}(m)}/t_{\rm
rel}^{\textsf{EX}(m)}\sim\log(Nm)\xrightarrow{N\to\infty}\infty.$ (24)

This asymptotic separation between the scaling of relaxation and mixing times
is a necessary condition for a cut-off phenomenon, reviewed in Sec. II.2. It
was conjectured that in sufficiently “generic” models, such separation should
also be a sufficient condition999 Note that the condition $t_{\rm
rel}=o(t_{\rm mix})$ is not sufficient to get a cutoff phenomenon in total
variation distance in general (on the other hand it is sufficient if the
distance from equilibrium is measured with a $L^{p}$ norm, $p>1$ Chen and
Saloff-Coste (2008)). It is possible to have a so-called pre-cutoff where
multiple step-like transitions arise in the total variation distance. However,
this usually requires some particular structure in the connectivity of the
underlying graph (e.g. bottlenecks separating highly connected regions). .
This cutoff effect would imply a sequence of dynamical transitions in our
dynamics: for any moment $m\leq{\rm poly}(N)$, there is a time scale at which
the ensemble of output states becomes indistinguishable from Haar-random at
the $m$-copy level; this $m$-copy pseudothermalization time in this model is
given by

$t_{\rm PTherm}^{(m)}\sim N2^{N}[\log(mN)\pm O(1)],$ (25)

with the $\pm$ term denoting the temporal width of the dynamical crossover to
equilibrium. This result is sketched in Fig. 4(b). We note that, while each
pseudothermalization transition is sharp (in the sense of $t_{\rm
rel}^{(m)}/t_{\rm mix}^{(m)}\to 0$ at large $N$), consecutive transitions
$t_{\rm PTherm}^{(m)}$, $t_{\rm PTherm}^{(m+1)}$ are not sharply separated.
This is due to the slow growth of $t_{\rm mix}^{(m)}$ with $m$, namely the
fact that $[t_{\rm mix}^{(m+1)}-t_{\rm mix}^{(m)}]/t_{\rm
rel}^{(m)}\sim\log\frac{m+1}{m}$ does not diverge (and in fact it goes to zero
at large $m$).

This exactly solvable model serves as a nice illustration of the idea of
pseudothermalization, including the possible emergence of a cutoff phenomenon.
At the same time, it has several drawbacks: the all-to-all interactions
require complex couplings or compilation into deep local circuits; even worse,
the resulting time scales for relaxation and mixing are exponentially long in
system size. However, these limitations are artifacts of the model, and not
intrinsic to the concept of pseudoentanglement or pseudothermalization. Below
we introduce a local automaton gate set that does not suffer from these
issues.

###  IV.3 Local circuit model

Figure 5: Local automaton circuit model for the spreading of
pseudoentanglement. (a) Circuit model, analogous to Fig. 4 but with local
$\textsf{CC}X$ gates. (b) Schematic summary of bounds on the
pseudothermalization time. Solid lines denote unconditional bounds while the
dot-dashed line depends on Conjecture 1. The yellow area indicates allowed
scalings. The lower bounds come from Eq. (30) (counting) and Eq. (31)
(locality). The dashed-line upper bound comes from Eq. (29); the solid-line
upper bound is Eq. (27). (c) Consecutive $m$-copy pseudothermalization
transitions in this model may be overlapping (left) or sharply separated
(right) in the thermodynamic limit: answering this question requires tighter
bounds on the mixing time.

In the interest of building a more physically realistic model of dynamics, we
seek to restrict the gate set $\mathcal{G}$ to geometrically local gates on a
$d$-dimensional array of qubits.

The simplest automaton gate set comprises just local bit flips:
$\mathcal{G}=\\{X_{i}\\}_{i=1}^{N}$. It is straightforward to see that the
action of $\mathcal{G}$ on $\Sigma_{1}$ is irreducible (we can map any
bitstring to any other via bitflips); however it is already not irreducible on
$\Sigma_{2}$: given a set $S=\\{\mathbf{z},\mathbf{z}^{\prime}\\}$, the
“relative coordinate” $\mathbf{z}\oplus\mathbf{z}^{\prime}$ is invariant under
bit flips. Moving up in complexity, one might try two-qubit automaton gates,
generated by NOTs and CNOTs; however, it is possible to show that these gates
are also insufficient to achieve irreducibility on $\Sigma_{m}$ above $m=3$.
Three-qubit automaton gates, on the other hand, can induce irreducible and
aperiodic Markov chains on $\Sigma_{m}$ for all $m\leq D/4$, as we show in
Appendix E.

In particular, for the case of a one-dimensional qubit array, we consider the
gate set
$\mathcal{G}=\\{u_{iab}=\textsf{C}_{i-1,a}\textsf{C}_{i+1,b}{X}_{i}\\}$ where
$a,b\in\\{0,1\\}$ label a control bitstring and $i\in[N]$ labels the target
site (we take periodic boundary conditions, identifying $i=0$ and $i=N$), as
shown in Fig. 5(a). The $u_{iab}$ gate flips qubit $i$ if and only if its left
neighbor $i-1$ is in state $a$ and its right neighbor $i+1$ is in state $b$.
There are $N$ choices for the target site $i$ and 4 choices for the control
string $(a,b)$, giving a total of $4N$ gates in $\mathcal{G}$. Up to bit
flips, each $u_{iab}$ is a Toffoli gate. It is easy to see that the gate set
generates $\textsf{C}X$ gates, SWAP gates, and bit flip gates as well—it in
fact generates arbitrary even permutations of the computational basiss, see
Appendix E.

####  IV.3.1 First moment

We start our analysis of this dynamics from the first moment, $m=1$. The
Markov chain on $\Sigma_{1}\equiv\mathbb{Z}_{2}^{N}$ corresponds once again to
the simple random walk on the hypercube, see Appendix D.1. With each gate the
walker, a bitstring $\mathbf{z}$, stays in place with probability 3/4 (if the
control string $a,b$ of the randomly-sampled $u_{iab}$ gate does not match
$\mathbf{z}$), otherwise it hops to a neighbor (a bitstring
$\mathbf{z}^{\prime}$ differing only at one bit) chosen uniformly at random
based on the value of $i\in[N]$. Using again the known results for the simple
random walk on the hypercube, Appendix D.1, we have

$t_{\rm rel}^{(1)}=2N,\qquad t_{\rm mix}^{(1)}\sim N\log(N),$ (26)

with the only difference from the case of Sec. IV.2 (with $m=1$) being the
different idling time between hops (4 instead of $D/2$). Just like in Sec.
IV.2, the random walk exhibits a cutoff Diaconis  _et al._ (1990), with the
approach to equilibrium taking place in a relatively narrow time window,
$t\sim N[\log N\pm O(1)]$. Intuitively, to achieve equilibrium the walker has
to hop in all $N$ directions with high probability (equivalently, all $N$
qubits need to be flipped with high probability). This is an instance of the
so-called “coupon collector’s problem”, which gives a scaling $\sim N\log(N)$.
Note that time $t$ is defined as total number of gates acting on the system;
defining a rescaled time $\tau\equiv t/N$ such that each qubit flips of order
once per unit $\tau$, equilibration happens at $\tau\sim\log(N)\pm O(1)$.

####  IV.3.2 Higher moments: Relaxation time

For higher moments $m>1$, unlike the all-to-all model of Sec. IV.2, the Markov
chain on $\Sigma_{m}$ is not exactly solvable (to our knowledge). To make
progress we resort to a combination of rigorous bounds and numerical
calculations.

Figure 6: (a) Largest 20 eigenvalues of the Markov chain transition matrix
$\Gamma_{S,S^{\prime}}$ for different system sizes $N=5,\dots 10$ (indicated
at the top of each panel) and subset sizes $m$, obtained from exact
diagonalization (Lanczos). We consider all pairs $(N,m)$ such that the
dimension $|\Sigma_{m}|=\binom{2^{N}}{m}$ of the transition matrix is below
$2\times 10^{7}$. For $m=1$ we recover the spectrum of the single-particle
random walk, $1-\frac{c}{2N}$, $c\in[N]$; a new leading eigenvalue
$\lambda_{1}>1-\frac{1}{2N}$ appears for $m=2$ and persists for all accessible
values of $m$. (b) Relaxation time $t_{\rm rel}=1/(1-\lambda_{1})$ as a
function of system size $N$ for two Markov chains: the full Markov chain over
pairs of bitstrings, $\Gamma$ ($\lambda_{1}$ taken from (a), $m=2$), and the
random walk of the “relative coordinate” $\Gamma|_{R}$ (see main text). The
two agree to numerical precision up to the system sizes accessible to both
($N\leq 10$); the data suggests an asymptotically linear scaling of $t_{\rm
rel}^{(2)}$ with $N$.

We begin by analyzing the relaxation time $t_{\rm rel}^{(m)}$ of the Markov
chain on $\Sigma_{m}$, for $m>1$, through exact numerical diagonalization of
the Markov chain transition matrix $\Gamma$ for accessible values of the
system size $N$ and subset size $m$. Exact diagonalization is limited by the
dimensionality of the problem, given by the size of the subset space101010The
effective dimensionality is lower due to the problem’s translation and
reflection invariance, but only by a factor of $\approx N$.
$|\Sigma_{m}|=\binom{D}{m}=\binom{2^{N}}{m}$. This grows explosively with $m$
and $N$, limiting us to $mN\lesssim 30$.

Results for the spectra of $\Gamma$ for different system sizes $N$ and moments
$m$ are shown in Fig. 6(a). Recall that the spectrum of $\Gamma$ has a leading
eigenvalue $\lambda_{0}=1$, associated to the uniform distribution $\pi_{m}$
(its unique steady state), and then subleading eigenvalues
$1>\lambda_{1}\geq\lambda_{2}\geq\dots$; the inverse spectral gap defines the
relaxation time: $t_{\rm rel}^{(m)}=1/(1-\lambda_{1})$. For $m=1$ the problem
reduces once again to the simple random walk on the hypercube
$\mathbb{Z}_{2}^{N}$, so the leading eigenvalue $\lambda_{0}=1$ is followed by
a $N$-fold degenerate multiplet of eigenvalues $\lambda_{1}=1-\frac{1}{2N}$,
giving the relaxation time $t_{\rm rel}^{(1)}=2N$. These eigenvalues appear in
the spectra for all values of $m$, where they are associated to the $N$
components of the dipole moment of the $m$-particle distribution on the
hypercube, see Appendix D.2.

However, for $m\geq 2$, we observe a new second-largest eigenvalue in the
spectrum of $\Gamma$, $1-1/2N<\lambda_{1}<1$, giving longer relaxation times
$t_{\rm rel}^{(m)}>t_{\rm rel}^{(1)}$. This is unlike the all-to-all model of
Sec. IV.2, where Aldous’ spectral gap conjecture Eq. (10) (see Sec. II.3)
forces $t_{\rm rel}^{(m)}$ to be $m$-independent. Nonetheless, Fig. 6(a) shows
no further $m$-dependence in $\lambda_{1}$ beyond $m=2$. This suggests that
$m=1$ is special in this problem.

By numerical inspection, we find that the new slowest mode of $\Gamma$ for
$m=2$, the function $f_{1}(S)$ such that $\Gamma f_{1}=\lambda_{1}f_{1}$,
depends on the subset $S=\\{\mathbf{z},\mathbf{z}^{\prime}\\}$ only through
the “relative coordinate”
$\mathbf{r}\equiv\mathbf{z}\oplus\mathbf{z}^{\prime}$. In Appendix D.3 we show
that the subspace $R$ spanned by functions that depend on $S$ only through the
relative coordinate $\mathbf{r}$ is invariant under the Markov chain. Thus we
can block-diagonalize $\Gamma=\Gamma|_{R}\oplus\Gamma|_{R^{\perp}}$, where
$\Gamma|_{R}$ describes a random walk of $\mathbf{r}$. Numerical
diagonalization confirms that the leading eigenvalue $\lambda_{1}$ of $\Gamma$
indeed comes from the $\Gamma|_{R}$ block: results for $\lambda_{1}$ agree to
numerical precision between the two Markov chains. Furthermore, due to the
reduced dimensionality ($\sim D$ for $\Gamma|_{R}$ compared to $\sim D^{2}$
for the full $\Gamma$), we are able to numerically diagonalize the restricted
Markov chain up to larger values of system size $N$. Results, shown in Fig.
6(b), indicate a clear linear scaling $t_{\rm rel}^{(2)}=\Theta(N)$.

These observations ($t_{\rm rel}^{(2)}$ scaling linearly in $N$ and $t_{\rm
rel}^{(m)}$ being independent of $m\geq 2$ for all numerically accessible
values) lead us to formulate the following conjecture:

######  Conjecture 1 (scaling of relaxation time).

For all $m\geq 2$, the Markov chain on $\Sigma_{m}$ induced by the local
automaton gate set $\mathcal{G}$ has relaxation time $t_{\rm
rel}^{(m)}=\Theta(N)$, independent of $m$.

We cannot rule out the possibility of further changes in the relaxation times
$t_{\rm rel}^{(m)}$ at larger values of $m$ beyond those accessible in
numerics. In Appendix F we show an unconditional upper bound

$t_{\rm rel}^{(m)}\leq O(mN^{5}),$ (27)

which we believe is very loose. Proving (or disproving) our Conjecture 1,
and/or tightening bounds around the relaxation times $t_{\rm rel}^{(m)}$, are
left as goals for future research. However, since a precise determination of
relaxation times in the arguably simpler models of Ref. Gowers (1996); Hoory
_et al._ (2005); Brodsky and Hoory (2008) is still lacking, we suspect these
may be hard problems.

####  IV.3.3 Higher moments: Mixing time

Figure 7: Numerical simulation of local random automaton circuits on a chain
of $N=64$ qubits, initialized in the product state $\ket{0}^{\otimes
N-N_{A}}\otimes\ket{+}^{\otimes N_{A}}$ with $N_{A}=4$ (we shift the $\ket{+}$
qubits to the center of the plot for clarity). Data averaged over $10^{4}$
random circuit realizations. Left: local expectation values averaged over
circuit realizations, $\overline{\langle Z_{i}\rangle}$, can relax under local
bit flips and thus equilibrate quickly. Right: higher-moment quantities, such
as $\overline{\langle Z_{i}\rangle^{2}}$ (related to the variance of $\langle
Z_{i}\rangle$ across circuit realizations), require the propagation of
coherence across the system. Their equilibration is thus constrained by a
light cone, $t/N\geq\Omega(N)$. Note $t$ is the total number of applied gates.

Assuming the validity of Conjecture 1, $t_{\rm rel}^{(m)}=\Theta(N)$, we move
on to the question of mixing times. A general two-sided bound on the mixing
time in terms of the relaxation time is given by Jerison (2013); Levin and
Peres (2017)

$(t^{(m)}_{\rm rel}-1)\log(\frac{1}{2\epsilon})\leq t_{\rm
mix}^{(m)}(\epsilon)\leq t_{\rm rel}^{(m)}\log(|\Sigma_{m}|/\epsilon),$ (28)

where $\epsilon$ is the accuracy threshold on total variation distance used to
define mixing; see Sec. II.2. We take $\epsilon=1/4$ in the following.
Plugging Conjecture 1 in Eq. (28) yields the asymptotic bounds

$\Omega(N)\leq t_{\rm mix}^{(m)}\leq O(mN^{2}),$ (29)

where we use $\log|\Sigma_{m}|=\log\binom{D}{m}\leq mN$.

The lower bound Eq. (28) clearly cannot be used to gap $t_{\rm mix}^{(m)}$
away from $t_{\rm rel}^{(m)}$, which is the necessary condition for a cut-off
phenomenon. To study the possible presence of a cut-off we must find stronger
lower bounds on the mixing time. We show two such bounds next.

######  Fact IV.1 (counting bound).

For all $m$, the mixing for the Markov chain induced by the local gate set
$\mathcal{G}$ obeys

$t_{\rm mix}^{(m)}\geq\Omega(mN/\log N).$ (30)

###### Proof.

Let us count the number of states that can be reached from an initial state
$S_{0}\in\Sigma_{m}$ in $t$ time steps. As each state is connected to at most
$|\mathcal{G}|=4N$ different states, the volume of a ball of radius $t$ in
$\Sigma_{m}$ (viewed here as a graph with edges given by the action of
$\mathcal{G}$) is at most $(4N)^{t}$. To achieve mixing, the distribution
$p_{t}$ needs to be supported on, say, at least $1/2$ of the states. Thus we
must have $(4N)^{t}\geq|\Sigma_{m}|/2$, giving
$t\log(4N)\geq\log\binom{D}{m}\geq mN$. ∎

######  Fact IV.2 (locality bound).

For all $m\geq 2$, the mixing for the Markov chain induced by the local gate
set $\mathcal{G}$ on a 1D qubit chain obeys

$t_{\rm mix}^{(m)}\geq\Omega(N^{2}).$ (31)

###### Proof.

It is easy to see that, starting from a subset $S_{0}$ where all bitstrings
agree outside the subsystem $A$ ($z_{j}=0$ for all sites $j=1,\dots N-N_{A}$,
for all bitstrings $\mathbf{z}\in S_{0}$), the region outside a “light cone”
emanating from $A$ remains in a product state: local automaton gates cannot
create superpositions, so coherence must spread outward from $A$ one gate at a
time. As the gate sequence is random, the light cone fluctuates instance by
instance, however its width is bounded above111111 We can bound $\ell(t)$
above by noting that $\ell(t+1)=\ell(t)+1$ with probabiliy at most $2/N$ (the
$\textsf{CC}X$ gate must have the target qubit $i$ at either light cone
front), and $\ell(t+1)\leq\ell(t)$ otherwise. Then $\overline{\delta\ell}\leq
2/N$ and thus $\ell(t)\leq\ell(0)+2t/N$ with finite probability. This gives an
upper bound on the butterfly velocity $v\leq 1$.  by $\ell(t)\leq
N_{A}+2vt/N$, for some constant ‘butterfly velocity’ $v$, with finite
probability. Whenever $\ell(t)<N$, the distribution $p_{t}(S)$ is supported
only on a small fraction of the state space $\Sigma_{m}$ (all subsets $S$ that
correspond to a product state outside the light cone), and thus must be far
from uniform121212 Another way to argue this is by the distinguishing
statistics method in Ref. Levin and Peres (2017) applied to, e.g.,
$\overline{\langle Z_{i}\rangle^{2}}$ for a qubit $i$ far from $A$. .
Equilibrium can only be reached after the light cone has spread across the
whole system, giving $\ell(t)\geq\Omega(N)$ or $t\geq\Omega(N^{2})$. ∎

We note that, while the counting bound (Eq 30) applies to all $m$, the
locality bound Eq. 31) does not apply to $m=1$ (indeed we have $t_{\rm
mix}^{(1)}\sim N\log N$, Eq. (26)). This is because “subsets” with $m=1$ are
just individual bitstrings, so equilibration does not require coherence (i.e.
superpositions) to spread across the system. This physical difference in the
effects of locality is illustrated in Fig. 7. We show results of numerical
simulations for circuit-averaged quantities $\overline{\langle Z_{i}\rangle}$
(a first-moment quantity that can equilibrate under local bitflips) in Fig.
7(a) and $\overline{\langle Z_{i}\rangle^{2}}$ (a second-moment quantity that
requires superpositions to equilibrate) in Fig. 7(b). Again we recall that
“time” $t$ here is measured as total number of gates in the circuit, so $t\sim
N^{2}$ here denotes ballistic, not diffusive, spreading of correlations.

Putting Facts IV.1 and IV.2 together with the upper bound Eq. (29), we arrive
at

$\Omega(N^{2},mN/\log(N))\leq t_{\rm mix}^{(m)}\leq O(mN^{2}).$ (32)

These bounds are summarized in Fig. 5(b). We remark again that the upper bound
Eq. (29) assumes Conjecture 1; an unconditional, likely loose bound $t_{\rm
mix}^{(m)}\leq O(m^{2}N^{6})$ can be derived based on Eq. (27) (see Appendix
F).

Several interesting consequences follow form Eq. (32). First of all,
independently of $m$, we have

$\lim_{N\to\infty}t_{\rm mix}^{(m)}/t_{\rm rel}^{(m)}=\infty,$ (33)

which is a necessary condition for a cut-off, see Sec. II.2. We thus
conjecture the presence of asymptotically sharp pseudothermalization
transitions in this model. While numerical simulations are strongly limited by
the explosive growth of the state space $|\Sigma_{m}|=\binom{2^{N}}{m}$,
results of exact simulations of $N=5$ qubits, shown in Fig. 8, are consistent
with the emergence of a cutoff—we clearly see the development of a plateau,
followed by a decay $\sim\lambda_{1}^{t-\delta t}$ of the total variation
distance, with a temporal offset $\delta t\sim m$ in the accessible
range131313$N=5$, $m=8$ yields $|\Sigma_{8}|=\binom{32}{8}\gtrsim 10^{7}$
states. $1\leq m\leq 8$.

Figure 8: Exact numerical simulation of the Markov chain induced by the local
automaton gate set on $N=5$ qubits for subset sizes $m=1,\dots 8$. The initial
quantum state is taken to be $\ket{00+++}$, a subset state with $K=8$. The
total variation distance from uniformity decays as $\|p_{t}-\pi_{m}\|_{\rm
tv}\sim\lambda^{t-\delta t}$ at late times (bottom panel, dashed lines). Fit
results for $\lambda$ (not shown) are consistent with the results of Fig. 6,
with $\lambda\simeq 0.90$ for $m=1$ and $\lambda\simeq 0.92$ for all $m\geq
2$; the temporal offset $\delta t$, a proxy for the mixing time (up to
additive constants), is found to scale linearly in $m$ (bottom panel, inset).

We can then break down the behavior of different moments $m$, based on their
scaling with $N$, to obtain an interesting (though still incomplete) picture
of the dynamics:

  1. (i)

For constant $m$, we have $t_{\rm mix}^{(m)}=\Theta(N^{2})$, and thus the
pseudothermalization time

$t_{\rm PTherm}^{(m)}\sim N^{2}[1\pm O(1/N)],$ (34)

with the $\pm$ term denoting the width of the time window in which the
transition to equilibrium happens. This suggests a picture based on a
‘pseudothermalization front’ spreading ballistically through the system, with
fast equilibration behind the front, so that $m$-copy pseudothermalization
happens suddenly, and simultaneously for all constant $m$, once the front has
swept across the system.

  2. (ii)

In the intermediate range $\omega(1)\leq m\leq O(N\log N)$, our bounds are
insufficient to pin down the scaling of $t_{\rm PTherm}^{(m)}$; we have

$t_{\rm PTherm}^{(m)}\sim N^{2}[f(m)\pm O(1/N)],$ (35)

for some function $\Omega(1)\leq f(m)\leq O(m)$. Depending on the scaling of
$f(m)$, it is possible that the ‘pseudothermalization front’ picture ceases to
apply, with a parametrically slower equilibration behind the ballistic front
and a separation of consecutive pseudothermalization transitions $t_{\rm
PTherm}^{(m)}$, $t_{\rm PTherm}^{(m+1)}$.

  3. (iii)

At even higher moments $\omega(N\log N)\leq m\leq O({\rm poly}(N))$, the
counting bound becomes dominant over the locality bound, so the ‘font’ picture
is not relevant anymore. We have

$t_{\rm PTherm}^{(m)}\sim N[mg(N)\pm O(1)],$ (36)

for some function $\Omega(1/\log N)\leq g(N)\leq O(N)$. Depending on the
scaling of $g(N)$, consecutive pseudothermalization times $t_{\rm
PTherm}^{(m)}$, $t_{\rm PTherm}^{(m+1)}$ may be sharply separated, or the two
crossover windows may orverlap. These two options are sketched in Fig. 5(c).

  4. (iv)

Finally, for $m\geq\omega({\rm poly}(N))$, we have $t_{\rm mix}^{(m)}\sim m$
up to ${\rm poly}(N)$ factors.

####  IV.3.4 Higher dimension

The analysis above could be straightforwardly adapted to higher-dimensional
qubit arrays. For example one can choose a gate set
$\mathcal{G}=\\{u_{i,\hat{v},\hat{w},a,b}\\}$, with $i\in[N]$ a site,
$\hat{v},\hat{w}$ two different directions on the lattice, and
$a,b\in\\{0,1\\}$ control bits, where the gate $u$ flips qubit $i$ iff its
neighbors $i+\hat{v}$, $i+\hat{w}$ are in states $a,b$ respectively.
$\mathcal{G}$ contaits $4\binom{2d}{2}N$ gates, with $d$ the spatial
dimension. In $d=1$ one recovers the gate set used earlier (in that case the
choice of directions $\hat{v}=+\hat{x}$, $\hat{w}=-\hat{x}$ is unique).

For constant $d$, all the results of this section proceed unchanged, except
for the locality bound, Eq. (31); there, we must replace Eq. (31) with $t_{\rm
mix}^{(m)}\geq\Omega(N^{1+1/d})$, with $d$ the spatial dimension. This still
meets the requirement for a cut-off, Eq. (33), but is not enough to pin down
the scaling of pseudothermalization time at constant $m$: we have instead
$\Omega(N^{1+1/d})\leq t_{\rm mix}^{(m)}\leq O(N^{2})$.

Qualitatively, this opens the question of whether, in dimension $d\geq 2$, the
‘pseudothermalization front’ picture still holds; i.e., if the state behind
the light cone equilibrates quickly (which would give $t_{\rm mix}^{(m)}\sim
N^{1+1/d}$), or whether it equilibrates more slowly, in such a way that even
after the light cone has swept through the whole system it is still possible
to distinguish the states from Haar-random ones with a constant number of
copies $m$.

##  V Generation of Pseudoentanglement

Figure 9: Local automaton circuit model for the generation of
pseudoentanglement from product states. (a) Circuit model: the same as Fig.
5(a) but with a product initial state $\ket{0}^{\otimes
N-N_{A}}\otimes\ket{+}^{\otimes N_{A}}$ (no random phases). (b) Schematic
summary of bounds for the pseudothermalization time. Solid lines indicate
unconditional bounds while the dot-dashed line depends on Conjecture 1.
Results are the same as in Fig. 5(b) but with a modified upper bound
$\tilde{O}(m^{2}N^{2})$, Eq. (44). (c) Given the more complex mapping between
quantum ensemble distinguishability and the mixing of classical Markov chains,
Eq. (38-39), the precise nature of pseudothermalization (e.g., whether there
is a cut-off) remains unclear in this case.

Thus far we have discussed the propagation of pseudoentanglement: we have
taken initial states that already contain random subset-phase ensembles on
subsystems, and studied how these turn into random subset-phase states on the
whole system under automaton circuit dynamics. This approach was motivated
partly by analogy with the problem of entanglement spreading in quantum
quenches, and partly by technical simplifications that arise when random
phases are already provided by the input states. In this section, we address
the question of generating pseudoentanglement from initial product states. To
avoid complications associated with random phases, we choose to focus on the
random subset states introduced in Ref. Giurgica-Tiron and Bouland (2023) and
reviewed in Sec. II.1. We use only the local gate set of Sec. IV.3.

Let us take an initial subset state

$\ket{\psi_{S_{0}}}=\ket{0}^{\otimes N-k}\otimes\ket{+}^{\otimes k},$ (37)

corresponding to subset $S_{0}=\\{\mathbf{z}:\ z_{1}=\dots=z_{N-k}=0\\}$ of
cardinality $K=2^{k}$, shown by Fig. 9(a). To be pseudoentangled, random
subset ensembles need to have $\omega({\rm poly}(N))\leq K\leq o(2^{N}/{\rm
poly}(N))$, therefore we take the number of qubits initially in superposition,
$k$, to be $\omega(\log N)\leq k<cN$ with $0<c<1$ a constant. After $t$ time
steps, we have a subset state $\ket{\psi_{S_{t}}}$ where $S_{t}$ is a random
variable on $\Sigma_{K}$ distributed according to a time-evolved distribution
$p_{t}$. To characterize the generation of pseudoentanglement, we look at the
trace distance between $m$-th moment operators, $\|\rho_{t}^{(m)}-\rho_{{\rm
Haar}}^{(m)}\|_{{\rm tr}}$, with
$\rho_{t}^{(m)}=\sum_{S}p_{t}(S)\ket{\psi_{S}}\bra{\psi_{S}}^{\otimes m}$. As
long as $m\leq O({\rm poly}(N))$, the application of Proposition 2 yields

$\displaystyle\|\rho_{t}^{(m)}-\rho_{{\rm Haar}}^{(m)}\|_{\rm tr}$
$\displaystyle\simeq\left\|\mathcal{M}^{(m)}_{p_{t}}\right\|_{\rm tr},$ (38)
$\displaystyle(\mathcal{M}^{(m)}_{p_{t}})_{S^{\prime},S^{\prime\prime}}$
$\displaystyle=\binom{K}{m}^{-1}\binom{K}{m+\delta}f^{(m+\delta)}_{p_{t}}(S^{\prime}\cup
S^{\prime\prime}),$ (39)

where $\delta\equiv|S^{\prime}\setminus S^{\prime\prime}|$ is the difference
between subsets $S^{\prime},S^{\prime\prime}\in\Sigma_{m}$ (so $S^{\prime}\cup
S^{\prime\prime}\in\Sigma_{m+\delta}$), and

$f^{(m+\delta)}_{p_{t}}=\Phi_{m+\delta\leftarrow K}[p_{t}]-\pi_{m+\delta}$
(40)

is the deviation from uniformity of a Markov chain on $\Sigma_{m+\delta}$
induced by the random automaton circuit, as in Sec. IV.

In this case, the pseudothermalization problem is not directly mapped onto the
mixing of a single classical random walk. Rather, it is mapped in a nontrivial
way, via Eq. (39), to the mixing of $m+1$ distinct random walks, indexed by
$\delta=0,1,\dots m$, and distributed according to $\Phi_{m+\delta\leftarrow
K}[p_{t}]$, defined by Eq. (16).

More specifically, the diagonal terms in $\mathcal{M}^{(m)}_{p_{t}}$
correspond to $\delta=0$, i.e., to the walk on $\Sigma_{m}$ distributed as
$\Phi_{m\leftarrow K}[p_{t}]$. This is the same walk that controlled the
spreading of pseudoentanglement in Sec. IV. However, the off-diagonal terms,
corresponding to $\delta>0$, make this problem more complex. While we do not
have an exact solution in this case, one can still obtain nontrivial bounds on
the pseudothermalization time $t_{\rm PTherm}^{(m)}$, as we discuss next.

First, a lower bound on $t_{\rm PTherm}^{(m)}$ is obtained by using the
monotonicity of trace norm under quantum channels:
$\|\mathcal{D}(\mathcal{M})\|_{\rm tr}\leq\|\mathcal{M}\|_{\rm tr}$, where
$\mathcal{D}$ is a dephasing super-operator that annihilates off-diagonal
matrix elements:
$[\mathcal{D}(\mathcal{M})]_{S^{\prime},S^{\prime\prime}}=\delta_{S^{\prime},S^{\prime\prime}}\mathcal{M}_{S^{\prime},S^{\prime}}$.
This gives

$\|\rho_{t}^{(m)}-\rho_{{\rm Haar}}^{(m)}\|_{\rm tr}\geq 2\|\Phi_{m\leftarrow
K}[p_{t}]-\pi_{m}\|_{\rm tv}$ (41)

(up to negligible terms $O(m/\sqrt{K})$), which in turn implies

$t_{\rm PTherm}^{(m)}\geq t_{\rm mix}^{(m)}.$ (42)

Here $t_{\rm mix}^{(m)}$ is the mixing time of the Markov chain
$\Phi_{m\leftarrow K}[p_{t}]$, studied in Sec. IV; the lower bounds we derived
there thus apply directly to the pseudothermalization time in the present
problem (in one spatial dimension):

$t_{\rm PTherm}^{(m)}\geq\Omega(N^{2},mN/\log(N)).$ (43)

Furthermore, since the matrix $\mathcal{M}^{(m)}_{p_{t}}$ that controls
pseudothermalization, Eq. (39), consists of differences
$\Phi_{m+\delta\leftarrow K}[p_{t}]-\pi_{m+\delta}$, we can guarantee
pseudothermalization once all the walks ($\delta=0,1,\dots m$) have come
sufficiently close to equilibrium. By properly taking into account the
multiplicity of each $\delta$, we derive the following upper bound on the
$m$-th pseudothermalization time (see Appendix G):

$t_{\rm PTherm}^{(m)}\leq\tilde{O}(m^{2}N^{2}),$ (44)

where $\tilde{O}$ indicates an upper bound up to logarithmic factors. Again,
this is conditional on Conjecture 1.

Based on these bounds, summarized in Fig. 9(b), we can conclude that

  1. (i)

at early times $t\leq\tilde{O}(N^{2},mN)$, the system can be efficiently
distinguished from Haar-random at the $m$-copy level.

  2. (ii)

at late times $t\geq\tilde{\Omega}(m^{2}N^{2})$, the system cannot be
efficiently distinguished from Haar-random at the $m$-copy level.

The nature of the crossover between these two regimes is an interesting open
question, see Fig. 9(c). Under Conjecture 1, we have $t^{(m)}_{\rm
mix}/t^{(m)}_{\rm rel}\to\infty$ in this limit for all $m$, suggestive of a
cutoff phenomenon in the associated Markov chains $\Phi_{m\leftarrow
K}[p_{t}]$. How would such cutoff phenomena be reflected in the quantum
distinguishability measure Eq. (38)? Would the latter undergo an abrupt
transition, would it decay gradually, or would it exhibit more complex
behavior (e.g. a sequence of pre-cutoffs associated to each $\delta$)? We
leave this interesting question as a goal for future work.

##  VI Discussion

###  VI.1 Summary

In this work, we discussed the dynamical generation and propagation of
pseudoentanglement. We considered automaton circuits as models of dynamics
that preserve the relevant ensembles of states: both the subset-phase states
and the subset states, which have been introduced as “standard models” of
pseudoentanglement in previous work Aaronson  _et al._ (2023); Giurgica-Tiron
and Bouland (2023). Specifically, we considered two models of automaton
dynamics: an all-to-all interacting circuit, and a circuit of geometrically-
local, 3-body, Toffoli-like gates on a $d$-dimensional qubit lattice (focusing
mostly on $d=1$). The all-to-all interacting model is exactly solvable, but
the time scales for generation and propagation of pseudoentanglement are
unphysically large (exponential in system size); the local circuit model does
not suffer from these exponential time scales, but is less analytically
tractable, so that our results are limited to a combination of rigorous bounds
and conjectures supported by numerics.

We separately analyzed two related problems: the spreading of
pseudoentanglement when a small pseudoentangled subsystem is placed in contact
with a large trivial system, and the generation of pseudoentanglement starting
from a product state. In both cases the key technical step involves mapping
the problem onto the equilibration of certain classical random walks over
subsets of the computational basis. This creates a connection to an
interesting feature of Markov chains on highly-connected spaces: the cutoff
phenomenon, whereby mixing happens in a step-like fashion in the limit of
large system size, giving asymptotically sharp transitions to equilibrium. We
argued that our Markov chains in many cases exhibit a cutoff phenomenon,
giving sharp dynamical transitions in the distinguishability of our output
states from Haar-random ones.

We named the dynamical approach toward pseudoentangled ensembles
pseudothermalization, motivated by the analogy with the role of entanglement
in quantum thermalization. This new concept adds to a growing body of research
on the nuances of thermalization, chaos, ergodicity, scrambling, etc, that can
be revealed only by higher-moment statistical quantities, thus going beyond
the first- or second-moment quantities (e.g., local expectation values
Srednicki (1994), out-of-time-ordered correlators von Keyserlingk  _et al._
(2018)) that are traditionally the basis of such notions.

Other recently introduced ideas in this direction include deep thermalization
Cotler  _et al._ (2023); Ho and Choi (2022); Ippoliti and Ho (2023, 2022);
Lucas  _et al._ (2023); Bhore  _et al._ (2023); Shrotriya and Ho (2023); Chan
and De Luca (2024), complete Hilbert-space ergodicity Pilatowsky-Cameo  _et
al._ (2023, 2024), and thermalization via free probability or cumulants Fava
_et al._ (2023); Pappalardi  _et al._ (2023). However, those are stronger
notions, based on information-theoretic (pseudo-)randomness of certain state
ensembles generated for example by quantum measurements or by quasiperiodic
drive sequences. The randomness can be quantified through the notion
$\epsilon$-approximate quantum state $m$-design, corresponding to the trace
distance between moment operators being smaller than a tolerance $\epsilon$
(as in Eq. (2)); in those stronger versions of thermalization, one generally
requires or expects an exponentially small tolerance $\epsilon\sim 1/\exp(N)$,
or possibly even exact state designs ($\epsilon=0$). For $m\geq 2$, these
conditions imply high entanglement. Pseudothermalization in contrast is based
on a computational notion of (pseudo-)randomness, where the observer is
limited to polynomial resources (e.g. number of state copies) in their attempt
to tell the states apart from truly random ones. Again in the language of
state designs, this corresponds to a tolerance $\epsilon=o(1/{\rm poly}(N))$
that is superpolynomially small (in order to trick the polynomially-bound
observer), but not exponentially small, thus allowing for weakly entangled
states.

###  VI.2 Outlook

Our work opens a number of directions for future research. First of all, it
would be desirable to improve and/or rigorize our results in several ways: by
tightening bounds, especially on the mixing times in the local circuit model
and on pseudothermalization without random phases; by proving or disproving
our Conjecture 1 on the relaxation times in the local circuit model; and by
proving the existence of cutoff phenomena in our Markov chains, for which we
only showed a necessary condition (divergent ratio of mixing time and
relaxation time). These improvements could have significant qualitative
consequences on the physics. For example, one may hope to show that, in the
local circuit model for pseudoentanglement spreading, Sec. IV.3, the $m$-copy
pseudothermalization transitions are sharply separated from one another, in
the sense that the difference of consecutive mixing times $t_{\rm
mix}^{(m+1)}-t_{\rm mix}^{(m)}$ is sufficiently larger than the relaxation
time $t_{\rm rel}^{(m)}$.

Another direction for future work is to identify physically realistic models
of dynamics where pseudoentanglement is achieved in polynomial time. In the
all-to-all model of Sec. IV.2 that time scale is exponential for trivial
reasons; in the local model of Sec. IV.3, the $m$-copy pseudothermalization
time is polynomial in system size, but grows at least linearly with $m$, so
that convergence to a pseudoentangled ensemble, indistinguishabile from Haar-
random with any polynomial number of copies $m=O({\rm poly}(N))$, requires
super-polynomial time. This can be a mild super-polynomial scaling, e.g.
$N^{\log N}$. Nonetheless it would be interesting to identify models of
dynamics that generate pseudoentanglement in polynomial time. Such models
cannot be random local automaton circuits, due to the counting bound, Eq.
(30). Pseudorandom subset(-phase) state ensembles can be prepared efficiently
by using quantum-secure pseudorandom permutations and phase functions Aaronson
_et al._ (2023). In this work we have opted for truly random subsets, as these
make for simpler models of circuit dynamics and are more in keeping with the
spirit of quantum thermalization; it would be interesting to identify
pseudorandom models of dynamics that can produce pseudoentanglement in
polynomial time while maintaining a simple, “random-looking”, physically
reasonable structure. A related construction of pseudorandom unitaries was
recently shown to not be robust to perturbations Haug  _et al._ (2023), which
might suggest intrinsic difficulties for this direction.

More broadly, subset(-phase) states are only one class of pseudoentangled
state ensembles. Due to their structure, these ensembles suggest automaton
circuits as a natural models of dynamics. It would be interesting to identify
different ensembles of pseudoentangled states and the types of quantum
dynamics that might naturally produce them.

Recent works have introduced related notions of ‘pseudocoherence’ Haug  _et
al._ (2023), ‘pseudomagic’ Gu  _et al._ (2023), etc, which are analogous to
pseudoentanglement but based on different resources (coherence, magic, etc). A
natural extensions of our work would be to study the dynamics of other such
pseudoresources in appropriate circuit models.

It is also interesting to ask whether, in analogy with measurement-induced
entanglement phase transitions Li  _et al._ (2018); Skinner  _et al._ (2019);
Gullans and Huse (2020); Iaconis  _et al._ (2020), one might obtain phase
transitions in the dynamics of pseudoentanglement by enriching our automaton
circuits with suitable projective measurements and other elements. For
example, while the size of subsets is invariant under automaton gates, it
generally decreases under $Z$ measurements, and can increase under Hadamard
gates and $X$ measurements. Can the competition of these elements give rise to
nontrivial phase transitions in the pseudoentanglement properties of late-time
states?

Finally, the possible emergence of cutoff phenomena in pseudothermalization
raises more general questions about the sharpness of thermalization in generic
quantum dynamics. Are sharp transitions in distinguishability a generic
feature of chaotic quantum dynamics Kastoryano  _et al._ (2012); Oh and Kais
(2023)? What are the right tools to analyze them? Can there particularly
structured classes of time evolutions (e.g., with symmetry) exhibit different
levels of sharpness? We leave these questions as an interesting target for
future research.

###### Acknowledgements.

We thank Nick Hunter-Jones, Ethan Lake, and Scott Aaronson for stimulating
discussions. We are especially grateful to Nick Hunter-Jones for pointing us
to Ref. Gowers (1996). Some numerical simulations were carried out on HPC
resources provided by the Texas Advanced Computing Center (TACC) at the
University of Texas at Austin.

## References

  * Horodecki and Horodecki (2009) Ryszard Horodecki and Pawe Horodecki, “Quantum entanglement,” Rev. Mod. Phys. 81, 865–942 (2009).
  * Chitambar and Gour (2019) Eric Chitambar and Gilad Gour, “Quantum resource theories,” Reviews of Modern Physics 91, 025001 (2019).
  * Kitaev and Preskill (2006) Alexei Kitaev and John Preskill, “Topological Entanglement Entropy,” Physical Review Letters 96, 110404 (2006).
  * Li and Haldane (2008) Hui Li and F. D. M. Haldane, “Entanglement Spectrum as a Generalization of Entanglement Entropy: Identification of Topological Order in Non-Abelian Fractional Quantum Hall Effect States,” Physical Review Letters 101, 010504 (2008).
  * Jiang  _et al._ (2012) Hong-Chen Jiang, Zhenghan Wang, and Leon Balents, “Identifying topological order by entanglement entropy,” Nature Physics 8, 902–905 (2012).
  * Rigol  _et al._ (2008) Marcos Rigol, Vanja Dunjko, and Maxim Olshanii, “Thermalization and its mechanism for generic isolated quantum systems,” Nature 452, 854–858 (2008).
  * Nandkishore and Huse (2015) Rahul Nandkishore and David A. Huse, “Many-Body Localization and Thermalization in Quantum Statistical Mechanics,” Annual Review of Condensed Matter Physics 6, 15–38 (2015).
  * Kaufman  _et al._ (2016) Adam M. Kaufman, M. Eric Tai, Alexander Lukin, Matthew Rispoli, Robert Schittko, Philipp M. Preiss, and Markus Greiner, “Quantum thermalization through entanglement in an isolated many-body system,” Science 353, 794–800 (2016).
  * Mezei and van der Schee (2020) Mark Mezei and Wilke van der Schee, “Black Holes Often Saturate Entanglement Entropy the Fastest,” Physical Review Letters 124, 201601 (2020).
  * Alhambra (2022) Alvaro M. Alhambra, “Quantum many-body systems in thermal equilibrium,”  (2022), 10.48550/arXiv.2204.08349.
  * Vidmar and Rigol (2016) Lev Vidmar and Marcos Rigol, “Generalized Gibbs ensemble in integrable lattice models,” Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
  * Nahum  _et al._ (2017) Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah, “Quantum Entanglement Growth under Random Unitary Dynamics,” Physical Review X 7, 031016 (2017).
  * Bertini  _et al._ (2019) Bruno Bertini, Pavel Kos, and Tomaz Prosen, “Entanglement Spreading in a Minimal Model of Maximal Many-Body Quantum Chaos,” Physical Review X 9, 021033 (2019).
  * Zhou and Nahum (2020) Tianci Zhou and Adam Nahum, “Entanglement Membrane in Chaotic Many-Body Systems,” Phys. Rev. X 10, 031066 (2020).
  * Brandao  _et al._ (2021) Fernando G. S. L. Brandao, Wissam Chemissany, Nicholas Hunter-Jones, Richard Kueng, and John Preskill, “Models of quantum complexity growth,” PRX Quantum 2, 030316 (2021).
  * Fisher  _et al._ (2022) Matthew P. A. Fisher, Vedika Khemani, Adam Nahum, and Sagar Vijay, “Random Quantum Circuits,” arXiv e-prints (2022), 10.48550/arXiv.2207.14280.
  * Potter and Vasseur (2022) Andrew C. Potter and Romain Vasseur, “Entanglement Dynamics in Hybrid Quantum Circuits,” in  _Entanglement in Spin Chains: From Theory to Quantum Technology Applications_, Quantum Science and Technology, edited by Abolfazl Bayat, Sougato Bose, and Henrik Johannesson (Cham, 2022) pp. 211–249.
  * Mittal and Hunter-Jones (2023) Shivan Mittal and Nicholas Hunter-Jones, “Local random quantum circuits form approximate designs on arbitrary architectures,”  (2023), 10.48550/arXiv.2310.19355.
  * Pilatowsky-Cameo  _et al._ (2023) Saul Pilatowsky-Cameo, Ceren B. Dag, Wen Wei Ho, and Soonwon Choi, “Complete Hilbert-Space Ergodicity in Quantum Dynamics of Generalized Fibonacci Drives,” Physical Review Letters 131, 250401 (2023).
  * Pilatowsky-Cameo  _et al._ (2024) Saul Pilatowsky-Cameo, Iman Marvian, Soonwon Choi, and Wen Wei Ho, “Hilbert-Space Ergodicity in Driven Quantum Systems: Obstructions and Designs,”  (2024), 10.48550/arXiv.2402.06720.
  * Page (1993) Don N. Page, “Average entropy of a subsystem,” Phys. Rev. Lett. 71, 1291–1294 (1993).
  * Amico  _et al._ (2008) Luigi Amico, Rosario Fazio, Andreas Osterloh, and Vlatko Vedral, “Entanglement in many-body systems,” Reviews of Modern Physics 80, 517–576 (2008).
  * Islam  _et al._ (2015) Rajibul Islam, Ruichao Ma, Philipp M. Preiss, M. Eric Tai, Alexander Lukin, Matthew Rispoli, and Markus Greiner, “Measuring entanglement entropy in a quantum many-body system,” Nature 528, 77–83 (2015).
  * Elben  _et al._ (2018) A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, “R\’enyi Entropies from Random Quenches in Atomic Hubbard and Spin Models,” Physical Review Letters 120, 050406 (2018).
  * Brydges  _et al._ (2019) Tiff Brydges, Andreas Elben, Petar Jurcevic, Benoit Vermersch, Christine Maier, Ben P. Lanyon, Peter Zoller, Rainer Blatt, _et al._ , “Probing Renyi entanglement entropy via randomized measurements,” Science 364, 260–263 (2019).
  * Huang  _et al._ (2020) Hsin-Yuan Huang, Richard Kueng, and John Preskill, “Predicting many properties of a quantum system from very few measurements,” Nature Physics 16, 1050–1057 (2020).
  * Elben  _et al._ (2023) Andreas Elben, Steven T. Flammia, Hsin-Yuan Huang, Richard Kueng, John Preskill, Benoit Vermersch, and Peter Zoller, “The randomized measurement toolbox,” Nature Reviews Physics 5, 9–24 (2023).
  * Aaronson  _et al._ (2023) Scott Aaronson, Adam Bouland, Bill Fefferman, Soumik Ghosh, Umesh Vazirani, Chenyi Zhang, and Zixin Zhou, “Quantum Pseudoentanglement,”  (2023), 10.48550/arXiv.2211.00747.
  * Giurgica-Tiron and Bouland (2023) Tudor Giurgica-Tiron and Adam Bouland, “Pseudorandomness from Subset States,”  (2023), 10.48550/arXiv.2312.09206.
  * Ji  _et al._ (2018) Zhengfeng Ji, Yi-Kai Liu, and Fang Song, “Pseudorandom Quantum States,” in  _Advances in Cryptology - CRYPTO 2018_, edited by Hovav Shacham and Alexandra Boldyreva (Cham, 2018) pp. 126–152.
  * Lu  _et al._ (2023) Chuhan Lu, Minglong Qin, Fang Song, Penghui Yao, and Mingnan Zhao, “Quantum Pseudorandom Scramblers,”  (2023), 10.48550/arXiv.2309.08941.
  * Zhandry (2016) Mark Zhandry, “A Note on Quantum-Secure PRPs,”  (2016), 10.48550/arXiv.1611.05564.
  * Gopalakrishnan (2018) Sarang Gopalakrishnan, “Operator growth and eigenstate entanglement in an interacting integrable Floquet system,” Physical Review B 98, 060302 (2018).
  * Iaconis  _et al._ (2019) Jason Iaconis, Sagar Vijay, and Rahul Nandkishore, “Anomalous subdiffusion from subsystem symmetries,” Physical Review B 100, 214301 (2019).
  * Iaconis (2021) Jason Iaconis, “Quantum State Complexity in Computationally Tractable Quantum Circuits,” PRX Quantum 2, 010329 (2021).
  * Klobas  _et al._ (2021) Katja Klobas, Bruno Bertini, and Lorenzo Piroli, “Exact Thermalization Dynamics in the “Rule 54” Quantum Cellular Automaton,” Physical Review Letters 126, 160602 (2021).
  * Diaconis and Shahshahani (1981) Persi Diaconis and Mehrdad Shahshahani, “Generating a random permutation with random transpositions,” Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 57, 159–179 (1981).
  * Diaconis (1996) P Diaconis, “The cutoff phenomenon in finite Markov chains.” Proceedings of the National Academy of Sciences 93, 1659–1664 (1996).
  * Chen and Saloff-Coste (2008) Guan-Yu Chen and Laurent Saloff-Coste, “The cutoff phenomenon for ergodic Markov processes.” Electronic Journal of Probability [electronic only] 13, 26–78 (2008).
  * Lacoin and Leblond (2011) Hubert Lacoin and Remi Leblond, “Cutoff phenomenon for the simple exclusion process on the complete graph,”  (2011), 10.48550/arXiv.1010.4866.
  * Chen (2021) Joe P. Chen, “The cutoff profile for exclusion processes in any dimension,”  (2021), 10.48550/arXiv.2106.03685.
  * Aldous and Diaconis (1986) David Aldous and Persi Diaconis, “Shuffling Cards and Stopping Times,” The American Mathematical Monthly 93, 333–348 (1986).
  * Roberts and Yoshida (2017) Daniel A. Roberts and Beni Yoshida, “Chaos and complexity by design,” Journal of High Energy Physics 2017, 121 (2017).
  * Hunter-Jones (2019) Nicholas Hunter-Jones, “Unitary designs from statistical mechanics in random quantum circuits,”  (2019), 10.48550/arXiv.1905.12053.
  * Cotler  _et al._ (2023) Jordan S. Cotler, Daniel K. Mark, Hsin-Yuan Huang, Felipe Hernandez, Joonhee Choi, Adam L. Shaw, Manuel Endres, and Soonwon Choi, “Emergent Quantum State Designs from Individual Many-Body Wave Functions,” PRX Quantum 4, 010311 (2023).
  * Ippoliti and Ho (2022) Matteo Ippoliti and Wen Wei Ho, “Solvable model of deep thermalization with distinct design times,” Quantum 6, 886 (2022).
  * Fava  _et al._ (2023) Michele Fava, Jorge Kurchan, and Silvia Pappalardi, “Designs via Free Probability,”  (2023), 10.48550/arXiv.2308.06200.
  * Nielsen and Chuang (2010) Michael A Nielsen and Isaac L Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Levin and Peres (2017) David A Levin and Yuval Peres,  _Markov chains and mixing times_ , Vol. 107 (American Mathematical Soc., 2017).
  * Hoory  _et al._ (2006) Shlomo Hoory, Nathan Linial, and Avi Wigderson, “Expander graphs and their applications,” Bulletin of the American Mathematical Society 43, 439–561 (2006).
  * Peres (2004) Yuval Peres, “ARCC workshop: Sharp thresholds for mixing times,”  (2004), online notes from a workshop of the American Institute of Mathematics, transcribed by S. Goel, A. Nachmias, P. Ralph and D. Shah.
  * Caputo  _et al._ (2010) Pietro Caputo, Thomas M. Liggett, and Thomas Richthammer, “Proof of Aldous’ spectral gap conjecture,” Journal of the American Mathematical Society 23, 831–851 (2010).
  * Calabrese and Cardy (2005) Pasquale Calabrese and John Cardy, “Evolution of entanglement entropy in one-dimensional systems,” Journal of Statistical Mechanics: Theory and Experiment 2005, P04010 (2005).
  * Calabrese and Cardy (2009) Pasquale Calabrese and John Cardy, “Entanglement entropy and conformal field theory,” Journal of Physics A: Mathematical and Theoretical 42, 504005 (2009).
  * Prosen and Žnidarič (2009) Tomaž Prosen and Marko Žnidarič, “Matrix product simulations of non-equilibrium steady states of quantum spin chains,” Journal of Statistical Mechanics: Theory and Experiment 2009, P02035 (2009).
  * Znidaric  _et al._ (2011) Marko Znidaric, Bojan Zunkovic, and Tomaz Prosen, “Transport properties of a boundary-driven one-dimensional gas of spinless fermions,” Physical Review E 84, 051115 (2011).
  * Morningstar  _et al._ (2022) Alan Morningstar, Luis Colmenarez, Vedika Khemani, David J. Luitz, and David A. Huse, “Avalanches and many-body resonances in many-body localized systems,” Physical Review B 105, 174205 (2022).
  * Sels (2022) Dries Sels, “Bath-induced delocalization in interacting disordered spin chains,” Physical Review B 106, L020202 (2022).
  * Hermon and Pymar (2020) Jonathan Hermon and Richard Pymar, “The exclusion process mixes (almost) faster than independent particles,”  (2020), 10.48550/arXiv.1808.10846.
  * Diaconis  _et al._ (1990) Persi Diaconis, R. L. Graham, and J. A. Morrison, “Asymptotic analysis of a random walk on a hypercube with many dimensions,” Random Structures & Algorithms 1, 51–72 (1990).
  * Gowers (1996) W. T. Gowers, “An Almost m-wise Independent Random Permutation of the Cube,” Combinatorics, Probability and Computing 5, 119–130 (1996).
  * Hoory  _et al._ (2005) Shlomo Hoory, Avner Magen, Steven Myers, and Charles Rackoff, “Simple permutations mix well,” Theoretical Computer Science Automata, Languages and Programming: Algorithms and Complexity (ICALP-A 2004), 348, 251–261 (2005).
  * Brodsky and Hoory (2008) Alex Brodsky and Shlomo Hoory, “Simple permutations mix even better,” Random Structures & Algorithms 32, 274–289 (2008).
  * Jerison (2013) Daniel Jerison, “General mixing time bounds for finite Markov chains via the absolute spectral gap,”  (2013), 10.48550/arXiv.1310.8021.
  * Srednicki (1994) Mark Srednicki, “Chaos and quantum thermalization,” Physical Review E 50, 888–901 (1994).
  * von Keyserlingk  _et al._ (2018) C. W. von Keyserlingk, Tibor Rakovszky, Frank Pollmann, and S. L. Sondhi, “Operator Hydrodynamics, OTOCs, and Entanglement Growth in Systems without Conservation Laws,” Phys. Rev. X 8, 021013 (2018).
  * Ho and Choi (2022) Wen Wei Ho and Soonwon Choi, “Exact Emergent Quantum State Designs from Quantum Chaotic Dynamics,” Physical Review Letters 128, 060601 (2022).
  * Ippoliti and Ho (2023) Matteo Ippoliti and Wen Wei Ho, “Dynamical Purification and the Emergence of Quantum State Designs from the Projected Ensemble,” PRX Quantum 4, 030322 (2023).
  * Lucas  _et al._ (2023) Maxime Lucas, Lorenzo Piroli, Jacopo De Nardis, and Andrea De Luca, “Generalized deep thermalization for free fermions,” Physical Review A 107, 032215 (2023).
  * Bhore  _et al._ (2023) Tanmay Bhore, Jean-Yves Desaules, and Zlatko Papic, “Deep thermalization in constrained quantum systems,”  (2023), 10.48550/arXiv.2307.03769.
  * Shrotriya and Ho (2023) Harshank Shrotriya and Wen Wei Ho, “Nonlocality of Deep Thermalization,”  (2023), 10.48550/arXiv.2305.08437.
  * Chan and De Luca (2024) Amos Chan and Andrea De Luca, “Projected state ensemble of a generic model of many-body quantum chaos,”  (2024).
  * Pappalardi  _et al._ (2023) Silvia Pappalardi, Felix Fritzsch, and Tomaz Prosen, “General Eigenstate Thermalization via Free Cumulants in Quantum Lattice Systems,”  (2023), 10.48550/arXiv.2303.00713.
  * Haug  _et al._ (2023) Tobias Haug, Kishor Bharti, and Dax Enshan Koh, “Pseudorandom unitaries are neither real nor sparse nor noise-robust,”  (2023), 10.48550/arXiv.2306.11677.
  * Gu  _et al._ (2023) Andi Gu, Lorenzo Leone, Soumik Ghosh, Jens Eisert, Susanne Yelin, and Yihui Quek, “A little magic means a lot,”  (2023), 10.48550/arXiv.2308.16228.
  * Li  _et al._ (2018) Yaodong Li, Xiao Chen, and Matthew P. A. Fisher, “Quantum Zeno effect and the many-body entanglement transition,” Physical Review B 98, 205136 (2018).
  * Skinner  _et al._ (2019) Brian Skinner, Jonathan Ruhman, and Adam Nahum, “Measurement-Induced Phase Transitions in the Dynamics of Entanglement,” Physical Review X 9, 031009 (2019).
  * Gullans and Huse (2020) Michael J. Gullans and David A. Huse, “Dynamical Purification Phase Transition Induced by Quantum Measurements,” Physical Review X 10, 041020 (2020).
  * Iaconis  _et al._ (2020) Jason Iaconis, Andrew Lucas, and Xiao Chen, “Measurement-induced phase transitions in quantum automaton circuits,” Physical Review B 102, 224311 (2020).
  * Kastoryano  _et al._ (2012) Michael J. Kastoryano, David Reeb, and Michael M. Wolf, “A Cutoff Phenomenon for Quantum Markov Chains,” Journal of Physics A: Mathematical and Theoretical 45, 075307 (2012).
  * Oh and Kais (2023) Sangchul Oh and Sabre Kais, “Cutoff phenomenon and entropic uncertainty for random quantum circuits,” Electronic Structure 5, 035004 (2023).
  * Shende  _et al._ (2003) V.V. Shende, A.K. Prasad, I.L. Markov, and J.P. Hayes, “Synthesis of reversible logic circuits,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 22, 710–722 (2003).
  * Aaronson  _et al._ (2015) Scott Aaronson, Daniel Grier, and Luke Schaeffer, “The Classification of Reversible Bit Operations,”  (2015), 10.48550/arXiv.1504.05155.
  * Grewal  _et al._ (2023) Sabee Grewal, Vishnu Iyer, William Kretschmer, and Daniel Liang, “Low-Stabilizer-Complexity Quantum States Are Not Pseudorandom,” in  _14th Innovations in Theoretical Computer Science Conference (ITCS 2023)_, Leibniz International Proceedings in Informatics (LIPIcs), Vol. 251, edited by Yael Tauman Kalai (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2023) pp. 64:1–64:20.

##  Appendix A Proof of Proposition 1

Here we prove Proposition 1 by analyzing the trace distance

$\left\|\rho^{(m)}_{\mathcal{E}_{p}}-\rho^{(m)}_{\rm Haar}\right\|_{\rm tr}$
(45)

between the $m$-th moment operators for the subset-phase ensemble
$\mathcal{E}_{p}$ (with uniformly random phases but arbitrarily-distributed
subsets, $S\sim p$) and the Haar ensemble.

###  A.1 Preliminaries

In this Appendix we will denote computational basis states by integers
$z\in[D]=\\{0,1,\dots D-1\\}$ rather than by bitstrings
$\mathbf{z}\in\\{0,1\\}^{N}$. We will reserve the latter notation for tuples
of computational basis states, which are used later.

Eq. (45) takes place within the symmetric sector of the $m$-fold replicated
Hilbert space $\mathcal{H}^{\otimes m}$. It is thus helpful to introduce an
orthonormal basis of that subspace:

$\ket{\mathbf{T}}\equiv\binom{m}{\mathbf{T}}^{-1/2}\sum_{\mathbf{z}:\,{\rm
type}(\mathbf{z})=\mathbf{T}}\ket{\mathbf{z}},$ (46)

where $\mathbf{z}\in[D]^{m}$ is an $m$-tuple of basis states and ${\rm
type}(\mathbf{z})\in[m+1]^{D}$ is a vector that counts the multiplicity of
each basis state $z$ in the $m$-tuple $\mathbf{z}$ (the notation is borrowed
from Ref. Aaronson  _et al._ (2023)): i.e., if ${\rm
type}(\mathbf{z})=\mathbf{T}$, then $T_{i}$ is the number of times the basis
state $i\in[D]$ appears in $\mathbf{z}$. As a consequence we have
$\sum_{i}T_{i}=m$. Finally, we used a shorthand for the multinomial
coefficient

$\binom{m}{\mathbf{T}}=\frac{m!}{T_{0}!\cdots T_{D-1}!}.$

Types that consist only of 0’s and 1’s are called “unique types”. These
describe non-degenerate $m$-tuples $\mathbf{z}$, where all elements are
distinct. As such, unique types are equivalent to subsets of cardinality $m$:
$\mathbf{T}\leftrightarrow S=\\{z:\,T_{z}=1\\}$. With slight abuse of notation
we will denote a unique-type state $\ket{\mathbf{T}}$ by the corresponding
subset, $\ket{S}$.

We may write the Haar ensemble’s moment in terms of type states as

$\rho^{(m)}_{\rm Haar}=\binom{D+m-1}{m}^{-1}\sum_{\mathbf{T}\in{\rm
types}}\outerproduct{\mathbf{T}}{\mathbf{T}}.$ (47)

Further, it is useful to introduce the “unique-type state”

$\displaystyle\rho_{{\rm unique}}^{(m)}$
$\displaystyle=\binom{D}{m}^{-1}\sum_{\mathbf{T}\in\text{unique
types}}\outerproduct{\mathbf{T}}{\mathbf{T}}$
$\displaystyle=\binom{D}{m}^{-1}\sum_{S^{\prime}\in\Sigma_{m}}\outerproduct{S^{\prime}}{S^{\prime}}.$
(48)

It is easy to verify that, as $D\to\infty$ with finite $m$, most types become
unique, and $\rho_{\rm Haar}^{(m)}\simeq\rho_{\rm unique}^{(m)}$. We also
introduce the same object restricted to basis states within a particular
subset $S\in\Sigma_{K}$:

$\displaystyle\rho_{S,{\rm unique}}^{(m)}$
$\displaystyle=\binom{K}{m}^{-1}\sum_{S^{\prime}\in\Sigma_{m}:\,S^{\prime}\subset
S}\outerproduct{S^{\prime}}{S^{\prime}}.$ (49)

###  A.2 Proof

Our proof builds on Ref. Aaronson  _et al._ (2023), specifically the two facts
below:

######  Fact A.1.

We have

$\left\|\rho_{\rm Haar}^{(m)}-\sum_{S\in\Sigma_{K}}\pi_{K}(S)\rho_{S,{\rm
unique}}^{(m)}\right\|_{\rm tr}\leq O\left(\frac{m^{2}}{K}\right),$

with $\pi_{K}(S)=\binom{D}{K}^{-1}$ the uniform distribution on $\Sigma_{K}$.

###### Proof.

See Proposition 2.3 in Ref. Aaronson  _et al._ (2023). ∎

######  Fact A.2.

We have

$\left\|\rho_{\mathcal{E}_{p}}^{(m)}-\sum_{S\in\Sigma_{K}}p(S)\rho_{S,{\rm
unique}}^{(m)}\right\|_{\rm tr}\leq O\left(\frac{m^{2}}{K}\right),$

###### Proof.

See Proposition 2.4 in Ref. Aaronson  _et al._ (2023). ∎

By using these two facts and the triangle inequality, we can rephrase the LHS
of Eq. (45) as

$\displaystyle\left\|\rho^{(m)}_{\mathcal{E}_{p}}-\rho^{(m)}_{\rm
Haar}\right\|_{\rm tr}$
$\displaystyle=\left\|\sum_{S\in\Sigma_{K}}(p(S)-\pi_{K}(S))\rho_{S,{\rm
unique}}^{(m)}\right\|_{\rm tr}$
$\displaystyle\qquad+O\left(\frac{m^{2}}{K}\right).$ (50)

Now, using the definition of unique-type state Eq. (49) and orthonormality of
the unique-type states, we arrive at

$\displaystyle\left\|\rho^{(m)}_{\mathcal{E}_{p}}-\rho^{(m)}_{\rm
Haar}\right\|_{\rm tr}$
$\displaystyle=\sum_{S^{\prime}\in\Sigma_{m}}|f(S^{\prime})|+O\left(\frac{m^{2}}{K}\right),$
(51) $\displaystyle f(S^{\prime})$
$\displaystyle=\binom{K}{m}^{-1}\sum_{\begin{subarray}{c}S\in\Sigma_{K}:\\\
S^{\prime}\subset S\end{subarray}}[p(S)-\pi_{K}(S)]$ (52)

Finally, we obtain the statement of Proposition 1 by introducing the map
$\Phi_{K\leftarrow m}$ of Eq. (16).

##  Appendix B Properties of the $\Phi$ maps

Here we briefly analyze and prove some properties of the $\Phi$ maps.

######  Fact B.1.

$\Phi_{m\leftarrow K}$ maps probability distributions on $\Sigma_{K}$ to
probability distributions on $\Sigma_{m}$.

###### Proof.

It is straightforward to verify that $\Phi_{m\leftarrow K}[p]$ corresponds to
the distribution on $\Sigma_{m}$ obtained by composing these two processes:
(i) draw $S\in\Sigma_{K}$ according to $p$; (ii) draw (without replacement)
$m$ elements from $S$ uniformly at random. ∎

######  Fact B.2.

The maps $\Phi$ define a semigroup:

$\Phi_{c\leftarrow b}\circ\Phi_{b\leftarrow a}=\Phi_{c\leftarrow a}$ (53)

for all $D\geq a>b>c\geq 1$.

###### Proof.

We have

$\displaystyle\Phi_{c\leftarrow b}\circ\Phi_{b\leftarrow
a}[p](S^{\prime\prime})$
$\displaystyle=\binom{a}{b}^{-1}\binom{b}{c}^{-1}\sum_{\begin{subarray}{c}S,S^{\prime}:\\\
S^{\prime\prime}\subset S^{\prime}\subset S\end{subarray}}p(S)$
$\displaystyle=\binom{a}{b}^{-1}\binom{b}{c}^{-1}\binom{a-c}{b-c}\binom{a}{c}$
$\displaystyle\qquad\times\Phi_{c\leftarrow a}[p](S^{\prime\prime}),$ (54)

where in the second line we carried out the summation over $S^{\prime}$ (the
$c$ elements of $S^{\prime}\cap S^{\prime\prime}$ are fixed, whereas the $b-c$
elements of $S^{\prime}\setminus S^{\prime\prime}$ must be chosen from
$S\setminus S^{\prime\prime}$, which has $a-c$ elements). Straightforward
algebra shows the product of binomial coefficients is 1. Alternatively,
without algebra, one can note that, given a set $S$ with $|S|=a$, the
following two processes are equivalent: (i) draw a uniformly-random subset
$S^{\prime}\subset S$ with $|S^{\prime}|=b$, and then draw a uniformly-random
subset $S^{\prime\prime}\subset S^{\prime}$ with $|S^{\prime\prime}|=c$; (ii)
draw a uniformly-random subset $S^{\prime\prime}\subset S$ with
$|S^{\prime\prime}|=c$. ∎

Next we focus on the relationship between the distributions $\Phi_{m\leftarrow
K}[p]$ with different values of $m$. Let us consider the relevant situation of
a random walk on the permutation group $\mathcal{S}_{D}$, represented by a
probability distribution $\mathbb{P}_{t}(\sigma)$. This induces random walks
on the spaces $\Sigma_{m}$ for all $m$: picking an initial subset
$S_{0}\in\Sigma_{m}$, we can define

$p_{t,S_{0}}^{(m)}(S)=\sum_{\sigma\in\mathcal{S}_{D}:\,\sigma(S_{0})=S}\mathbb{P}_{t}(\sigma),$
(55)

with $\sigma(S_{0})$ representing the action of permutations on subsets.

######  Fact B.3.

Given a Markov chain $p_{t,S_{0}}^{(K)}(S)$ on $\Sigma_{K}$ induced by a
random walk on permutations as in Eq. (55), we have

$\Phi_{m\leftarrow
K}\left[p_{t,S_{0}}^{(K)}\right](S^{\prime})=\mathbb{E}_{S_{0}^{\prime}}\left[p^{(m)}_{t,S_{0}^{\prime}}(S^{\prime})\right],$
(56)

where the average is taken over subsets $S_{0}^{\prime}\in\Sigma_{m}$,
$S_{0}^{\prime}\subset S_{0}$, according to the uniform measure.

###### Proof.

Plugging in the definition of $\Phi$, the LHS of Eq. (56) reads

$\displaystyle\Phi_{m\leftarrow K}\left[p_{t,S_{0}}^{(K)}\right](S^{\prime})$
$\displaystyle=\binom{K}{m}^{-1}\sum_{S:\,S^{\prime}\subset
S}p_{t,S_{0}}^{(K)}(S)$
$\displaystyle=\binom{K}{m}^{-1}\sum_{S:\,S^{\prime}\subset
S}\sum_{\begin{subarray}{c}\sigma\in\mathcal{S}_{D}:\\\
\sigma(S_{0})=S\end{subarray}}\mathbb{P}_{t}(\sigma)$
$\displaystyle=\binom{K}{m}^{-1}\sum_{\begin{subarray}{c}\sigma\in\mathcal{S}_{D}:\\\
S^{\prime}\subset\sigma(S_{0})\end{subarray}}\mathbb{P}_{t}(\sigma).$ (57)

But the condition $S^{\prime}\subset\sigma(S_{0})$ is verified if and only if
$S^{\prime}=\sigma(S_{0}^{\prime})$ for some subset $S_{0}^{\prime}\subset
S_{0}$. We therefore obtain

$\displaystyle\Phi_{m\leftarrow K}\left[p_{t,S_{0}}^{(K)}\right](S^{\prime})$
$\displaystyle=\binom{K}{m}^{-1}\sum_{S_{0}^{\prime}\subset
S_{0}}\sum_{\begin{subarray}{c}\sigma\in\mathcal{S}_{D}:\\\
S^{\prime}=\sigma(S_{0}^{\prime})\end{subarray}}\mathbb{P}_{t}(\sigma)$
$\displaystyle=\binom{K}{m}^{-1}\sum_{S_{0}^{\prime}\subset
S_{0}}p^{(m)}_{t,S_{0}^{\prime}}(S^{\prime})$ (58)

where we have recognized the definition of
$p_{t,S_{0}^{\prime}}^{(m)}(S^{\prime})$, the Markov chain induced on
$\Sigma_{m}$, with initial state $S_{0}^{\prime}$ averaged uniformly over all
subsets of $S_{0}$ of the appopriate cardinality, as in Eq. (56). ∎

######  Fact B.4 (monotonicity).

The total variation distances

$\Delta_{m}\equiv\|\Phi_{m\leftarrow K}[p]-\pi_{m}\|_{\rm tv}$ (59)

are non-increasing in $m$.

###### Proof.

We have

$\displaystyle\Delta_{m-1}$ $\displaystyle=\|\Phi_{m-1\leftarrow
K}[p-\pi_{K}]\|_{\rm tv}$ $\displaystyle=\|\Phi_{m-1\leftarrow
m}\circ\Phi_{m\leftarrow K}[p-\pi_{K}]\|_{\rm tv}$
$\displaystyle\leq\|\Phi_{m\leftarrow K}[p-\pi_{K}]\|_{\rm tv}=\Delta_{m}.$
(60)

The first line uses the fact that $\Phi_{m\leftarrow K}[\pi_{K}]=\pi_{m}$ for
all $m<K$; the second uses the semigroup law of $\Phi$ maps, Fact B.2; the
last one uses monotonicity of total variation norm under stochastic maps:
$\|\Phi[\cdots]\|_{\rm tv}\leq\|\cdots\|_{\rm tv}$. ∎

##  Appendix C Proof of Proposition 2

We begin by recalling a useful fact proven in Ref. Giurgica-Tiron and Bouland
(2023):

######  Fact C.1.

For $\omega({\rm poly}(N))<K<o(2^{N})$, we have

$\left\|\sum_{S\in\Sigma_{K}}\pi_{K}(S)\outerproduct{\psi_{S}}{\psi_{S}}^{\otimes
m}-\rho_{\rm Haar}^{(m)}\right\|_{\rm tr}\leq
O\left(\frac{mK}{D},\frac{m^{2}}{K}\right).$ (61)

###### Proof.

See Theorem 1 in Ref. Giurgica-Tiron and Bouland (2023). ∎

Using this fact, in the relevant range of $K$, we can recast the LHS of
Proposition 2 in the more convenient form

$\left\|\sum_{S}[p(S)-\pi_{K}(S)]\outerproduct{\psi_{S}}{\psi_{S}}^{\otimes
m}\right\|_{{\rm tr}}$ (62)

up to small error $O(mK/D,m^{2}/K)$.

Next, using the type-state basis introduced in Appendix A, Eq. (46), we can
express the tensor power $\ket{\psi_{S}}^{\otimes m}$ as

$\ket{\psi_{S}}^{\otimes
m}=\frac{1}{K^{m/2}}\sum_{\mathbf{T}\in\mathcal{T}^{m}_{S}}\binom{m}{\mathbf{T}}^{1/2}\ket{\mathbf{T}}.$
(63)

Here we use $\mathcal{T}_{S}^{m}$ to represent all types of $m$-tuples that
have nonzero components only for states in $S$. We also define pure unique-
type states by

$\displaystyle\ket{\psi_{S,{\rm
uni}}^{(m)}}=\binom{K}{m}^{-1/2}\sum_{{S^{\prime}\in\Sigma_{m}:\,S^{\prime}\subset
S}}\ket{S^{\prime}},$ (64)

where $\ket{S^{\prime}}$ represents the unique-type basis state
$\ket{\mathbf{T}}$ with 1’s only on states in $S^{\prime}$. As also noted in
Appendix A, in the limit of $D\gg m$ most types are unique. This is helpful as
it allows us to replace the state $\ket{\psi_{S}}$ in Eq. (62) by its unique-
type counterpart, $\ket{\psi_{S,{\rm uni}}}$, with a small error. To show
this, we first need a simple fact about the trace distance between pure
states:

######  Fact C.2.

For any two pure states $\ket{\phi}$, $\ket{\chi}$ we have

$\|\outerproduct{\phi}{\phi}-\outerproduct{\chi}{\chi}\|_{\rm
tr}=2\sqrt{1-|\innerproduct{\phi}{\chi}|^{2}}.$ (65)

###### Proof.

See Fact 2 in Ref. Aaronson  _et al._ (2023) (set
$\rho_{1}=\outerproduct{\phi}{\phi}$, $\rho_{2}=\outerproduct{\chi}{\chi}$). ∎

We can now show the following:

######  Lemma C.1 (unique types suffice).

We have

$\left\|\outerproduct{\psi_{S}}{\psi_{S}}^{\otimes
m}-\outerproduct{\psi_{S,{\rm uni}}^{(m)}}{\psi_{S,{\rm
uni}}^{(m)}}\right\|_{{\rm tr}}\leq O\left(\frac{m}{\sqrt{K}}\right).$ (66)

###### Proof.

Using Fact C.2, the LHS reads

$\displaystyle 2\sqrt{1-\left|\innerproduct{\psi_{S}^{\otimes m}}{\psi_{S,{\rm
uni}}^{(m)}}\right|^{2}}.$

Writing both $\ket{\psi_{S}^{\otimes m}}$ and $\ket{\psi_{S,{\rm uni}}^{(m)}}$
in their type basis representation, we get the inner product

$\displaystyle\innerproduct{\psi_{S}^{\otimes m}}{\psi_{S,{\rm uni}}^{(m)}}$
$\displaystyle=K^{-m/2}\binom{K}{m}^{-1/2}$
$\displaystyle\qquad\times\sum_{\mathbf{T}\in\mathcal{T}_{S}^{m}}\binom{m}{\mathbf{T}}^{1/2}\delta(\mathbf{T}\text{
is unique})$ $\displaystyle=(m!)^{1/2}K^{-m/2}\binom{K}{m}^{1/2}.$ (67)

Using some straightforward algebra and the ‘birthday asymptotics’
$\frac{K!}{(K-m)!}=K^{m}\left(1+O(m^{2}/K)\right)$ (see e.g. Fact 1 in Ref.
Giurgica-Tiron and Bouland (2023)) we finally get

$\displaystyle\sqrt{1-\left|\innerproduct{\psi_{S}^{\otimes m}}{\psi_{S,{\rm
uni}}^{(m)}}\right|^{2}}$ $\displaystyle=O\left(\frac{m}{\sqrt{K}}\right).$
(68)

∎

Putting it all together, we can rewrite Eq. (62) as

$\displaystyle\left\|\sum_{S\in\Sigma_{K}}(p(S)-\pi_{K}(S))\outerproduct{\psi_{S,{\rm
uni}}^{(m)}}{\psi_{S,{\rm uni}}^{(m)}}\right\|_{\rm tr}$ (69)

now up to error $O(m/\sqrt{K})$. Finally, plugging in the definition of
unique-type pure states, Eq. (64), yields

$\displaystyle\left\|\sum_{S\in\Sigma_{K}}(p(S)-\pi_{K}(S))\binom{K}{m}^{-1}\sum_{\begin{subarray}{c}S^{\prime},S^{\prime\prime}\in\Sigma_{m}:\\\
S^{\prime},S^{\prime\prime}\subset
S\end{subarray}}\outerproduct{S^{\prime}}{S^{\prime\prime}}\right\|_{\rm tr}$
$\displaystyle=$
$\displaystyle\left\|\sum_{S^{\prime},S^{\prime\prime}\in\Sigma_{m}}\outerproduct{S^{\prime}}{S^{\prime\prime}}\left.\binom{K}{m}^{-1}\sum_{\begin{subarray}{c}S\in\Sigma_{K}:\\\
S^{\prime}\cup S^{\prime\prime}\subset
S\end{subarray}}(p(S)-\pi_{K}(S))\right.\right\|_{\rm tr}.$ (70)

We can recognize the coefficient of
$\outerproduct{S^{\prime}}{S^{\prime\prime}}$ as

$\binom{K}{m}^{-1}\binom{K}{m+\delta}\Phi_{m+\delta\leftarrow
K}[p-\pi](S^{\prime}\cup S^{\prime\prime})$ (71)

which is the matrix $\mathcal{M}_{p}^{(m)}$ in Proposition 2.

##  Appendix D Random walk on the hypercube

###  D.1 Single particle

In this Appendix we review, for completeness, the solution of the simple
random walk on the hypercube $\mathbb{Z}_{2}^{N}$ Diaconis  _et al._ (1990).
We label each vertex of the hypercube by a bitstring
$\mathbf{z}\in\\{0,1\\}^{N}$; two vertices $\mathbf{x},\mathbf{y}$ are
neighbors iff the corresponding bitstrings have Hamming distance 1:
$|\mathbf{x}\oplus\mathbf{y}^{\prime}|=1$ ($\oplus$ is sum modulo 2). In each
time step the walker hops randomly in any one of the $N$ directions,
represented by the Markov chain transition matrix

$\Gamma_{\mathbf{x},\mathbf{y}}=\frac{1}{N}\delta_{|\mathbf{x}\oplus\mathbf{y}|,1}.$
(72)

Note that as written, the walk is not aperiodic: the parity of the Hamming
weight alternates in each time step (i.e., the hypercube is bipartite). This
problem can be eliminated by allowing idling, i.e. replacing
$\Gamma_{\mathbf{x},y}\mapsto
p\Gamma_{\mathbf{x},y}+(1-p)\delta_{\mathbf{x},\mathbf{y}}$ (the walker takes
a step with probability $p$, remains in place otherwise). This does not change
the eigenfunctions and changes the eigenvalues in a trivial way,
$\lambda\mapsto p\lambda+(1-p)$.

The transition matrix in Eq. (72) is diagonalized via the Hadamard transform,
$H_{\mathbf{a},\mathbf{x}}\equiv\frac{1}{\sqrt{D}}(-1)^{\mathbf{a}\cdot\mathbf{x}}$:

$\displaystyle\tilde{\Gamma}_{\mathbf{a},\mathbf{b}}$
$\displaystyle\equiv\frac{1}{D}\sum_{\mathbf{x},\mathbf{y}}\Gamma_{\mathbf{x},\mathbf{y}}(-1)^{\mathbf{a}\cdot\mathbf{x}+\mathbf{b}\cdot\mathbf{y}}$
$\displaystyle=\frac{1}{ND}\sum_{\mathbf{x},\mathbf{r}}\delta_{|\mathbf{r}|-1}(-1)^{(\mathbf{a}\oplus\mathbf{b})\cdot\mathbf{x}+\mathbf{b}\cdot\mathbf{r}}$
$\displaystyle=\delta_{\mathbf{a},\mathbf{b}}\frac{1}{N}\sum_{i=1}^{N}(-1)^{b_{i}}\equiv\delta_{\mathbf{a},\mathbf{b}}\lambda_{\mathbf{a}}.$
(73)

Here we have introduced $\mathbf{r}\equiv\mathbf{x}\oplus\mathbf{y}$ and used
$\sum_{\mathbf{x}}(-1)^{\mathbf{x}\cdot\mathbf{c}}=D\delta_{\mathbf{c},\boldsymbol{0}}$.
The eigenvalue associated to “wave vector” $\mathbf{a}\in\\{0,1\\}^{N}$ is

$\lambda_{\mathbf{a}}=\frac{1}{N}\sum_{i=1}^{N}(1-2a_{i})=1-2\frac{|\mathbf{a}|}{N},$
(74)

which is only a function of Hamming weight $|\mathbf{a}|$.

There is a non-degenerate $\lambda_{0}=1$ eigenvalue, associated to
$\mathbf{a}=\boldsymbol{0}$ (the uniform state), then a gap of magnitude
$2/N$, and then an $N$-fold degenerate eigenvalue $\lambda_{1}=1-2/N$. The
relaxation time is therefore $t_{\rm rel}=1/(1-\lambda_{1})=N/2$. One can also
identify the mixing time exactly and prove that the problem (with idling)
features a cutoff Diaconis  _et al._ (1990). A simple way to estimate the
mixing time is to require that $N\lambda_{1}^{t}$ be small, which gives
$t_{\rm mix}\sim N\log(N)$ (the correct scaling up to constant factors).

The gate set of Sec. IV.3 induces a simple random walk on the hypercube with
idling parameter $p=1/4$, giving rescaled eigenvalues
$\lambda_{1}=1-\frac{2}{N}\mapsto\frac{1}{4}\lambda_{1}+\frac{3}{4}=1-\frac{1}{2N}$,
thus the relaxation time $t_{\rm rel}=2N$.

###  D.2 Many particles: multipole eigenmodes

The eigenvalues $\lambda_{\mathbf{a}}$, Eq. (74), also appear in the spectra
of the Markov chains with $m>1$. We can construct the corresponding eigenmodes
directly. Let us introduce the functions

$\displaystyle f_{\mathbf{a}}(S)$
$\displaystyle=\frac{n_{S}(Z_{\mathbf{a}}=+1)-n_{S}(Z_{\mathbf{a}}=-1)}{m},$
(75)

where for each bitstring $\mathbf{a}\in\mathbb{Z}_{2}^{N}$ we defined
$Z_{\mathbf{a}}=\bigotimes_{i=1}^{N}Z_{i}^{a_{i}}$, and
$n_{S}(Z_{\mathbf{a}}=\pm 1)$ is the number of elements $\mathbf{z}\in S$ such
that $\bra{\mathbf{z}}Z_{\mathbf{a}}\ket{\mathbf{z}}=\pm 1$. The action of the
transition matrix $\Gamma_{S,S^{\prime}}$ on these functions is:

$\displaystyle\sum_{S^{\prime}}\Gamma_{S,S^{\prime}}f_{\mathbf{a}}(S^{\prime})$
$\displaystyle=\sum_{S^{\prime}}\frac{|\\{u\in\mathcal{G}:uS=S^{\prime}\\}|}{|\mathcal{G}|}f_{\mathbf{a}}(S^{\prime})$
$\displaystyle=\sum_{u\in\mathcal{G}}\frac{1}{|\mathcal{G}|}f_{\mathbf{a}}(uS)$
(76)

Here we have used the definition of the transition probability
$\Gamma_{S,S^{\prime}}$ as the fraction of gates $u$ in the gate set
$\mathcal{G}$ that maps subset $S$ to $S^{\prime}$, and then changed variables
via $S^{\prime}\equiv uS$. Summing over $u$ over-counts each subset
$S^{\prime}$ in the original sum by a factor of
$|\\{u^{\prime}\in\mathcal{G}:u^{\prime}S=uS\\}|$, which precisely cancels the
numerator in $\Gamma_{S,S^{\prime}}$. Now we note that

$f_{\mathbf{a}}(uS)\propto
1-2n_{uS}(Z_{\mathbf{a}}=-1)=1-2n_{S}(uZ_{\mathbf{a}}u=-1).$ (77)

If $u$ is the $\textsf{CC}X$ gate that flips qubit $i$ controlled by its
neighbors $i\pm 1$ being in states $c_{\pm}\in\\{-1,1\\}$, we have

$uZ_{i}u=\left[I-2\frac{I+c_{+}Z_{i+1}}{2}\frac{I+c_{-}Z_{i-1}}{2}\right]Z_{i}$
(78)

and $uZ_{j}u=Z_{j}$ for all $j\neq i$. The operator $uZ_{i}u$ has eigenvalues
$\pm 1$ in equal numbers. The $+1$ eigenspace is specified by $Z_{i}=-1$ and
$c_{-}Z_{i-1}=c_{+}Z_{i+1}=+1$ ($D/8$ states), or $Z_{i}=+1$ and
$c_{\xi}Z_{i+\xi}=-1$ for at least one $\xi\in\\{-1,+1\\}$ ($3D/8$ states).
The remaining $D/2$ computational basis states form the $-1$ eigenspace. It is
easy to see that each computational basis state with $Z_{i}=+1$ is in the $+1$
eigenspace of $uZ_{i}u$ for 3 choices of control bits, and in the $-1$
eigenspace for one choice ($c_{+}=c_{-}=+1$). Therefore, when we sum over the
control bits $c_{\pm}\in\\{0,1\\}$ (as part of our sum over $u\in\mathcal{G}$)
we find
$\sum_{c_{\pm}}n_{S}(u_{i,c_{\pm}}Z_{\mathbf{a}}u_{i,c_{\pm}}=+1)=(4-2\delta_{a_{i},1})n_{S}(Z_{i}=+1)$.
On net, we find

$\displaystyle\frac{1}{4N}\sum_{u\in\mathcal{G}}f_{\mathbf{a}}(uS)$
$\displaystyle=\frac{1}{4N}\sum_{i=1}^{N}\sum_{c_{\pm}=\pm
1}f_{\mathbf{a}}(u_{i,c_{-},c_{+}}S)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{4-2\delta_{a_{i},1}}{4}f_{\mathbf{a}}(S)$
$\displaystyle=\left(1-\frac{|\mathbf{a}|}{2N}\right)f_{\mathbf{a}}(S)=\lambda_{\mathbf{a}}f_{\mathbf{a}}(S),$
(79)

where we recognize the eigenvalue $\lambda_{\mathbf{a}}$ of the simple random
walk (with $3/4$ idling probability).

These eigenmodes for the $m$-particle problem have an intuitive
interpretation. $|\mathbf{a}|=0$ represents a monopole moment (i.e. total
particle number), which is conserved; $|\mathbf{a}|=1$ yields the $N$
components of the dipole moment, or center of mass position, which are the
slowest decaying modes in the problem; $|\mathbf{a}|=2$ yields the
$\binom{N}{2}$ components of the quadrupole moment, which decay more quickly;
and so on for higher multipole moments.

For $m=1$ the $2^{N}$ $\\{f_{\mathbf{a}}\\}$ modes exhaust the spectrum,
whereas for $m\geq 2$ they make up only a small part. Still, they are very
important, since they describe the decay of circuit-averaged expectation
values $\overline{\langle Z_{\mathbf{a}}\rangle}$. Indeed, we have

$\displaystyle\overline{\langle Z_{\mathbf{a}}(t)\rangle}$
$\displaystyle=\sum_{S}p_{t}(S)\bra{\psi_{S}}Z_{\mathbf{a}}\ket{\psi_{S}}=(p_{t}|f_{\mathbf{a}})$
$\displaystyle=(p_{0}|\Gamma^{t}|f_{\mathbf{a}})=\left(1-\frac{|\mathbf{a}|}{2N}\right)^{t}(p_{0}|f_{\mathbf{a}}),$
(80)

where we used the definition
$\bra{\psi_{S}}Z_{\mathbf{a}}\ket{\psi_{S}}=f_{\mathbf{a}}(S)$, introduced the
inner product $(a|b)=\sum_{S}a(S)b(S)$ between real functions on $\Sigma_{m}$,
and used the Hermiticity of the transition matrix $\Gamma$ (reversibility of
the Markov chain) to express the time-evolved probability distribution
$(p_{t}|=(p_{0}|\Gamma^{t}$.

It follows that any eigenmode $h(S)$ of the $m$-particle random walk whose
eigenvalue is not in the single-particle spectrum does not show up in circuit-
averaged expectation values $\overline{\langle Z_{\mathbf{a}}\rangle}$:

$\displaystyle\sum_{S}h(S)\bra{\psi_{S}}Z_{\mathbf{a}}\ket{\psi_{S}}=\sum_{S}h(S)f_{\mathbf{a}}(S)=(h|f_{\mathbf{a}})=0,$
(81)

due to orthogonality of the eigenmodes (Hermiticity of $\Gamma$). Such
eigenmodes instead can appear in correlation functions across circuit
realizations, e.g. $\overline{\langle Z_{\mathbf{a}}\rangle\langle
Z_{\mathbf{b}}\rangle}$.

###  D.3 Two particles: relative coordinate eigenmodes

Here we provide more details on eigenmodes of the two-particle problem that
depend subset $S=\\{\mathbf{z},\mathbf{z}^{\prime}\\}$ only through the
relative coordinate $\mathbf{r}=\mathbf{z}\oplus\mathbf{z}^{\prime}$. Let us
consider how $\mathbf{r}$ might change under the action of our local automaton
gate set $\mathcal{G}=\\{u_{iab}\\}$. There are two cases to consider:

  * •

If $r_{i-1}=r_{i+1}=0$ (i.e., $z_{i-1}=z^{\prime}_{i-1}$ and
$z_{i+1}=z^{\prime}_{i+1}$), then either both $z_{i}$ and $z_{i}^{\prime}$
flip, or neither one flips; either way, $\mathbf{r}$ is unchanged.

  * •

Otherwise, there is a $1/4$ probability (based on randomly drawn control bits
$a,b$) that $u_{iab}$ flips $z_{i}$ but not $z_{i}^{\prime}$, a $1/4$
probability that it flips $z_{i}^{\prime}$ but not $z_{i}$, and a $1/2$
probability that it flips neither. In all, there is a $1/2$ probability that
$r_{i}$ flips.

The rules above describe a random walk on $\mathbb{Z}_{2}^{N}$ based on
“activated” bit flips, where $r_{i}$ can flip only if one of its neighbors
$r_{i\pm 1}$ is 1. Clearly $\mathbf{r}=\boldsymbol{0}$ is a steady state, but
it does not describe a valid subset (it gives
$\mathbf{z}=\mathbf{z}^{\prime}$). One can check that the restriction to
$\mathbb{Z}_{2}^{N}\setminus\\{\boldsymbol{0}\\}$ is irreducible, aperiodic
(since lazy moves are allowed) and reversible, thus it has the uniform
distribution as unique steady state.

One can write down the transition matrix $\Gamma|_{R}$ for the relative
coordinate as a quantum Hamiltonian:

$\displaystyle\Gamma|_{R}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left[\left(I-\Pi_{i}\right)\frac{I+X_{i}}{2}+\Pi_{i}\right]\equiv
I-\frac{H}{N},$ (82) $\displaystyle\Pi_{i}$
$\displaystyle=\outerproduct{0}{0}_{i-1}\otimes\outerproduct{0}{0}_{i+1}\otimes
I_{\text{not }i\pm 1}.$ (83)

The linear-in-$N$ relaxation time corresponds to the Hamiltonian $\mathcal{H}$
being gapped.

##  Appendix E Properties of the local automaton gate set

In this Appendix we give more details on the local automaton gate set
$\mathcal{G}$ used in Sec. IV.3. We recall that the gate set is defined as
$\mathcal{G}=\\{u_{iab}=\textsf{C}_{i-1,a}\textsf{C}_{i+1,b}{X}_{i}\\}$, with
$i\in[N]$, $a,b\in\\{0,1\\}$, and periodic boundary conditions (i.e., $i\pm 1$
are taken modulo $N$). Each time step consists of drawing a gate
$u_{iab}\in\mathcal{G}$ uniformly at random and applying it to the system.
This defines a random walk on the permutation group $\mathcal{S}_{D}$ with
transition matrix

$\Gamma_{\sigma,\sigma^{\prime}}=\frac{1}{|\mathcal{G}|}\sum_{u\in\mathcal{G}}\delta_{\sigma^{\prime},u\circ\sigma},$
(84)

where with slight abuse of notation we treat the gate $u$ as a permutation of
the computational basis. Through the natural action of permutations on
subsets,

$S\xrightarrow{\sigma}\\{\sigma(i):\,i\in S\\}\equiv\sigma(S),$

this also defines a Markov chain on the subset spaces $\Sigma_{m}$ for each
$m$, with transition matrices

$\Gamma_{S,S^{\prime}}=\frac{1}{|\mathcal{G}|}\sum_{u\in\mathcal{G}}\delta_{S^{\prime},u(S)}$
(85)

(we use $\Gamma$ for all these transition matrices and rely on the indices to
distinguish them).

To show that $\Gamma_{S,S^{\prime}}$ equilibrates to the unique steady state
$\pi_{m}$ (the uniform distribution on $\Sigma_{m}$), we need the following
three facts to hold:

######  Lemma E.1 (reversibility).

The Markov chain $\Gamma_{S,S^{\prime}}$ on $\Sigma_{m}$ is reversible for all
$m$.

###### Proof.

Since the gate set $\mathcal{G}$ is such that $u^{2}=I$ for all
$u\in\mathcal{G}$, we have $u(S)=S^{\prime}$ if and only if $u(S^{\prime})=S$.
Thus $\Gamma_{S,S^{\prime}}=\Gamma_{S^{\prime},S}$ for all
$S,S^{\prime}\in\Sigma_{m}$. ∎

######  Lemma E.2 (irreducibility).

The Markov chain $\Gamma_{S,S^{\prime}}$ on $\Sigma_{m}$ is irreducible for
all $m\leq D-2$.

###### Proof.

It is well known Shende  _et al._ (2003); Aaronson  _et al._ (2015) that the
Toffoli ($\textsf{CC}X$) and NOT ($X$) gates are “almost” universal for
classical computation, in the sense that they generate the subgroup of
$\mathcal{S}_{D}$ made of even permutations (permutations that can be
decomposed into an even number of transpositions). Our gate set includes
Toffoli gates between consecutive qubits ($u_{i11}$); further, it can generate
NOT gates ($X_{i}=u_{i00}u_{i01}u_{i10}u_{i11}$) and CNOT gates between
neighboring qubits, with either choice of control and target:

$\textsf{C}_{i-1,a}X_{i}=u_{ia0}u_{ia1},\qquad\textsf{C}_{i+1,a}X_{i}=u_{i0a}u_{i1a}.$

Witht the local CNOTs we can make a local SWAP gate:

$\textsf{SWAP}_{i,i+1}=\textsf{C}_{i,0}{X}_{i+1}\circ\textsf{C}_{i+1,0}{X}_{i}\circ\textsf{C}_{i,0}{X}_{i+1}.$

Finally, by composing local SWAPs we can make arbitrary permutations of the
qubits, and thus generate the whole (non-local) Toffoli + NOT gate set.
Therefore $\mathcal{G}$ generates all even permutations. Now, given two sets
$S,S^{\prime}\in\Sigma_{m}$, we can find a permutation
$\sigma\in\mathcal{S}_{D}$ that maps one to the other: $S^{\prime}=\sigma(S)$.
If $\sigma$ is even, we are done. If $\sigma$ is odd, then let $i,j\in[D]$ be
two basis states that are not in subset $S$ (we have assumed $|S|=m\leq D-2$
so that two such states exist), and let $\tau_{ij}\in\mathcal{S}_{D}$ be the
transposition $i\leftrightarrow j$. Then $u\circ\tau_{ij}$ is an even
permutation that maps $S$ to $S^{\prime}$. Since any two sets can be connected
by an even permutation and any even permutation can be decomposed into a
product of $u\in\mathcal{G}$, the Markov chain can connect any two subsets
$S,S^{\prime}\in\Sigma_{m}$ and is thus irreducible. ∎

######  Lemma E.3 (aperiodicity).