# Solvable model of deep thermalization with distinct design times

Matteo Ippoliti  Department of Physics, Stanford University, Stanford, CA
94305, USA<EMAIL_ADDRESS>Wen Wei Ho  Department of Physics, Stanford
University, Stanford, CA 94305, USA  Department of Physics, National
University of Singapore, Singapore 117542<EMAIL_ADDRESS>
###### Abstract

We study the emergence over time of a universal, uniform distribution of
quantum states supported on a finite subsystem, induced by projectively
measuring the rest of the system. Dubbed deep thermalization, this phenomenon
represents a form of equilibration in quantum many-body systems stronger than
regular thermalization, which only constrains the ensemble-averaged values of
observables. While there exist quantum circuit models of dynamics in one
dimension where this phenomenon can be shown to arise exactly, these are
special in that deep thermalization occurs at precisely the same time as
regular thermalization. Here, we present an exactly-solvable model of chaotic
dynamics where the two processes can be shown to occur over different time
scales. The model is composed of a finite subsystem coupled to an infinite
random-matrix bath through a small constriction, and highlights the role of
locality and imperfect thermalization in constraining the formation of such
universal wavefunction distributions. We test our analytical predictions
against exact numerical simulations, finding excellent agreement.

##  1 Introduction

The advent of quantum simulators—systems of many quantum particles whose
interactions are finely programmable [1]—is beginning to allow the exploration
of novel, complex many-body phenomena, particularly out of equilibrium [2, 3,
4, 5, 6]. Beyond their advanced control capabilities, one unique aspect of
quantum simulators is the nature of the information they provide on the
underlying quantum state. In many physical realizations (e.g., ultracold atoms
with quantum gas microscopes, individually trapped ions or Rydberg atoms,
superconducting qubits etc.), measurements involve the simultaneous read-out
of the state of individual degrees of freedom, such as the occupation number
of particles at a site in an optical lattice, or whether an atom is in its
ground or excited state. In such cases, we can think of a measurement in a
given run of the experiment as obtaining a classical ‘snapshot’ of the
configuration of the entire system. Running an experiment multiple times, we
then get a collection of snapshots, which constitutes microscopically-resolved
information about the state. This represents a fundamentally new way of
interrogating quantum many-body systems (as compared to more traditional
macroscopic, coarse-grained probes, such as linear response to external
perturbations or expectation values of order parameters). The vast amount of
information that is generated by this type of probe immediately raises a
number of questions. For example, what universal phenomena can we uncover
based on the information content of these bit-strings? Can such information be
used to characterize novel phases of matter, whether in or out of equilibrium?
And what sorts of applications, especially for quantum information science,
can we develop that leverage such large amounts of data?

There have been many recent developments along these fronts. Measuring local
degrees of freedom during an otherwise-unitary many-body evolution was found
to give rise to a novel kind of nonequilibrium, ‘measurement-induced’ phases
and associated phase transitions, driven by a competition between entanglement
creation and destruction [7, 8, 9, 10, 11, 12, 13, 14, 15]. Advanced data
processing methods have been developed which allow for the efficient
extraction of many properties of a quantum system by manipulating only
compact, classical descriptions of its state, dubbed ‘classical shadows’ [16,
17, 18], which are constructed from snapshots in locally-randomized
measurement bases. In the context of digital quantum simulators, the
statistical analysis of classical snapshots is at the heart of the random
circuit sampling problem, an ideal test of quantum computational advantage
with noisy, intermediate-scale quantum hardware, which has been the subject of
intense theoretical and experimental study in recent years [19, 20, 21, 22,
23, 24, 25]. Relatedly, methods based on the statistics of such classical
snapshots have been proposed to quantify the many-body fidelity between two
quantum states, allowing the benchmarking of quantum simulator performance
[26, 27].

In this work, we consider a novel nonequilibrium universality enabled by this
measurement paradigm, which can be understood as a generalization of the
concept of quantum thermalization in isolated quantum many-body systems—the
process by which local subsystems approach thermal equilibrium despite the
global dynamics being unitary and thus reversible [28, 29, 30, 31].
Conventionally, such a phenomenon is explained by the fact that a local
subsystem $A$ typically becomes entangled with its surroundings $B$ (the
‘bath’, which is taken to be large), over time. Ignoring the precise state of
the bath, the state on $A$ is then described by a statistical mixture of pure
quantum states which replicates thermal equilibrium predictions on average
(i.e., expectation values in an appropriate Gibbs ensemble), in accordance
with statistical mechanical principles. However, the measurement capabilities
of quantum simulators, as described above, motivate a refinement of the above
concept: it is not necessary to assume that the bath is inaccessible; instead,
one can study properties of the local subsystem $A$ conditioned upon
particular configurations of the bath $B$, information which is accessible
from global measurement snapshots. In particular, one can imagine obtaining
(classical) information of the state of $B$, which can then be used to index a
pure quantum state on $A$ arising from the collapse of the global wavefunction
after the measurement; the collection of all such quantum states is called the
projected ensemble [32, 26].

Much of the interest in studying these ensembles lies in the possibility that
they might exhibit a universal limiting form at late times. Indeed, it has
been conjectured [32, 26] that for ensembles built from quantum many-body
states undergoing chaotic quench dynamics, and in the absence of conservation
laws or at very high energy density (corresponding to infinite temperature),
the limiting form should be the unitarily-invariant (i.e., Haar) distribution
over the Hilbert space—i.e., any pure state should be as likely as any other.
One can view this as a generalization of the fundamental statistical
mechanical principle that a local subsystem tends to maximize its entropy
subject to global conservation laws (none, in this case). For the projected
ensemble, the entropy in question is that of the distribution of pure states
over the Hilbert space. This condition can therefore be referred to as deep
thermalization, as opposed to regular thermalization, which only constrains
the behavior of ensemble-averaged observables, i.e., the distribution’s mean.

That such universal distributions emerge was recently rigorously proven in a
large family of models of quantum dynamics in one dimension [33, 34].
Formally, the convergence to the uniform ensemble can be systematically probed
by studying how well each statistical moment of the ensemble replicates that
of the Haar distribution, a condition known in quantum information theory as
forming an approximate quantum state $k$-design [35, 36, 37], where integer
$k\geq 1$ enumerates the statistical moments. It is interesting, then, to ask
about the timescales $t_{k}$ for which a projected ensemble forms a good
$k$-design in quench dynamics (within some fixed accuracy, captured by an
appropriate distance measure), and the physical processes which govern these
time scales. Note that by construction, the design times $t_{k}$ are
monotonically increasing ($t_{k+1}\geq t_{k}$), and that $t_{1}$ is the
regular thermalization time. Ref. [33] found that $t_{k}=t_{1}$ for all $k$ in
a Floquet Ising model tuned to a self-dual point, a result which was later
extended to dual-unitary circuits111These are special circuits that are
unitary not only in time, but also in space—a condition which makes these
models analytically tractable in many cases [38, 39, 40, 41, 42]. under more
general conditions [43, 34]. Further, breaking of dual-unitarity in
$(1+1)$-dimensional unitary circuits was shown to gap $t_{\infty}$ away from
$t_{1}$ in general, via a mapping to measurement-induced dynamical
purification [34], which limits the transfer of information over long spatial
distances in the system. Aside from these results, our understanding of how
$t_{k}$ scales with $k$ and with other parameters of dynamics in general is
still rather limited.

Here, we consider a distinct and complementary mechanism that can gap
$t_{\infty}$ away from $t_{1}$, having to do with local bottlenecks to the
transfer of information in the system. To do so, we will focus on a model of
chaotic quantum dynamics with minimal structure, meant to capture the basic
feature of a spatially-local system: that the subsystem of interest $A$
interacts with the rest of the universe (an unstructured ‘bath’) through the
mediation of another subsystem (informally, its boundary). This minimal
requirement is enough to give rise to an interesting separation between
regular and deep thermalization times, with the model remaining simple enough
to allow for exact analytical results using random-matrix methods.

We find that the projected ensemble on the subsystem of interest $A$, in the
limit of an infinite random-matrix bath, almost always realizes the so-called
Scrooge [44, 26] or Gaussian Adjusted Projected (GAP) [45, 46] ensemble—a
maximally-entropic distribution of pure states consistent with a given density
matrix, which is of independent interest in quantum information theory. We
further compute statistical moments of this distribution and derive an exact
form for the design times $t_{k}$ in our model. As a consequence of this
result we show that (i) $t_{k}$ grows as $\log(k)$ for sufficiently small
moments $k$, and (ii) this growth asymptotes at large $k$ to
$t_{\infty}=2t_{1}$. This represents the second analytically-solved model of
deep thermalization, after the dual-unitary circuit model of Ref. [33], and
the first to exhibit a nontrivial separation between distinct design times
$t_{k}$. We conjecture that the latter phenomenology should also be present in
systems in any dimension as long as they harbor local interactions, and
whenever thermalization is approximate, i.e., when $\rho_{A}$ is close to, but
not exactly equal to, a thermal equilibrium state (dual-unitary circuits are
exceptional in this sense as thermalization there is exact [47]). Thus, our
findings provide a better understanding of the physical phenomenon of deep
thermalization beyond the specific setting of the model studied in this work.

The paper is structured as follows. In Sec. 2 we briefly review key concepts
and introduce our model of dynamics. An analytical solution for the design
times is derived in Sec. 3 and is tested against exact numerical simulations
in Sec. 4. Finally, we summarize our results and discuss directions for future
research in Sec. 5.

##  2 Setup

In this section, we first present a brief review of the object of interest,
the projected ensemble, as well as the quantum information-theoretic formalism
used to characterize it, namely quantum state designs (Sec. 2.1). For a more
detailed review, see Refs. [32, 34]. We then define the specific model of
dynamics considered in this work (Sec. 2.2).

###  2.1 The projected ensemble and quantum state designs

Consider a pure quantum state $\ket{\Psi}_{AB}$ on a bipartite system $AB$. An
ensemble $\mathcal{E}$ of pure states on $A$, called the projected
ensemble222Such an ensemble is also an instance of a geometric quantum state
[48, 49]. [32], can be naturally defined from it via the outcomes of
projective measurements on subsystem $B$:

$\mathcal{E}=\\{(p(z),\ket{\psi_{z}}_{A}):\ z=1,\dots d_{B}\\}.$ (1)

Here $d_{B}$ is the Hilbert space dimension of $B$, $z$ labels an orthonormal
basis333This is typically taken to be the computational basis in a system made
of qubits. We do not make use of this structure in this paper. of $B$, $p(z)$
is the probability of obtaining the outcome $z$ upon measuring $B$, and
$\ket{\psi_{z}}_{A}$ is the post-measurement pure state of subsystem $A$,
conditional on the outcome $z$. Explicitly, these are given by

$p(z)=\|\,{}_{B}\\!\innerproduct{z}{\Psi}_{AB}\|^{2},\qquad\ket{\psi_{z}}_{A}={}_{B}\\!\innerproduct{z}{\Psi}_{AB}/\sqrt{p(z)}.$
(2)

If $\ket{\Psi}_{AB}$ is a typical state in the many-body Hilbert space
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, for sufficiently large $d_{B}$ the
projected ensemble can be shown to acquire a universal form with high
probability, namely its states are distributed according to the Haar (i.e.,
unitarily-invariant) measure on $A$ [32]. The same behavior was conjectured to
emerge at late times in quantum chaotic many-body dynamics from simple initial
states without conservation laws or at infinite temperature, hence
representing a novel kind of nonequilibrium universality. This was supported
by extensive numerical results [32]; later, rigorous exact results were
obtained for dynamics in $(1+1)$-dimensional non-integrable dual-unitary
circuits [33, 43, 34]. We note the projected ensemble is not only a
theoretical construct, but is in fact also accessible in many present-day
experimental platforms of programmable quantum simulators, owing to their
ability to obtain site-resolved measurement data of the global system: by
classically post-processing the data (for example, analyzing the distribution
of bit-strings corresponding to measurement outcomes of $A$, correlated to a
particular measurement outcome of $B$), one can probe the statistics of the
projected ensemble [26].

As mentioned in the introduction, the approach of the distribution of states
in the projected ensemble to the uniform distribution can be characterized in
a systematic fashion: one can ask about the extent to which the projected
ensemble $\mathcal{E}$ reproduces the $k$th moment of the an ensemble of
states drawn from the Haar measure. If the agreement is perfect, $\mathcal{E}$
is said to be a quantum state $k$-design, a concept from quantum information
theory which has significant relevance in applications like randomized
benchmarking and tomography [50, 35, 51, 52, 17, 18]. By allowing for a finite
error $\epsilon>0$ in agreement (quantified by a distance measure, introduced
below), this can be extended to the notion of an approximate quantum state
$k$-design.

Concretely, one can construct “moment operators” from $\mathcal{E}$, which are
density matrices on a $k$-fold replicated Hilbert space
$\mathcal{H}_{A}^{\otimes k}$, defined for all integer values of $k\geq 1$:

$\rho^{(k)}=\sum_{z}p(z)(\ket{\psi_{z}}\bra{\psi_{z}})^{\otimes k}$ (3)

(we omit the subscript $A$ for brevity). It is worth emphasizing that
$\rho^{(k)}\neq\rho^{\otimes k}$ (i.e. the $k$-th moment operator for the
ensemble $\mathcal{E}$ is _not_ equal to $k$ copies of the density matrix
$\rho=\sum_{z}p(z)\outerproduct{\psi_{z}}{\psi_{z}}$), and in fact
$\rho^{(k)}$ is generally not a function of $\rho$ alone; in other words, the
ensemble $\mathcal{E}$ genuinely contains more information than its associated
density matrix $\rho$.

For the particular ensemble of states distributed according to the Haar
measure, the moment operators read

$\rho^{(k)}_{H}=\int{\rm d}\psi\ (\ket{\psi}\bra{\psi})^{\otimes
k}=\binom{k+d_{A}-1}{k}^{-1}{\hat{\Pi}}_{\text{symm}}=\frac{(d_{A}-1)!}{(k+d_{A}-1)!}\sum_{\sigma\in
S_{k}}\hat{\sigma}$ (4)

where $\hat{\sigma}$ is an operator that shuffles the $k$ replicas of the
Hilbert space according to a permutation $\sigma\in S_{k}$, and
${\hat{\Pi}}_{\text{symm}}=\frac{1}{k!}\sum_{\sigma\in S_{k}}\hat{\sigma}$ is
the projector on the permutation-symmetric sector of $\mathcal{H}_{A}^{\otimes
k}$, which has dimension $\binom{k+d_{A}-1}{k}$. The condition of
$\mathcal{E}$ being a quantum state $k$-design is equivalent to the statement
that $\rho^{(k)}=\rho_{H}^{(k)}$.

As the moment operator is a density matrix, we can quantify the distance
between two moment operators $\rho^{(k)}_{1},\rho^{(k)}_{2}$ via standard
distance measures for density matrices, like the trace norm
$\frac{1}{2}\|\rho^{(k)}_{1}-\rho^{(k)}_{2}\|_{1}$. In this work, we will find
it convenient to consider a distance measure based on the Frobenius norm

$\Delta^{(k)}:=\frac{\|\rho^{(k)}-\rho_{H}^{(k)}\|_{2}}{\|\rho_{H}^{(k)}\|_{2}};$
(5)

$\Delta^{(k)}=0$ is then equivalent to $\mathcal{E}$ being a quantum state
$k$-design, while $\Delta^{(k)}\leq\epsilon$ is then equivalent to
$\mathcal{E}$ being an $\epsilon$-approximate quantum state $k$-design. A
closely related quantity is the frame potential of the ensemble, defined from
the purity of the moment operators:

$F^{(k)}:=\|\rho^{(k)}\|_{2}^{2}=\text{Tr}[(\rho^{(k)})^{2}]=\sum_{z_{1},z_{2}}p(z_{1})p(z_{2})\left|\innerproduct{\psi_{z_{1}}}{\psi_{z_{2}}}\right|^{2k}.$
(6)

For the Haar ensemble, the frame potential is
$F_{H}^{(k)}=\binom{k+d_{A}-1}{k}^{-1}$. The distance measure $\Delta^{(k)}$,
Eq. (5), may be rewritten in terms of the frame potentials of the two
ensembles (projected and Haar) as

$\Delta^{(k)}=\left(\frac{F^{(k)}}{F_{H}^{(k)}}-1\right)^{1/2}.$ (7)

This directly implies $F^{(k)}\geq F_{H}^{(k)}$, with equality if and only if
the ensemble forms an exact state $k$-design: $\rho^{(k)}=\rho_{H}^{(k)}$.

Now, consider a quantum many-body state $|\Psi\rangle_{AB}$ which arises in
quench dynamics, beginning from a disentangled state. We can define the time
taken for its projected ensemble to form an $\epsilon$-approximate state
$k$-design, i.e., the minimum $t$ such that $\Delta^{(k)}(t)\leq\epsilon$ for
some arbitrary small threshold $\epsilon$. This defines the design times
$t_{k}$ for each $k$. Due to the fact that $\Delta^{(k+1)}\geq\Delta^{(k)}$ as
shown in Ref. [34], these times are non-decreasing in $k$, $t_{k+1}\geq
t_{k}$. For $k=1$ one recovers the ‘regular thermalization’ time, since the
first moment operator is simply the reduced density matrix of $A$:
$\rho^{(1)}=\text{Tr}_{B}\outerproduct{\Psi}{\Psi}=\rho_{A}$; thus $t_{1}$ is
the time taken to achieve $\rho_{A}\simeq\rho_{H}^{(1)}=\mathbb{I}/d_{A}$,
i.e., thermalization to infinite temperature444Here we implicitly assume
either absence of any conservation laws, or an initial state whose expectation
values of all conserved quantities match the infinite-temperature ones [28,
29].. The design times with $k>1$ may be viewed as probing timescales
associated with stronger forms of thermalization: equilibration not just of
(averaged) expectation values, but also of statistical distributions of the
conditional wavefunctions $|\psi_{z}\rangle$ themselves, which we may call
deep thermalization. In particular, $t_{\infty}\equiv\lim_{k\to\infty}t_{k}$
(which exists due to monotonicity of $k$, though a priori may be infinite) is
the time taken to completely reproduce the Haar measure. Such deep
thermalization is genuinely distinct from (and stronger than) conventional
thermalization: it is fully possible to achieve the latter but not the former.
As a paradigmatic example, a maximally-entangled “EPR pair” state
$\ket{\Psi}_{AB}=\frac{1}{\sqrt{d_{A}}}\sum_{i=1}^{d_{A}}\ket{i}_{A}\otimes\ket{i}_{B}$
has an infinite-temperature reduced density matrix
$\rho_{A}=\frac{1}{d_{A}}\sum_{i=1}^{d_{A}}\outerproduct{i}{i}$, but its
projected ensemble $\mathcal{E}=\\{(p(z)=d_{A}^{-1},\ket{z}_{A}):\ z=1,\dots
d_{A}\\}$ is very far from the Haar measure555It is easy to see that the
ensemble’s frame potential is $d_{A}^{-1}$, independent of $k$. Thus
$\Delta^{(1)}=0$ but e.g. $\Delta^{(2)}=\sqrt{(d_{A}-1)/2}\geq 1$..

###  2.2 Model

A key feature of any geometrically local, short-range interacting system
(which many physical systems are) is that a contiguous subsystem $A$ exchanges
information with the rest of the system through the mediation of other local
subsystems. For example, a ball of radius $\xi$ in a $d$-dimensional lattice,
$\\{\mathbf{r}\in\mathbb{Z}^{d}:\ |\mathbf{r}|<\xi\\}$, is surrounded by a
shell $\\{\mathbf{r}\in\mathbb{Z}^{d}:\ \xi<|\mathbf{r}|<\xi^{\prime}\\}$
($\xi^{\prime}$ being another arbitrary radius $>\xi$) which contains a finite
number of degrees of freedom666In fact, there are infinitely many such shells,
with increasing radii $\xi^{\prime},\xi^{\prime\prime},\dots$, giving a
hierarchy of potential bottlenecks that may prevent or slow down the
transmission of information over long distances; Ref. [34] studies an example
of this mechanism in $(1+1)$-dimensional quantum circuits, where the loss of
information about measurement outcomes over long distances limits the
effective size of the projected ensemble and slows down deep thermalization.
Here instead we consider the effect of a single local bottleneck.; with short-
range interactions, any exchange of information with the outside world
($|\mathbf{r}|>\xi^{\prime}$) must necessarily be mediated by this shell,
which may act as a bottleneck and slow down the dynamics (compared e.g. to a
non-local, all-to-all interacting system).

[width=]fig_schematic.pdf

Figure 1: (a) A schematic of the model: a tripartite system $AB_{1}B_{2}$,
with interactions between $A$ and $B_{1}$ (blue dashed lines) and between
$B_{1}$ and $B_{2}$ (green dashed lines). The Hilbert space dimensions are
also listed. (b) Quantum circuit representing the dynamics of the model: an
initial state $\ket{0}_{AB_{1}B_{2}}$ evolves under unitary gates $V_{T}\cdots
U_{0}V_{0}$ (here $T=3$), then subsystem $B_{1}B_{2}$ is projectively measured
in a fixed basis, yielding an outcome $z$ and a pure state $\ket{\psi_{z}}$ on
$A$. (c) The quantum circuit can be split into two blocks separated by a time-
like cut that intersects the worldline of subsystem $B_{1}$ $2T$ times. The
bond connecting the two blocks has dimension $d_{C}=d_{B_{1}}^{2T}$. One
tensor can be viewed as a state $\ket{\phi_{z}}$ of subsystem $C$ (green box,
right), while the other is a map $E$ from $C$ to $A$ (blue box, left). We
analyze the two parts in Sec. 3.1 and 3.2 respectively.

Here we consider a minimal example of a system featuring this mechanism, with
a single “local bottleneck” mediating interactions between a finite subsystem
of interest and the rest of the system, with the latter harboring nonlocal,
all-to-all interactions. Our aim is to understand the effect of such locality
constraints on the projected ensemble and on the emergence of higher state
designs, i.e. deep thermalization. Concretely, we consider a bipartite system
$AB$, where $B$ is the “bath” to be measured in order to generate the
projected ensemble at $A$. In addition, we partition $B$ into two subsystems
$B_{1}$ and $B_{2}$. $A$, $B_{1}$ and $B_{2}$ have Hilbert spaces of dimension
$d_{A}$, $d_{B_{1}}$ and $q$, respectively. Subsystem $A$ interacts only with
degrees of freedom in $B_{1}$, while $B_{1}$ and $B_{2}$ interact strongly and
without spatial structure. This is shown schematically in Fig. 1(a).

The system is initialized in a disentangled state
$\ket{0}_{A}\ket{0}_{B_{1}}\ket{0}_{B_{2}}$ and subsequently evolves under a
unitary circuit comprised of alternating gates on $AB_{1}$, $U_{t}\in
U(d_{A}d_{B_{1}})$, and on $B_{1}B_{2}$, $V_{t}\in U(d_{B_{1}}q)$. The gates
$U_{t}$, $V_{t}$ are sampled from the Haar measure on the respective unitary
groups (this means there is no additional locality structure inside the
subsystems). We define a time step to consist of a unitary
$(U_{t}\otimes\mathbb{I}_{B_{2}})(\mathbb{I}_{A}\otimes V_{t})$; the system is
evolved for $T$ time steps ($t=0,\dots T-1$) and also by an additional gate
$(\mathbb{I}_{A}\otimes V_{T})$, following which $B$ is measured in an
orthonormal basis $\\{\ket{z}:\ z=1,\dots d_{B_{1}}q\\}$, inducing a projected
ensemble at $A$. This is illustrated in Fig. 1(b).

At this point it is worth remarking on different sources of randomness in the
problem. We wish to characterize the universal distributions on $A$ induced
solely by randomness in the measurement outcomes at $B$ (that is, for a fixed
realization of random gates $\\{U_{t},V_{t}\\}$). However, in our calculations
we will utilize averages over the gates $\\{V_{t}\\}$. Note though that this
is purely a technical tool to simplify calculations: we will find that the
behavior of the (gate-averaged) projected ensemble is typical, allowing us to
make statements about the behavior of the projected ensemble for almost all
circuit realizations.

Finally we remark on the structure of our toy model and how it relates to more
realistic chaotic dynamics. Informally, subsystem $B_{1}$ plays the role of a
local bottleneck that mediates correlations between $A$ and the rest of the
universe, $B_{2}$, which is a large, unstructured random-matrix environment
whose dimension $q$ will be taken to infinity. Clearly, this is a very
favorable scenario for the generation of randomness on $A$: realistic local
models in $d$-dimensional space will generically feature additional
“bottlenecks” at all length scales, which we disregard here by taking a single
bottleneck and an unstructured random environment. For this reason, this
serves as a useful minimal example of the effects of locality on deep
thermalization: it is natural to conjecture that the exact results derivable
in this model may turn into bounds for more realistic and structured models of
dynamics.

##  3 Analytical solution

In this Section we derive the main results of our work: an analytical
computation of the distance measures $\Delta^{(k)}$ for the model of dynamics
defined in Sec. 2.2, leading to exact expressions for all design times and
thus for deep thermalization, in the $q\to\infty$ limit. Our strategy is to
divide the problem in two parts, as sketched in Fig. 1(c): first a treatment
of the environment $B_{2}$ in the large-$q$ limit (Sec. 3.1) and then an
analysis of subsystem $A$ through a replica method (Sec. 3.2). The results are
then unpacked in Sec. 3.4.

###  3.1 Environment in the large-$q$ limit

Here we show that, in the large-$q$ limit, subsystem $B$ (the “environment”
with which subsystem $A$ interacts, and which is measured at the final time)
can be replaced by a Haar-random state in the ‘bond space’, i.e. the temporal
cut777States defined on such temporal cuts in quantum circuits have been
recently studied for various purposes, including simulation of quantum
dynamics [53, 54, 55, 56, 57], realizations of monitored quantum circuits [58,
59, 60], measurement-based quantum computation [61] and quantum process
tomography [62]. $C$ separating $A$ from $B$.

[width=1]fig_env.pdf

Figure 2:  (a) Section of the quantum circuit in Fig. 1(c) representing the
environment. The tensor enclosed in the green dashed box can be viewed as a
state $\ket{\Phi}$ defined on systems $B$ (top legs) and $C$ (left legs). Upon
projecting system $B$ on a basis state $\ket{z}$, one obtains a state
$\ket{\phi_{z}}$ on the bond space $C$, which has dimension
$d_{C}=d_{B_{1}}^{2T}$ ($T=3$ in this example). (b) Tensor network contraction
representing the pseudo-frame potential $F^{(n,k)}$, Eq. (10). $R=2(n+k)$ is
the total number of replicas. Semicircles represent operators acting on $R$
replicas of the Hilbert space, $e$ denotes the identity and $\mathcal{Q}$ is
as in Eq. (11). (c) The same tensor network contraction upon averaging
$V_{0},\dots V_{T}$ over the Haar measure on $U(d_{B_{1}}q)$, to leading order
in $q\to\infty$. Averaging $(V_{t}\otimes V_{t}^{\ast})^{\otimes R}$ yields a
sum over permutations $\sigma_{t}\in S_{R}$, see Eq. (12); the contractions
$(\sigma_{t}|\sigma_{t+1})\propto\delta_{\sigma_{t},\sigma_{t+1}}$ (to leading
order in $q$) force all permutations $\sigma_{t}$ to coincide. (d) Schematic
of the inner product $(\sigma|\mathcal{Q})$. (e) For $z_{1}\neq z_{2}$ (as in
Eq. (14)), the inner product vanishes unless the permutation $\sigma$ takes
the form $\tau\pi_{1}\pi_{2}$, with $\tau$ a transposition of the central $2k$
replicas as shown and $\pi_{1,2}\in S_{n+k}$ arbitrary.

The circuit setup is illustrated in Fig. 2(a), where the environment ($B$,
comprised of $B_{1}$ and $B_{2}$, with dimensions $d_{B_{1}}$ and $q$
respectively, see Sec. 2) is initialized in $\ket{0}$ and measured at the end
to find computational basis state $\ket{z}$ ($z\in\\{1,\dots d_{B_{1}}q\\}$).
This gives a state $\ket{\phi_{z}}$ at the timelike subsystem $C$, which is
made of $2T$ qudits of dimension $d_{B_{1}}$ each ($T$ is the circuit depth as
defined earlier), giving a total dimension $d_{C}=d_{B_{1}}^{2T}$. We may view
this state as belonging to a projected ensemble on $C$ from a bipartite state
$\ket{\Phi}_{BC}$, where $B$ is being measured. Concretely, let us introduce
the unnormalized states
$|\tilde{\phi}_{z}\rangle_{C}\equiv{}_{B}\\!\innerproduct{z}{\Phi}_{BC}$, with
$\ket{\Phi}$ the tensor in Fig. 2(a), which obeys888This is seen by
contracting the corresponding tensor network from the final time backwards,
using unitarity of the $V_{t}$ gates; this produces a factor of $d_{B_{1}}$
per time step, so overall $d_{B_{1}}^{T}=d_{C}^{1/2}$.
$\innerproduct{\Phi}{\Phi}=d_{C}^{1/2}$. Then, the ensemble at $C$ is given by

$\mathcal{E}_{C}=\left\\{\left(p(z)=\frac{\langle\tilde{\phi}_{z}|\tilde{\phi}_{z}\rangle}{d_{C}^{1/2}},\ket{\phi_{z}}=\frac{|\tilde{\phi}_{z}\rangle}{\||\tilde{\phi}_{z}\rangle\|}\right):\
z=1,\dots d_{B_{1}}q\right\\}.$ (8)

Note the normalization of $p(z)$ follows from
$\sum_{z}p(z)=d_{C}^{-1/2}\innerproduct{\Phi}{\Phi}=1$.

In what follows, we are going to show that for $q\to\infty$, $\mathcal{E}_{C}$
approaches an ensemble of states distributed according to the Haar measure on
$\mathcal{H}_{C}$. To do so, we aim to compute the frame potential:

$F^{(k)}=\sum_{z_{1},z_{2}}p(z_{1})p(z_{2})\left|\innerproduct{\phi_{z_{1}}}{\phi_{z_{2}}}\right|^{2k}=\frac{1}{d_{C}}\sum_{z_{1},z_{2}}\innerproduct{\tilde{\phi}_{z_{1}}}{\tilde{\phi}_{z_{1}}}^{1-k}\left|\innerproduct{\tilde{\phi}_{z_{1}}}{\tilde{\phi}_{z_{2}}}\right|^{2k}\innerproduct{\tilde{\phi}_{z_{2}}}{\tilde{\phi}_{z_{2}}}^{1-k}.$
(9)

To make progress, we use a replica technique [43] to define a quantity
$F^{(n,k)}$ which reduces to the frame potential $F^{(k)}$ upon taking the
“replica limit” $n\to 1-k$:

$\displaystyle F^{(n,k)}$
$\displaystyle=\frac{1}{d_{C}}\sum_{z_{1},z_{2}}\innerproduct{\tilde{\phi}_{z_{1}}}{\tilde{\phi}_{z_{1}}}^{n}\left|\innerproduct{\tilde{\phi}_{z_{1}}}{\tilde{\phi}_{z_{2}}}\right|^{2k}\innerproduct{\tilde{\phi}_{z_{2}}}{\tilde{\phi}_{z_{2}}}^{n}=\frac{1}{d_{C}}\Tr\left(\outerproduct{\Phi}{\Phi}^{\otimes
R}_{BC}\mathcal{Q}_{B}\right)$ (10) $\displaystyle\mathcal{Q}_{B}$
$\displaystyle=\sum_{z_{1},z_{2}}\ket{z_{1}^{\otimes n}z_{1}^{\otimes
k}z_{2}^{\otimes k}z_{2}^{\otimes n}}_{B}\bra{z_{1}^{\otimes n}z_{2}^{\otimes
k}z_{1}^{\otimes k}z_{2}^{\otimes n}}_{B}$ (11)

where $R=2(n+k)$ denotes the total number of replicas. The $F^{(n,k)}$
quantity, represented as a tensor network contraction in Fig. 2(b), involves
replicated Haar-random unitaries $(V_{t}\otimes V_{t}^{\ast})^{\otimes R}$.
These tensor products can be averaged exactly [63], and to leading order in
$q\to\infty$ one gets

$\mathbb{E}_{V_{t}\sim\text{Haar}}[(V_{t}\otimes V_{t}^{\ast})^{\otimes
R}]=\sum_{\sigma,\pi\in
S_{R}}\text{Wg}_{d_{B_{1}}q}(\sigma^{-1}\pi,R)|\hat{\sigma})(\hat{\pi}|\simeq\text{Wg}_{d_{B_{1}}q}(e,R)\sum_{\sigma\in
S_{R}}|\hat{\sigma})(\hat{\sigma}|$ (12)

where $\text{Wg}_{d}(\sigma,R)$ is the Weingarten function [63] for $R$
replicas of a $d$-dimensional Hilbert space evaluated on a permutation
$\sigma\in S_{R}$, $e\in S_{R}$ is the identity permutation, and
$\hat{\sigma}$ represents an operator that shuffles the $R$ replicas according
to a permutation $\sigma$: $\hat{\sigma}\ket{i_{1},\dots
i_{R}}=\ket{i_{\sigma^{-1}(1)},\dots i_{\sigma^{-1}(R)}}$. Additionally, here
and in the following we use a standard operator-state mapping notation where
an operator $O=\sum_{a,b}O_{ab}\outerproduct{a}{b}$ defines a state
$|O)=\sum_{ab}O_{ab}\ket{a}\otimes\ket{b}$ on a doubled Hilbert space, with
inner product $(O|O^{\prime})=\text{Tr}(O^{\dagger}O^{\prime})$. Thus e.g.
$|\hat{\sigma})(\hat{\sigma}|$ in Eq. (12) is a super-operator whose action is
defined by $O\mapsto\hat{\sigma}{\rm Tr}(\hat{\sigma}^{\dagger}O)$.

The computation of $F^{(n,k)}$, Eq. (10), features $T+1$ uncorrelated gates
$V_{t}$, all of which are independently averaged according to Eq. (12). We
have that $\text{Wg}_{d_{B_{1}}q}(e,R)=(d_{B_{1}}q)^{-R}$ to leading order in
large $q$, contributing an overall prefactor of $(d_{B_{1}}q)^{-R(T+1)}$. In
addition, the averaging produces $T$ inner products between permutations
(acting on $B_{2}$ only), as sketched in Fig. 2(c). We have
$(\hat{\sigma}|\hat{\pi})=q^{|\sigma^{-1}\pi|}$ where $|\cdots|$ denotes the
number of cycles in a permutation. To leading order in large $q$, this forces
all permutations to be the same: if $\sigma=\pi$ we have
$q^{|\sigma^{-1}\pi|}=q^{|e|}=q^{R}$; otherwise $\sigma^{-1}\pi$ has strictly
fewer than $R$ cycles giving a subleading contribution. Combining all these
inner products yields a factor of $q^{RT}$.

In all, we obtain

$\displaystyle\mathbb{E}[F^{(n,k)}]$ $\displaystyle\simeq
d_{C}^{-1}d_{B_{1}}^{-R(T+1)}q^{-R}\sum_{\sigma\in
S_{R}}\sum_{z_{1},z_{2}=1}^{d_{B_{1}}q}\bra{0^{\otimes
R}}\hat{\sigma}\ket{0^{\otimes
R}}_{B}\text{Tr}_{C}(\hat{\sigma})\text{Tr}_{B}(\hat{\sigma}\mathcal{Q})$ (13)

where the three factors inside the summation arise from contraction of
$\hat{\sigma}$ with the circuit’s boundary conditions at the initial time
($\ket{0}$), final time (operator $\mathcal{Q}$ from Eq. (11), representing
projection on all possible $z_{1}$, $z_{2}$ outcomes on various replicas), and
spatial edge (tracing out $C$) respectively, as illustrated in Fig. 2(c). For
the first factor, we note that (as $\ket{0}^{\otimes R}$ is permutation-
invariant) $\bra{0^{\otimes R}}\hat{\sigma}\ket{0^{\otimes
R}}=\innerproduct{0}{0}^{R}=1$. The second factor is sketched diagrammatically
in Fig. 2(d). The summation over $z_{1,2}$ (see Eq. (11)) breaks up into two
pieces: $d_{B_{1}}q$ terms with $z_{1}=z_{2}$, and
$(d_{B_{1}}q)^{2}-d_{B_{1}}q$ terms with $z_{1}\neq z_{2}$. There is no
further dependence on $q$, so in the $q\to\infty$ limit we can focus on the
latter terms. We get

$\mathbb{E}[F^{(n,k)}]\simeq d_{C}^{-1-R/2}(d_{B_{1}}q)^{2-R}\sum_{\sigma\in
S_{R}}\bra{0^{\otimes n}0^{\otimes k}1^{\otimes k}1^{\otimes
n}}\hat{\sigma}\ket{0^{\otimes n}1^{\otimes k}0^{\otimes k}1^{\otimes
n}}_{B}\text{Tr}_{C}\hat{\sigma}$ (14)

where we have inserted 0 and 1 as arbitrary (orthogonal) computational basis
states on $B$. The term $\bra{0^{\otimes n}1^{\otimes k}1^{\otimes
k}0^{\otimes n}}\hat{\sigma}\ket{0^{\otimes n}1^{\otimes k}0^{\otimes
k}1^{\otimes n}}$ serves to restrict the summation to permutations in the form
$\sigma=\tau\pi_{1}\pi_{2}$, where $\tau$ is a transposition of the central
$2k$ replicas, while $\pi_{1,2}$ act on the first and last $n+k$ replicas,
respectively, as illustrated in Fig. 2(e). One can see that any permutation
not in this form would produce an inner product $\innerproduct{0}{1}=0$.

In conclusion we arrive at

$\displaystyle\mathbb{E}[F^{(n,k)}]$ $\displaystyle\simeq
d_{C}^{-1-R/2}(d_{B_{1}}q)^{2-R}\sum_{\pi_{1,2}\in
S_{n+k}}\text{Tr}[\hat{\tau}\cdot(\hat{\pi}_{1}\otimes\hat{\pi}_{2})]$
$\displaystyle=d_{C}^{-1-R/2}(d_{B_{1}}q)^{2-R}\left[\frac{(R/2+d_{C}-1)!}{(d_{C}-1)!}\right]^{2}\text{Tr}[\hat{\tau}\cdot(\rho_{H,C}^{(n+k)}\otimes\rho_{H,C}^{(n+k)})]$
$\displaystyle=d_{C}^{-1-R/2}(d_{B_{1}}q)^{2-R}\left[\frac{(R/2+d_{C}-1)!}{(d_{C}-1)!}\right]^{2}F_{H}^{(k)}$
(15)

where in the last line we have traced out the $2n$ replicas (on which the
transposition $\hat{\tau}$ does not act) obtaining
$(\rho_{H,C}^{(k)})^{\otimes 2}$, and then noted that
$\text{Tr}(\hat{\tau}A^{\otimes 2})=\text{Tr}(A^{2})$. Finally, having
obtained an analytical function of $R$, we continue it to $R=2(n+k)=2$ (i.e.
$n=1-k$), which gives the desired result:

$\lim_{q\to\infty}\mathbb{E}[F^{(k)}]=F_{H}^{(k)}.$ (16)

Given that $F^{(k)}\geq F_{H}^{(k)}$ in each realization (following Eq. (7)),
we can conclude that $F^{(k)}=F_{H}^{(k)}$ almost always (over the Haar
measure for gates $V_{t}$). This means that statistically speaking, we can
replace the entire subsystem $B$ with a Haar-random state on the temporal cut
$C$ in order to analyze the formation of state designs on $A$.

Some comments are in order here. First, we note that the same result holds in
dual-unitary circuits with suitable initial states and final measurement bases
in 1+1 dimensions [33, 43, 34]: there, too, one obtains the Haar distribution
in the bond space $C$ in the limit of an infinite bath, $N_{B}\to\infty$ where
$N_{B}$ is the number of degrees of freedom in $B$ (the corresponding
thermodynamic limit in our case is $q\to\infty$). In that case, the Haar
distribution emerges essentially from a deep random unitary circuit running in
the space direction; here instead it emerges from measuring part of a given
highly-entangled state $\ket{\Phi}_{BC}$, similarly to the case of a global
Haar-random state $\ket{\Psi}_{AB}$ studied in Ref. [32]. Secondly, we note
that realistic models of dynamics may fail to produce Haar-randomness in the
bond space ($C$ in our case); this is indeed the subject of Ref. [34], where
it was shown that $(1+1)$-dimensional unitary circuits may fail to achieve
this due to dynamical purification in the space direction, i.e., the loss of
information about measurement outcomes over very long distances. Here, by
taking an all-to-all interacting random-matrix environment, we have avoided
any such issues. Therefore this is a favorable scenario for the formation of
state designs at $A$, and we expect that results obtained in this model may
provide bounds for more realistic short-range-interacting models.

###  3.2 Projected ensemble on $A$

With the above result we can greatly simplify the setup, by discarding the
thermodynamically-large environment $B_{2}$ and focusing only on subsystem $A$
and the temporal subsystem (or bond space) $C$, as sketched in Fig. 1(c).
After observing outcome $z$ on $B$, we have a state $\ket{{\phi}_{z}}$ at $C$
which gets mapped to $A$ by a linear map999The map $E$ is obtained by
combining the gates $\\{U_{t}\\}$ as shown diagrammatically in Fig. 1(c).
Formally, one has
$E_{a,\mathbf{b}}=\sum_{\alpha_{0},\dots\alpha_{T}}\delta_{\alpha_{0},0}\delta_{\alpha_{T},a}\prod_{t=0}^{T-1}\bra{\alpha_{t+1},b_{2t+1}}U_{t}\ket{\alpha_{t},b_{2t}}$,
where $a$ and each $\alpha_{t}$ run over states of $A$, each $b_{t}$ runs over
states of $B_{1}$, and the multi-index $\mathbf{b}=(b_{0},\dots b_{2T-1})$
runs over states of $C$. $E$. As we just saw, in the large-$q$ (thermodynamic)
limit the input state $\ket{\phi_{z}}$ is effectively Haar-random on $C$, thus
the question becomes: how well is this randomness propagated from $C$ to $A$?

In the case of dual-unitary circuits in $1+1$ dimensions, the randomness is
propagated optimally: as long as $d_{C}\geq d_{A}$, the Haar measure at $C$
induces the Haar measure at $A$. This follows from the perfect thermalization
of dual-unitary circuits [47]: because $\rho_{A}=\mathbb{I}_{A}/d_{A}$, $E$
must act as an isometry on a subspace of $\mathcal{H}_{C}$ (of the appropriate
dimension $d_{A}$), which is enough to produce the Haar measure on
$\mathcal{H}_{A}$ [33]. However, for generic interactions (i.e. non-dual-
unitary, or in higher dimension, etc), this is not the case: $\rho_{A}$ equals
$\mathbb{I}/d_{A}$ only approximately, and it is not a priori clear what the
induced distribution of states at $A$ looks like. This is the question we
address next.

In terms of the states
$|\tilde{\phi}_{z}\rangle_{C}={}_{B}\\!\innerproduct{z}{\Phi}_{BC}$ studied in
Sec. 3.1, the projected ensemble at $A$ takes the form

$\mathcal{E}=\left\\{\left(p(z)=\langle\tilde{\phi}_{z}|E^{\dagger}E|\tilde{\phi}_{z}\rangle,\
\ket{\psi_{z}}=\frac{E|\tilde{\phi}_{z}\rangle}{\|E|\tilde{\phi}_{z}\rangle\|}\right):\
z=1,\dots d_{B_{1}}q\right\\}$ (17)

It follows that the expectation of a function $f(\ket{\psi})$ on the ensemble
$\mathcal{E}$ is

$\mathbb{E}f=\sum_{z}p(z)f(\ket{\psi_{z}})=\sum_{z}\innerproduct{\tilde{\phi}_{z}}{\tilde{\phi}_{z}}\bra{\phi_{z}}E^{\dagger}E\ket{\phi_{z}}f\left(\frac{E\ket{\phi_{z}}}{\|E\ket{\phi_{z}}\|}\right)$
(18)

where we have introduced the normalized vectors
$\ket{\phi_{z}}=|\tilde{\phi}_{z}\rangle/\||\tilde{\phi}_{z}\rangle\|$. The
result of Sec. 3.1 is that the ensemble in Eq. (8) is Haar-distributed in the
thermodynamic limit $q\to\infty$, i.e.,

$\lim_{q\to\infty}\sum_{z}d_{C}^{-1/2}\innerproduct{\tilde{\phi}_{z}}{\tilde{\phi}_{z}}g(\ket{\phi_{z}})=\int{\rm
d}\phi\ g(\ket{\phi})$ (19)

with ${\rm d}\phi$ the Haar measure and $g$ any sufficiently regular function.
With this substitution, Eq. (18) becomes

$\mathbb{E}f=d_{C}^{1/2}\int{\rm d}\phi\
\bra{\phi}E^{\dagger}E\ket{\phi}f\left(\frac{E\ket{\phi}}{\|E\ket{\phi}\|}\right)$
(20)

which defines the limiting form of the projected ensemble at $A$:

$\lim_{q\to\infty}\mathcal{E}=\left\\{\left({\rm d}\psi={\rm d}\phi\
d_{C}^{1/2}\bra{\phi}E^{\dagger}E\ket{\phi},\
\ket{\psi}=\frac{E\ket{\phi}}{\|E\ket{\phi}\|}\right)\right\\},$ (21)

where again ${\rm d}\phi$ is the Haar measure.

The resulting ensemble, Eq. (21), is known as the Scrooge ensemble [44, 26] or
Gaussian Adjusted Projected (GAP) ensemble [45, 46], which is the maximally
entropic101010More precisely, it minimizes a measure of accessible information
[44], thus it is “maximally stingy” with its information, hence the name.
ensemble of pure states compatible with a given reduced density matrix—in this
case,

$\rho_{A}=\int{\rm d}\psi\ \outerproduct{\psi}{\psi}=\int{\rm d}\phi\
d_{C}^{1/2}E\outerproduct{\phi}{\phi}E^{\dagger}=d_{C}^{-1/2}EE^{\dagger}.$
(22)

As a technical remark, we note that Eq. (21) is seemingly different from the
conventional construction of the Scrooge/GAP ensemble (as presented in [44,
26]), as the former involves objects acting on two Hilbert spaces of different
dimensions ($A$ and $C$) while the latter only one. However, it can be
straightforwardly shown that the two definitions are equivalent, see Appendix
A.

###  3.3 Moments of the Scrooge/GAP ensemble near infinite temperature

To recapitulate, the analysis of Sec. 3.1 and 3.2 has shown that, in the
large-$q$ (thermodynamic) limit, the projected ensemble constructed from our
model of dynamics realizes (almost always) the Scrooge/GAP ensemble, Eq. (21),
induced by the appropriate density matrix, Eq. (22). With this result in hand,
we are now in a position to study the formation of higher state designs at
$A$. However, we remark that the derivation which follows is more general: it
may be applied to the Scrooge/GAP ensemble built from any density matrix
$\rho$ on any system, as long as $\rho$ is very close to thermalization at
infinite temperature. Such an ensemble may potentially be obtained through
other means (i.e., not by the projected ensemble on our model) and is an
object of independent interest.

In order to characterize the emergence of quantum state designs, we aim to
compute the $k$th moment operators for the ensemble in Eq. (21):

$\displaystyle\rho_{A}^{(k)}$ $\displaystyle=\int{\rm
d}\psi(\outerproduct{\psi}{\psi})^{\otimes k}=d_{C}^{1/2}\int{\rm d}\phi\
\|E\ket{\phi}\|^{2}\left(\frac{E\outerproduct{\phi}{\phi}E^{\dagger}}{\|E\ket{\phi}\|^{2}}\right)^{\otimes
k}$ $\displaystyle=d_{C}^{1/2}\int{\rm d}\phi\
\bra{\phi}E^{\dagger}E\ket{\phi}^{1-k}(E\outerproduct{\phi}{\phi}E^{\dagger})^{\otimes
k}.$ (23)

At this point, for $k>1$ the presence of a negative power $1-k$ prevents us
from directly performing the integral over the Haar measure ${\rm d}\phi$.
This issue can be overcome by using a replica trick [43] analogous to the one
in Eq. (10): we define a “pseudo-moment operator”

$\rho^{(n,k)}_{A}=d_{C}^{1/2}\int{\rm d}\phi\
\bra{\phi}E^{\dagger}E\ket{\phi}^{n}(E\ket{\phi}\bra{\phi}E^{\dagger})^{\otimes
k}$ (24)

which reduces to $\rho^{(k)}_{A}$ for $n=1-k$. We take $n>0$ in order to
perform the calculation; at the end we will take the “replica limit” $n\to
1-k$ and obtain the desired result. By using the fact that
$\innerproduct{\phi}{\phi}=\text{Tr}\outerproduct{\phi}{\phi}$, we may re-
write this as

$\displaystyle\rho^{(n,k)}_{A}$ $\displaystyle=d_{C}^{1/2}\Tr_{A_{1},\dots
A_{n}}\int{\rm d}\phi\ (E\ket{\phi}\bra{\phi}E^{\dagger})^{\otimes(n+k)}$ (25)

where we have introduced $n$ new replicas of $A$ in addition to the $k$
replicas that go into the definition of the $k$-th moment operator, and the
trace is taken over these new replicas $A_{1},\dots A_{n}$.

We can now perform the integral over the Haar measure ${\rm d}\phi$ exactly:

$\displaystyle\rho^{(n,k)}_{A}$ $\displaystyle=d_{C}^{1/2}\Tr_{A_{1},\dots
A_{n}}\left[E^{\otimes(n+k)}\rho_{H,C}^{(n+k)}(E^{\dagger})^{\otimes(n+k)}\right]$
$\displaystyle=d_{C}^{1/2}\frac{(d_{C}-1)!}{(d_{A}-1)!}\frac{(n+k+d_{A}-1)!}{(n+k+d_{C}-1)!}\Tr_{A_{1},\dots
A_{n}}\left[\rho_{H,A}^{(n+k)}(EE^{\dagger})^{\otimes(n+k)}\right]$ (26)

where $\rho_{H,A}^{(n+k)}$ and $\rho_{H,C}^{(n+k)}$ are the moment operators
of the Haar ensemble on $A$ and $C$, respectively. In the second line we have
commuted $E^{\otimes(n+k)}$ past $\rho_{H,C}^{(n+k)}$ by using the
decomposition of the Haar moment operators into permutations and the fact that
$E^{\otimes(n+k)}$ is permutation-symmetric: for any integer $m>0$, we have

$\displaystyle E^{\otimes m}\rho_{H,C}^{(m)}$
$\displaystyle=\frac{(d_{C}-1)!}{(m+d_{C}-1)!}\sum_{\sigma\in S_{m}}E^{\otimes
m}\hat{\sigma}_{C}=\frac{(d_{C}-1)!}{(m+d_{C}-1)!}\sum_{\sigma\in
S_{m}}\hat{\sigma}_{A}(\hat{\sigma}_{A}^{-1}E^{\otimes m}\hat{\sigma}_{C})$
$\displaystyle=\frac{(d_{C}-1)!}{(m+d_{C}-1)!}\frac{(m+d_{A}-1)!}{(d_{A}-1)!}\rho_{H,A}^{(m)}E^{\otimes
m}.$ (27)

Next, we use Eq. (22) to eliminate $EE^{\dagger}$ from Eq. (26), obtaining

$\displaystyle\rho^{(n,k)}_{A}$
$\displaystyle=g(n,k,d_{A},d_{C})\Tr_{A_{1},\dots
A_{n}}[\rho_{H,A}^{(n+k)}(d_{A}\rho_{A})^{\otimes(n+k)}]\;,$ (28)

where we have collected the awkward prefactors into a function $g$ which can
be seen to equal 1 for all $d_{A}$, $d_{C}$ when $n+k=1$, and can thus be
safely dropped for our purpose.

The expression Eq. (28) is noteworthy as it only involves the reduced density
matrix $\rho_{A}$ and the moment operator of the Haar ensemble. It is clear
that if $\rho_{A}=\mathbb{I}/d_{A}$ (perfect thermalization), then
$\rho_{A}^{(n,k)}=\text{Tr}_{A_{1},\dots
A_{n}}[\rho_{H,A}^{(n+k)}]=\rho_{H,A}^{(k)}$, independent of $n$, and we have
an exact quantum state $k$-design for all $k$, i.e., the exact Haar measure at
$A$. Conversely, any deviation $\rho_{A}\neq\mathbb{I}/d_{A}$ (imperfect
thermalization) directly maps onto a deviation between the $k$th moment
operators. This deviation is what we aim to quantify next.

Let us assume that $\rho_{A}^{(1)}=(\mathbb{I}+\mu)/d_{A}$, where $\mu$ is a
Hermitian operator that represents the slowest observable in the process of
relaxation towards equilibrium at infinite temperature. By definition, we have

$\Delta^{(1)}=\frac{\|\rho_{A}-\mathbb{I}/d_{A}\|_{2}}{\|\mathbb{I}/d_{A}\|_{2}}=\frac{\|\mu\|_{2}}{\|\mathbb{I}\|_{2}}.$
(29)

At late times, when $\|\mu\|_{2}$ is small, we can expand Eq. (28) in powers
of $\mu$:

$\displaystyle\rho^{(n,k)}_{A}$
$\displaystyle=\rho_{A,H}^{(k)}+\Tr_{A_{1},\dots
A_{n}}\left[\rho_{H,A}^{(n+k)}\sum_{j=1}^{n+k}\mu_{j}\right]+O(\mu^{2}),$ (30)

where $\mu_{j}$ denotes a $\mu$ operator acting on replica $A_{j}$ only (with
$\mathbb{I}$ on all other replicas). Focusing on the linear term, we can split
the sum into two parts:

  * •

Terms with $j>n$: tracing out $A_{1},\dots A_{n}$ directly gives
$\rho_{H,A}^{(k)}\sum_{j=1}^{k}\mu_{j}$;

  * •

Terms with $1\leq j\leq n$: due to permutation symmetry these $n$ terms are
all identical, giving $n$ contributions of
$\text{Tr}_{A_{k+1}}[\rho_{H,A}^{(k+1)}\mu_{k+1}]$.

The latter can be computed as follows. Decomposing $\rho_{H,A}^{(k+1)}$ into
permutations $\sigma\in S_{k+1}$, we find that permutations with
$\sigma(k+1)=k+1$ give a contribution $\propto{\rm Tr}(\mu)=0$ ($\mu$ is
traceless by definition). On the other hand, permutations with
$\sigma(k+1)=j\leq k$ give a contribution $\propto\mu_{j}\hat{\pi}$, where
$\pi\in S_{k}$ is a permutation induced on the remaining $k$
replicas111111There is a one-to-one correspondence between $\pi\in S_{k}$ and
$\sigma\in S_{k+1}$ with the constraint $\sigma(k+1)=j$.. We have

$\displaystyle{\rm Tr}_{A_{k+1}}\rho_{H,A}^{(k+1)}\mu_{k+1}$
$\displaystyle=\frac{(d_{A}-1)!}{(k+d_{A})!}\sum_{j=1}^{k}\sum_{\begin{subarray}{c}\sigma\in
S_{k+1}:\\\ \sigma(k+1)=j\end{subarray}}{\rm
Tr}_{A_{k+1}}\hat{\sigma}\mu_{k+1}$
$\displaystyle=\frac{(d_{A}-1)!}{(k+d_{A})!}\sum_{j=1}^{k}\mu_{j}\sum_{\pi\in
S_{k}}\hat{\pi}=\frac{1}{k+d_{A}}\sum_{j=1}^{k}\mu_{j}\rho_{H,A}^{(k)}.$ (31)

Adding both contributions, we conclude

$\displaystyle\rho_{A}^{(n,k)}-\rho_{H,A}^{(k)}$
$\displaystyle=\left(1+\frac{n}{d_{A}+k}\right)\rho_{H,A}^{(k)}\sum_{j=1}^{k}\mu_{j}+O(\mu^{2}).$
(32)

At this point, all the $n$ auxiliary replicas have been traced out, and the
parameter $n$ appears only as a prefactor. It is possible to take the replica
limit121212Strictly speaking, this assumes that the operations of replica
limit and Taylor expansion to $O(\mu)$ commute. As we will see in Sec. 4, this
assumption is justified by excellent agreement with numerical simulations.
$n\to 1-k$, which gives

$\rho_{A}^{(k)}-\rho_{H,A}^{(k)}=\frac{d_{A}+1}{d_{A}+k}\rho_{H,A}^{(k)}\sum_{j=1}^{k}\mu_{j}+O(\mu^{2}).$
(33)

This is one of the central results of our work, directly connecting imperfect
thermalization (quantified by $\mu$) to the discrepancy between $k$-th moments
of the projected and Haar ensembles.

With this result in hand, we may now evaluate the distance measure Eq. (5):

$\displaystyle(\Delta^{(k)})^{2}$
$\displaystyle=\frac{\text{Tr}[(\rho_{A}^{(k)}-\rho_{A.H}^{(k)})^{2}]}{\text{Tr}[(\rho_{H,A}^{(k)})^{2}]}=\left(\frac{d_{A}+1}{d_{A}+k}\right)^{2}\Tr[\rho_{H,A}^{(k)}\sum_{i,j=1}^{k}\mu_{i}\mu_{j}]$
(34)

where we used the facts that $\sum_{i}\mu_{i}$ commutes with
$\rho_{H,A}^{(k)}$ (being permutation-symmetric) and that
$(\rho_{H,A}^{(k)})^{2}/\text{Tr}[(\rho_{H,A}^{(k)})^{2}]=\rho_{H,A}^{(k)}$ in
order to simplify the result. The summation over replicas $i$, $j$ breaks up
into two parts:

  * •

Terms $i=j$ give $k$ identical contributions of
$\text{Tr}(\rho_{H,A}^{(1)}\mu_{1}^{2})=\|\mu\|_{2}^{2}/d_{A}$.

  * •

Terms $i\neq j$ give $k(k-1)$ identical contributions of
$\text{Tr}(\rho_{H,A}^{(2)}\mu_{1}\mu_{2})$; using
$\rho_{H,A}^{(2)}=\frac{\hat{e}+\hat{\tau}}{d_{A}(d_{A}+1)}$ (where $e,\tau\in
S_{2}$ are the identity permutation and the transposition, respectively), and
the fact that $\mu$ is traceless, we get $\|\mu\|_{2}^{2}/[d_{A}(d_{A}+1)]$.

Putting it all together, and using Eq. (29) to eliminate $\|\mu\|_{2}$, we
arrive at our main result:

$\displaystyle\Delta^{(k)}$
$\displaystyle=f(k,d_{A})\Delta^{(1)}+O((\Delta^{(1)})^{2}),$ (35)
$\displaystyle f(k,d_{A})$
$\displaystyle=\left(\frac{1+d_{A}}{1+d_{A}/k}\right)^{1/2}.$ (36)

###  3.4 Scaling of design times

With the result above, Eq. (35-36), we are now in a position to discuss the
scaling of design times in the projected ensemble for our model. We begin by
observing that, by the definition Eq. (5), we have
$\Delta^{(1)}=\sqrt{2^{N_{A}-S_{2}(A)}-1}$, where
$S_{2}(A)=-\log_{2}\text{Tr}\rho_{A}^{2}$ is the second Renyi entropy of
subsystem $A$, measured in bits. In chaotic dynamics starting from a
disentangled state, we expect the entropy to grow linearly in time [64],
$S_{2}(A)\sim v_{E}t$ (the proportionality constant $v_{E}$ is known as the
entanglement velocity), until near saturation to the maximal value
$\log_{2}(d_{A})=N_{A}$. At that point, achieved when $t\approx
t^{\ast}=N_{A}/v_{E}$, we expect a crossover and an exponential convergence to
the asymptote: $N_{A}-S_{2}(A)\sim 2^{-v_{E}^{\prime}(t-t^{\ast})}$.

In models described by the statistical mechanics of an entanglement membrane
[65, 66, 67, 68, 69], one expects that the two velocities governing the early-
time ballistic growth and the late-time exponential convergence should in fact
be the same, $v_{E}\equiv v_{E}^{\prime}$, related to the surface tension of
the entanglement membrane131313We note however that a “two-step relaxation
process” (corresponding to two different velocities for the early-time linear
growth and the late-time exponential saturation) has been observed in one-
dimensional random unitary circuits, depending nontrivially on circuit
architecture and boundary conditions [70, 71].. Our model also displays this
behavior in the $q\to\infty$ limit: as shown in Appendix B, we have

$S_{2}(A)=\left\\{\begin{aligned} &v_{E}t&(t\ll t^{\ast})\\\
&N_{A}-\frac{1}{\ln 2}2^{-v_{E}(t-t^{\ast})}&(t\gg
t^{\ast})\end{aligned}\right.,\qquad v_{E}=2\log_{2}(d_{B_{1}}),\qquad
t^{\ast}=N_{A}/v_{E}.$ (37)

We note that $\log_{2}(d_{B_{1}})$ plays the role of $|\partial A|$, the size
of the boundary of subsystem $A$ (i.e., the number of qubits with which $A$
interacts directly). In a spatially local system, one indeed expects
$v_{E}\propto|\partial A|$, and $t^{\ast}\propto|A|/|\partial A|\sim\xi$ if
$A$ is a ball of radius $\xi$ in a $d$-dimensional lattice. Motivated by these
considerations, in the following we use $v_{E}$ for the late-time rate of
convergence of $S_{2}$ towards $N_{A}$; however this is a slight abuse of
notation, and one should be mindful that this identification need not work in
general, e.g. in integrable systems, dual-unitary circuits, etc.

As a consequence of Eq. (37) we have, at late times,

$\Delta^{(1)}=\sqrt{2^{N_{A}-S_{2}(A)}-1}\sim 2^{-v_{E}(t-t^{\ast})/2}$ (38)

and thus the thermalization time, defined with respect to an arbitrary
threshold $\epsilon\ll 1$:

$\Delta^{(1)}<\epsilon\quad\iff\quad
t>t_{1}\equiv\frac{N_{A}}{v_{E}}+\frac{2\log_{2}(1/\epsilon)}{v_{E}}.$ (39)

By using our main result Eq. (35-36) we can carry out the same reasoning for
the $k$-th moment, which gives

$\Delta^{(k)}=f(k,d_{A})\Delta^{(1)}<\epsilon\quad\iff\quad t>t_{k}\equiv
t_{1}+\frac{2\log_{2}f(k,d_{A})}{v_{E}}.$ (40)

We have thus obtained the exact design times for this problem, for all $k\geq
1$, in the $q\to\infty$ limit:

$t_{k}=t_{1}+\frac{1}{v_{E}}\log_{2}\frac{1+d_{A}}{1+d_{A}/k}.$ (41)

This is the second exact calculation of design times in a model of quantum
dynamics (after the self-dual kicked Ising model studied in Ref. [33]), and
the first to show a separation between distinct design times. In the following
we comment on the scaling of design times with $k$ and $d_{A}$ in two regimes
of interest.

#### Large $d_{A}$.

If $k\ll d_{A}$, the design times from Eq. (41) approximately reduce to

$t_{k}\approx t_{1}+\frac{1}{v_{E}}\log_{2}k,$ (42)

showing a nontrivial growth with $k$. However this growth is only logarithmic
in $k$ (in other words, the successive $k$-designs are formed extremely
quickly in dynamics, with $k(t)\sim\exp(t)$). It is interesting to ask how
this compares to other models. Numerical results on chaotic Hamiltonian
dynamics at infinite temperature [32], while in a limited dynamic range
($1\leq k\leq 4$), are suggestive of a stronger, possibly algebraic growth
$t_{k}\sim k^{\alpha}$. Such a difference in scaling might originate from
features of dynamics that are absent in our toy model: geometric locality and
energy conservation, both of which could slow down deep thermalization. The
latter in particular is known to parametrically change the scaling of the
thermalization time (due to energy diffusion) and it is interesting to ask how
it might affect deep thermalization.

#### Large $k$.

The logarithmic growth of $t_{k}$ in Eq. (42) eventually arrests altogether
when $k\approx d_{A}$. This is also the regime in which state designs for
different $k$ cease to be fully independent (i.e., for such a $k$, a
$k$-design is automatically an approximate $(k+1)$-design [52]). At large $k$
we have

$t_{\infty}=\lim_{k\to\infty}t_{k}=t_{1}+\frac{\log_{2}(d_{A}+1)}{v_{E}}\simeq\frac{2N_{A}}{v_{E}}+\frac{2\log_{2}(1/\epsilon)}{v_{E}},$
(43)

where the approximation holds for $d_{A}\gg 1$. The time scale $t_{\infty}$
governs the emergence of the Haar measure at $A$ (fulfilling the state
$k$-design condition for all $k$’s). Remarkably, we see that (i) a uniform
distribution of states (to within fixed accuracy $\epsilon$) is generated in
finite time in this model, and (ii) in the limit $N_{A}\to\infty$ (taken after
$k\to\infty$) the emergence of the Haar distribution is governed by a
renormalized entanglement velocity, $v_{\text{Haar}}\equiv v_{E}/2$.
Appropriately, this takes longer than thermalization; precisely a factor of
$2$ longer in this model.

##  4 Numerical simulations

To test our analytical results, we perform exact numerical simulations of the
circuit model. We aim to compute the ratio of design distances
$\Delta^{(k)}/\Delta^{(1)}$ in order to confirm the relationship in Eq.
(35-36); the scaling of design times follows from the behavior of this ratio.

In order to maximize the accessible environment size, we choose subsystems $A$
and $B_{1}$ of minimal size, $d_{A}=d_{B_{1}}=2$. We also set $q=2^{L-2}$ for
integer $L$ so that the total Hilbert space corresponds to $L$ qubits. We
initialize a state $\ket{0}^{\otimes L}$ and evolve it under the circuit in
Fig. 1(b), where the gates $U_{t}$ and $V_{t}$ are sampled from the Haar
measure on $U(4)$ and $U(2^{L-1})$ respectively. After each time step
$T=1,\dots t_{\text{max}}=20$ we construct the exact moment operators
$\rho^{(k)}$ for $k=1,\dots k_{\text{max}}=7$ and evaluate the distance
measures $\Delta^{(k)}(t)$. We repeat this process for many realizations of
the Haar-random gates and average the results. See Appendix C for further
details on the implementation and efficiency of this algorithm.

[width=]fig_numerics.pdf

Figure 3: Results of numerical simulations. (a-c) $k$-design distance measure
$\Delta^{(k)}$ as a function of circuit depth $T$, computed as described in
Sec. 4 (see also Appendix C), averaged over 4000 realizations ($L=12$, 20) and
1000 realizations ($L=28$) of the Haar-random unitary evolution. (d-f) Ratios
$\Delta^{(k)}/\Delta^{(1)}$ at various times $T$ (marked by vertical dashed
lines of the corresponding color in (a-c)), compared to the analytical
prediction $f(k,d_{A}=2)=\sqrt{3k/(k+2)}$ (dotted line), Eq. (36).

Numerical results for the gate-averaged $\Delta^{(k)}(t)$ are shown in Fig.
3(a-c), for systems of size up to $L=28$ (i.e. $q\simeq 6.7\times 10^{7}$).
Focusing first on $k=1$, i.e. regular thermalization, we see an exponential
decay followed by saturation to a minimum. This final plateau is seen to be a
finite-size effect: it corresponds to the finite value one would obtain from a
global Haar-random state on $AB$, which is $\sim 1/\sqrt{q}$, as we explain in
Appendix D. For $q\to\infty$ the exponential decay would continue
indefinitely.

Moving now to $k>1$, we see that the design distances $\Delta^{(k)}$, while
qualitatively replicating the behavior of $\Delta^{(1)}$, weakly increase with
$k$. The $k$-dependence is visibly slower in the intermediate-time exponential
regime and faster in the late-time plateau. In Fig. 3(d-f) we analyze this
$k$-dependence in detail, at different system sizes and times. For times
corresponding to the exponential-decay regime, we see an excellent agreement
with the analytical prediction

$\frac{\Delta^{(k)}}{\Delta^{(1)}}\simeq
f(k,d_{A})=\left(\frac{1+d_{A}}{1+d_{A}/k}\right)^{1/2}$ (44)

(dashed line). The quality of the agreement is limited by two factors:

  * •

At early time, $\Delta^{(1)}(T)$ is not yet very small, so the
$O((\Delta^{(1)})^{2})$ corrections in Eq. (35) are significant. Their effect
can be seen clearly e.g. in Fig. 3(f), with curves for $T=1$, $2$, $3$ showing
progressively better agreement.

  * •

At later time, finite-$q$ effects become dominant. As soon as the
$\Delta^{(k)}(T)$ curves start peeling off from the exponential decay to
settle onto the finite-size plateau value $\sim q^{-1/2}$, the agreement with
the prediction gets worse. The $k$-dependence crosses over to that of a Haar-
random state, which in this case ($d_{A}=d_{B_{1}}=2$) is $\sim\sqrt{k}$, as
we show in Appendix D.

As a consequence of these two effects, the quality of the agreement with the
analytical prediction is non-monotonic: it improves with time at first, and
then becomes worse once finite-size effects become substantial.

However, in the $q\to\infty$ limit, finite-size effects would disappear and we
would recover the predicted behavior for arbitrarily late times. This is
corroborated by the finite-size trends visible in Fig. 3(d-f): with increasing
$q$, the time window of validity of our solution increases, and the agreement
itself improves with time within this window. This represents strong evidence
in support of our main analytical result, Eq. (35-36).

##  5 Discussion and Conclusion

We have presented a toy model of chaotic quantum dynamics where the process of
deep thermalization can be studied exactly in a suitable thermodynamic limit.
The simple structure of the model (a subsystem of interest $A$, a “bottleneck”
$B_{1}$, and an infinite bath $B_{2}$) enables this solvability while
capturing a basic feature of spatially-local, short-range interacting systems
in any dimension: that information may only be exchanged in and out of a
subsystem via the mediation of other subsystems. This simple fact has
important consequences on the process of thermalization and deep
thermalization alike.

Our key result, in Eq. (35, 36), shows how approximate thermalization,
described by a nonzero distance $\Delta^{(1)}$ between the reduced density
matrix $\rho_{A}$ and the infinite-temperature thermal state
$\mathbb{I}/d_{A}$, translates to a larger distance $\Delta^{(k)}$ between the
higher moments of the projected and Haar ensembles. In other words, the two
ensembles become increasingly more distinguishable as one looks at higher
moments of their distributions. A nontrivial growth of design times $t_{k}$ vs
$k$ follows: a logarithmic growth for $k\ll d_{A}$ ($d_{A}$ being the Hilbert
space dimension of $A$), saturating at a finite value $t_{\infty}=2t_{1}$.

This is the second known exact solution for the design times and deep
thermalization in models of chaotic dynamics (following the results of Ref.
[33] on the self-dual kicked Ising chain and its generalization to dual-
unitary circuits in one dimension [43, 34]) and it is the first to exhibit
distinct time scales $t_{k}$. This crucially depends on imperfect (or
asymptotic) thermalization, which is the generic fate of chaotic many-body
systems. Dual-unitary circuits (with suitable initial states) are exceptional
in this regard in that they instead thermalize perfectly in a finite time: one
has $\rho_{A}=\mathbb{I}/d_{A}$ exactly after a finite time $t_{1}=N_{A}$ in
quench dynamics [47]. Thus they evade our present result and show a collapse
of design times, $t_{k}=t_{1}$ for all $k$.

Further, it is worth commenting on how our results compare to other bounds on
deep thermalization that have been obtained in two settings. First, Ref. [72]
proved that, if one chooses the measurement basis in the bath $B$ Haar-
randomly, then the projected ensemble yields an approximate state design with
high probability whenever the reduced density matrix is close to
thermalization at infinite temperature: $\rho_{A}\approx\mathbb{I}_{A}/d_{A}$.
More precisely, given a deviation from ideal thermalization
$\delta=\frac{1}{2}\|\rho_{A}-\mathbb{I}_{A}/d_{A}\|_{1}\leq\frac{1}{2d_{A}}$,
one gets an $\epsilon$-approximate state $k$-design with

$\epsilon\leq 2k\sqrt{d_{A}\delta}+d_{A}\delta.$ (45)

This bound is directly applicable to our model, as the final unitary gate
$V_{T}$ Haar-randomizes the basis on the bath. We have found that our model
realizes an $\epsilon$-approximate state $k$-design with
$\epsilon=f(k,d_{A})\Delta^{(1)}$, where $f$ is as in Eq. (36). This
implies141414We note that Ref. [72] uses the trace norm, while we use the
Frobenius norm; however the norm inequality $\|A\|_{2}\leq\|A\|_{1}$ (for all
$A$) gives $\Delta^{(1)}\leq 2\delta\sqrt{d_{A}}$.
$\epsilon<(2k+1)d_{A}\delta$; using the assumption that $d_{A}\delta<1$, and
this is seen to be consistent with the bound from Ref. [72], Eq. (45). On the
other hand, the techniques we used here may be directly adapted to analyze the
setting considered in Ref. [72], yielding $\epsilon=\Theta(\delta)$. This
suggests that the bound in Eq. (45), which obeys $\epsilon=O(\sqrt{\delta})$
for $\delta\to 0$, may not be tight.

Secondly, Ref. [34] derived a bound on $t_{\infty}$ in one-dimensional unitary
circuits: $t_{\infty}>(v_{E}/v_{p})t_{1}$, where $v_{p}$ is an additional
velocity scale (the “purification velocity”). The derivation of this bound is
based on loss of information about measurement outcomes over long distances,
and generally only applies in one spatial dimension. On the contrary, the
exact results about our toy model characterize the effects of imperfect
thermalization ($\rho_{A}\neq\mathbb{I}_{A}/d_{A}$, a strictly local
condition) on the formation of higher state designs, in a situation that is
arguably a best-case scenario for deep thermalization (a large, all-to-all
interacting random-matrix bath, as opposed to a structured, geometrically-
local one). Thus it is natural to conjecture that the results for this model
might instead hold as bounds for more general models: e.g., we may conjecture
that $t_{\infty}>2t_{1}$ in geometrically-local models in any dimension. As
the arguments leading to the two bounds are a priori independent, one would
expect both of them to apply jointly in one-dimensional chaotic unitary
circuits, giving $t_{\infty}>\max\\{2,v_{E}/v_{p}\\}t_{1}$. It would be
interesting to further investigate the interplay between the two bounds.

Further, it is tempting to speculate that our result for the ratio
$\Delta^{(k)}/\Delta^{(1)}$ may be optimal, i.e., that no ensemble consistent
with a given density matrix $\rho$ can approximate a $k$-design better than
the Scrooge/GAP ensemble; this idea seems supported by the numerical results
of Fig. 3(d-f), where the data never dips below the Scrooge/GAP prediction.
Proving this statement is an interesting question for future work.

By a similar reasoning, our results for the design times $t_{k}$ may also
represent bounds in more general settings—e.g., one may conjecture that the
growth of $t_{k}$ with $k$ should be at least logarithmic (for $k\ll d_{A}$)
in any local, non-perfectly-thermalizing model. This appears consistent with
numerical results on Hamiltonian dynamics in spin chains [32], where in
addition to locality one also has energy diffusion as a possible bottleneck.
It would be interesting to understand the impact of conservation laws on deep
thermalization, starting e.g. from quantum circuit models with $U(1)$ charge
conservation [73, 74, 75, 76] or integrable circuit models [77, 78, 79, 80]. A
particularly tractable setting may be fermionic Gaussian dynamics, for which
the existence of a limiting ensemble was recently shown [81].

## Acknowledgements

We thank S. Choi, T. Rakovszky and V. Khemani for discussions and for previous
collaborations on related topics. M.I. thanks S. L. Sondhi for insightful
discussions. M. I. is supported by the Gordon and Betty Moore Foundation’s
grant GBMF8686. W. W. H. is supported in part by the Stanford Institute of
Theoretical Physics, and by the National University of Singapore (NUS) start-
up grants A-8000599-00-00 and A-8000599-01-00. Numerical simulations were
performed on Stanford Research Computing Center’s Sherlock cluster.

##  Appendix A Equivalence between Eq. (21) and the Scrooge/GAP ensemble

The Scrooge/GAP ensemble [44, 45, 46] is conventionally defined via
“$\rho$-distortion” of the Haar measure on $\mathcal{H}_{A}$: i.e., states
$\ket{\psi}=\sqrt{\rho_{A}}\ket{\phi}/\bra{\phi}\rho_{A}\ket{\phi}$ are drawn
with weights $d_{A}\langle\phi|\rho_{A}|\phi\rangle d\phi$, where $d\phi$ is
the Haar measure on $\mathcal{H}_{A}$. Here, $\sqrt{\rho_{A}}$ is the unique
positive Hermitian square root of $\rho_{A}$. Our ensemble, Eq. (21), is
instead defined from the Haar measure on a larger Hilbert space
$\mathcal{H}_{C}$ (with $d_{C}>d_{A}$) via a rectangular square root of
$\rho_{A}$: the $d_{A}\times d_{C}$ matrix $E$, which obeys
$EE^{\dagger}\propto\rho_{A}$. However, the two definitions are equivalent (in
the sense that they yield precisely the same ensemble on $\mathcal{H}_{A}$).
This can be seen as follows. From a singular-value decomposition of $E$, we
have $E=USPV^{\dagger}$, where $U$ and $V$ are unitary matrices on
$\mathcal{H}_{A}$ and $\mathcal{H}_{C}$ respectively, $S$ is a positive
diagonal $d_{A}\times d_{A}$ matrix of singular values, and $P$ is a
rectangular $d_{A}\times d_{C}$ matrix that projects onto the first $d_{A}$
entries (i.e., $P_{ij}=\delta_{ij}$ if $j\leq d_{A}$, 0 otherwise). From
$\rho_{A}\propto EE^{\dagger}=US^{2}U^{\dagger}$ we have
$E\propto\sqrt{\rho_{A}}UPV^{\dagger}$; further noting that
$UP=P(U\otimes\mathbb{I})$, where the identity operator has size
$d_{C}/d_{A}$, and absorbing $(U\otimes\mathbb{I})V^{\dagger}$ in the Haar
measure on $\mathcal{H}_{C}$, we see that in the ensemble of Eq. (21) we may
replace $E$ by $\sqrt{\rho_{A}}P$. Finally, we see that
$P\ket{\phi}/\|P\ket{\phi}\|$ is a Haar-random vector of size $d_{A}$ (i.e., a
projector acting on a Haar-random state produces another Haar-random state of
the appropriate size, up to normalization). We can thus replace the Haar
measure ${\rm d}\phi$ on $C$ with the analogous measure on $A$, obtaining the
usual $\rho$-distortion construction of the Scrooge/GAP ensemble.

##  Appendix B Thermalization of subsystem $A$

Here we derive the exact form of the purity of subsystem $A$ as a function of
time $T$, in the limit $q\to\infty$, averaged over Haar-random gates $U_{t}$.

Starting from the reduced density matrix $\rho_{A}$ as written in Eq. (22), we
have

$P_{A}(T)\equiv{\rm Tr}\rho_{A}^{2}(T)=d_{B_{1}}^{-2T}{\rm
Tr}(EE^{\dagger})^{2}.$ (46)

This is represented diagrammatically in Fig. 4(a), where the boundary
conditions are permutations of the two replicas, $e,\tau\in S_{2}$: the
identity $e$ before and after each gate on $B_{1}$, and the transposition
$\tau$ at the final time on $A$. We note that $d_{B_{1}}^{-2}|e)(e|_{B_{1}}$
is a super-operator implementing the fully depolarizing (or erasure) channel
on each replica of $B_{1}$—i.e., it is the partial trace on $B_{1}$ followed
by introducing the fully-mixed state $\mathbb{I}_{B_{1}}/d_{B_{1}}$. This
Markovian evolution is a consequence of taking the $q\to\infty$ limit on
$B_{2}$, which introduces irreversibility in the local dynamics of $A$ and
$B_{1}$.

We compute the purity, Eq. (46), averaged over the choice of gates
$\\{U_{t}:t=0,\dots T-1\\}$ which make up $E$. Each $U_{t}$ is sampled
independently and identically form the Haar measure. We compute the average
iteratively, starting from the last gate, $U_{T-1}$. Acting on the operator
$|\hat{\tau})_{A}\otimes|\hat{e})_{B_{1}}$ above it, the Haar-averaged gate
yields

$[\mathbb{E}(U\otimes U^{\ast})^{\otimes
2}]|\hat{\tau})_{A}\otimes|\hat{e})_{B_{1}}=\frac{(d_{A}^{2}-1)d_{B_{1}}}{d_{A}^{2}d_{B_{1}}^{2}-1}|\hat{\tau})_{A}\otimes|\hat{\tau})_{B_{1}}+\frac{d_{A}(d_{B_{1}}^{2}-1)}{d_{A}^{2}d_{B_{1}}^{2}-1}|\hat{e})_{A}\otimes|\hat{e})_{B_{1}}.$
(47)

This result admits a simple interpretation. At the final time $T$ we have a
“domain wall” boundary condition, with permutation $\tau$ in $A$ and $e$ in
$B_{1}$. This configuration is evolved backward in time and gives two possible
configurations: one where the whole system is polarized along $e$, and one
where it is polarized along $\tau$. This structure of the Haar average for
second-moment quantities was used previously to obtain exact results on
entanglement growth, by mapping to a simple random walk for an extended
$1+1$-dimensional quantum circuit [65, 82, 59].

[width=0.7]fig_therm.pdf

Figure 4: Purity of subsystem $A$. (a) Tensor network contraction expressing
${\rm Tr}(\rho_{A}^{2})$ after evolution for time $T=3$. The $|\cdots)$
symbols denote operators on two replicas of the appropriate Hilbert space,
with $|0)=\outerproduct{0}{0}^{\otimes 2}$ (initial state on $A$), and $|e)$
and $|\tau)$ are replica permutations. Note that $d_{B_{1}}^{-2}|e)$
corresponds to two replicas of the completely-mixed state on $B_{1}$,
$(\mathbb{I}_{B_{1}}/d_{B_{1}})^{\otimes 2}$. Upon averaging $U_{T-1}$ over
the Haar measure one obtains two terms: (b) the same diagram as in (a), but
with $T\mapsto T-1$; (c) a diagram with a uniform $|e)$ boundary condition at
the final time, which simply computes the trace of the state and is trivial
due to trace-preservation. This gives the recursion in Eq. (49).

Moving down along the tensor network, we must contract the result of Eq. (47)
with $d_{B_{1}}^{-2}|\hat{e})(\hat{e}|_{B_{1}}$ (i.e. the erasure channel on
$B_{1}$), obtaining

$\frac{1}{d_{B_{1}}^{2}}|\hat{e})(\hat{e}|_{B_{1}}[\mathbb{E}(U\otimes
U^{\ast})^{\otimes
2}]|\hat{\tau})_{A}\otimes|\hat{e})_{B_{1}}=\frac{d_{A}^{2}-1}{d_{A}^{2}d_{B_{1}}^{2}-1}|\hat{\tau})_{A}\otimes|\hat{e})_{B_{1}}+\frac{d_{A}(d_{B_{1}}^{2}-1)}{d_{A}^{2}d_{B_{1}}^{2}-1}|\hat{e})_{A}\otimes|\hat{e})_{B_{1}}$
(48)

(we use the fact that $(\hat{e}|\hat{e})_{B_{1}}=d_{B_{1}}^{2}$ and
$(\hat{e}|\hat{\tau})_{B_{1}}=d_{B_{1}}$). The first term, shown in Fig. 4(b),
is identical to the initial boundary condition (a domain wall between $\tau$
and $e$) up to a prefactor; the second term, shown in Fig. 4(c), is a uniform
$|e)_{AB_{1}}$ state, which is invariant under the remaining backward
evolution151515This follows from trace preservation, and thus holds under any
quantum channels, including in particular the unitary and depolarizing
channels considered in this calculation.. This is enough to derive the
following recursion:

$\mathbb{E}P_{A}(T)=\frac{(d_{A}^{2}-1)\mathbb{E}P_{A}(T-1)+d_{A}(d_{B_{1}}^{2}-1)}{d_{A}^{2}d_{B_{1}}^{2}-1}.$
(49)

The recursion is solved exactly (with initial condition $P_{A}(0)=1$) by

$\mathbb{E}P_{A}(T)=\frac{1}{d_{A}}+\left(1-\frac{1}{d_{A}}\right)\left[\frac{d_{A}^{2}-1}{d_{A}^{2}d_{B_{1}}^{2}-1}\right]^{T}$
(50)

This obeys several sanity checks: $d_{A}^{-1}<P_{A}(T)\leq 1$ (corresponding
to a valid purity), $P_{A}(\infty)=d_{A}^{-1}$ (corresponding to the
completely mixed state), and $P_{A}(T)\equiv 1$ for $d_{B_{1}}=1$ (trivial
subsystem $B_{1}$, so that $A$ is in fact isolated and remains pure).

For a large subsystem $d_{A}\gg 1$, Eq. (50) reduces to $P_{A}(T)\simeq
d_{A}^{-1}+d_{B}^{-2T}$. This has two limiting behaviors for the annealed
average of the second Renyi entropy:

$S_{2}(A,T)=-\log_{2}\mathbb{E}P_{A}(T)=\left\\{\begin{aligned}
&2\log_{2}(d_{B_{1}})T&\text{ if }d_{B_{1}}^{2T}\ll d_{A}\\\
&N_{A}-\frac{1}{\ln 2}d_{A}d_{B_{1}}^{-2T}&\text{ if }d_{B_{1}}^{2T}\gg
d_{A}\end{aligned}\right.$ (51)

which yields Eq. (37), proving our claim that the early-time slope of
$S_{2}(A,T)$ coincides with the time constant of the exponential approach to
the asymptote $S_{2}(A)=N_{A}$.

This result is also expected to hold more generally in extended chaotic
systems where entanglement is described by a ‘minimal membrane’ picture [65,
67, 68, 66, 69]. There, the entropy is given by the free energy of a classical
statistical system with a partition function
$\mathcal{Z}=\sum_{\gamma}e^{-E(\gamma)}$. Here $\gamma$ labels membranes that
partition the spacetime volume into regions whose boundaries are $A$ and its
complement, and $E(\gamma)$ is a suitable cost functional akin to a surface
tension. In a local quantum circuit of depth $T$ in arbitrary spatial
dimension $d$, for a contiguous region $A$, there are two distinct kinds of
membranes to consider: temporal ones, $E(\gamma)\sim|\partial A|T$, and
spatial ones, $E(\gamma)\sim N_{A}$. The former are dominant at small $T$,
giving ballistic growth of entanglement, motivating the definition of an
entanglement velocity $v_{E}$ via $E(\gamma)=v_{E}T$. The latter become
dominant at large $T$, giving saturation to a volume-law entangled state. In
general we have $S\simeq-\log_{2}(2^{-N_{A}}+2^{-v_{E}T})$, which, expanded at
large $T$, gives

$S\simeq N_{A}-\log_{2}(1+2^{N_{A}-v_{E}T})\simeq N_{A}-\frac{1}{\ln
2}2^{-v_{E}(T-t^{\ast})},\quad t^{\ast}=N_{A}/v_{E}.$ (52)

Note however that exceptions to the above behavior have been recently
discovered in finite one-dimensional random-unitary circuits with certain
architectures and boundary conditions [70, 71].

##  Appendix C Details on computational method

Here we present some details on the numerical method used to produce the data
in Fig. 3. As explained in Sec. 4, we use $d_{A}=d_{B_{1}}=2$ and $q=2^{L-2}$,
and we begin with the all-zero state $\ket{\Psi(0)}=\ket{0}^{\otimes L}$.
There are two computational tasks to discuss: how to implement time evolution,
and how to build the moment operator $\rho^{(k)}$ from a given state.

#### Time evolution.

The action of gates $U_{t}\in U(4)$ is implemented via dense matrix-matrix
multiplication ($\Psi$ being viewed as a $4\times q$ matrix) which takes time
$O(q)$. Doing the same thing for $V_{t}$ would take time $O(q^{2})$; however
this is not necessary. Instead we use the following algorithm:

  1. 1.

decompose

$\ket{\Psi(t)}=c_{0}\ket{0}_{A}\ket{\chi_{0}}_{B}+c_{1}\ket{1}_{A}\ket{\chi_{1}}_{B},$
(53)

where $\ket{\chi_{0}}$, $\ket{\chi_{1}}$ are normalized;

  2. 2.

compute the overlap $\innerproduct{\chi_{0}}{\chi_{1}}\equiv\cos\theta
e^{i\varphi}$;

  3. 3.

generate two orthonormal Haar-random vectors $\ket{\eta_{0}}$,
$\ket{\eta_{1}}$;

  4. 4.

update $\ket{\Psi(t)}\mapsto
c_{0}\ket{0}_{A}\ket{\eta_{0}}_{B}+c_{1}\ket{1}_{A}(\cos\theta
e^{i\varphi}\ket{\eta_{0}}+\sin\theta\ket{\eta_{1}})_{B}$.

This method only deals with the action of $V_{t}$ on the subspace spanned by
$\ket{\chi_{0}}$ and $\ket{\chi_{1}}$, which is all that is needed to evolve
$\ket{\Psi}$, and requires time $O(q)$.

#### Moment operator.

For a fixed value of $k$, $\rho^{(k)}$ lives in a replicated Hilbert space
$\mathcal{H}_{A}^{\otimes k}$ of dimension $d_{A}^{k}$; however it is
supported only in the symmetric sector, of dimension $\binom{k+d_{A}-1}{k}$.
In the case of interest, $d_{A}=2$, this dimension reduces to $k+1$: a basis
is given by symmetrized versions of $\ket{0^{\otimes l}1^{\otimes k-l}}$ for
$l=0,\dots k$. All the information in $\rho^{(k)}$ is thus obtained from the
matrix elements

$\bra{0^{\otimes l_{1}}1^{\otimes k-l_{1}}}\rho^{(k)}\ket{0^{\otimes
l_{2}}1^{\otimes
k-l_{2}}}=\sum_{z}p(z)\innerproduct{0}{\psi_{z}}^{l_{1}}\innerproduct{1}{\psi_{z}}^{k-l_{1}}\innerproduct{\psi_{z}}{0}^{l_{2}}\innerproduct{\psi_{z}}{1}^{k-l_{2}}$
(54)

for $l_{1}\leq l_{2}$ (the remaining terms are given by Hermiticity). For
fixed $l_{1}$, $l_{2}$ the sum is computed in time $O(q)$ and thus the overall
time cost to get $\rho^{(k)}$ is $O(qk^{2})$. Iterating this process for all
$k$ up to $k_{\text{max}}$ costs $O(qk_{\text{max}}^{3})$.

Note that for general $d_{A}$ one has $\binom{d_{A}-1+k}{k}\sim k^{d_{A}-1}$
(at large $k$), giving total cost $O(qk_{\text{max}}^{2d_{A}-1})$; thus, in
order to analyze the largest possible range in $k$, we set $d_{A}=2$.

##  Appendix D $\Delta^{(k)}$ for Haar-random states

Here we derive the form of the design distance measures $\Delta^{(k)}$ for
Haar-random states on a bipartite system $AB$ and show that it matches the
post-saturation behavior observed in our numerical simulations, Fig. 3. We
employ the same replica trick as in Sec. 3.1 and compute the Haar-average of
the pseudo-frame potential

$F^{(n,k)}=\sum_{z_{1},z_{2}=1}^{d_{B}}\innerproduct{\tilde{\psi}_{z_{1}}}{\tilde{\psi}_{z_{1}}}^{n}|\innerproduct{\tilde{\psi}_{z_{1}}}{\tilde{\psi}_{z_{2}}}|^{2k}\innerproduct{\tilde{\psi}_{z_{2}}}{\tilde{\psi}_{z_{2}}}^{n}={\rm
Tr}(\mathcal{Q}_{B}\outerproduct{\Psi}{\Psi}_{AB}^{\otimes R})$ (55)

where $\ket{\Psi}_{AB}$ is the state of $AB$ (to be averaged over the Haar
measure) and $\ket{\tilde{\psi}_{z}}_{A}={}_{B}\\!\innerproduct{z}{\Psi}_{AB}$
is the unnormalized post-measurement state of $A$. The operator
$\mathcal{Q}_{B}$ is as in Eq. (11) and $R=2(n+k)$ is the total number of
replicas. On average we have

$\mathbb{E}F^{(n,k)}={\rm
Tr}(\mathcal{Q}_{B}\rho_{H,AB}^{(R)})=\frac{(d_{A}d_{B}-1)!}{(R+d_{A}d_{B}-1)!}\sum_{\sigma\in
S_{R}}{\rm Tr}(\mathcal{Q}_{B}\hat{\sigma}_{AB}).$ (56)

Next, we have
$(\mathcal{Q}|\hat{\sigma})_{B}=d_{B}+d_{B}(d_{B}-1)\delta_{\sigma,\tau\pi_{1}\pi_{2}}$,
where $\tau$ and $\pi_{1,2}$ are as in Fig. 2(e) and the two terms come from
expanding $\mathcal{Q}$ into $z_{1}=z_{2}$ and $z_{1}\neq z_{2}$ terms.
Following the same manipulations as in Sec. 3.1, one gets

$\sum_{\sigma\in S_{R}}{\rm
Tr}(\mathcal{Q}_{B}\hat{\sigma}_{AB})=d_{B}\frac{(R+d_{A}-1)!}{(d_{A}-1)!}+d_{B}(d_{B}-1)\left[\frac{(R/2+d_{A}-1)!}{(d_{A}-1)!}\right]^{2}F_{H}^{(k)}$
(57)

and thus, taking $n=1-k$ (i.e. $R=2$),

$\mathbb{E}F^{(k)}=\frac{d_{A}+1+d_{A}(d_{B}-1)F_{H}^{(k)}}{d_{A}d_{B}-1}$
(58)

Finally, we obtain

$\mathbb{E}(\Delta^{(k)})^{2}=\frac{\mathbb{E}F^{(k)}}{F_{H}^{(k)}}-1=\frac{d_{A}+1}{d_{A}d_{B}-1}\left[\binom{k+d_{A}-1}{k}-1\right].$
(59)

Incidentally, we note that this identity may serve as an alternative approach
to derive the results of Ref. [32] on obtaining an $\epsilon$-approximate
$k$-design from the projected ensemble on a Haar-random state with a given
probability.

Thus far we have made no assumptions on $d_{A,B}$. Next, we plug in the values
used in our numerical simulations, see Sec. 4: $d_{A}=2$, $d_{B}=2q$. We find

$(\Delta^{(k)})_{\text{r.m.s.}}=[\mathbb{E}(\Delta^{(k)})^{2}]^{1/2}=\left(\frac{3k}{4q+1}\right)^{1/2}.$
(60)

This shows at once that we have a finite-size “floor” $\Delta^{(1)}\sim
q^{-1/2}$, attained when the dynamics has successfully Haar-randomized the
entire system $AB$, and that in this regime the $k$-dependence is
$\Delta^{(k)}=\sqrt{k}\Delta^{(1)}$. Notably the ratio
$\Delta^{(k)}/\Delta^{(1)}$ is unbounded, in contrast with our result Eq. (35)
which for the same value of $d_{A}$ gives
$\Delta^{(k)}=\sqrt{3k/(k+2)}\Delta^{(1)}\leq\sqrt{3}\Delta^{(1)}$. This
accounts for the distinct behaviors visible in Fig. 3(a-c) in the two regimes.

## References

  * [1] Ehud Altman, Kenneth R. Brown, Giuseppe Carleo, Lincoln D. Carr, Eugene Demler, Cheng Chin, et al.  “Quantum Simulators: Architectures and Opportunities”.  PRX Quantum 2, 017003 (2021).
  * [2] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S. Zibrov, Manuel Endres, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin.  “Probing many-body dynamics on a 51-atom quantum simulator”.  Nature 551, 579–584 (2017).
  * [3] Soonwon Choi, Joonhee Choi, Renate Landig, Georg Kucsko, Hengyun Zhou, Junichi Isoya, Fedor Jelezko, Shinobu Onoda, Hitoshi Sumiya, Vedika Khemani, Curt von Keyserlingk, Norman Y. Yao, Eugene Demler, and Mikhail D. Lukin.  “Observation of discrete time-crystalline order in a disordered dipolar many-body system”.  Nature 543, 221–225 (2017).
  * [4] Xiao Mi, Matteo Ippoliti, Chris Quintana, Ami Greene, Zijun Chen, Jonathan Gross, et al.  “Time-crystalline eigenstate order on a quantum processor”.  Nature 601, 531–536 (2022).
  * [5] J. Randall, C. E. Bradley, F. V. van der Gronden, A. Galicia, M. H. Abobeih, M. Markham, D. J. Twitchen, F. Machado, N. Y. Yao, and T. H. Taminiau.  “Many-body–localized discrete time crystal with a programmable spin-based quantum simulator”.  Science 374, 1474–1478 (2021).
  * [6] Philipp T. Dumitrescu, Justin G. Bohnet, John P. Gaebler, Aaron Hankin, David Hayes, Ajesh Kumar, Brian Neyenhuis, Romain Vasseur, and Andrew C. Potter.  “Dynamical topological phase realized in a trapped-ion quantum simulator”.  Nature 607, 463–467 (2022).
  * [7] Brian Skinner, Jonathan Ruhman, and Adam Nahum.  “Measurement-Induced Phase Transitions in the Dynamics of Entanglement”.  Physical Review X 9, 031009 (2019).
  * [8] Yaodong Li, Xiao Chen, and Matthew P. A. Fisher.  “Quantum Zeno effect and the many-body entanglement transition”.  Physical Review B 98, 205136 (2018).
  * [9] Soonwon Choi, Yimu Bao, Xiao-Liang Qi, and Ehud Altman.  “Quantum Error Correction in Scrambling Dynamics and Measurement-Induced Phase Transition”.  Physical Review Letters 125, 030505 (2020).
  * [10] Michael J. Gullans and David A. Huse.  “Dynamical Purification Phase Transition Induced by Quantum Measurements”.  Physical Review X 10, 041020 (2020).
  * [11] Matteo Ippoliti, Michael J. Gullans, Sarang Gopalakrishnan, David A. Huse, and Vedika Khemani.  “Entanglement Phase Transitions in Measurement-Only Dynamics”.  Physical Review X 11, 011030 (2021).
  * [12] Andrew C. Potter and Romain Vasseur.  “Entanglement dynamics in hybrid quantum circuits” (2021).  arXiv:2111.08018.
  * [13] Matthew P. A. Fisher, Vedika Khemani, Adam Nahum, and Sagar Vijay.  “Random Quantum Circuits” (2022).  arXiv:2207.14280.
  * [14] Crystal Noel, Pradeep Niroula, Daiwei Zhu, Andrew Risinger, Laird Egan, Debopriyo Biswas, Marko Cetina, Alexey V. Gorshkov, Michael J. Gullans, David A. Huse, and Christopher Monroe.  “Measurement-induced quantum phases realized in a trapped-ion quantum computer”.  Nature Physics 18, 760–764 (2022).
  * [15] Jin Ming Koh, Shi-Ning Sun, Mario Motta, and Austin J. Minnich.  “Experimental Realization of a Measurement-Induced Entanglement Phase Transition on a Superconducting Quantum Processor” (2022).  arXiv:2203.04338.
  * [16] Scott Aaronson.  “Shadow Tomography of Quantum States” (2018).  arXiv:1711.01053.
  * [17] 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).
  * [18] Andreas Elben, Steven T. Flammia, Hsin-Yuan Huang, Richard Kueng, John Preskill, Benoît Vermersch, and Peter Zoller.  “The randomized measurement toolbox” (2022).  arXiv:2203.11374.
  * [19] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, et al.  “Quantum supremacy using a programmable superconducting processor”.  Nature 574, 505–510 (2019).
  * [20] Yulin Wu, Wan-Su Bao, Sirui Cao, Fusheng Chen, Ming-Cheng Chen, Xiawei Chen, et al.  “Strong Quantum Computational Advantage Using a Superconducting Quantum Processor”.  Physical Review Letters 127, 180501 (2021).
  * [21] Scott Aaronson and Lijie Chen.  “Complexity-theoretic foundations of quantum supremacy experiments”.  In Proceedings of the 32nd Computational Complexity Conference.  Pages 1–67.  CCC ’17Dagstuhl, DEU (2017).
  * [22] Yiqing Zhou, E. Miles Stoudenmire, and Xavier Waintal.  “What Limits the Simulation of Quantum Computers?”.  Physical Review X 10, 041038 (2020).
  * [23] Xun Gao, Marcin Kalinowski, Chi-Ning Chou, Mikhail D. Lukin, Boaz Barak, and Soonwon Choi.  “Limitations of Linear Cross-Entropy as a Measure for Quantum Advantage” (2021).  arXiv:2112.01657.
  * [24] Abhinav Deshpande, Pradeep Niroula, Oles Shtanko, Alexey V. Gorshkov, Bill Fefferman, and Michael J. Gullans.  “Tight bounds on the convergence of noisy random circuits to the uniform distribution” (2021).  arXiv:2112.00716.
  * [25] Alexander M. Dalzell, Nicholas Hunter-Jones, and Fernando G. S. L. Brandão.  “Random quantum circuits transform local noise into global white noise” (2021).  arXiv:2111.14907.
  * [26] Joonhee Choi, Adam L. Shaw, Ivaylo S. Madjarov, Xin Xie, Jacob P. Covey, Jordan S. Cotler, Daniel K. Mark, Hsin-Yuan Huang, Anant Kale, Hannes Pichler, Fernando G. S. L. Brandão, Soonwon Choi, and Manuel Endres.  “Emergent Randomness and Benchmarking from Many-Body Quantum Chaos” (2021).  arXiv:2103.03535.
  * [27] Daniel K. Mark, Joonhee Choi, Adam L. Shaw, Manuel Endres, and Soonwon Choi.  “Benchmarking Quantum Simulators using Quantum Chaos” (2022).  arXiv:2205.12211.
  * [28] Mark Srednicki.  “Chaos and quantum thermalization”.  Physical Review E 50, 888–901 (1994).
  * [29] Marcos Rigol, Vanja Dunjko, and Maxim Olshanii.  “Thermalization and its mechanism for generic isolated quantum systems”.  Nature 452, 854–858 (2008).
  * [30] Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn.  “Colloquium: Many-body localization, thermalization, and entanglement”.  Reviews of Modern Physics 91, 021001 (2019).
  * [31] Álvaro M. Alhambra.  “Quantum many-body systems in thermal equilibrium” (2022).  arXiv:2204.08349.
  * [32] 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 wavefunctions” (2021).  arXiv:2103.03536.
  * [33] Wen Wei Ho and Soonwon Choi.  “Exact Emergent Quantum State Designs from Quantum Chaotic Dynamics”.  Physical Review Letters 128, 060601 (2022).
  * [34] Matteo Ippoliti and Wen Wei Ho.  “Dynamical purification and the emergence of quantum state designs from the projected ensemble” (2022).  arXiv:2204.13657.
  * [35] Andris Ambainis and Joseph Emerson.  “Quantum t-designs: t-wise Independence in the Quantum World”.  In Twenty-Second Annual IEEE Conference on Computational Complexity (CCC’07).  Pages 129–140.  (2007).
  * [36] D. Gross, K. Audenaert, and J. Eisert.  “Evenly distributed unitaries: On the structure of unitary designs”.  Journal of Mathematical Physics 48, 052104 (2007).
  * [37] Richard A. Low.  “Pseudo-randomness and Learning in Quantum Computation” (2010).  arXiv:1006.5227.
  * [38] M. Akila, D. Waltner, B. Gutkin, and T. Guhr.  “Particle-time duality in the kicked Ising spin chain”.  Journal of Physics A: Mathematical and Theoretical 49, 375101 (2016).
  * [39] Bruno Bertini, Pavel Kos, and Tomaž Prosen.  “Exact Spectral Form Factor in a Minimal Model of Many-Body Quantum Chaos”.  Physical Review Letters 121, 264101 (2018).
  * [40] Bruno Bertini, Pavel Kos, and Tomaž Prosen.  “Exact Correlation Functions for Dual-Unitary Lattice Models in 1+1 Dimensions”.  Phys. Rev. Lett. 123, 210601 (2019).
  * [41] Sarang Gopalakrishnan and Austen Lamacraft.  “Unitary circuits of finite depth and infinite width from quantum channels”.  Physical Review B 100, 064309 (2019).
  * [42] Pieter W. Claeys and Austen Lamacraft.  “Maximum velocity quantum circuits”.  Physical Review Research 2, 033032 (2020).
  * [43] Pieter W. Claeys and Austen Lamacraft.  “Emergent quantum state designs and biunitarity in dual-unitary circuit dynamics”.  Quantum 6, 738 (2022).
  * [44] Richard Jozsa, Daniel Robb, and William K. Wootters.  “Lower bound for accessible information in quantum mechanics”.  Physical Review A 49, 668–677 (1994).
  * [45] Sheldon Goldstein, Joel L. Lebowitz, Roderich Tumulka, and Nino Zanghì.  “On the Distribution of the Wave Function for Systems in Thermal Equilibrium”.  Journal of Statistical Physics 125, 1193–1221 (2006).
  * [46] Sheldon Goldstein, Joel L. Lebowitz, Christian Mastrodonato, Roderich Tumulka, and Nino Zanghì.  “Universal Probability Distribution for the Wave Function of a Quantum System Entangled with its Environment”.  Communications in Mathematical Physics 342, 965–988 (2016).
  * [47] Lorenzo Piroli, Bruno Bertini, J. Ignacio Cirac, and Tomaž Prosen.  “Exact dynamics in dual-unitary quantum circuits”.  Physical Review B 101, 094304 (2020).
  * [48] Fabio Anza and James P. Crutchfield.  “Beyond density matrices: Geometric quantum states”.  Physical Review A 103, 062218 (2021).
  * [49] Fabio Anza and James P. Crutchfield.  “Quantum Information Dimension and Geometric Entropy”.  PRX Quantum 3, 020355 (2022).
  * [50] Joseph M. Renes, Robin Blume-Kohout, A. J. Scott, and Carlton M. Caves.  “Symmetric informationally complete quantum measurements”.  Journal of Mathematical Physics 45, 2171–2180 (2004).
  * [51] E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland.  “Randomized benchmarking of quantum gates”.  Physical Review A 77, 012307 (2008).
  * [52] Daniel A. Roberts and Beni Yoshida.  “Chaos and complexity by design”.  Journal of High Energy Physics 2017, 121 (2017).
  * [53] M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac.  “Matrix Product States for Dynamical Simulation of Infinite Chains”.  Phys. Rev. Lett. 102, 240603 (2009).
  * [54] M. B. Hastings and R. Mahajan.  “Connecting entanglement in time and space: Improving the folding algorithm”.  Phys. Rev. A 91, 032306 (2015).
  * [55] Alessio Lerose, Michael Sonner, and Dmitry A. Abanin.  “Influence Matrix Approach to Many-Body Floquet Dynamics”.  Physical Review X 11, 021040 (2021).
  * [56] Michael Sonner, Alessio Lerose, and Dmitry A. Abanin.  “Influence functional of many-body systems: Temporal entanglement and matrix-product state representation”.  Annals of Physics 435, 168677 (2021).
  * [57] Giacomo Giudice, Giuliano Giudici, Michael Sonner, Julian Thoenniss, Alessio Lerose, Dmitry A. Abanin, and Lorenzo Piroli.  “Temporal Entanglement, Quasiparticles, and the Role of Interactions”.  Physical Review Letters 128, 220401 (2022).
  * [58] Matteo Ippoliti and Vedika Khemani.  “Postselection-Free Entanglement Dynamics via Spacetime Duality”.  Physical Review Letters 126, 060501 (2021).
  * [59] Matteo Ippoliti, Tibor Rakovszky, and Vedika Khemani.  “Fractal, Logarithmic, and Volume-Law Entangled Nonthermal Steady States via Spacetime Duality”.  Physical Review X 12, 011045 (2022).
  * [60] Tsung-Cheng Lu and Tarun Grover.  “Spacetime duality between localization transitions and measurement-induced transitions”.  PRX Quantum 2, 040319 (2021).
  * [61] David T. Stephen, Wen Wei Ho, Tzu-Chieh Wei, Robert Raussendorf, and Ruben Verresen.  “Universal measurement-based quantum computation in a one-dimensional architecture enabled by dual-unitary circuits” (2022).  arXiv:2209.06191.
  * [62] Felix A. Pollock, Cesar Rodriguez-Rosario, Thomas Frauenheim, Mauro Paternostro, and Kavan Modi.  “Non-Markovian quantum processes: Complete framework and efficient characterization”.  Physical Review A 97, 012127 (2018).
  * [63] Georg Köstenberger.  “Weingarten Calculus” (2021).  arXiv:2101.00921.
  * [64] Hyungwon Kim and David A. Huse.  “Ballistic Spreading of Entanglement in a Diffusive Nonintegrable System”.  Physical Review Letters 111, 127205 (2013).
  * [65] Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah.  “Quantum Entanglement Growth under Random Unitary Dynamics”.  Physical Review X 7, 031016 (2017).
  * [66] Cheryne Jonay, David A. Huse, and Adam Nahum.  “Coarse-grained dynamics of operator and state entanglement” (2018).  arXiv:1803.00089.
  * [67] Tianci Zhou and Adam Nahum.  “Emergent statistical mechanics of entanglement in random unitary circuits”.  Phys. Rev. B 99, 174205 (2019).
  * [68] Tianci Zhou and Adam Nahum.  “Entanglement Membrane in Chaotic Many-Body Systems”.  Phys. Rev. X 10, 031066 (2020).
  * [69] Yaodong Li and Matthew P. A. Fisher.  “Statistical mechanics of quantum error correcting codes”.  Physical Review B 103, 104306 (2021).
  * [70] Jac Bensa and Marko Znidaric.  “Fastest Local Entanglement Scrambler, Multistage Thermalization, and a Non-Hermitian Phantom”.  Physical Review X 11, 031019 (2021).
  * [71] Jac Bensa and Marko Znidaric.  “Two-step phantom relaxation of out-of-time-ordered correlations in random circuits”.  Physical Review Research 4, 013228 (2022).
  * [72] Henrik Wilming and Ingo Roth.  “High-temperature thermalization implies the emergence of quantum state designs” (2022).  arXiv:2202.01669.
  * [73] Tibor Rakovszky, Frank Pollmann, and C. W. von Keyserlingk.  “Diffusive Hydrodynamics of Out-of-Time-Ordered Correlators with Charge Conservation”.  Phys. Rev. X 8, 031058 (2018).
  * [74] Vedika Khemani, Ashvin Vishwanath, and David A. Huse.  “Operator Spreading and the Emergence of Dissipative Hydrodynamics under Unitary Evolution with Conservation Laws”.  Phys. Rev. X 8, 031057 (2018).
  * [75] Nicholas Hunter-Jones.  “Operator growth in random quantum circuits with symmetry” (2018).  arXiv:1812.08219.
  * [76] Utkarsh Agrawal, Aidan Zabalo, Kun Chen, Justin H. Wilson, Andrew C. Potter, J. H. Pixley, Sarang Gopalakrishnan, and Romain Vasseur.  “Entanglement and charge-sharpening transitions in U(1) symmetric monitored quantum circuits” (2021).  arXiv:2107.10279.
  * [77] Sarang Gopalakrishnan.  “Operator growth and eigenstate entanglement in an interacting integrable Floquet system”.  Physical Review B 98, 060302 (2018).
  * [78] Katja Klobas, Bruno Bertini, and Lorenzo Piroli.  “Exact Thermalization Dynamics in the “Rule 54” Quantum Cellular Automaton”.  Physical Review Letters 126, 160602 (2021).
  * [79] Berislav Buča, Katja Klobas, and Tomaž Prosen.  “Rule 54: exactly solvable model of nonequilibrium statistical mechanics”.  Journal of Statistical Mechanics: Theory and Experiment 2021, 074001 (2021).
  * [80] Hansveer Singh, Romain Vasseur, and Sarang Gopalakrishnan.  “The Fredkin staircase: An integrable system with a finite-frequency Drude peak” (2022).  arXiv:2205.08542.
  * [81] Maxime Lucas, Lorenzo Piroli, Jacopo De Nardis, and Andrea De Luca.  “Generalized Deep Thermalization for Free Fermions” (2022).  arXiv:2207.13628.
  * [82] C. W. von Keyserlingk, Tibor Rakovszky, Frank Pollmann, and S. L. Sondhi.  “Operator Hydrodynamics, OTOCs, and Entanglement Growth in Systems without Conservation Laws”.  Physical Review X 8, 021013 (2018).