# Designs via Free Probability

Michele Fava, Jorge Kurchan, Silvia Pappalardi

► To cite this version:

Michele Fava, Jorge Kurchan, Silvia Pappalardi. Designs via Free Probability. Physical Review X, 2025, 15 (1), pp.011031. 10.1103/PhysRevX.15.011031. hal-04186851

HAL Id: hal-04186851

https://hal.science/hal-04186851v1

Submitted on 17 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Designs via Free Probability

Michele Fava,¹ Jorge Kurchan,² and Silvia Pappalardi ³,*,

¹Philippe Meyer Institute, Laboratoire de Physique de l'École Normale Supérieure (ENS),
Université PSL, 24 rue Lhomond, F-75231 Paris, France

²Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS,
Sorbonne Université, Université de Paris, F-75005 Paris, France

³Institut für Theoretische Physik, Universität zu Köln, Zülpicher Straße 77, 50937 Köln, Germany

<figure><img src="image_2.png" /></figure>

(Received 27 October 2023; revised 16 May 2024; accepted 24 December 2024; published 12 February 2025)

Unitary designs have become a vital tool for investigating pseudorandomness, since they approximate the statistics of the uniform Haar ensemble. Despite their central role in quantum information, their relations to quantum chaotic evolution and, in particular, to the eigenstate thermalization hypothesis (ETH) are still largely debated issues. This work provides a bridge between the latter and *k* designs through free probability theory. First, by introducing the more general notion of *k*-freeness, we show that it can be used as an alternative probe to designs. In turn, free probability theory comes with several tools, useful, for instance, for the calculation of mixed moments or the so-called *k*-fold quantum channels. Our second result is the connection to quantum dynamics. Quantum ergodicity and, correspondingly, ETH apply to a restricted class of physical observables, as already discussed in the literature. In this spirit, we show that unitary evolution with generic Hamiltonians always leads to freeness at sufficiently long times but only when the operators considered are restricted within the ETH class. Our results provide a direct link between unitary designs, quantum chaos, and the eigenstate thermalization hypothesis and shed new light on the universality of late-time quantum dynamics.

DOI: 10.1103/PhysRevX.15.011031

Subject Areas: Quantum Physics

I. INTRODUCTION

The purpose of this paper is to make explicit the
connections between three lines of thought: *k* designs,
the “full” version of the eigenstate thermalization hypoth-
esis, and free probability.

In recent years, there has been a growing interest in
the study of random unitaries and their uniform distribu-
tion, the Haar ensemble. These concepts are essential in
quantum information, with applications ranging from
cryptography [1] and quantum key distribution [2–4] to
randomized benchmarking and measurements [5–14],
as well as many others [15]. However, sampling unitaries
from the Haar ensemble in practical setups is challenging,
especially for large Hilbert spaces. To address this issue,
*k* designs have been introduced. These correspond to
ensembles of unitaries which approximate the Haar ensem-
ble by reproducing its correlations up to a certain order
*k* [16–18]. Random unitaries play a fundamental role also

in the study of the chaotic evolution of many-body systems.
The reasons are multiple: On one side, many universal
properties of chaotic many-body systems are understood
from random matrix theory; on the other, chaotic evolution
is expected to lead the system to uniformly sample Hilbert
space at sufficiently long times, as some form of ergodicity.
Several works in these directions have studied the emer-
gence of design properties in different models, such as
random unitary circuits [15,19–29] and many-body models
with Brownian couplings [30,31]. Other approaches have
focused on the random properties of ensembles of states
[26,32]. In this case, a popular setup is to consider the
ensemble as arising from projective measurements on a
subsystem, the so-called deep thermalization [33–42].

A framework for understanding quantum statistical mechanics is provided by the eigenstate thermalization hypothesis (ETH), which describes the statistical properties of a class of physical observables $\mathcal{A}$ when represented in the basis of generic many-body Hamiltonians [43,44]. This approach has a well-defined classical limit in Khichin’s statistical mechanics [45], that, rather than on ergodicity, focuses on observables and their properties in the thermodynamic limit. The validity of ETH in local many-body systems has been verified in a plethora of numerical experiments [46–59], making ETH one of the most established paradigms for thermalization in quantum systems.

<footnote>*Contact author: pappalardi@thp.uni-koeln.de</footnote>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.</footnote>
Recently, the introduction of out-of-time order correlators (OTOCs) as a probe of quantum chaos and information scrambling [60–62] has emphasized the importance of higher-order correlations to understand the connection between these two concepts. On the one hand, higher-order OTOCs encoding correlations of the order of *k* and naturally raised the question of whether evolution may give rise to *k* designs, as discussed by Roberts and Yoshida in Ref. [63]. On the other hand, the nontrivial behavior of OTOC motivated the development of a “full” version of ETH [64] which encompasses all correlations between matrix elements of observables. This perspective has attracted a lot of interest from the many-body [65–72] to the high-energy communities [73–75]. More recently, a direct connection between the full ETH and multitime correlation functions has been achieved using the language of free probability [76]. This is a generalization of classical probability theory to noncommuting variables, where “independence” between random variables becomes “free independence”—or “freeness”—standard cumulants become “free cumulants” and many other concepts are revised to be applied to noncommuting objects [77–79]. When dealing with large $D \times D$ (random) matrices, asymptotic freeness applies, valid when the size of the matrix *D* diverges.

In the language of free probability, the full eigenstate thermalization hypothesis takes a very simple form, since the statistics of matrix elements is directly encoded in terms of the thermal free cumulants [76]—a property that was later confirmed numerically [80].

Interestingly, free probability was initially anticipated in planar field theory [81–83], and it has afterward emerged in several other branches of quantum mechanics at both a technical and a conceptual level [84–88]. Lately, it has also appeared in the description of mesoscopic models of quantum transport [89–92] and for studying the thermalization of Wigner matrices [93]. All these findings suggest that free probability, which focuses on noncommuting objects, could have a wider range of applicability in the study of many-body quantum systems.

In this work, we use free probability to study unitary *k* designs and to elucidate the connections to the long-time evolution of large chaotic quantum systems, as summarized in Fig. 1. Our starting point is the definition of *asymptotic k-freeness* which relaxes the one of freeness, restricting it to the order of *k*. Building on this, we show that *k* designs imply asymptotic *k*-freeness for every set of operators, in the limit of large matrix size. As a consequence, our result provides many useful tools for the study of designs in quantum information. These include an asymptotic expansion for quantum channels in terms of free cumulants or a closed expression to compute OTOC.

Second, we consider Hamiltonian time evolution, and we restrict ourselves to physical ETH observables $\mathcal{A}$ based on local operators and, in particular, those having a good

<figure><img src="image_3.png" /><figcaption>FIG. 1. Bridge between *k* designs (on the left) and the eigenstate thermalization hypothesis (on the right) provided by free probability. The connection steps on the definition of asymptotic *k*-freeness valid in the limit of large Hilbert space dimension $\dim \mathcal{H} = D \to \infty$. As a first result, we show that *k* designs imply *k*-freeness asymptotically, i.e., in the limit of large *D*, for every observable. Second, we show that generic unitary time evolution implies *k*-freeness at sufficiently large times but only for observables restricted to the subclass $\mathcal{A}$ of physical observables for which ETH applies.</figcaption></figure>

semiclassical limit. A relation between time evolution and designs restricted to $\mathcal{A}$ was already discussed under the name of “partial unitary design” by Kaneko, Iyoda, and Sagawa in Ref. [94]. In the same ETH spirit, we employ a free probability approach which shifts the attention from designs to the relation between different observables. Importantly, we shall see that the set of unitaries obtained by time evolution implies arbitrary asymptotic *k*-freeness, restricted to $\mathcal{A}$. This result holds at any finite temperature and is mathematically rigorous upon infinite time average. In physical chaotic many-body Hamiltonians, we shall argue that *k*-freeness is obtained at a finite free-*k* time $t_k$ which depends on the Hamiltonian microscopic and on the observables.

The results of this work provide a novel perspective on the relationship between unitary design and the long-time universality of chaotic evolution. Our findings show that at sufficiently long times Hamiltonian dynamics are as good as designs in making physical observables free.

The rest of the paper is organized as follows. In Sec. II, we present an overview of the content of the paper, introducing the concepts and summarizing our new results, with a comment on the differences with the other works in the literature. To make the paper self-contained, in Sec. III, we provide a pedagogical overview of the notions and notations that are studied in this paper. The next two sections contain the bulk of new results: In Sec. IV, we discuss *k* designs via free probability, and, in Sec. V, we provide the bridge between ETH and long-time dynamics. Finally, Sec. VI is devoted to the conclusions and the open questions.

## II. SUMMARY

In this section, we aim to give the reader a summary of the content and results presented in this paper. Our main
findings are summarized and illustrated in Fig. 1 and discussed below.

### A. Full eigenstate thermalization hypothesis

In recent years, there has been a renewed interest in quantum chaos, spurred by new ideas coming from the black-hole and holography community on the one side and from quantum information on the other. Unlike the approaches in the 1970s and 1980s, which were mostly concentrated on semiclassical tools, the interest now is directed on the thermodynamic limit.

In a very influential paper, Maldacena, Shenker, and Stanford derived a temperature-dependent “bound to chaos” [61], formulated in terms of a Lyapunov exponent expressed as a function of out-of-order correlators containing correlations between four matrix elements.

Considering the matrix elements of an operator in the eigenbasis of a chaotic Hamiltonian as a random variable is the subject of the eigenstate thermalization hypothesis, developed by Deutsch and Srednicki [43,95] in the 1990s. When these ideas were put to test in the context of the OTOC, it became clear that ETH as was known was incomplete, because it dealt only with two-element correlations and neglected all those correlations that were essential for the computation of four-point functions.

The full form of ETH, afterward developed in Ref. [64], is an ansatz for the k-point correlations:

$$ \overline{A_{i_1 i_2} A_{i_2 i_3} \dots A_{i_k i_1}} = e^{(1-k)S(e_+)} F_{e^+}^{(k)}(\omega_{i_1 i_2}, \dots, \omega_{i_{k-1} i_k}) \quad \text{for different indices } i_1, i_2, \dots, i_k, \quad (1) $$

where $e_+ = (E_{i_1} + \dots + E_{i_k})/k$, $\omega_{i_k i_{k+1}} = E_k - E_{k+1}$ are the energy differences and S is the thermodynamic entropy. The smooth functions $F^{(k)}$ define the operator. This expression is for different indices, but what about the host of possibilities with different index repetitions that we are now forced to take into account?

### B. Free probability: The framework

It turns out that the treatment of all this proliferation of possibilities, plus the development of the combinatorics that ensues, had already been developed in the field of free probability introduced by Voiculescu [96] and greatly developed by Speicher [78] and many others. Free probability may be seen as generalizing the notion of independence for noncommutative variables. For our purposes, the main ingredient is the definition of “free cumulants” of operators, for example, $\kappa_{2k}(A_1, A_2, A_3, \dots, A_{2k})$, a set of a fixed combination of traces of products, in the spirit of (and taking the place of) ordinary cumulants. The full ETH is expressed by saying that expectations of products with different indices are directly connected with the free cumulants [76].

An important notion we use here is the following—Two operators B and A are said to be free (asymptotically, for large size) if their free cumulants vanish as follows:

$$ \kappa_{2l}(A, B, A, B, \dots, A, B) = 0 \quad \forall \ l \le k. \quad (2) $$

Actually, the definition of freeness involves arbitrary mixed sequences of A and B; for simplicity, we illustrate our arguments using alternating sequences. Importantly, B and a randomly rotated matrix $A = UA'U^\dagger$ are asymptotically free when U is taken to be uniformly distributed over the full unitary group. It is common knowledge that, for large matrices, certain physical transformations bring the eigenvectors of A' and B to become “generic” or “typical” with respect to each other. The definition of freeness replaces this vague notion and provides the tools for computing mixed correlations from the individual asymptotic properties of A and B. In this work, we relax this condition by defining k-freeness when Eq. (2) is true only up to some k. A comment is in order. In most practical settings, rather than requiring the free cumulant Eq. (2) to vanish exactly, it is more convenient to require that its absolute value is smaller than some arbitrary precision e. For simplicity, in what follows, we always write “= 0” as in Eq. (2), although it is intended that there are small corrections, which vanish only asymptotically in the appropriate limit.

### C. k designs

A third development coming from quantum information is the problem of computing flat, structureless averages over the set of unitary matrices, which formally corresponds to the integration over the unitary group with the Haar measure. As already mentioned in the introduction, this notion has many practical applications in quantum information, but sampling unitaries directly from the Haar measure is unpractical, and it is convenient to find simpler ensembles that can approximate the uniform measure to some degree.

The basic idea had been introduced in the context of integrals over the d-dimensional sphere $S^d$ [16]: It is to find a set of points (k design) such that for all polynomials of degree k the averages over the set coincide with the full flat average. Reference [18] extended the notion of k design to the unitary group U(D) (D being the dimension) by introducing unitary designs. A set $\mathcal{E}$ of unitaries U is a unitary k design if the averages of all polynomials of degree k of the matrix elements of U taken over the set $\mathcal{E}$ and over all the Haar ensemble coincide. This is a measure of how good our set $\mathcal{E}$ of unitaries is as a proxy for the whole group. The same definition can be given in terms of a different object: the superoperator known as k-fold quantum channel

$$ \Phi_{\mathcal{E}}^{(k)}(A_1 \otimes \dots \otimes A_k) = \mathbb{E}_{U \sim \mathcal{E}} [U^{\dagger \otimes k} A_1 \otimes \dots \otimes A_k U^{\otimes k}], $$
which quantifies how the combined action of $U$ and $U^†$ acts on observables defined on $k$ copies of the Hilbert space, i.e., the $k$ replicas in the left side in Fig. 1. Or equivalently this can be written as, with summation conversion:

$$ \left\langle a_1, \dots, a_k \left| \Phi_{\mathcal{E}}^{(k)} (A^1 \otimes \dots \otimes A^k) \right| b_1, \dots, b_k \right\rangle = \mathbb{E} [U_{a_1 \bar{a}_1}^\dagger \dots U_{a_k \bar{a}_k}^\dagger U_{\bar{b}_1 b_1} \dots U_{\bar{b}_k b_k}] A_{\bar{a}_1 \bar{b}_1}^1 \dots A_{\bar{a}_k \bar{b}_k}^k . \quad (3) $$

In this language, $\mathcal{E}$ is a $k$ design if $\Phi_{\mathcal{E}}^{(k)}$ coincides with the Haar $k$-fold channel $\Phi_{\text{Haar}}^{(k)}$.

In this work, we look at the properties of $k$-fold quantum channels asymptotically, in the limit of large dimension $D$.

### D. Main results

In this work, we show how free probability provides new insights into the study of unitary designs and their connection to the eigenstate thermalization hypothesis. This is based on the definition of $k$-freeness which is obtained by asking that the condition (2) is valid only up to some finite $k$. The main results are pictorially illustrated in Fig. 1 and can be summarized as follows:

(1) Unitary $k$ designs can be rephrased as asymptotic $k$-freeness for any operator. The relation between freeness and designs is the following: Two operators transformed with respect to one another with a randomly chosen element of a $k$ design are asymptotically $k$-free. Specifically, in the limit of large $D$

$$ k \text{ design} \Rightarrow \begin{cases} \kappa_{2k}(A_1^U, B_1, \dots, A_k^U, B_k) = 0 \\ \forall A_i, B_i, \text{ up to } k, \end{cases} \quad (4) $$

where $A_i^U = UA_iU^†$ and $U$ is sampled from a $k$ design. As a consequence, we may use all the tools and combinatorics developed in the free probability literature to study designs. Particularly useful for quantum information are

(1.1) formulas for computing $2k$-OTOC in terms of individual moments of $A_1, \dots, A_k$ and $B_1, \dots, B_k$ separately;

(1.2) expressions for the $k$-fold quantum channel as an expansion for large $D$ in terms of free cumulants.

(2) ETH implies that unitary time evolution $U(t) = e^{-iHt}$ leads at infinite times to freeness for all finite $k$ but only when restricted to ETH physical operators $\mathcal{A}$. This means that, in the limit of long times, physical observables $A(t) = U^†(t)AU(t)$ and $B$ become $k$-free, for all finite $k$, if $A$ and $B$ obey ETH. Specifically,

$$ \text{ETH}_{\mathcal{A}} \underset{\text{large times}}{\longrightarrow} \begin{cases} \kappa_{2k}^{\beta}[A_1(t), B_1, \dots, A_k(t), B_k] = 0 \\ \forall A_i, B_i \in \mathcal{A}, \quad \forall k. \end{cases} \quad (5) $$

This applies to the canonical version of free cumulants at any inverse temperature $\beta$.

(2.1) This result is rigorous upon infinite time average $\lim_{t_{\text{max}} \to \infty} (1/t_{\text{max}}) \int_0^{t_{\text{max}}} f(t) dt$.

(2.2) We define a physical operator-dependent time-scale $t_k$, named free-$k$ time, for which the $2k$th thermal-free cumulant vanishes.

Analogously to the Haar case, more general statements with more times are straightforward. For instance, the limit of large time of the OTOC $\langle A_1(t)B_1 \dots A_k(t)B_k \rangle_{\beta}$ has a simple expression in terms of thermal moments of $A_1, \dots, A_k$ and $B_1, \dots, B_k$ separately. The definition of $(k)$-freeness requires the vanishing of all mixed cumulants. However, in this context, the vanishing of alternating free cumulants is enough. This is known for Haar rotated matrices [97]; for ETH, it follows from the fact that small powers of ETH observables obey ETH.

### III. NOTIONS AND NOTATIONS

This section covers the notation and the known results that serve as the basis for this work. We start in Sec. III A by introducing $k$ channels, and we discuss Haar averages and $k$ designs in Sec. III B. In Sec. III C, we provide a small introduction to free probability. Finally, in Sec. III D, we discuss the full ETH and its relation with free probability.

#### A. $k$-fold quantum channels

Let $\mathcal{H}$ be a $D$-dimensional Hilbert space and $\mathcal{U}(\mathcal{H})$ the set of unitary operators in $\mathcal{H}$. We consider an arbitrary ensemble of unitaries $\mathcal{E} = \{p_j, U_j\}$, where $p_j$ are probability distributions $(\sum_i p_i = 1)$ and $U_j$ are unitary operators. The ensemble average of a function $f(U)$ over $\mathcal{E}$ is given by

$$ \mathbb{E}_{U \sim \mathcal{E}}[f(U)] := \sum_{i} p_i f(U_i). \quad (6) $$

One may work as well with continuous ensembles $\mathcal{E}$, specified by a probability distribution $d\mu(U)$ over the set of unitaries. In this case, it is intended that the definition above should be replaced by

$$ \mathbb{E}_{U \sim \mathcal{E}}[f(U)] := \int d\mu(U)f(U). \quad (7) $$

A natural quantity to study is how a given ensemble acts on the density matrices of the system. An equivalent point of view we adopt focuses instead on the observables, in which case we think of the unitaries as acting on the observables in the Heisenberg picture. Adopting this second point of view, we define the quantum channel $\Phi_{\mathcal{E}}$ associated with the ensemble $\mathcal{E}$ through its action on observables $O \in \mathcal{H}$:

$$ \Phi_{\mathcal{E}}(O) = \mathbb{E}_{U \sim \mathcal{E}}[U^{\dagger} O U]. \quad (8) $$
To further investigate the correlation property of the ensemble $\mathcal{E}$, we introduce $k$ copies of $\mathcal{H}$—often referred to as replicas. In the replicated space, we can think of a more general class of observables $O$, defined as generic operators acting on the replicated space $\mathcal{H}^{\otimes k}$. Given an ensemble $\mathcal{E}$, we can then introduce the $k$-fold quantum channel $\Phi_{\mathcal{E}}^{(k)}(O)$ defined as

$$ \Phi_{\mathcal{E}}^{(k)}(O) = \mathbb{E}_{U \sim \mathcal{E}} [U^{\dagger \otimes k} O U^{\otimes k}]; \quad (9) $$

i.e., first, each replica is evolved with the same unitary $U$ drawn from the ensemble; afterward, the average over $U$ is performed.

The notion of $k$-fold quantum channel is extremely useful from a physical perspective in its own right. An example of a setting where the $k$-fold channel plays an important role is when $U$ describes the time evolution of a noisy system for a given disorder realization. In this case, we might be interested in computing nonlinear functions of the evolved system’s density matrix—e.g., the exponential of its Rényi entropies—for a given disorder realization and average over disorder only afterward, as commonly done in the context of random unitary circuits [98–101] (see also, e.g., Ref. [102] for a review), noisy fermion transport [89,90,92,103–109], or the Sachdev-Ye-Kitaev model [110–113]. Another application closer to the following discussion is the computation of the average $2k$-OTOC [later defined in Eq. (48)], which can be simply computed in terms of the $k$-fold channel [63,80,114].

## B. The Haar channel and $k$ designs

Among all choices of ensembles, one of the simplest ones, particularly relevant to quantum information, is the uniform ensemble. This is the ensemble induced by the Haar measure $d\mu_{\text{Haar}}$ over the unitary group. The Haar measure is uniquely characterized by the fact that it is left- and right-invariant under multiplication by a unitary and that it is normalized. In other words, for any fixed $V \in \mathcal{U}(\mathcal{H})$, the measure over $U$ satisfies

$$ d\mu_{\text{Haar}}(VU) = d\mu_{\text{Haar}}(UV) = d\mu_{\text{Haar}}(U), \quad (10) $$

$$ \int d\mu_{\text{Haar}}(U) = 1, \quad (11) $$

respectively, due to left-and-right invariance and normalization. The Haar measure can then be used to define the Haar $k$-fold channel combining Eqs. (7) and (9):

$$ \Phi_{\text{Haar}}^{(k)}(O) = \int d\mu_{\text{Haar}}(U) U^{\dagger \otimes k} O U^{\otimes k}. \quad (12) $$

Because of its properties, the Haar channel is useful in many applications in quantum information [1–15]. In practice, however, many applications utilize the Haar

channels only up to a finite value of $k$. Furthermore, in many practical setups, it is difficult to directly sample unitary gates from the Haar ensemble at large $D$. For these reasons, it is convenient to find different ensembles that can reproduce the properties of the Haar $k$-fold channel at least up to some finite $k$.

This has inspired the notion of $k$ design. An ensemble $\mathcal{E}$ is said to be a $k$ design if

$$ \Phi_{\mathcal{E}}^{(k)}(O) = \Phi_{\text{Haar}}^{(k)}(O). \quad (13) $$

Intuitively, a unitary $k$ design is indistinguishable from the Haar ensemble by studying up to the $k$th moment of $\mathcal{E}$. That is, for any polynomial $p_k$ containing up to $k$th powers of $U$ and $U^\dagger$, then

$$ \mathbb{E}_{U \sim \text{Haar}}[p_k(U)] = \mathbb{E}_{U \sim \mathcal{E}}[p_k(U)]. \quad (14) $$

By this definition, if an ensemble is a $k$ design, then it is also a $(k-1)$ design and, in general, a $k'$ design for all $k' < k$.

Furthermore, the action of the Haar $k$ channel can be computed explicitly using the Schur-Weyl duality (see, e.g., Sec. II in Ref. [63] for an elementary introduction). The underlying idea is that, for any $O$, $\Phi_{\text{Haar}}^{k}(O)$ commutes with any unitary operator of the form $V^{\otimes k}$ with $V \in \mathcal{U}(\mathcal{H})$. From this, it can be shown that it needs to be a linear superposition of operators acting as permutations of the $k$ replicas in $\mathcal{H}^{\otimes k}$.

To be more explicit, we introduce a representation of the permutation group $S_k$ onto the $k$-replica Hilbert space $\mathcal{H}^{\otimes k}$. A permutation $\alpha \in S_k$, specified as $\alpha = [\alpha(1), \dots, \alpha(k)]$, acts on $\mathcal{H}^{\otimes k}$ through the permutation operator $P_\alpha$ as follows:

$$ \begin{aligned} P_{\alpha}(|i_1\rangle \otimes |i_2\rangle \otimes \dots \otimes |i_k\rangle) \\ = |\alpha(i_1)\rangle \otimes |\alpha(i_2)\rangle \otimes \dots \otimes |\alpha(i_k)\rangle \end{aligned} \quad (15) $$

with $|i_l\rangle \in \mathcal{H}$ for $1 \le l \le k$. From Schur-Weyl duality, it can be shown that

$$ \Phi_{\text{Haar}}^{(k)}(O) = \sum_{\alpha \beta \in S_k} \mathrm{W}_{\alpha, \beta}(D) \mathrm{Tr}(P_{\beta}O) P_{\alpha^{-1}}, \quad (16) $$

where $\mathrm{W}_{\alpha, \beta}(D) = \mathrm{W}_{\alpha^{-1}, \beta}(D)$ is the unitary Weingarten matrix, given by the inverse of the Gram matrix [115] [116,117]

$$ Q_{\alpha, \beta} = \mathrm{Tr}(P_{\alpha}^{-1} P_{\beta}) = D^{\#(\alpha^{-1}\beta)}, \quad (17) $$

and $\#\alpha$ counts the number of cycles in the permutation $\alpha$ [118]. While the above result provides a direct way to compute $\Phi_{\text{Haar}}^{(k)}(O)$ for finite $k$, the general expressions obtained in this way quickly become unwieldy as $k$ increases. One first simplification can be obtained by
considering the large-$D$ limit. For example, in this limit, the first three orders for $O = A^{\otimes k}$ read

$$ \Phi_{\text{Haar}}^{(1)}(A) = \langle A \rangle P_{(1)}, \quad (18a) $$

$$ \Phi_{\text{Haar}}^{(2)}(A^{\otimes 2}) = (\langle A \rangle^2 + o(1)) P_{(1,2)} + (\langle A^2 \rangle - \langle A \rangle^2 + o(1)) \frac{P_{(2,1)}}{D}, \quad (18b) $$

$$ \begin{aligned} \Phi_{\text{Haar}}^{(3)}(A^{\otimes 3}) &= (\langle A \rangle^3 + o(1)) P_{(1,2,3)} \\ &\quad + [\langle A \rangle (\langle A^2 \rangle - \langle A \rangle^2) + o(1)] \frac{P_{(2,1,3)} + P_{(1,3,2)} + P_{(3,2,1)}}{D} \\ &\quad + [\langle A^3 \rangle - 3\langle A \rangle \langle A^2 \rangle + 2\langle A \rangle^3 + o(1)] \frac{P_{(2,3,1)} + P_{(3,1,2)}}{D^2}, \end{aligned} \quad (18c) $$

where $o(1)$ denotes a vanishing number for $D \to \infty$.

We have also defined the normalized expectation [119] for $A \in \mathcal{H}$:

$$ \langle A \rangle = \frac{\text{Tr}(A)}{D}, \quad (19) $$

i.e., the expectation value of A in the fully mixed state.

Before proceeding, we remark that, in the expansion in Eq. (18) above, we neglect subleading corrections to the coefficient of each permutation operator but retain all permutation operators, independently of the power of $D$ suppressing their effect. As we see in the next section, this turns out to be the correct way of expanding the output of the channel for our purposes, as each permutation operator ends up contributing an $O(1)$ amount to, e.g., the $k$-OTOC.

From the above expression, we also begin to see that each $k$-fold channel contains information about $A$ not present in the $(k-1)$-fold channel. This new information appears in the last term of each of the expressions above, i.e., as the coefficient of the nontrivial cyclic permutations. To better understand this structure for arbitrary $k$, we turn to free probability.

### C. Free probability

Free probability—or “noncommutative probability”— can be thought of as the generalization of classical probability to noncommutative random variables, where the concept of “freeness” extends the one of “independence.” Let us discuss some basic properties; for all the details, we refer to Refs. [79,97,120].

For this purpose, we consider two random matrices $A$ and $B$, which could be drawn from different ensembles. We then introduce the expectation value

$$ \langle \bullet \rangle_{\mathcal{E}} = \lim_{D \to \infty} \frac{1}{D} \mathbb{E}_{\mathcal{E}}[\text{Tr}(\bullet)], \quad (20) $$

where here $\mathcal{E}$ represents a generic random matrix ensemble [121]. Note the slight abuse of notation, since we also denote $\langle \bullet \rangle = \text{Tr}(\bullet)/D$ when the argument does not involve any random matrix. Before proceeding, we comment on a technical assumption: The ensembles of $A$ and $B$ must be defined as a function of $D$; this can be done naturally in the context of many-body physics, where local (or few-body) operators can be naturally extended to act on a larger Hilbert space. In free probability, it is further required that the moments $\langle A^k \rangle_{\mathcal{E}}$ exist for all $k$ (and similarly for $B$). Again, this assumption is naturally satisfied for local (or few-body) operators embedded in a many-body system. Finally, an operator with zero expectation value is said to be centered.

In terms of the expectation value above, we can define the notion of freeness for $A$ and $B$. These are said to be free if [77,96]

$$ \left\langle (A^{n_1} - \langle A^{n_1} \rangle_{\mathcal{E}})(B^{m_1} - \langle B^{m_1} \rangle_{\mathcal{E}}) \cdots \times (A^{n_k} - \langle A^{n_k} \rangle_{\mathcal{E}})(B^{m_k} - \langle B^{m_k} \rangle_{\mathcal{E}}) \right\rangle_{\mathcal{E}} = 0 \quad (21) $$

for all $n_1, \dots, m_k \ge 1$.

The product of the terms above is said to be alternating, since consecutive powers of neither $A$ nor $B$ appear. The definition of freeness can then be enunciated as “$A$ and $B$ are free if the alternating product of centered elements is centered.” The definition can be extended to different sequences of matrices $\{A_1, A_2, \dots\}$ and $\{B_1, B_2, \dots\}$. Let us remark that the definition of freeness always depends on the expectation value $\langle \bullet \rangle_{\mathcal{E}}$; cf. Eq. (20).

This definition is not particularly transparent, but it has to be understood as a rule for computing mixed moments of noncommuting variables from knowledge of the moments of individual variables. For instance, if $A$ and $B$ are two free random matrices, we can simplify the expectation value of the product

$$ \langle ABAB \rangle_{\mathcal{E}} = \langle A^2 \rangle_{\mathcal{E}} \langle B \rangle_{\mathcal{E}}^2 + \langle B^2 \rangle_{\mathcal{E}} \langle A \rangle_{\mathcal{E}}^2 - \langle B \rangle_{\mathcal{E}}^2 \langle A \rangle_{\mathcal{E}}^2, \quad (22) $$

as can be seen by expanding the products and setting to zero all alternating products, e.g., using that

$$ \left\langle (A - \langle A \rangle_{\mathcal{E}})(B - \langle B \rangle_{\mathcal{E}})(A - \langle A \rangle_{\mathcal{E}})(B - \langle B \rangle_{\mathcal{E}}) \right\rangle_{\mathcal{E}} = 0. \quad (23) $$

The definition of freeness can be rephrased in a compact way using the concept of free cumulants and the combinatorial theory of freeness, based on noncrossing partitions, as developed by Speicher [78]. In the following, we first define noncrossing partitions. In terms of these, we introduce the notion of free cumulants, and finally we show how expressing expectation values in terms of free cumulants simplifies the definition of freeness (21).
<figure><img src="image_4.png" /><figcaption>FIG. 2. Graphical representation of six noncrossing partitions appearing in the expansion of $\langle A_1 A_2 A_3 A_4 \rangle$ and their relation to free cumulants through Eqs. (24) and (25). In every picture, each shaded area correspond to a block of the partition encompassing a subset of $\{1, 2, 3, 4\}$. In fact, every noncrossing partition of a set of four elements can be obtained by a rotation of one of the partitions drawn above. The number $[\times n]$ written below the partition indicates how many distinct noncrossing partitions can be obtained in this way. Note that, since the partitions are noncrossing, terms like $\{\{1, 3\}, \{2, 4\}\}$ are not included, as they would induce a crossing of the blocks in the plane. Through Eq. (25), each noncrossing partition $\pi$ is associated with $\kappa_\pi$, a product of free cumulants $\kappa_{l \in \mathbb{Z}}$ over the blocks. This is explicitly reported below the picture of every partition.</figcaption></figure>

A partition of a set $\{1, \dots, n\}$ is a decomposition in blocks that do not overlap and whose union is the whole set. To define noncrossing partitions, we can proceed as described in Fig. 2 for $n = 4$. We fix an order, by arranging $n$ points on a circle and depicting each block as a polygon (including the degenerate cases of one-vertex and two-vertex polygons) whose vertices are a subset of the $n$ points. Partitions in which “blocks do not cross”—i.e., polygons do not overlap—are called noncrossing partitions. The set of all noncrossing partitions of $\{1, \dots, n\}$ is denoted $NC(n)$. The free cumulants $\kappa_n$ are defined implicitly from the moment-cumulant formula, which states that the moments of variables $A_i$ read

$$ \langle A_1 \dots A_n \rangle_{\mathcal{E}} := \sum_{\pi \in NC(n)} \kappa_{\pi}(A_1, \dots, A_n), \quad (24) $$

$$ \kappa_{\pi}(A_1, \dots, A_n) := \prod_{\beta \in \pi} \kappa_{|\beta|}(A_{j_1}, \dots, A_{j_{|\beta|}}), \quad (25) $$

where $\beta = (j_1, \dots, j_\beta)$ denotes the element of a block of the noncrossing partition $\pi$, with $|\beta|$ numbers of elements. See Fig. 2 for a graphical illustration.

Note that the formula above shall indeed be used as a recursive definition of $\kappa_n$, since all the other terms in Eq. (24) are defined by imposing that the same definition holds for $n-1$. For example, in the simple case where $A_i = A$, the $\kappa_n$ for $n \le 4$ is given by

$$ \kappa_1(A) = \langle A \rangle_{\mathcal{E}}, \quad (26a) $$

$$ \kappa_2(A, A) = \langle A^2 \rangle_{\mathcal{E}} - \langle A \rangle_{\mathcal{E}} \langle A \rangle_{\mathcal{E}}, \quad (26b) $$

$$ \kappa_3(A, A, A) = \langle A^3 \rangle_{\mathcal{E}} + 2 \langle A \rangle_{\mathcal{E}}^3 - 3 \langle A \rangle \langle A^2 \rangle_{\mathcal{E}}, \quad (26c) $$

$$ \kappa_4(A, A, A, A) = \langle A^4 \rangle_{\mathcal{E}} - 2 \langle A^2 \rangle_{\mathcal{E}}^2 - 4 \langle A \rangle_{\mathcal{E}} \langle A^3 \rangle_{\mathcal{E}} \\ + 10 \langle A \rangle_{\mathcal{E}}^2 \langle A^2 \rangle_{\mathcal{E}} - 5 \langle A \rangle_{\mathcal{E}}^4. \quad (26d) $$

Let us stress that the definition of free cumulants depends on the expectation value $\langle \bullet \rangle_{\mathcal{E}}$, which in the case of random matrices involves an average over the ensemble. Before proceeding, we point out that, if the sum in Eq. (24) were to run over all partitions, both crossing and noncrossing, then we would recover the standard definition of cumulants. As a consequence of this, standard cumulants and free cumulants coincide for $n < 3$, as all partitions are noncrossing in this case.

The inversion of Eq. (24) can be made systematic using combinatorial tools [97], leading to the so-called Möbius inversion formula:

$$ \kappa_{\pi}(A_1, \dots, A_n) = \sum_{\sigma \in NC(n), \sigma \le \pi} \langle A_1, \dots, A_n \rangle_{\sigma} \mu(\sigma, \pi), \quad (27) $$

where $\langle \bullet \rangle_{\sigma}$ is the product of moments, one for each term of the partition $\sigma$, viz.

$$ \langle A_1, \dots, A_n \rangle_{\sigma} := \prod_{\beta \in \sigma} \langle A_{j_1} \dots A_{j_{\beta}} \rangle_{\mathcal{E}}, \quad (28) $$

$\mu(\sigma, \pi)$ is the so-called Möbius function, and $\sigma \le \pi$ indicates that the sum is restricted to the partitions where each block of $\sigma$ is contained in one of the blocks of $\pi$.

A salient property of free cumulants is that they characterize freeness by the vanishing of the mixed [122] free cumulants, as shown, e.g., in Ref. [120].

In the case of large random matrices, one has that two families $\{A_1, A_2, \dots\}$ and $\{B_1, B_2, \dots\}$ are asymptotically free if all mixed cumulants are zero:

$$ \kappa_n(a_1, a_2, \dots, a_n) = 0 \quad \forall n, \quad (29) $$

where $a_1, \dots, a_n$ are letters of $\{A_1, A_2, \dots\}$ and $\{B_1, B_2, \dots\}$, containing at least a pair $(A_i, B_j)$ for some
i, j. In particular, for any n and for any set of indices $i_1, ..., i_n$, and $j_1, ..., j_n$,

if $\{A_1, A_2, ...\}$ and $\{B_1, B_2, ...\}$ are asymptotically free, then $\kappa_{2n}(A_{i_1}, B_{j_1}, A_{i_2}, B_{j_2}, ..., A_{i_n}, B_{j_n}) = 0 \quad \forall n$. (30)

This works as a rule for the computation of mixed moments [as in Eq. (22)] in terms of moments in which $A_1, ..., A_n$ and $B_1, ..., B_n$ appear separately. Even more nicely, one can find an explicit expression that relates the two. Concretely, if $A_1, ..., A_n$ and $B_1, ..., B_n$ are asymptotically free, then the mixed moments are given by [120]

$$ \langle A_1 B_1 A_2 B_2 \dots A_n B_n \rangle_{\epsilon} = \sum_{\pi \in NC(n)} \kappa_{\pi}(A_1, \dots, A_n) \langle B_1, \dots, B_n \rangle_{\pi^{*}}, \quad (31) $$

where $\pi^*$ is the dual of the partition $\pi$ (also known as Kraweras complement), which is defined below. Graphically (see Fig. 3 for $n=4$), the blocks composing $\pi^*$ are the maximal blocks (polygons) with vertices on “B” that do not cross the blocks of $\pi$. As a simple check, one can verify that, for the simple case of $\langle ABAB \rangle_{\epsilon}$, the general result (31) immediately yields Eq. (22).

While the right-hand side of Eq. (31) is not manifestly symmetric under $A \leftrightarrow B$, one can also exchange the role of A and B. In this case, one obtains

$$ \langle A_1 B_1 A_2 B_2 \dots A_n B_n \rangle_{\epsilon} = \sum_{\pi \in NC(n)} \langle A_1, \dots, A_n \rangle_{K^{-1}(\pi)} \kappa_{\pi}(B_1, \dots, B_n), \quad (32) $$

where $K^{-1}$ is the inverse of the duality transformation $\bullet^*$, i.e., $K^{-1}(\pi^*) = \pi$ [123].

Thus, the definition of freeness in Eq. (21) is greatly simplified when expressed in terms of cumulants as in Eq. (30). First of all, one gets rid of the condition of “centerness” for the variables and relaxes the condition of

“alternating” to “mixed.” Second, there are several properties of free variables which become particularly easy. For instance, the free cumulants of the sum of two-free variables are given by the sum of the free cumulants, e.g.,

$$ \kappa_n(A + B) = \kappa_n(A) + \kappa_n(B), \quad (33) $$

where, for ease of notation, we denote $\kappa_n(X, ..., X)$ with $X$ appearing $n$ times as $\kappa_n(X)$. These nice properties translate over the generating functions, such as the $R$ transform or the $S$ transform; see, e.g., Ref. [120]. We, however, do not deal with these topics in the present manuscript.

While the notion of freeness might seem exotic, many well-known distributions of random matrices produce free matrices asymptotically in the $D \to \infty$ limit. For example, independent Gaussian random matrices [77] are free with respect to each other, and Wigner matrices are free with respect to deterministic ones [120]—in both cases, the statement holds asymptotically in the $D \to \infty$ limit. Another instance of free random matrices that is relevant for us is the following. Given two nonrandom (deterministic) matrices $A$ and $B$ (such that again $\langle A^n \rangle$ is finite for all $n$ and similarly for $B$) and given $U \sim \text{Haar}$, then [77]

$$ U^{\dagger} A U \quad \text{and} \quad B \text{ are asymptotically free as } D \to \infty. \quad (34) $$

This theorem says that unitarily invariant random matrix models are asymptotically free from deterministic matrices. As a consequence, all the relations of Eqs. (30) and (31) hold. The eigenvalue distribution is not changed by the random rotation that, however, changes the relationship between the eigenvectors of $A$ and $B$ which we would have called “generic” or “typical” and we can now refer to as “free.”

## D. The full eigenstate thermalization hypothesis

In a many-body system of $N$ particles described by a Hamiltonian $H$, the information about higher-order

<figure><img src="image_5.png" /><figcaption>FIG. 3. Graphical illustration of some of the terms appearing in Eq. (31) for $n=4$. The shaded blue areas denote noncrossing partitions of the set $\{A_1, A_2, A_3, A_4\}$, as in Fig. 2. To each noncrossing partition $\pi$, one can associate a dual noncrossing partition $\pi^*$ on the set $\{B_1, B_2, B_3, B_4\}$, denoted here as shaded red areas. The contributions of these terms to the sum in Eq. (31) is the product of two terms reported below each partition. With $[\times n]$ we indicate that there are $n$ arrangements of that noncrossing diagram.</figcaption></figure>
correlations is encoded in equilibrium correlators depending on $n$ times, i.e.,

$$ \langle A(t_1)A(t_2)...A(t_n) \rangle^{\beta} \quad \text{with} \quad \langle \cdot \rangle^{\beta} = \frac{1}{Z} \text{Tr}(e^{-\beta H} \cdot) \quad (35) $$

and $Z = \text{Tr}(e^{-\beta H})$. In principle, correlation functions at the order of $n$ contain new information that is not encoded in lower moments, which implies that matrix elements $A_{ij} = \langle E_i | A | E_j \rangle$ of physical observables in the energy eigenbasis need to be correlated. This lead Ref. [64] to introduce a full version of the eigenstate thermalization hypothesis as an ansatz on the statistical properties of the product of $n$ matrix elements $A_{ij}$. Specifically, averages of products with distinct indices $i_1 \neq i_2 \cdots \neq i_n$ read

$$ \overline{A_{i_1 i_2} A_{i_2 i_3} \cdots A_{i_n i_1}} = e^{-(n-1)S(E^+)} F_{e^+}^{(n)}(\omega_{i_1 i_2}, \dots, \omega_{i_{n-1} i_n}), \quad (36) $$

while, with repeated indices, it factorizes in the large-$D$ limit as

$$ \overline{A_{i_1 i_2} \cdots A_{i_{k-1} i_1} A_{i_1 i_{k+1}} \cdots A_{i_n i_1}} = \overline{A_{i_1 i_2} \cdots A_{i_{k-1} i_1} A_{i_1 i_{k+1}} \cdots A_{i_n i_1}}. \quad (37) $$

In Eq. (36), $Ne^+ = E^+ = (E_{i_1} + \cdots + E_{i_n})/q$ is the average energy, $\omega = (\omega_{i_1 i_2}, \dots, \omega_{i_{n-1} i_n})$ with $\omega_{ij} = E_i - E_j$ are $n-1$ energy differences, and $F_{e^+}^{(n)}(\omega)$ is a smooth function of the energy density $e^+ = E^+/N$ and $\omega$. Thanks to the explicit entropic factor, the functions $F_{e^+}^{(n)}(\omega)$ are of the order of one, and they contain all the physical information. For $n=1,2$, one retrieves the standard ETH ansatz [44]; i.e., $F_e^{(1)} = A(e)$ is the microcanonical diagonal value, while $F_e^{(2)} = |f(e, \omega)|^2$ encodes dynamical correlation functions.

The relation between the full ETH ansatz (36) and (37) and the equilibrium correlators (35) becomes more apparent by the use of free probability [76]. Revisiting the definition (24), thermal free cumulants $\kappa_n^\beta$ are defined implicitly by

$$ \langle A(t_1)A(t_2)...A(t_n) \rangle^{\beta} = \sum_{\pi \in NC(n)} \kappa_{\pi}^{\beta}(A(t_1), A(t_2), \dots, A(t_n)), \quad (38) $$

where $\kappa_\pi^\beta$ are products of thermal free cumulants, one for each block of $\pi$. Exactly as in Eq. (24), this is just an implicit definition of cumulants in terms of moments, but the ETH ansatz (36) and (37) implies a particularly simple form for them at the leading order in $N$. Namely, the thermal free cumulants of ETH-obeying systems are given only by summations with distinct indices [76]

$$ \kappa_n^\beta(A(t_1), A(t_2), \dots, A(t_n)) = \frac{1}{Z} \sum_{i_1 \neq i_2 \neq \cdots \neq i_n} e^{-\beta E_i} A_{i_1 i_2} A_{i_2 i_3} \cdots A_{i_n i_1} e^{i t_1 \omega_{i_1 i_2} + \cdots + i t_n \omega_{i_n i_1}}. \quad (39) $$

This result shows that all the non-Gaussianities of the full ETH (36) are encoded precisely in the thermal free cumulants. Explicitly, their Fourier transform $\text{FT}[\cdot] = \int d\omega e^{i\omega \cdot t} \cdot$ reads

$$ \kappa_n^\beta(A(t_1), A(t_2), \dots, A(t_n)) = \text{FT}\left[F_{e_\beta}^{(n)}(\omega)e^{-\beta\omega\cdot\mathcal{L}_n}\right], \quad (40) $$

where the thermal weight with $\mathcal{L}_n = [(n-1/n), \dots, (1/n), 0]$ corresponds to a generalization of the fluctuation-dissipation theorem [124]. The validity of the full ETH ansatz and its free cumulants decomposition (39) has been recently tested numerically in chaotic local many-body systems in Ref. [125]. Conversely, in Bethe-ansatz-integrable models, the structure of the statistics of the matrix elements of operators is more involved, as recently explored in Ref. [126].

## IV. FROM $k$ DESIGN TO $k$-FREENESS

We now present the results in the context of sets of $D \times D$ matrices and the full unitary group, where there is a solid mathematical basis. In the section below, we extrapolate them to the quantum context, for which the rotational invariance is only local.

Let us rephrase the properties of $k$ designs in terms of free probability. We start by considering the case of only two operators $A$ and $B$ on $\mathcal{H}$. As discussed in Sec. III C, if $U$ is sampled according to the Haar distribution, then $A^U = U^\dagger A U$ and $B$ are asymptotically free [77]. This means that $\forall k$ mixed cumulants vanish [79]:

$$ \kappa_{2k}(A^U, B, \dots, A^U, B) = 0, \quad (41) $$

where we recall the ensemble average over $U$ is already included in the definition of the free cumulants (24) via Eq. (20). It now follows that if $U$ is sampled from a $k$ design, then the mixed free cumulants vanish up to the order of $k$; in other words, $A^U$ and $B$ are asymptotically $k$-free.

As such, they are characterized by the properties that characterize free variables up to the order of $k$ but not above. In particular, if $\Phi_\epsilon(A^{\otimes k}) = \Phi_{\text{Haar}}(A^{\otimes k})$, one has (1) the vanishing of the $2k$th free cumulant

$$ \kappa_{2k}(A^U, B, \dots, A^U, B) = 0. \quad (42) $$

This should be seen as an alternative and straightforward way to check if an ensemble is a $k$ design directly from a connected correlation of some observables.

(2) Mixed $2k$-OTOC can be directly computed from the individual moments of $A$ and $B$ from

$$ \langle A^U B A^U B \cdots A^U B \rangle_{\text{Haar}} = \sum_{\pi \in NC(k)} \kappa_{\pi}(A^k) \langle B^k \rangle_{\pi^*}, \quad (43) $$
hence providing a direct generalization to arbitrary $k$ of the results of Ref. [63].

(3) In the large-$D$ limit, the $k$-fold Haar channel can be decomposed in terms of free cumulants:

$$ \Phi_{\text{Haar}}^{(k)}(A^{\otimes k}) = \sum_{\alpha \in S_k} \frac{P_{\alpha^{-1}}}{D^{k-\#\alpha}} (\kappa_{\alpha}(A, \dots, A) + o(1)), \quad (44) $$

where $\kappa_{\alpha}(A, \dots, A)$ are related to the free cumulants $\kappa_{\pi}$. This expression extends Eqs. (18) to general $k$, and it should be understood as an expansion in large $D$ for the coefficients of each permutation operator. More precisely, the subleading corrections to the coefficients, denoted as $o(1)$, scale like $O(D^{-2})$ for large $D$.

In what follows, we build some intuition about these results by discussing them in the general case of multiple operators $A_1, A_2, \dots, A_k$ and $B_1, B_2, \dots, B_k$. From the point of view of channels, this implies that the results apply to arbitrary $A_1 \otimes A_2 \cdots \otimes A_k$ and $B_1 \otimes B_2 \cdots \otimes B_k$.

### A. Haar $k$-fold channel and free cumulants

The structure of $\Phi_{\text{Haar}}^{(k)}$ can be highly simplified using free cumulants. At the leading order in $D$, we find

$$ \Phi_{\text{Haar}}^{(k)}(A_1 \otimes \cdots \otimes A_k) = \sum_{\alpha \in S_k} \frac{P_{\alpha^{-1}}}{D^{k-\#\alpha}} (\kappa_{\alpha}(A_1, \dots, A_k) + o(1)), \quad (45) $$

where we recall that $\#\alpha$ is the number of cycles of $\alpha$.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

$$ \kappa_{\alpha} = \kappa_3(A_1, A_2, A_3) \times \kappa_2(A_4, A_5)\kappa_1(A_6) $$

(a)

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

$$ \kappa_{\alpha} = \kappa_3(A_1, A_3, A_2) \times \kappa_2(A_4, A_5)\kappa_1(A_6) $$

(b)

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

$$ \kappa_{\alpha} = \kappa_3(A_1, A_2, A_6, A_4) \times \kappa_2(A_3, A_5) $$

(c)

FIG. 4. Graphical illustration of some of the coefficients $\kappa_{\alpha}$. Each panel reports a permutation of six elements. The top row reports $\alpha$ both in the form $\alpha = (\alpha(1) \alpha(2) \dots)$ and in a more graphical form as well. Below this, the permutation is drawn as a set of arrows on the circle (cf. Fig. 2). (a) reports a permutation that satisfies both conditions (i) and (ii)—see the main text—therefore, $\kappa_{\alpha}$ can be directly computed as the product of free cumulants over the orbits. (b) reports a permutation that satisfies condition (i) but not (ii), since the cycle in green is clockwise. The coefficient $\kappa_{\alpha}$ can be obtained as the product of free cumulants, but the ordering of matrices must be consistent with the order of the cycle. (c) reports a more general permutation. The coefficient $\kappa_{\alpha}$ in this case can be computed by reordering the elements in such a way that, after reordering, conditions (i) and (ii) are satisfied.

While the derivation of this expression can be found in Appendix B 1, let us discuss here that the coefficients $\kappa_{\alpha}$ can be most easily expressed in terms of free cumulants as follows.

A permutation $\alpha \in S_n$ naturally induces a partition $\pi$ of $n$ elements through the set of its orbits $\pi = \text{Orb}(\alpha)$. An orbit—also referred to as a cycle in the context of permutations—is a subset of the $n$ elements that can be obtained by repeatedly acting with $\alpha$ on a certain element; e.g., $\{1, \alpha(1), \alpha[\alpha(1)], \dots\}$ is an orbit. Varying the initial element, different subsets can be obtained, inducing a partition over $\{1, \dots, n\}$. We begin by considering a permutation that satisfies the following two conditions: (i) Within each of its orbits, the elements are permuted along the direction specified by the cyclic permutation $\gamma = (2, 3, \dots, n, 1)$; in practice, with respect to our graphical notation where the elements are ordered in the counterclockwise direction, this means that, within each cycle, the permutation acts in the counterclockwise direction and (ii) the associated partition $\pi$ formed by its orbits is noncrossing. In Figs. 4(a) and 4(b), respectively, we provide an example of a permutation satisfying and violating condition (i). If $\alpha$ satisfies both conditions, then

$$ \kappa_{\alpha}(A_1, \dots, A_k) = \kappa_{\pi}(A_1, \dots, A_k). \quad (46) $$

On the other hand, if $\alpha$ does not satisfy either of the conditions, this is simply due to the arbitrary ordering given to the operators $A_1, \dots, A_k$ in the tensor product. Therefore, one can simply reorder the operators $A_1, \dots, A_k$ in such a way that the partition $\alpha$ satisfies conditions (i) and (ii) after
reordering. Formally, this is due to the fact that for any permutation $\rho$

$$ \kappa_{\alpha}(A_1, \dots, A_k) = \kappa_{\rho^{-1}\alpha\rho}(A_{\rho(1)}, \dots, A_{\rho(k)}) \quad \forall \rho \in S_k, \quad (47) $$

and, for a given $\alpha$, $\rho$ can be chosen in such a way that, on the right-hand side, $\rho^{-1}\alpha\rho$ satisfies (i) and (ii). The

computation of $\kappa_{\alpha}$ can, therefore, always be reduced to the evaluation of free cumulants through Eqs. (46) and (47). See Fig. 4(c) for an example of such reordering.

We have, thus, expressed Haar channels as a sum of free cumulants as defined in Eq. (24). In the first few orders, this reads

$$ a. \quad \Phi^{(1)}(A_1) = \kappa_1(A_1) $$

$$ b. \quad \Phi^{(2)}(A_1 \otimes A_2) = (\kappa_1(A_1) \kappa_1(A_2) + o(1)) + \left( \times \frac{\kappa_2(A_1 A_2)}{D} + o(1) \right) $$

$$ c. \quad \Phi^{(3)}(A_1 \otimes A_2 \otimes A_3) = (\kappa_1(A_1) \kappa_1(A_2) \kappa_1(A_3) + o(1)) + \left( \times \frac{\kappa_2(A_1 A_2) \kappa_1(A_3)}{D} + \times \frac{\kappa_2(A_2 A_3) \kappa_1(A_1)}{D} + \times \frac{\kappa_2(A_1 A_3) \kappa_1(A_2)}{D} + o(1) \right) \\ + \left( \times \frac{\kappa_3(A_1 A_2 A_3)}{D^2} + \times \frac{\kappa_3(A_1 A_3 A_2)}{D^2} + o(1) \right) $$

where the colored lines indicate the permutation operators: For instance, parallel lines correspond to identities, while two lines that cross to the swap $P_{(2,1)}$ and so on. In the simple case where $A_i = A$, these expressions recover the expansions of $\Phi_{\text{Haar}}^{(k)}$ in Eq. (18) introduced in Sec. III B.

In summary, it is known that $k$th Haar channels are the sum of permutations. Here, we have identified coefficients at the leading order in $D$ as free cumulants. The subleading corrections of the order of $O(D^{-2})$ are irrelevant for the asymptotic properties of mixed OTOCs (as we discuss below).

$$ \langle A_1^U B_1 \dots A_k^U B_k \rangle_{\mathcal{E}} = \frac{1}{D} \text{Tr}(\Phi_{\mathcal{E}}(A_1 \otimes \dots \otimes A_k) B_1 \otimes \dots \otimes B_k P_{\gamma}), \quad (48) $$

where $P_{\gamma}$ is the cyclic permutation over $k$ replicas.

Using the result for the Haar $k$-fold channel (44), in Appendix B 2, we show that $2k$-OTOCs can be written as

$$ \langle A_1^U B_1 \dots A_k^U B_k \rangle_{\text{Haar}} = \sum_{\pi \in \text{NC}(k)} \kappa_{\pi}(A_1, \dots, A_k) \langle B_1 \dots B_k \rangle_{\pi}, \quad (49) $$

## B. Mixed OTOCs of $k$ designs

We are now in the position to review the derivation of Eq. (43) [97] using the above rewriting of the $k$-fold quantum channel. Multipoint $2k$-mixed OTOCs are known to contain information about $k$ channels [63], since they can be reexpressed as

which corresponds to the anticipated result of Eq. (43) generalized to different matrices. Thus, we recognize that this is the same formula of the product of mixed moments of free variables in Eq. (31). The example of an eight-OTOC shown in Fig. 3 reads

$$ \langle A_1^U B_1 A_2^U B_2 A_3^U B_3 A_4^U B_4 \rangle_{\text{Haar}} = \kappa_1(A_1)\kappa_1(A_2)\kappa_1(A_3)\kappa_1(A_4)\langle B_1 B_2 B_3 B_4 \rangle + \kappa_2(A_1, A_2)\kappa_1(A_3)\kappa_1(A_4)\langle B_1 B_2 B_3 \rangle\langle B_4 \rangle + [3] \\ + \kappa_2(A_1, A_2)\kappa_2(A_3, A_4)\langle B_1 B_3 \rangle\langle B_2 \rangle\langle B_4 \rangle + [1] + \kappa_3(A_1, A_2, A_3)\kappa_1(A_4)\langle B_1 B_4 \rangle\langle B_2 \rangle\langle B_3 \rangle + [3] \\ + \kappa_4(A_1, A_2, A_3, A_4)\langle B_1 \rangle\langle B_2 \rangle\langle B_3 \rangle\langle B_4 \rangle, \quad (50) $$

where $[n]$ represents the other combination of observables which can be obtained by the same noncrossing partition; i.e., see the $[\times n]$ in Fig. 3.

## V. DESIGNS FOR ETH CLASS OF OPERATORS

So far, we have introduced the notion of $k$-freeness and showed that for an ensemble $\mathcal{E}$ forming a $k$ design the $2k$-mixed OTOC can be written in terms of free

cumulants, yielding an expression identical to when $A^U$ and $B$ are $k$-free. Here, we show that the same applies to unitary evolution with chaotic Hamiltonians at sufficiently long times.

### A. Time-average ensemble of a Hamiltonian evolution

In this section, we consider the special ensemble $\mathcal{E}$ of unitaries generated by the time evolution under a fixed
chaotic Hamiltonian, i.e., for a fixed Hamiltonian $H$,

$$ \mathcal{E}_{H,t_{\text{Max}}} = \{U(t) = e^{-iHt}\} \quad \text{with} \quad t \in [0, t_{\text{Max}}], \quad (51) $$

with the time $t$ uniformly distributed in the interval. Hence, the average over the ensemble shall be thought of as the standard time average

$$ \mathbb{E}_{t_{\text{Max}}} \{f[U(t)]\} \equiv \frac{1}{t_{\text{Max}}} \int_{0}^{t_{\text{Max}}} f[U(t)] dt. \quad (52) $$

Even if $H$ is chaotic, $\mathcal{E}_{H,t_{\text{Max}}}$ does not approximate a $k$ design at long times. This was already pointed out in Ref. [63]. Formally, this can be seen by explicitly computing the distance induced by the Frobenius norm ($\|A\|_2 = \sqrt{A^\dagger A}$) from the Haar $k$-fold channel $\Phi_{\text{Haar}}^{(k)}$ to the channel associated to $\mathcal{E}_{H,t_{\text{Max}}}$. It turns out that the minimum distance is achieved for large times, at which

$$ \| \Phi_{\text{Haar}}^{(k)} - \Phi_{\mathcal{E}_{H,t_{\text{Max}}}}^{(k)} \|_2 \sim \sqrt{k!} D^{k/2}. \quad (53) $$

The explicit calculation is reported in Sec. III. 3 in Ref. [63]. Intuitively, it is expected that the Hamiltonian evolution is very different from a $k$ design, given that the Hamiltonian hosts $D$ conserved quantities, the projectors on the eigenstates.

Since one deals with Hamiltonian evolution, it is natural to consider the thermal average introduced above, i.e.,

$$ \langle \bullet \rangle^{\beta} = \frac{1}{Z} \text{Tr}(e^{-\beta H \bullet}), \quad (54) $$

which reduces to the usual definition for $\beta = 0$, i.e., $\langle \bullet \rangle^{\beta=0} = \langle \bullet \rangle$.

While the $k$-fold channel acting on a generic set of operators is not close to the Haar $k$ channel [cf. Eq. (53)], one could ask if its action on the restricted set of operators $\mathcal{A}$, namely, those satisfying full ETH, is similar to that of the Haar $k$ channel. We show in the following that this is the case and that Eq. (49) can be recovered from the full ETH formalism of Refs. [64,76,125].

Before discussing long times, let us remark that these arguments put us in a position to understand easily what becomes of these considerations in the *semiclassical limit*. In fact, all operators of which we normally take semiclassical expectations are by nature ETH class in $\mathcal{A}$; they have the requisites of smoothness in the semiclassical limit that guarantee this. What has become then, in the semiclassical limit, of the (vast majority of) operators that are not ETH? A moment's thought tells us that these are operators that have a “bad” semiclassical expression, such as tensor products of two pure states, or operators of the form $A^{c/\hbar}$.

## B. Long-time averages

Within this spirit, we restrict ourselves to $A, B \in \mathcal{A}$, and we show below that ETH implies arbitrary $k$-freeness upon time average, namely,

$$ \mathbb{E}_{t_{\text{Max}} \to \infty} \left[ \kappa_{2k}^{\beta}(A(t), B, A(t), B, \dots, A(t), B) \right] = 0. \quad (55) $$

This result holds for *any finite temperature* defining the thermal average in Eq. (54), and for infinite temperature ($\beta = 0$) it reduces to the factorization of a $k$ design in Eq. (49). As a consequence of Eq. (55), at infinite times $2k$-OTOCs factorize in terms of moments which depend only on $A$ or $B$ separately, namely,

$$ \mathbb{E}_{t_{\text{Max}} \to \infty} \langle A(t) B \dots A(t) B \rangle^{\beta} = \sum_{\pi \in \mathcal{NC}(k)} \kappa_{\pi}^{\beta}(A) \langle B \rangle_{\pi^{*}}^{\beta}, \quad (56) $$

analogously to the factorization property if $A(t)$ and $B$ were $k$-free; cf. Eq. (31). Note that, even if the average over the time ensemble $\mathbb{E}_{t_{\text{Max}} \to \infty}$ appears in a linear fashion, the thermal free cumulants defined here are sufficient to encode nonlinearly the time-evolution ensemble. This is due to the fact that—at leading order in $1/D$—the “ensemble” average of the product of moments is equal to the product of averaged moments for $D \to \infty$, as we further elaborate in Appendix C.

To derive Eq. (55), we use the result of Refs. [76,125], which showed that, when ETH holds, free cumulants can be rewritten in a simple form in terms of matrix elements of the operators $A$ and $B$ in the Hamiltonian eigenbasis:

$$ \begin{aligned} & \kappa_{2k}^{\beta}(A(t), B, A(t), B, \dots, A(t), B) \\ & = \sum_{\substack{i_1 \dots i_k \\ i_1 \dots i_k \text{ distinct}}} \frac{e^{-\beta E_{i_1}}}{Z} A_{i_1 i_1} B_{i_1 i_2} \dots A_{i_{k-1} i_k} B_{i_k i_1} e^{-i t \sum_{l=1}^{k} (E_{i_l} - E_{\tilde{i}_l})}. \end{aligned} \quad (57) $$

Here, the indices $i_l$ and $\tilde{i}_l$ label the eigenstates of the Hamiltonian and $E_i$ their energies—we label as $\tilde{i}$ the eigenstates that correspond to a backward evolution. Equation (57) shall be considered as a part of the ETH ansatz itself; it has been heuristically argued in Ref. [76]. If we were not restricting the sum to a pairwise distinct set of indices, the above would simply be an expansion of the $2k$-OTOC. The nontrivial aspect is that restricting the sum to pairwise distinct indices isolates the contribution of the $2k$th free cumulant to the OTOC [76,125].

When averaging over time in a sufficiently long time window (which we further discuss below), only terms satisfying

$$ \sum_{l=1}^{k} E_{i_l} = \sum_{l=1}^{k} E_{\tilde{i}_l} \quad (58) $$
survive. Furthermore, if the system is chaotic, we can assume that there will not be exact resonances in the spectrum, and the only way to achieve the equality above is to have the set {$i_l$} be a permutation of the set {$\tilde{i}_l$}. Therefore, under these assumptions, all the terms in the expansion (57) vanish upon time integration. Concurrently, the factorization (56) applies to the (2k)-OTOCs of ETH operators in A.

This is one of the main results of this paper: We have shown that ETH implies that chaotic time evolution generates arbitrary k-free observables for all finite k.

Even if these findings are here derived explicitly for OTOC of only two A and B, Eq. (55) is expected to hold also for distinct ETH observables $A_i, B_i \in A$, namely,

$$ \mathbb{E}_{t_{\text{Max}} \to \infty} \left[ \kappa_{2k}^{\beta}(A_1(t), B_1, A_2(t), B_2, \dots, A_k(t), B_k) \right] = 0. \quad (59) $$

In fact, the full ETH in Eq. (36) holds also for distinct operators [76], implying, therefore, the more general results in Eq. (59).

### C. Physical long times

The discussion above tells us that, as $t_{\text{Max}} \to \infty$, the ensemble $\mathcal{E}_{H,t_{\text{Max}}}$ reproduces the Haar ensemble as far as 2k-OTOCs of A ETH-class operators are concerned. Formally, this result is guaranteed to hold when averaging over an exponentially long timescales. Physically, the natural question is whether there is a finite timescale at which the statement above becomes approximately exact not only on average, but at each time. We, thus, define the free-k timescale after which alternating free cumulants vanish, [127] namely,

$$ \kappa_{2k}^{\beta}(A(t), B, A(t), B, \dots, A(t), B) \approx 0 \quad \text{for } t \gtrsim t_k. \quad (60) $$

Intuitively, this is a reasonable expectation in many cases, i.e., far from thermal phase transitions. In ETH, the magnitude of matrix elements is assumed to depend smoothly on the energy of the eigenvectors. Therefore, in Fourier space, $\kappa_{2k}^{\beta}(\omega)$ [cf. Eq. (40)] will be a smooth function of the frequency vectors, and the energy scale over which $\kappa_{2k}^{\beta}(\omega)$ behaves smoothly will determine $1/t_k$.

The idea that operators in chaotic dynamics are organized in a hierarchy of increasing complexity with k is well established; see, e.g., Ref. [63]. Here, we propose to address this complexity via the value of the kth free cumulants. The free-k time $t_k$ shall increase with k and may be compared to the design time. This picture is consistent with the recent numerical investigation of free cumulants in local many-body systems [125], which demonstrated that the $\kappa_4(t)$ is characterized by a much longer timescale to respect $\kappa_2(t)$; in other words, $t_2 \gg t_1$.

The physics is particularly rich in systems in finite spatial dimensions. Here, locality and the finite butterfly velocity introduce extra timescales that can diverge with the linear system size [i.e., scaling like $\ln(D)$], which will modify the smoothness properties of $\kappa_{2k}^{\beta}$ [66]. Understanding how locality enters in this hierarchy of timescales for generic k remains an interesting challenge.

### D. The postulates of Deutsch and Srednicki meet free probability

The basic principles postulated by Deutsch [43] and Srednicki [95] find a nice and compact expression in terms of free probability. Consider a generic Hamiltonian $H$ and perturb it with a noncommuting one $H'$ as

$$ H_{\lambda} = H + c\lambda H', \quad (61) $$

choosing $\lambda$ of $O(1)$, but $c$ as a constant that goes to zero in the thermodynamic limit, e.g., as a power law $N^{-a}$ with $a > 0$. The argument underlying ETH is that $H'$ does not change the macroscopic physics of $H$ but is strong enough with respect to the level spacing so as to act as a “blender” for many nearby levels. On the other hand, for finite $\beta$, we may approximate $e^{-\beta H_{\lambda}} \sim e^{-\beta H} = \rho$ for all $\lambda$ of $O(1)$.

Thus, one may consider the unitary $U_{\lambda}$ which changes the basis between the perturbed and the unperturbed Hamiltonian, namely,

$$ [U_{\lambda}]_{nm} = \langle E_n | E_m^{(\lambda)} \rangle, \quad (62) $$

where $E_n$ (and $E_m^{(\lambda)}$) are the eigenvectors of $H$ (and $H_{\lambda}$) in Eq. (61). The set of unitaries $U_{\lambda}$ could now play the role of our ETH ensemble. Differently from fully random unitaries, this matrix is structured in energy, with matrix elements that decay at large $\omega = E_n - E_m$, as depicted on the right-hand side in Fig. 1. Thus, rotational invariance (that characterizes Haar rotations) emerges here only close to the diagonal— for small frequencies—leading to the concept of local rotational invariance, which has been used to deduce the energy-dependent ETH correlations in Ref. [64]. See also Refs. [76,128] for other discussions on $U_{\lambda}$.

Even if the unitaries $U_{\lambda}$ are clearly not a group, they enjoy another important property. If we consider a set $\{H_{\lambda_1}, \dots, H_{\lambda_k}\}$ with different perturbations $\lambda_i$, the corresponding observables A written in the basis of macroscopically indistinguishable Hamiltonians, $\{A_{\lambda_1}, \dots, A_{\lambda_k}\}$, are mutually free, with measure $\rho$. In other words, different elements of the $U_{\lambda}$ ensemble (62) generate mutually free observables.

Understanding the impact of mutual freeness and other properties of the ETH ensemble on designs is an interesting question which is left as a direction for future work.
## VI. OUTLOOK: THERMALIZATION FOR ENSEMBLES OF STATES

We now discuss how the formalism introduced before may be used to address questions related to thermalization, such as if time-evolved states

$$|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle \quad (63)$$

can be considered as “pseudorandom.”

To give a precise meaning to the term pseudorandom, we consider ensemble of states $S_{|\psi\rangle} = \{p_j, |\psi_i\rangle\}$, where $|\psi_i\rangle$ are random states distributed with probability $p_i$. The information of the full distribution is encoded in the moments

$$\rho_S^{(k)} = \mathbb{E}_{\psi \sim S}[(|\psi\rangle\langle\psi|)^{\otimes k}] = \sum_i p_i (|\psi_i\rangle\langle\psi_i|)^{\otimes k}. \quad (64)$$

An ensemble $S$ is called a quantum state design if its $k$th moment (64) equals the $k$th moment of a random ensemble, i.e., $\rho_S^{(k)} = \rho_{Haar}^{(k)}$. The rationale for the name $\rho_{Haar}^{(k)}$ derives from the fact that random states can always be written as the application of random unitaries $U_i$ to a reference state $|0\rangle$, i.e., $|\psi_i\rangle = U_i|0\rangle$. It is then clear that information about the ensemble of states $S = \{p_i, U_i|0\rangle\}$ is included in the knowledge about the ensemble of unitaries $\mathcal{E} = \{p_i, U_i\}$ that we discuss above. For instance, Eq. (64) can be written in terms of the $k$-fold quantum channel as

$$\rho_{Haar}^{(k)} = \Phi_{Haar}^{(k)}(|0\rangle\langle0|^{\otimes k}). \quad (65)$$

As mentioned above, one would like to understand to what extent unitary evolution generated by a chaotic Hamiltonian can lead to a quantum state design. For this purpose, we define the ensemble of states

$$S_{H,t_{\text{Max}}} = \{|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle\} \text{ with } t \in [0, t_{\text{Max}}] \quad (66)$$

and ask whether at long times it shares properties of a quantum state design. Clearly, $S_{H,t_{\text{Max}}}$ cannot form a quantum-state design, not even for $k=1$, since the time evolution conserves many quantities, such as the projectors onto all its eigenstates, hence the energy.

In this outlook, we consider some simple instances (such as moments or Rényi entropies) to illustrate how free probability can be used to show that quantum state designs are asymptotically equivalent to time-evolved states at long times. For an exhaustive discussion on the topic, we refer to Refs. [32,94,129]. Finally, we conclude by discussing possible connections and differences between $k$-freeness and the recently introduced notion of “deep thermalization.”

## A. Moments of expectation values

We begin discussing moments of observable expectation values, whose information is encoded in the $k$th moment in Eq. (64), i.e.,

$$\mathbb{E}_{\psi \sim S} [\langle\psi|A|\psi\rangle^k] = \text{Tr}(\rho_S^{(k)} A^{\otimes k}). \quad (67)$$

In the case of a quantum-state design, we can now use the free-probability approach outlined in Sec. IV. A straightforward calculation (see Appendix D 1) at the leading order in $1/D$ leads to

$$\begin{aligned} \mathbb{E}_{\psi \sim \text{Haar}} [\langle \psi | A | \psi \rangle^k] &= \text{Tr}(|0\rangle\langle 0|^{\otimes k} \Phi_{\text{Haar}}^{(k)}(A^{\otimes k})) \\ &\approx \kappa_1(A) \dots \kappa_1(A) + \mathcal{O}(D^{-1}). \end{aligned} \quad (68)$$

Hence, the average of the $k$th power of an expectation value is dominated by the $k$th power of the average of the observable. Since the fluctuations are subleading at large $D$, they are not accounted for by our asymptotic approach.

The same happens for chaotic dynamics on $|\psi(t)\rangle$: The time evolution at infinite times leads to asymptotic state designs on the diagonal ensemble identified by the initial state $|\psi(0)\rangle$, i.e.,

$$\langle\bullet\rangle_{\text{diag}} = \text{Tr}(\rho_{\text{diag}}\bullet) \text{ with } \rho_{\text{diag}} = \sum_i |c_i|^2 |E_i\rangle\langle E_i|, \quad (69)$$

where $c_i = \langle\psi(0)|E_i\rangle$ are the overlaps with the initial state, which, for generic initial states, can be assumed to be well peaked around an extensive energy $E_0$, with subextensive energy fluctuations.

The $k$th moment of observables on the time-dependent state can be rewritten as

$$\begin{aligned} \langle\psi(t)|A|\psi(t)\rangle^k &= \sum_{\substack{i_1 \cdots i_k \\ i_1 \cdots i_k}} c_{i_1} c_{i_1}^* \cdots c_{i_k} c_{i_k}^* A_{i_1 i_1} \cdots A_{i_k i_k} e^{-it \sum_{l=1}^k (E_{i_l} - E_{i_{l+1}})}. \end{aligned} \quad (70)$$

Averaging over infinite times and using the nonresonance condition discussed in Sec. V, in Appendix D 2 we show that ETH implies

$$\mathbb{E}_{t_{\text{Max}} \to \infty} [\langle\psi(t)|A|\psi(t)\rangle^k] \approx \langle A \rangle_{\text{diag}}^k + \mathcal{O}(e^{-S}). \quad (71)$$

In summary, we have discussed how the $k$th moments of physical observables are asymptotically the same for quantum state designs and time-evolved states at large times: They are given by the $k$th power of the first cumulant, i.e., $[\kappa_1(A)]^k = \langle A \rangle^k$, the latter being evaluated over the diagonal ensemble for Hamiltonian time evolution.
Note that Eqs. (68) and (71) are immediately generalized to the fluctuations of the product of expectation values of different observables, such as $\langle A^{(1)} \rangle \dots \langle A^{(k)} \rangle$.

It is worth noting that, for different energy distributions, one might compare the time-average ensemble with more sophisticated choices of the random ensemble; see, e.g., Ref. [129].

## B. Rényi entropies

We now discuss the Rényi entropies across a bipartition between two complementary regions X and $\bar{X}$. The total Hilbert space is factorized $\mathcal{H} = \mathcal{H}_X \otimes \mathcal{H}_{\bar{X}}$ and has dimensions dim $\mathcal{H}_X = D_X$ and dim $\mathcal{H}_{\bar{X}} = D_{\bar{X}} = D/D_X$. Given a pure state $|\psi\rangle \in \mathcal{H} = \mathcal{H}_X \otimes \mathcal{H}_{\bar{X}}$ and its reduced density matrix $\rho_X = \text{Tr}_{\bar{X}}|\Psi\rangle\langle\Psi|$, we are interested in the higher-order purity

$$R^{(k)}(|\Psi\rangle) = \text{Tr}(\rho_X^k). \qquad (72)$$

Using the partial cycle permutation $S_X$, defined as the operator that cyclically permutes the k replicas over the subsystem X, i.e., $S_X|i_1, \bar{i}_1; i_2, \bar{i}_2; \dots i_k, \bar{i}_k\rangle = |i_2, \bar{i}_1; i_3, \bar{i}_2; \dots i_1, \bar{i}_k\rangle$, [130] Eq. (72) can be rewritten as

$$R^{(k)}(|\Psi\rangle) = \langle\Psi|^{\otimes k} S_X |\Psi\rangle^{\otimes k}; \qquad (73)$$

see, e.g., Ref. [94]. The kth Rényi entropy is defined from Eq. (72) as $S^k = (1/1 - k) \ln R^{(k)}$ [131].

Free probability, despite being asymptotic, allows one to show that quantum state designs and a time-average ensemble possess the same Rényi entropy at the leading order, including the so-called Page correction, also known as “residual entropy” [132].

The Rényi entropy of random Haar states has been intensively studied (see, e.g., Refs. [133–137]), including via free probability tools [135–137]. In Appendix D 3, we provide an alternative derivation using asymptotic free probability. The kth purity in Eq. (72) reads

$$\mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^k)] = \sum_{\pi \in NC(k)} \frac{1}{D_X^{k-\#\pi}} \frac{1}{D_{\bar{X}}^{k-\#\pi}} + \mathcal{O}(D^{-1}), \quad (74)$$

where $\#\pi$ and $\#\pi^*$ are the number of blocks of the partition $\pi$ and of its dual, respectively. This expression is symmetric upon exchanging X and $\bar{X}$. Moreover, when X and $\bar{X}$ become comparable, i.e., $D_X = D_{\bar{X}} = \sqrt{D}$, it gives

$$\mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^k)] = \frac{C_k}{D_X^{k-1}} + \mathcal{O}(D^{-1}), \quad (75)$$

where $C_k$ is the Catalan number which counts the non-crossing partitions of $NC(k)$. Upon taking the logarithm,

this contribution leads to the so-called residual entropy of the Rényi entropies: $S^{(k)} = \ln D_X - (\ln C_k/k - 1) + \mathcal{O}[(\ln D)^{-1}]$ [135–137].

Let us comment now on how one may apply these findings to the time-averaged entropy of the nonequilibrium state, defined as

$$\overline{R^{(k)}} = \frac{1}{t_{\text{Max}}} \int_{0}^{t_{\text{Max}}} dt R^{(k)}[|\psi(t)\rangle]. \qquad (76)$$

Using Eq. (73), we compute the long-time average for large enough $t_{\text{Max}}$, by expanding in the energy eigenbasis. With the assumption of the absence of exact resonances (cf. Sec. V), the average over $\exp[-it(E_{j_1} - E_{\bar{j}_1} + \dots + E_{j_k} - E_{\bar{j}_k})]$ selects only permutations of the energy indices. This leads to

$$\mathbb{E}_{t_{\text{Max}} \to \infty} \{ R^{(k)}[|\psi(t)\rangle] \} = \sum_{i_1, \dots, i_k} |c_{i_1}|^2 \dots |c_{i_k}|^2 \sum_{\beta \in S_k} \langle E_{i_1} \dots E_{i_k} | P_{\beta} S_X | E_{i_1} \dots E_{i_k} \rangle, \quad (77)$$

where $c_i = \langle E_i | \psi(0) \rangle$ is the diagonal ensemble in Eq. (69). To make further progress, one may consider the case where the coefficients $|c_j|^2 = 1/D$ for all j. Here, one has

$$\overline{R^{(k)}} = D^{-k} \sum_{\beta \in S_k} \text{Tr}(P_{\beta} S_X) \qquad (78)$$

[cf. Eq. (D9) in Appendix D 3], thus recovering the result (74) for a quantum state design.

This result is expected to apply to generic states, where $|c_j|^2 \sim 1/D$ averaged over small energy windows, also known as infinite-temperature initial states. In the spirit of ETH, we expect $\langle E_{i_1} \dots E_{i_k} | P_{\beta} S_X | E_{i_1} \dots E_{i_k} \rangle$ to be independent of small variations in the set of indices $i_1, \dots, i_k$ but rather to have only a smooth dependence on the energies $E_{i_1}, \dots, E_{i_k}$. In this way, small-scale fluctuations of $|c_j|^2$ do not modify the result in Eq. (78).

More generally, it would be interesting to better understand how this result generalizes to finite-temperature states. There, we expect energy conservation to affect the average value of $R^{(k)}$. In fact, it is known that the entropy of a finite-energy state has additional corrections due to the curvature of the density of states; cf. Refs. [138,139].

## C. k-freeness versus “deep thermalization”

The notion of k-freeness introduced in this paper is one way to extend the usual concept of thermalization. Indeed, the thermalization of the local density matrix corresponds to the case $k = 1$ discussed in this manuscript. This is, however, not the only possible generalization. For example,
recent advances in quantum simulators stimulated the idea of deep thermalization [33–42].

For concreteness, consider a lattice system, which is bipartite into two regions S and R. The regular notion of thermalization involves the reduced density matrix on S obtained after tracing over R, i.e., $\rho_S = \text{Tr}_R|\psi(t)\rangle\langle\psi(t)|$. To define deep thermalization, rather than simply tracing over R, we perform projective measurements onto each site of R. In this way, the state over R is collapsed onto a product state, specified by the array of measurement outcomes $\mathbf{m}$. In general, we then have a random state of the form $|\phi_{t,m}\rangle_S \otimes |\mathbf{m}\rangle_R$ selected with probability $P_{\mathbf{m}} = \langle\psi(t)|\mathbf{m}\rangle_{RR} |\mathbf{m}|\psi(t)\rangle$. The ensemble of states $|\phi_{t,m}\rangle_S$ and their probability $P_{\mathbf{m}}$ forms the so-called projected ensemble. Deep thermalization concerns whether this ensemble forms a quantum state design.

k-freeness and deep thermalization clearly involve different quantities and concepts. Nonetheless, one might ask if there is a common physical origin. While this is an important direction for future work, we here discuss an example based on the results of Ref. [140] that hints toward a negative answer: While k-freeness has the same properties of unitary dynamics, such as locality, deep thermalization is intrinsically nonlocal.

One way to address this question quantitatively is by studying if there is a relation between deep-thermalization times $\tilde{\tau}_k$ and free-k times $t_k$. Firstly, one must introduce a notion of locality within k-freeness. In the bipartition setting above, this can be most easily done by redefining $t_k$ as the maximum time overall set of possible operators A and B with support in S to achieve the condition (60). In the specific setting of Ref. [140], where S sits between two R regions, we consider how $t_k$ and $\tilde{\tau}_k$, obtained for open boundary conditions (BCs), can be modified upon adding a coupling between the first and final R sites. In this case, changing the BCs cannot affect $t_k$. The Lieb-Robinson bound [141] implies that $A(t)$ has support over a finite region and, thus, is independent of the BCs in the thermodynamic limit. Therefore, if $t_k$ is finite in the thermodynamic limit, it must be BC independent. On the other hand, this is not the case for deep thermalization at $k \ge 2$. Even in the thermodynamic limit and at finite times, the projected ensemble can drastically depend on the BCs. A remarkable example is the one presented in Ref. [140], where deep thermalization takes place at finite times, but nonetheless using periodic BCs halves $\tilde{\tau}_k$ with respect to the open BCs case. The origin of the nonlocality is the combination of the measurements and the sharing (classical communication) of the measurement outcomes over the whole system.

While this specific example indicates that, in general, there will not be a connection between k-freeness and deep thermalization, it would nonetheless be interesting to understand if there are limiting cases where these two phenomena share the same mechanism and timescales.

## VII. CONCLUSIONS AND PERSPECTIVES

In this paper, we propose to consider designs from the point of view of free probability. In this language, an operator A and $UAU^\dagger$, the latter rotated with an element of a k design, are k-free. This implies that the 2k-OTOCs factorize and the mixed free cumulants vanish.

A by-product of our work is the new notion of k-freeness, which, in turn, is accompanied by all the tools of free probability now available in connection with k designs.

The results of this work have broader consequences for the chaotic dynamics of many-body systems.

(i) ETH equipped with free probability allows one to identify the universality of late-time chaotic dynamics in the factorization of higher moments into products of free cumulants, each of them with observables evaluated at the same time. This acts as a sort of generalization of the mixing condition [142]. We have shown that unitary evolution always generates k free variables at k sufficiently long separated times. In other words, the notion of freeness gives a compact and operational way to understand the effect of chaotic time evolution by vanishing the appropriate free cumulants of all orders.

(ii) While other approaches to this issue are based on frame potentials or operator's norms, this point of view shifts the focus to moments and cumulants of physical observables. This must be considered as being in the original spirit of ETH, which takes us from postulating an ensemble of Hamiltonians having the same macroscopic physics, yet are random “in the small” (Berry and Deutsch), to the ensuing thermalization properties. This is a natural step, as when one goes from a classical description in terms of trajectories in phase space to one in terms of the evolution of physical observables.

(iii) This approach provides practical tools to verify when time evolution yields a k design, studying at what times the free cumulants vanish:

$$ \kappa_{2k}[A, A(t), \dots, A, A(t)] = 0. $$

This condition is easily implementable on classical simulations of quantum dynamics and may be experimentally achievable via some of the techniques to measure OTOCs; see, e.g., Refs. [143–147].

(iv) The definition of k-freeness does not rely on the existence of eigenstates, and it naturally applies to arbitrarily time-dependent and noisy evolutions (such as random circuits). As such, the vanishing-free cumulants above may be easily computed in a much larger class of quantum evolutions whose chaotic properties we may thus probe.

(v) Furthermore, this line of thought naturally includes the ETH ideas on equilibration at the level of
ensembles of states. Our free probability approach shows that the chaotic time evolution of pure states (equipped with the appropriate diagonal ensemble) is equivalent to quantum state designs from the point of view of observable fluctuations or entanglement entropies.

As mentioned at the outset, all the new notions introduced in this work may, and perhaps sometimes need to, be relaxed by replacing equalities with equalities up to a predefined precision $\epsilon$. Our approach is based on asymptotic definitions; however, a rigorous definition of an $\epsilon$-approximated $k$-freeness would require a thorough discussion of the finite $D$ corrections. While in the case of full matrices, further insights in this direction could be gained using the so-called “second-order freeness” [148], in the ETH case it remains an open question what to use as an effective dimension. Finite-size corrections may be relevant also in the discussion of fluctuations of multitime observables such as $\mathbb{E}[\langle A(t)A \rangle^k]$.

Physically, an important question concerns the implications of locality on the meaning and scaling of the free-$k$ times $t_k$. In this direction, it would be very interesting to understand what are the implications of $k$-freeness on the study of transport.

This problem is linked to the relationship between $k$-freeness and deep thermalization. As discussed in Sec. VI C, these two are generally distinct issues, but it would be worthwhile to investigate if there are limiting cases where the two phenomena share mechanisms or timescales.

Finally, an interesting but challenging question concerns the fate of late-time $k$-freeness when ETH does not apply, such as for different forms of ergodicity breaking or in the case of integrable systems [126].

## ACKNOWLEDGMENTS

We thank Laura Foini for insightful comments and collaboration on related works. S. P. thanks Xhek Turkeshi for the initial insights on the Haar ensemble and for useful discussions. We thank Takato Yoshimura for valuable feedback on the first version of Sec. IV. We also thank Andrea de Luca, Jacopo de Nardis, Denis Bernard, and Wen Wei Ho for valuable discussions on the topic. This work was completed during the workshop “Open QMBP” held at the Institut Pascal at Université Paris-Saclay, made possible with the support of the program “Investissements d’avenir” ANR-11-IDEX-0003-01. S. P. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Projektnummer 277101999—TRR 183 (Project No. B02) and under Germany’s Excellence Strategy—Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1-390534769.

## DATA AVAILABILITY

No data were created or analyzed in this study.

## APPENDIX A: EXAMPLES OF HAAR CHANNELS

The Haar quantum channel can be written explicitly by using the Schur-Weyl duality [149]:

$$ \Phi_{\text{Haar}}^{(k)}(O) = \sum_{\alpha\beta \in S_k} \text{Wg}_{\alpha,\beta}(D) \text{Tr}(P_{\beta}O) P_{\alpha^{-1}}, \quad (\text{A1}) $$

where $\text{Wg}_{\alpha,\beta}(D)$ is the Weingarten matrix, defined as the inverse of the Gram matrix [150] [116,117]:

$$ Q_{\alpha,\beta} = \text{Tr}(P_{\alpha}P_{\beta}) = D^{\#(\alpha^{-1}\beta)}, \quad (\text{A2}) $$

where $\#(\alpha)$ counts the number of cycles of the permutation $\alpha$. Let us show explicitly the first few $k$.

$$ (i) \quad k = 1. \quad S_1 = \{(1)\}, \quad P_{\alpha^{-1}} = \{\mathbb{I}\}, \quad Q_{11} = d, \\ \text{Wg}_{11}(D) = (1/D), \text{ and} $$

$$ \Phi_{\text{Haar}}^{(1)}(O) = \frac{\text{Tr}(O)}{D} \mathbb{I} = \langle A \rangle \mathbb{I}, \quad (\text{A3}) $$

where we define the normalized expectation; for $A \in \mathcal{H}$,

$$ \langle A \rangle = \frac{\text{Tr}(A)}{D}. \quad (\text{A4}) $$

$$ (ii) \quad k = 2. \quad S_2 = \{(12), (21)\}, \quad P_{\alpha^{-1}} = \{P_{12}, P_{21}\}, \\ Q = \begin{pmatrix} D^2 & D \\ D & D^2 \end{pmatrix}, \quad \text{Wg} = \frac{1}{D(D^2-1)} \begin{pmatrix} D & -1 \\ -1 & D \end{pmatrix}, $$

which leads to

$$ \Phi_{\text{Haar}}^{(2)}(O) = \frac{D\text{Tr}(O) - \text{Tr}(OP_{12})}{D(D^2-1)} P_{12} + \frac{D\text{Tr}(OP_{12}) - \text{Tr}(O)}{D(D^2-1)} P_{21}. \quad (\text{A5}) $$

If $O = A \otimes A$, then we have

$$ \begin{aligned} \text{Tr}(O) &= \text{Tr}(A^{\otimes 2}) = \text{Tr}(A)^2 = D^2 \langle A \rangle^2, \\ \text{Tr}(OP_{21}) &= \text{Tr}(A^{\otimes 2}P_{21}) = \text{Tr}(A^2) = D \langle A^2 \rangle, \end{aligned} $$

which leads to

$$ \begin{aligned} \Phi_{\text{Haar}}^{(2)}(A^{\otimes 2}) &= \frac{D^2 \langle A \rangle^2 - \langle A^2 \rangle}{(D^2 - 1)} P_{12} + \frac{D \langle A^2 \rangle - D \langle A \rangle^2}{(D^2 - 1)} P_{21} \\ &= \underset{D \to \infty}{\langle A \rangle^2} P_{12} + \frac{\langle A^2 \rangle - \langle A \rangle^2}{D} P_{21}. \end{aligned} \quad (\text{A6}) $$
(iii) $k=3$.—The permutation group is $S_3 = \{(123), (213), (321), (132), (231), (312)\}$,

$$Q = \begin{pmatrix} D^3 & D^2 & D^2 & D^2 & D & D \\ D^2 & D^3 & D & D & D^2 & D^2 \\ D^2 & D & D^3 & D & D^2 & D^2 \\ D^2 & D & D & D^2 & D^2 & D^2 \\ D & D^2 & D^2 & D^2 & D^3 & D \\ D & D^2 & D^2 & D^2 & D & D^3 \end{pmatrix}, \quad (A7)$$

and

$$W_g = \frac{1}{(D^2 - 1)(D^2 - 4)} \times \begin{pmatrix} 1 - \frac{2}{D} & -1 & -1 & -1 & \frac{2}{D} & \frac{2}{D} \\ -1 & 1 - \frac{2}{D} & \frac{2}{D} & \frac{2}{D} & -1 & -1 \\ -1 & \frac{2}{D} & 1 - \frac{2}{D} & \frac{2}{D} & -1 & -1 \\ -1 & \frac{2}{D} & \frac{2}{D} & 1 - \frac{2}{D} & -1 & -1 \\ \frac{2}{D} & -1 & -1 & -1 & 1 - \frac{2}{D} & \frac{2}{D} \\ \frac{2}{D} & -1 & -1 & -1 & \frac{2}{D} & 1 - \frac{2}{D} \end{pmatrix}, \quad (A8)$$

which for large $D$ results in

$$\begin{aligned} \Phi_{\text{Haar}}^{(3)}(A^{\otimes k}) &=_{D \to \infty} \langle A \rangle^3 P_{123} \\ &\quad + \langle A \rangle (\langle A^2 \rangle - \langle A \rangle^2) \frac{P_{213} + \text{perm}}{D} \\ &\quad + (\langle A^3 \rangle - 3\langle A \rangle\langle A^2 \rangle + 2\langle A \rangle^3) \frac{P_{231} + P_{312}}{D^2}. \end{aligned} \quad (A9)$$

It is clear that the new information contained in the $k$th channel that was not in the $k-1$ is associated with the cyclic permutation $P_{12...k}$ and has the form of a connected correlation function. To understand this structure, we need free probability.

# APPENDIX B: FREE PROBABILITY ON $k$ DESIGNS

## 1. Derivation of Eq. (45)

Let us here illustrate the main steps leading to the expansion of the $k$-fold Haar quantum channel in terms of free cumulants.

### a. Scaling with $D$

The starting point to show Eq. (45) is to expand the Weingarten matrix in Eq. (16) to leading order in $D$ [116]:

$$W_{g,\alpha,\beta}(D) \approx \frac{1}{D^{2k-\#(\beta^{-1}\alpha)}} [\mu(\beta, \alpha) + O(D^{-2})]. \quad (B1)$$

To extract the leading contributions, it is necessary to make the $D$ dependence explicit, i.e.,

$$\mathrm{Tr}(P_{\beta} A_1 \otimes A_2 \otimes \cdots \otimes A_k) = D^{\#\beta} \prod_{b \in \mathrm{Orb}(\beta)} \langle A_{j_1} \cdots A_{j_b} \rangle. \quad (B2)$$

In the expression above, $b$ runs over the orbits of $\beta$—in other words, the cycles of $\beta$. It is further intended that $j_1, \dots, j_b$ enumerate the element of the orbit in the appropriate order, i.e., $A_{j_2} = A_{\beta(j_1)}$, $A_{j_3} = A_{\beta(j_2)}$, and so on, until $A_{j_1} = A_{\beta(j_{|b|})}$.

The two equations lead to

$$\Phi_{\text{Haar}}^{(k)}(A_1 \otimes \cdots \otimes A_k) = \sum_{\alpha \in S_k} P_{\alpha^{-1}} \sum_{\beta \in S_k} \frac{\mu(\beta, \alpha)}{D^{2k-\#(\beta^{-1}\alpha)-\#\beta}} \times \prod_{b \in \text{Orb}(\beta)} \langle A_{j_1} \cdots A_{j_b} \rangle + O(D^{-2}). \quad (B3)$$

Therefore, focusing on a given $P_{\alpha^{-1}}$, its coefficient $\kappa_{\alpha}$ at leading order is determined by the permutations $\beta$ that minimize $2k - \#(\beta^{-1}\alpha) - \#\beta$. To understand which permutations $\beta$ give the leading contributions, we recall the notion of length $\ell(\sigma)$ of a permutation $\sigma$, i.e., the least number of two-element transpositions necessary to factorize $\sigma$ [97]. This is useful for our purposes, since the length is directly related to the number of cycles through $\ell(\sigma) = k - \#\sigma$. Furthermore, the length can be used to define a distance on $S_k$ $d(\alpha, \beta) = \ell(\alpha^{-1}\beta)$. In particular, we use that $d$ satisfies the triangular inequality. In our case, considering the distance between $\alpha$ and the identical permutation $\mathbb{I}$, we have that $\ell$ satisfies

$$\ell(\beta) + \ell(\beta^{-1}\alpha) \geq \ell(\alpha). \quad (B4)$$

From this, it follows that the permutations $\beta$ giving the leading contributions are the ones saturating the triangular inequality, i.e.,

$$\ell(\beta) + \ell(\beta^{-1}\alpha) = \ell(\alpha), \quad (B5)$$

whereas all the other permutations can be neglected. With reference to the distance $d$ induced by $\ell$, the set of permutations $\beta$ satisfying Eq. (B5) is often called geodesic; in this case, the geodesic between $\alpha$ and the identity $\mathbb{I}$. Equation (B5) constrains the sum over $\beta$ in Eq. (B3) and leads to a coefficient proportional to $1/D^{\ell(\alpha)} = 1/D^{k-\#\alpha}$ with corrections $O(D^{-2})$. Putting everything together, we then obtain Eqs. (45) with

$$\begin{aligned} \kappa_{\alpha}(A_1, A_2, \dots, A_k) = & \sum_{\beta: \ell(\beta) + \ell(\beta^{-1}\alpha) = \ell(\alpha)} \mu(\beta, \alpha) \\ & \times \prod_{b \in \text{Orb}(\beta)} \langle A_{j_1} \cdots A_{j_b} \rangle. \end{aligned} \quad (B6)$$
## b. The coefficients $\kappa_\alpha$ and noncrossing partitions

The next question we want to ask is whether there is a more explicit way of writing $\kappa_\alpha$. In fact, we now show that $\kappa_\alpha$ is simply related to the free cumulants $\kappa_\pi$ previously introduced.

This connection hinges on an important combinatorial property of the distance $d(\alpha, \beta) = \ell(\alpha^{-1}\beta)$, leading to an embedding of noncrossing partitions into the permutation group. We begin by restricting ourselves to the case where the permutation $\alpha$ satisfies conditions (i) and (ii) listed above Eq. (46). Namely, the cycles of $\alpha$ induce a non-crossing partition $\pi$ and permute the elements of each cycle counterclockwise with respect to our drawing convention. Then, the permutations $\beta$ contributing to $\kappa_\alpha$, i.e., those lying on the geodesic between the identity $I$ and $\alpha$, can be characterized as the set of permutations $\beta$ such that the orbits of $\beta$ induce a noncrossing partition $\sigma$, $\sigma \le \pi$ (each block of $\sigma$ is contained in a block of $\pi$), and all cycles of $\beta$ permute their elements counterclockwise with respect to our convention. For further details and a proof, we refer to Lecture 23 in Ref. [97]. This immediately leads to the result of Eq. (45) and to identification of $\kappa_\pi(A_1, A_2, \dots, A_k)$ with the free cumulants.

The previous statement can be extended to permutations which have a different cyclic order. This can be done by noting through Eq. (47), which we also report below for the reader's convenience:

$$ \kappa_{\alpha}(A_1 \dots A_k) = \kappa_{\rho^{-1}\alpha\rho}(A_{\rho(1)} \dots A_{\rho(k)}) \quad \forall \rho \in S_k. \quad (B7) $$

To show this, we observe that

$$ \Phi_{\text{Haar}}^{(k)}(A_1 \otimes \dots \otimes A_k) = P_{\rho}^{\dagger} \Phi_{\text{Haar}}^{(k)}(A_{\rho(1)} \otimes \dots \otimes A_{\rho(k)}) P_{\rho}. \quad (B8) $$

Expanding both sides using Eq. (45) then yields Eq. (47). Equation (47) is particularly useful, since for any $\alpha$ we can always choose $\rho$ such that $\alpha' = \rho^{-1}\alpha\rho$ has noncrossing cycles that permute their elements counterclockwise [see, e.g., Fig. 4(c) for an example]. [In fact two permutations, $\alpha$ and $\alpha'$, can always be related by conjugation $\rho^{-1}\alpha\rho$ if the lengths of their cycles are the same.] Hence, with an appropriate choice of $\alpha'$, we can reduce ourselves to the case studied in the previous paragraph. Correspondingly, the coefficient $\kappa_\alpha$ in Eq. (45) is given by the free cumulant of the associated permutation:

$$ \begin{aligned} \kappa_{\alpha}(A_1, A_2, \dots, A_k) &= \sum_{\sigma \le \pi} \mu(\pi, \sigma) \langle A_{\rho(1)}, A_{\rho(2)}, \dots, A_{\rho(k)} \rangle_{\sigma} \\ &= \kappa_{\pi}(A_{\rho(1)}, A_{\rho(2)}, \dots, A_{\rho(k)}). \end{aligned} \quad (B9) $$

## 2. Derivation of Eq. (49)

We now want to show how 2k-OTOCs can be written in the form (43) using the result for the Haar k-fold channel (44). For this purpose, we start from the 2k-OTOC in Eq. (48), and we plug in the result of the Haar kth channel of Eqs. (44) and (B6) as

$$ \begin{aligned} & \langle A_1^U B_1 \dots A_k^U B_k \rangle_{\text{Haar}} \\ &= \frac{1}{D} \mathrm{Tr}(\Phi_{\text{Haar}}(A_1 \otimes \dots \otimes A_k) B_1 \otimes \dots \otimes B_k P_{\gamma}) \end{aligned} \quad (B10a) $$

$$ = \frac{1}{D} \sum_{\alpha} \kappa_{\alpha}(A_1, \dots, A_k) \frac{\mathrm{Tr}(P_{\gamma} P_{\alpha^{-1}} B_1 \otimes \dots \otimes B_k)}{D^{k-\#\alpha}} \quad (B10b) $$

$$ = \sum_{\alpha} \frac{\kappa_{\alpha}(A_1, \dots, A_k)}{D^{k+1-\#\alpha-\#(\gamma\alpha^{-1})}} \prod_{b \in \mathrm{Orb}(\gamma\alpha^{-1})} \langle B_{j_1} \dots B_{j_b} \rangle, \quad (B10c) $$

where from the first to the second line we use Eq. (45) and then we make the D dependence explicit using Eq. (B2) for $\gamma\alpha^{-1}$. For large D, the leading contribution to this summation is given by the permutations $\alpha$ which minimize $k+1-\#\alpha-\#(\gamma\alpha^{-1})$. Using that the length of the cyclic permutation is $\ell(\gamma) = k-1$, it follows that the leading contribution satisfy

$$ \ell(\alpha) + \ell(\gamma\alpha^{-1}) = \ell(\gamma), \quad (B11) $$

which implies that, in the large-D limit, Eq. (B10c) can be rewritten as

$$ \begin{aligned} & \langle A_1^U B_1 \dots A_k^U B_k \rangle_{\text{Haar}} \\ &= \sum_{\alpha: \ell(\alpha) + \ell(\gamma\alpha^{-1}) = \ell(\gamma)} \kappa_{\alpha}(A_1, \dots, A_k) \langle B_1 \dots B_k \rangle_{\gamma\alpha^{-1}}, \end{aligned} \quad (B12) $$

where we use the definition of $\langle B_1 \dots B_k \rangle_{\gamma\alpha^{-1}}$ as the products of moments, one for each block of the permutation, e.g., given $\langle B_1 \dots B_k \rangle_{\gamma\alpha^{-1}} = \prod_{d \in \mathrm{Orb}(\gamma\alpha^{-1})} \langle \prod_{j \in d} B_j \rangle$.

The condition (B11) implies that $\alpha$ lies on the geodesic between the identity and the cyclic permutation $\gamma$. Hence, due to the embedding of the noncrossing partitions discussed above, $\mathrm{Orb}(\alpha) = \pi$ is a noncrossing partition. Furthermore, one can show that $\gamma\pi^{-1} = \pi^*$ corresponds to the dual of noncrossing partition of $\pi$ [97]. Note that, from this discussion, it is clear that all permutations $\alpha$ contributing to OTOCs satisfy the conditions (i) and (ii) discussed above in Sec. IVA. Hence, together with Eq. (46), this leads to Eq. (49).
## APPENDIX C: TIME AVERAGE OF THE CUMULANT VERSUS CUMULANT OF THE TIME-AVERAGED ENSEMBLE

In Sec. V, we show that the average of the free cumulant defined according to Eq. (39) is zero in the long-time limit; cf. Eq. (55). However, the time average of the cumulant $\kappa_{2k}^{\beta}$ is distinct, in principle, from the cumulant of the time average ensemble $\mathcal{E}_{H,t_{\text{Max}}}$, which in this appendix we simply denote as $\mathcal{E}$. For the purpose of this appendix, we explicitly denote the ensemble in the latter with $\kappa_{2k}^{\beta,\mathcal{E}}$. Furthermore, to show that $\mathcal{E}_{H,t_{\text{Max}}}$ acts as a k design on ETH-class operators, we need to show that $\kappa_{2k}^{\beta,\mathcal{E}}$ vanishes in the long-time limit. In this appendix, we demonstrate that, if the average of $\kappa_{2k}^{\beta}$ is zero (55), then also $\kappa_{2k}^{\beta,\mathcal{E}}$ is zero, up to corrections of size $1/D$.

We first begin by showing that $\kappa_{2k}^{\beta,\mathcal{E}}$ and the time average of $\kappa_{2k}^{\beta}$ do not coincide, in general. The simplest example of this is obtained for $k=2$. Writing explicitly $\kappa_{2k}^{\beta,\mathcal{E}}$ using Eq. (26), we have

$$ \begin{aligned} & \kappa_4^{\beta,\mathcal{E}}[A(t)BA(t)B] \\ &= \mathbb{E}_{t_{\text{Max}}} [\langle A(t)BA(t)B \rangle^{\beta}] - 2(\mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}])^2 \\ & \quad - 2\mathbb{E}_{t_{\text{Max}}} [\langle A(t)BA(t) \rangle^{\beta}] \langle B \rangle^{\beta} - 2\mathbb{E}_{t_{\text{Max}}} [\langle BA(t)B \rangle^{\beta}] \langle A \rangle^{\beta} \\ & \quad + 10\mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}] \langle A \rangle^{\beta} \langle B \rangle^{\beta} - 5(\langle A \rangle^{\beta})^2 (\langle B \rangle^{\beta})^2. \end{aligned} \quad (C1) $$

These are the same terms that contribute to $\mathbb{E}_{t_{\text{Max}}} [\kappa_4^{\beta}[A(t)BA(t)B]]$, except that the term in the third line of Eq. (C1) would be replaced by

$$ \mathbb{E}_{t_{\text{Max}}} [(\langle A(t)B \rangle^{\beta})^2] \neq (\mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}])^2. \quad (C2) $$

Nonetheless, at long times and at leading order in $1/D$, we argue that the time average can be split as

$$ \begin{aligned} & \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \dots \rangle^{\beta} \dots \langle A(t)B \dots \rangle^{\beta}] \\ &= \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \dots \rangle^{\beta}] \dots \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \dots \rangle^{\beta}] + O(1/D), \end{aligned} \quad (C3) $$

where each expectation value might contain an arbitrary number of terms.

For concreteness, we show this explicitly in the simplest possible case, i.e.,

$$ \begin{aligned} & \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta} \langle A(t)B \rangle^{\beta}] \\ &= \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}] \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}] + O(1/D), \end{aligned} \quad (C4) $$

but the argument extends straightforwardly to the most general case. We expand the rhs in the eigenbasis of the Hamiltonian as

$$ \begin{aligned} & \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta} \langle A(t)B \rangle^{\beta}] \\ &= \sum_{\substack{i,l \\ i,j}} \frac{e^{-\beta(E_i+E_j)}}{Z^2} A_{\bar{i}l} B_{\bar{i}\bar{l}} A_{\bar{j}l} B_{\bar{j}\bar{l}} \mathbb{E}_{t_{\text{Max}}} [e^{-it(E_i-E_i+E_j-E_j)}]. \end{aligned} \quad (C5) $$

Proceeding as in Sec. V, at long times, assuming that there are not exact resonances in the spectrum,

$$ \mathbb{E}_{t_{\text{Max}}} [e^{-it(E_i - E_i + E_j - E_j)}] = \delta_{i,\bar{i}}\delta_{j,\bar{j}} + (\delta_{i,\bar{j}}\delta_{j,\bar{i}} - \delta_{i,\bar{i}}\delta_{j,\bar{j}}\delta_{i,j}). \quad (C6) $$

Plugging the first of the two terms in Eq. (C5), we obtain exactly the factorized term $\mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}] \mathbb{E}_{t_{\text{Max}}} [\langle A(t)B \rangle^{\beta}]$. Instead, plugging the second term in Eq. (C5), we obtain a “crossing” partition contribution, namely,

$$ \sum_{i \neq \bar{i}} \frac{e^{-\beta(E_i + E_j)}}{Z^2} (A_{\bar{i}i} B_{\bar{i}\bar{i}})(A_{\bar{i}\bar{i}} B_{\bar{i}\bar{i}}) \sim 1/D, \quad (C7) $$

which is suppressed as $1/D$ [76], as can be seen using the standard ETH scaling $A_{\bar{i}i}B_{\bar{i}\bar{i}} \sim e^{-S(E_+)}$, with $E_+ = (E_i + E_{\bar{i}})/2$. Therefore, the approximation in Eq. (C4) holds.

The same arguments work, in general, at all orders. Expanding the rhs of Eq. (C3) in terms of eigenstates, there is the average of oscillating phases as in Eq. (C6). The average of oscillating phases always yields pairing schemes where the energy levels are equal within a given expectation value, but otherwise unconstrained, and terms constraining levels from different expectation values coincide. The former terms give the factorized contribution in the rhs of Eq. (C3). The latter is instead always to be suppressed by factors of $1/D$ or higher.

Finally, let us remark that the “asymptotic” factorization between ensemble averages of products and products of ensemble averages [cf. Eq. (C4)] holds in the same way for the Haar ensemble.

## APPENDIX D: FREE PROBABILITY ON ENSEMBLES OF STATES

In this appendix, we report the main steps in the application of free probability theory to ensembles of states discussed in Sec. VI.

### 1. Calculation of Eq. (68)

We derive the leading contribution to the moments of expectation values using free probability. First of all,
we rewrite Eqs. (65) and (67) using the cyclicity of the trace and invariance of the Haar distribution under $U \mapsto U^\dagger$. In this way, we obtain [151]

$$ 
\begin{aligned}
\mathbb{E}_{\psi \sim \text{Haar}} [\langle \psi | A | \psi \rangle^k] \\
&= \text{Tr}(|0\rangle\langle 0|^{\otimes k} \Phi_{\text{Haar}}^{(k)}(A^{\otimes k})) \\
&\approx \sum_{\alpha \in S_k} \frac{\kappa_{\alpha}(A, \dots, A)}{D^{k-\#\alpha}} \text{Tr}(P_{\alpha^{-1}}|0\rangle\langle 0|^{\otimes k}) + \mathcal{O}(D^{-2}) \\
&= \sum_{\alpha \in S_k} \frac{\kappa_{\alpha}(A, \dots, A)}{D^{k-\#\alpha}} + \mathcal{O}(D^{-2}) \\
&\approx \kappa_1(A) \dots \kappa_k(A) + \mathcal{O}(D^{-1}),
\end{aligned}
\quad \text{(D1)}
$$

with $\kappa_1(A) = \text{Tr}(A)/D$. Here, from the first to the second line we use Eq. (44) and to the third $\text{Tr}(P_\alpha|0\rangle\langle 0|^{\otimes k}) = 1 \forall \alpha$. Finally, we note that the leading contribution comes from the identity permutation $P_\alpha = \mathbb{I}$, while the other terms are suppressed as $D^{-1}$.

## 2. Derivation of Eq. (70)

We here derive the infinite-time limit of the $k$th moment of observables in Eq. (70). As discussed in Sec. V, when averaging over infinite times, i.e., using Eq. (52), only terms satisfying $\sum_{l=1}^k E_{\bar{i}_l} = \sum_{l=1}^k E_{\bar{i}_l}$ survive, and, using the nonresonance condition, this equality above is to have the set $\{i_l\}$ be a permutation of the set $\{\bar{i}_l\}$. This leads to

$$ 
\begin{aligned}
\mathbb{E}_{t_{\text{Max}} \to \infty} [\langle \psi(t) | A | \psi(t) \rangle^k] \quad & \text{(D2)} \\
&= \sum_{\substack{i_1 \dots i_k \\ \alpha \in S_k}} |c_{i_1}|^2 \dots |c_{i_k}|^2 A_{i_1 \alpha(i_1)} \dots A_{i_k \alpha(i_k)} \\
&\approx \langle A \rangle_{\text{diag}}^k + \mathcal{O}(e^{-S}).
\end{aligned}
\quad \text{(D3)}
$$

To obtain the last line above, one recognizes that the leading contribution to Eq. (D3) is given by the identity permutation:

$$ \left( \sum_{i_1} |c_{i_1}|^2 A_{i_1 i_1} \right)^k = \langle A \rangle_{\text{diag}}^k, \quad \text{(D4)} $$

while, for $\alpha$ different from the identity, one can use ETH to bound the other terms as

$$ 
\begin{aligned}
& \sum_{\substack{i_1, i_k \text{ dist} \\ \alpha \in S_k}} |c_{i_1}|^2 \dots |c_{i_k}|^2 A_{i_1 \alpha(i_1)} \dots A_{i_k \alpha(i_k)} \\
& \le \max A_{i_1 \alpha(i_1)} \dots A_{i_k \alpha(i_k)} \quad \text{(D5)}
\end{aligned}
$$

$$ \propto e^{-\#(\text{OD})S}, \quad \text{(D6)} $$

where $\#(\text{OD})$ counts the number of off-diagonal matrix elements included in the products, viz. $\#(\text{OD}) = \#\{j|\alpha(j) \neq j\}$. Here, $S$ is the minimum of the thermodynamic entropy where $|c_j|$ are non-negligible.

## 3. Rényi entropy for state designs using free probability

We here provide a derivation of Eq. (74) of the main text; see also Refs. [135–137]. The average of the purity $R^{(k)}$ over a state design (65) reads

$$ 
\begin{aligned}
\mathbb{E}_{|\psi\rangle \sim \text{Haar}} [\text{Tr}(\rho_X^k)] \\
&= \text{Tr}(|0\rangle\langle 0|^{\otimes k} \Phi_{\text{Haar}}^{(k)}(S_X)) \\
&= \sum_{\alpha \beta \in S_k} \text{Wg}(\alpha, \beta^{-1}) \text{Tr}(P_{\alpha^{-1}}|0\rangle\langle 0|^{\otimes k}) \text{Tr}(P_\beta S_X).
\end{aligned}
\quad \text{(D7)}
$$

In the second line, we use the exact result for the quantum channel in Eq. (16). We recall that the starting point of our free probability approach was the asymptotic expansion for the Weingarten matrix in Eq. (B1):

$$ \text{Wg}(\alpha, \beta^{-1}) \approx \frac{\mu(\alpha\beta^{-1})}{D^{2k-\#\alpha\beta^{-1}}} \approx \frac{\delta_{\alpha\beta}}{D^k} + \mathcal{O}(D^{-(k+1)}), \quad \text{(D8)} $$

where in the second line we use that the leading contribution to this matrix is diagonal in the replicas. While for observables the off-diagonal terms lead to finite results, we now show that the diagonal approximation is sufficient to get the Page correction to the entanglement [132]. We plug Eq. (D8) into Eq. (D7), together with $\text{Tr}(P_{\alpha^{-1}}|0\rangle\langle 0|^{\otimes k}) = 1 \forall \alpha$, and this leads to

$$ 
\begin{aligned}
\mathbb{E}_{\psi \sim \text{Haar}} [\text{Tr}(\rho_X^k)] &\approx \frac{1}{D^k} \sum_{\beta \in S_k} \text{Tr}(P_\beta S_X) \quad \text{(D9)} \\
&= \sum_{\beta \in S_k} \frac{D_X^{\#\gamma\beta^{-1}} D_{\bar{X}}^{\#\beta}}{D^k}, \quad \text{(D10)}
\end{aligned}
$$

where from the first to the second line we use $\text{Tr}(P_\beta S_X) = D_X^{\#\gamma\beta^{-1}} D_{\bar{X}}^{\#\beta}$, with $\gamma$ is the cycle permutation over the $k$ replicas with the same cyclicity as $S_X$. We now substitute in the denominator $D = D_X D_{\bar{X}}$ and recognize that the exponents of $D_{X/\bar{X}}$ are given by the permutation length function $\ell(\sigma) = k - \#\sigma$ (introduced above in Appendix B 1). This leads to

$$ \mathbb{E}_{\psi \sim \text{Haar}} [\text{Tr}(\rho_X^k)] \approx \sum_{\beta \in S_k} \frac{1}{D_X^{\ell(\gamma\beta^{-1})} D_{\bar{X}}^{\ell(\beta)}} \quad \text{(D11)} $$

This expression is maximal when $D_X = D_{\bar{X}} = \sqrt{D}$. Hence, in the large-$D$ limit, the leading contribution to the
summation is given by permutations that saturate the triangular inequality for the length function in Eq. (B4), which in this case reads

$$ \ell(\beta) + \ell(\gamma\beta^{-1}) = \ell(\gamma). $$

As discussed above, this implies that the permutations $\beta$ lie on the geodesic between the identity and the cyclic permutation. Thus, they correspond to noncrossing partitions, and one arrives to the desired result

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^k)] \approx \sum_{\pi \in \text{NC}(k)} \frac{1}{D_X^{k-\#\pi}} \frac{1}{D_{\bar{X}}^{k-\#\pi}}. \quad (\text{D12}) $$

For instance, at $k = 2, 3, 4$ Eq. (D12) reads

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^2)] = \frac{1}{D_X} + \frac{1}{D_{\bar{X}}}, \quad (\text{D13a}) $$

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^3)] = \frac{1}{D_{\bar{X}}^2} + \frac{3}{D_X D_{\bar{X}}} + \frac{1}{D_{\bar{X}}^2}, \quad (\text{D13b}) $$

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^4)] = \frac{1}{D_{\bar{X}}^3} + \frac{6}{D_{\bar{X}}^2 D_{\bar{X}}} + \frac{6}{D_X D_{\bar{X}}^2} + \frac{1}{D_{\bar{X}}^3}. \quad (\text{D13c}) $$

As expected, these expressions are symmetric upon exchanging $X$ and $\bar{X}$. Moreover, when $X$ and $\bar{X}$ are the same ($D_X = D_{\bar{X}} = \sqrt{D}$), it is well known that the Rényi entropies of random states read

$$ S^{(k)} = \frac{1}{1-k} \ln \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^k)] \\ = \ln D_X + c_{\text{Page}}^{(k)} + O[(\ln D)^{-1}], $$

where $c_{\text{Page}}^{(k)}$ is a number of the order of one, known as the Page correction or residual entropy; see, e.g., Refs. [132–137]. In this case, Eq. (D12) can be rewritten as

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^k)] \approx \frac{1}{\sqrt{D}^{k-1}} \sum_{\pi \in \text{NC}(k)} \frac{C_k}{\sqrt{D}^{k-1}}, \quad (\text{D14}) $$

where one uses that $\#\pi + \#\pi^* = k + 1$ and $C_k$ denotes the number of noncrossing partitions of $\text{NC}(k)$, known as the Catalan number. In the first few orders, this reads

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^2)] = \frac{2}{\sqrt{D}}, \quad (\text{D15a}) $$

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^3)] = \frac{5}{\sqrt{D^2}}, \quad (\text{D15b}) $$

$$ \mathbb{E}_{\psi \sim \text{Haar}}[\text{Tr}(\rho_X^4)] = \frac{14}{\sqrt{D^3}}. \quad (\text{D15c}) $$

Finally, Eq. (D14) allows us to identify the Page correction:

$$ c_{\text{Page}}^{(k)} = \frac{\ln C_k}{1-k}. \quad (\text{D16}) $$

[1] Hoi Fung Chau, *Unconditionally secure key distribution in higher dimensions by depolarization*, IEEE Trans. Inf. Theory **51**, 1451 (2005).

[2] Joseph M. Renes, Robin Blume-Kohout, Andrew J. Scott, and Carlton M. Caves, *Symmetric informationally complete quantum measurements*, J. Math. Phys. (N.Y.) **45**, 2171 (2004).

[3] A. Hayashi, T. Hashimoto, and M. Horibe, *Reexamination of optimal quantum state estimation of pure states*, Phys. Rev. A **72**, 032325 (2005).

[4] Sofyan Iblisdir and Jérémie Roland, *Optimal finite measurements and gauss quadratures*, Phys. Lett. A **358**, 368 (2006).

[5] Joseph Emerson, Robert Alicki, and Karol Życzkowski, *Scalable noise estimation with random unitary operators*, J. Opt. B **7**, S347 (2005).

[6] S. J. van Enk and C. W. J. Beenakker, *Measuring Trρ$^n$ on single copies of ρ using random measurements*, Phys. Rev. Lett. **108**, 110503 (2012).

[7] A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, *Rényi entropies from random quenches in atomic Hubbard and spin models*, Phys. Rev. Lett. **120**, 050406 (2018).

[8] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, *Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states*, Phys. Rev. A **99**, 052323 (2019).

[9] Hsin-Yuan Huang, Richard Kueng, and John Preskill, *Predicting many properties of a quantum system from very few measurements*, Nat. Phys. **16**, 1050 (2020).

[10] Marco Paini and Amir Kalev, *An approximate description of quantum states*, arXiv:1910.10543.

[11] Joshua Morris and Borivoje Dakić, *Selective quantum state tomography*, arXiv:1909.05880.

[12] Lukas Knips, Jan Dziewior, Waldemar Kłobus, Wiesław Laskowski, Tomasz Paterek, Peter J. Shadbolt, Harald Weinfurter, and Jasmin D. A. Meinecke, *Multipartite entanglement analysis from random correlations*, npj Quantum Inf. **6**, 51 (2020).

[13] Andreas Ketterer, Nikolai Wyderka, and Otfried Gühne, *Characterizing multipartite entanglement with moments of random correlations*, Phys. Rev. Lett. **122**, 120505 (2019).

[14] Andreas Elben, Steven T. Flammia, Hsin-Yuan Huang, Richard Kueng, John Preskill, Benoît Vermersch, and Peter Zoller, *The randomized measurement toolbox*, Nat. Rev. Phys. **5**, 9 (2022).

[15] Richard A. Low, *Pseudo-randomness and learning in quantum computation*, arXiv:1006.5227.

[16] Philippe Delsarte, Jean-Marie Goethals, and Johan Jacob Seidel, *Spherical codes and designs*, in Geometry and Combinatorics (Elsevier, New York, 1991), pp. 68–93.
[17] Christoph Dankert, *Efficient simulation of random quantum states and operators*, arXiv:quant-ph/0512217.

[18] Christoph Dankert, Richard Cleve, Joseph Emerson, and Etera Livine, *Exact and approximate unitary 2-designs and their application to fidelity estimation*, Phys. Rev. A **80**, 012304 (2009).

[19] Aram W. Harrow and Richard A. Low, *Random quantum circuits are approximate 2-designs*, Commun. Math. Phys. **291**, 257 (2009).

[20] Igor Tuche Diniz and Daniel Jonathan, *Comment on "Random quantum circuits are approximate 2-designs" by A. W. Harrow and R. A. Low (Commun. Math. Phys. 291, 257-302 (2009))*, Commun. Math. Phys. **304**, 281 (2011).

[21] Fernando G. S. L. Brandão, Aram W. Harrow, and Michał Horodecki, *Efficient quantum pseudorandomness*, Phys. Rev. Lett. **116**, 170502 (2016).

[22] Jonas Haferkamp, *Random quantum circuits are approximate unitary t-designs in depth $o(n^{5+o(1)})$, Quantum **6**, 795 (2022).

[23] Yoshifumi Nakata, Christoph Hirche, Masato Koashi, and Andreas Winter, *Efficient quantum pseudorandomness with nearly time-independent Hamiltonian dynamics*, Phys. Rev. X **7**, 021006 (2017).

[24] Aram W. Harrow and Saeed Mehraban, *Approximate unitary t-designs by short random quantum circuits using nearest-neighbor and long-range gates*, Commun. Math. Phys. **401**, 1531 (2023).

[25] Nicholas Hunter-Jones, *Unitary designs from statistical mechanics in random quantum circuits*, arXiv:1905.12053.

[26] Fernando G. S. L. Brandão, Wissam Chemissany, Nicholas Hunter-Jones, Richard Kueng, and John Preskill, *Models of quantum complexity growth*, PRX Quantum **2**, 030316 (2021).

[27] Jonas Haferkamp, Philippe Faist, Naga B. T. Kothakonda, Jens Eisert, and Nicole Yunger Halpern, *Linear growth of quantum circuit complexity*, Nat. Phys. **18**, 528 (2022).

[28] Joseph Emerson, Etera Livine, and Seth Lloyd, *Convergence conditions for random quantum circuits*, Phys. Rev. A **72**, 060302(R) (2005).

[29] Sumner N. Hearth, Michael O. Flynn, Anushya Chandran, and Chris R. Laumann, *Efficient local classical shadow tomography with number conservation*, Phys. Rev. Lett. **133**, 060802 (2024).

[30] Shao-Kai Jian, Gregory Bentsen, and Brian Swingle, *Linear growth of circuit complexity from Brownian dynamics*, arXiv:2206.14205.

[31] Anastasiia Tiutiakina, Andrea De Luca, and Jacopo De Nardis, *Frame potential of Brownian SYK model of Majorana and Dirac fermions*, arXiv:2306.11160.

[32] Saúl Pilatowsky-Cameo, Ceren B. Dag, Wen Wei Ho, and Soonwon Choi, *Complete Hilbert-space ergodicity in quantum dynamics of generalized Fibonacci drives*, Phys. Rev. Lett. **131**, 250401 (2023).

[33] Jordan S. Cotler, Daniel K. Mark, Hsin-Yuan Huang, Felipe Hernandez, Joonhee Choi, Adam L. Shaw, Manuel Endres, and Soonwon Choi, *Emergent quantum state designs from individual many-body wave functions*, PRX Quantum **4**, 010311 (2023).

[34] Joonhee Choi, Adam L. Shaw, Ivaylo S. Madjarov, Xin Xie, Ran Finkelstein, 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, *Preparing random states and benchmarking with many-body quantum chaos*, Nature (London) **613**, 468 (2023).

[35] Wen Wei Ho and Soonwon Choi, *Exact emergent quantum state designs from quantum chaotic dynamics*, Phys. Rev. Lett. **128**, 060601 (2022).

[36] Pieter W. Claeys and Austen Lamacraft, *Emergent quantum state designs and biunitarity in dual-unitary circuit dynamics*, Quantum **6**, 738 (2022).

[37] Matteo Ippoliti and Wen Wei Ho, *Solvable model of deep thermalization with distinct design times*, Quantum **6**, 886 (2022).

[38] Matteo Ippoliti and Wen Wei Ho, *Dynamical purification and the emergence of quantum state designs from the projected ensemble*, PRX Quantum **4**, 030322 (2023).

[39] Maxime Lucas, Lorenzo Piroli, Jacopo De Nardis, and Andrea De Luca, *Generalized deep thermalization for free fermions*, Phys. Rev. A **107**, 032215 (2023).

[40] Tanmay Bhore, Jean-Yves Desaules, and Zlatko Papić, *Deep thermalization in constrained quantum systems*, Phys. Rev. B **108**, 104317 (2023).

[41] Max McGinley and Michele Fava, *Shadow tomography from emergent state designs in analog quantum simulators*, Phys. Rev. Lett. **131**, 160601 (2023).

[42] Matteo Ippoliti and Wen Wei Ho, *Dynamical purification and the emergence of quantum state designs from the projected ensemble*, PRX Quantum **4**, 030322 (2023).

[43] J. M. Deutsch, *Quantum statistical mechanics in a closed system*, Phys. Rev. A **43**, 2046 (1991).

[44] Mark Srednicki, *The approach to thermal equilibrium in quantized chaotic systems*, J. Phys. A **32**, 1163 (1999).

[45] I. Aleksandr and Akovlevich Khinchin, *Mathematical Foundations of Statistical Mechanics* (Courier Corporation, New York, 1949).

[46] Tomaz Prosen, *Ergodic properties of a generic nonintegrable quantum many-body system in the thermodynamic limit*, Phys. Rev. E **60**, 3949 (1999).

[47] Marcos Rigol, Vanja Dunjko, and Maxim Olshanii, *Thermalization and its mechanism for generic isolated quantum systems*, Nature (London) **452**, 854 (2008).

[48] Giulio Biroli, Corinna Kollath, and Andreas M. Läuchli, *Effect of rare fluctuations on the thermalization of isolated quantum systems*, Phys. Rev. Lett. **105**, 250401 (2010).

[49] Anatoli Polkovnikov, Krishnendu Sengupta, Alessandro Silva, and Mukund Vengalattore, *"Colloquium": Non-equilibrium dynamics of closed interacting quantum systems*, Rev. Mod. Phys. **83**, 863 (2011).

[50] Tatsuhiko N. Ikeda, Yu Watanabe, and Masahito Ueda, *Finite-size scaling analysis of the eigenstate thermalization hypothesis in a one-dimensional interacting Bose gas*, Phys. Rev. E **87**, 012125 (2013).

[51] Robin Steinigeweg, Jacek Herbrych, and Peter Prelovšek, *Eigenstate thermalization within isolated spin-chain systems*, Phys. Rev. E **87**, 012118 (2013).
[52] Vincenzo Alba, *Eigenstate thermalization hypothesis and integrability in quantum spin chains*, Phys. Rev. B **91**, 155123 (2015).

[53] Wouter Beugeling, Roderich Moessner, and Masudul Haque, *Off-diagonal matrix elements of local operators in many-body quantum systems*, Phys. Rev. E **91**, 012144 (2015).

[54] David J. Luitz, *Long tail distributions near the many-body localization transition*, Phys. Rev. B **93**, 134201 (2016).

[55] David J. Luitz and Yevgeny Bar Lev, *Anomalous thermalization in ergodic systems*, Phys. Rev. Lett. **117**, 170404 (2016).

[56] Tyler LeBlond and Marcos Rigol, *Eigenstate thermalization for observables that break Hamiltonian symmetries and its counterpart in interacting integrable systems*, Phys. Rev. E **102**, 062113 (2020).

[57] Marlon Brenes, Tyler LeBlond, John Goold, and Marcos Rigol, *Eigenstate thermalization in a locally perturbed integrable system*, Phys. Rev. Lett. **125**, 070605 (2020).

[58] Felix Fritzsch and Tomaz Prosen, *Eigenstate thermalization in dual-unitary quantum circuits: Asymptotics of spectral functions*, Phys. Rev. E **103**, 062133 (2021).

[59] S. J. Garratt and J. T. Chalker, *Local pairing of Feynman histories in many-body Floquet models*, Phys. Rev. X **11**, 021051 (2021).

[60] A. Kitaev, *talk given at the Fundamental Physics Prize Symposium (2014)*, https://www.youtube.com/watch?v=OQ9qN8j7EZI.

[61] Juan Maldacena, Stephen H. Shenker, and Douglas Stanford, *A bound on chaos*, J. High Energy Phys. **08** (2016) 106.

[62] Pavan Hosur, Xiao-Liang Qi, Daniel A. Roberts, and Beni Yoshida, *Chaos in quantum channels*, J. High Energy Phys. **02** (2016) 004.

[63] Daniel A. Roberts and Beni Yoshida, *Chaos and complexity by design*, J. High Energy Phys. **04** (2017) 121.

[64] Laura Foini and Jorge Kurchan, *Eigenstate thermalization hypothesis and out of time order correlators*, Phys. Rev. E **99**, 042139 (2019).

[65] Laura Foini and Jorge Kurchan, *Eigenstate thermalization and rotational invariance in ergodic quantum systems*, Phys. Rev. Lett. **123**, 260601 (2019).

[66] Amos Chan, Andrea De Luca, and J. T. Chalker, *Eigenstate correlations, thermalization, and the butterfly effect*, Phys. Rev. Lett. **122**, 220601 (2019).

[67] Chaitanya Murthy and Mark Srednicki, *Bounds on chaos from the eigenstate thermalization hypothesis*, Phys. Rev. Lett. **123**, 230606 (2019).

[68] Jonas Richter, Anatoly Dymarsky, Robin Steinigeweg, and Jochen Gemmer, *Eigenstate thermalization hypothesis beyond standard indicators: Emergence of random-matrix behavior at small frequencies*, Phys. Rev. E **102**, 042127 (2020).

[69] Jiaozhi Wang, Mats H. Lamann, Jonas Richter, Robin Steinigeweg, Anatoly Dymarsky, and Jochen Gemmer, *Eigenstate thermalization hypothesis and its deviations from random-matrix theory beyond the thermalization time*, Phys. Rev. Lett. **128**, 180601 (2022).

[70] Marlon Brenes, Silvia Pappalardi, Mark T. Mitchison, John Goold, and Alessandro Silva, *Out-of-time-order correlations and the fine structure of eigenstate thermalization*, Phys. Rev. E **104**, 034120 (2021).

[71] Anatoly Dymarsky, *Bound on eigenstate thermalization from transport*, Phys. Rev. Lett. **128**, 190601 (2022).

[72] Zohar Nussinov and Saurish Chakrabarty, *Exact universal chaos, speed limit, acceleration, Planckian transport coefficient, “collapse” to equilibrium, and other bounds in thermal quantum systems*, Ann. Phys. (Amsterdam) **443**, 168970 (2022).

[73] Julian Sonner and Manuel Vielma, *Eigenstate thermalization in the Sachdev-Ye-Kitaev model*, J. High Energy Phys. **11** (2017) 149.

[74] Daniel Louis Jafferis, David K. Kolchmeyer, Baur Mukhametzhanov, and Julian Sonner, *Matrix models for eigenstate thermalization*, Phys. Rev. X **13**, 031033 (2023).

[75] Daniel Louis Jafferis, David K. Kolchmeyer, Baur Mukhametzhanov, and Julian Sonner, *JT gravity with matter, generalized ETH, and random matrices*, arXiv: 2209.02131.

[76] Silvia Pappalardi, Laura Foini, and Jorge Kurchan, *Eigenstate thermalization hypothesis and free probability*, Phys. Rev. Lett. **129**, 170603 (2022).

[77] Dan Voiculescu, *Limit laws for random matrices and free products*, Inventiones Mathematicae **104**, 201 (1991).

[78] Roland Speicher, *Free probability theory and non-crossing partitions*, Sém. Lotharingien Combinat. **39**, B39c (1997), https://www.mat.univie.ac.at/~slc/s/s39speicher.pdf.

[79] Roland Speicher, *Free probability theory*, arXiv:0911.0087.

[80] Silvia Pappalardi and Jorge Kurchan, *Quantum bounds on the generalized Lyapunov exponents*, Entropy **25**, 246 (2023).

[81] Edouard Brézin, Claude Itzykson, Giorgio Parisi, and Jean-Bernard Zuber, *Planar diagrams*, Commun. Math. Phys. **59**, 35 (1978).

[82] Predrag Cvitanović, *Planar perturbation expansion*, Phys. Lett. B **99**, 49 (1981).

[83] Predrag Cvitanović, P. G. Lauwers, and P. N. Scharbach, *The planar sector of field theories*, Nucl. Phys. **B203**, 385 (1982).

[84] Kurusch Ebrahimi-Fard and Frédéric Patras, *The combinatorics of Green's functions in planar field theories*, Front. Phys. **11**, 1 (2016).

[85] Benoît Collins and Ion Nechita, *Random matrix techniques in quantum information theory*, J. Math. Phys. (N.Y.) **57**, 015215 (2016).

[86] Matteo Bellitti, Siddhardh Morampudi, and Chris R. Laumann, *Hamiltonian dynamics of a sum of interacting random matrices*, Phys. Rev. B **100**, 184201 (2019).

[87] Hassan Shapourian, Shang Liu, Jonah Kudler-Flam, and Ashvin Vishwanath, *Entanglement negativity spectrum of random mixed states: A diagrammatic approach*, PRX Quantum **2**, 030347 (2021).

[88] Jinzhao Wang, *Beyond islands: A free probabilistic approach*, J. High Energy Phys. **10** (2023) 040.

[89] Michel Bauer, Denis Bernard, and Tony Jin, *Stochastic dissipative quantum spin chains (I): Quantum fluctuating discrete hydrodynamics*, SciPost Phys. **3**, 033 (2017).

[90] Ludwig Hruza and Denis Bernard, *Coherent fluctuations in noisy mesoscopic systems, the open quantum SSEP, and free probability*, Phys. Rev. X **13**, 011045 (2023).
[91] Michel Bauer, Denis Bernard, Philippe Biane, and Ludwig Hruza, *Bernoulli variables, classical exclusion processes and free probability*, Ann. Henri Poincare **25**, 125 (2024).

[92] Denis Bernard and Ludwig Hruza, *Exact entanglement in the driven quantum symmetric simple exclusion process*, SciPost Phys. **15**, 175 (2023).

[93] Giorgio Cipolloni, László Erdős, and Dominik Schröder, *Thermalisation for Wigner matrices*, J. Funct. Anal. **282**, 109394 (2022).

[94] Kazuya Kaneko, Eiki Iyoda, and Takahiro Sagawa, *Characterizing complexity of many-body quantum dynamics by higher-order eigenstate thermalization*, Phys. Rev. A **101**, 042126 (2020).

[95] Mark Srednicki, *Chaos and quantum thermalization*, Phys. Rev. E **50**, 888 (1994).

[96] Dan V. Voiculescu, Ken J. Dykema, and Alexandru Nica, *Free Random Variables* (American Mathematical Society, Providence, 1992), Vol. 1.

[97] Alexandru Nica and Roland Speicher, *Lectures on the Combinatorics of Free Probability* (Cambridge University Press, Cambridge, England, 2006), Vol. 13.

[98] Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah, *Quantum entanglement growth under random unitary dynamics*, Phys. Rev. X **7**, 031016 (2017).

[99] Adam Nahum, Sagar Vijay, and Jeongwan Haah, *Operator spreading in random unitary circuits*, Phys. Rev. X **8**, 021014 (2018).

[100] C. W. von Keyserlingk, Tibor Rakovszky, Frank Pollmann, and S. L. Sondhi, *Operator hydrodynamics, OTOCs, and entanglement growth in systems without conservation laws*, Phys. Rev. X **8**, 021013 (2018).

[101] Tianci Zhou and Adam Nahum, *Entanglement membrane in chaotic many-body systems*, Phys. Rev. X **10**, 031066 (2020).

[102] Matthew P. A. Fisher, Vedika Khemani, Adam Nahum, and Sagar Vijay, *Random quantum circuits*, Annu. Rev. Condens. Matter Phys. **14**, 335 (2023).

[103] Michael J. Gullans and David A. Huse, *Localization as an entanglement phase transition in boundary-driven Anderson models*, Phys. Rev. Lett. **123**, 110601 (2019).

[104] Michael J. Gullans and David A. Huse, *Entanglement structure of current-driven diffusive fermion systems*, Phys. Rev. X **9**, 021007 (2019).

[105] Denis Bernard and Tony Jin, *Open quantum symmetric simple exclusion process*, Phys. Rev. Lett. **123**, 080601 (2019).

[106] Denis Bernard and Tony Jin, *Solution to the quantum symmetric simple exclusion process: The continuous case*, Commun. Math. Phys. **384**, 1141 (2021).

[107] Denis Bernard and Lorenzo Piroli, *Entanglement distribution in the quantum symmetric simple exclusion process*, Phys. Rev. E **104**, 014146 (2021).

[108] Denis Bernard, Fabian Essler, Ludwig Hruza, and Marko Medenjak, *Dynamics of fluctuations in quantum simple exclusion processes*, SciPost Phys. **12**, 042 (2022).

[109] Tobias Swann, Denis Bernard, and Adam Nahum, *Space-time picture for entanglement generation in noisy fermion chains*, arXiv:2302.12212.

[110] Shao-Kai Jian and Brian Swingle, *Note on entropy dynamics in the Brownian SYK model*, J. High Energy Phys. 03 (2021) 042.

[111] Christoph Sünderhauf, Lorenzo Piroli, Xiao-Liang Qi, Norbert Schuch, and J. Ignacio Cirac, *Quantum chaos in the Brownian SYK model with large finite n: OTOCs and tripartite information*, J. High Energy Phys. 11 (2019) 038.

[112] Lakshya Agarwal and Shenglong Xu, *Emergent symmetry in Brownian SYK models and charge dependent scrambling*, J. High Energy Phys. 02 (2022) 1–57.

[113] Douglas Stanford, Zhenbin Yang, and Shunyu Yao, *Sub-leading weingartens*, J. High Energy Phys. 02 (2022) 045.

[114] Silvia Pappalardi, Laura Foini, and Jorge Kurchan, *Quantum bounds and fluctuation-dissipation relations*, SciPost Phys. **12**, 130 (2022).

[115] We assume the existence of the inverse $Q^{-1}$, valid for $k \le N$, which is the case under consideration.

[116] Benoît Collins, *Moments and cumulants of polynomial random variables on unitary groups, the Itzykson-Zuber integral, and free probability*, Int. Math. Res. Not. **2003**, 953 (2003).

[117] Yinzheng Gu, *Moments of random matrices and Weingarten functions*, Ph.D. thesis, Queen's University, 2013.

[118] Note that this expression is the same as in Ref. [63] upon substituting $\alpha^{-1}$ with $\alpha$. Here, we use $P_{\alpha^{-1}}$ for consistency with the standard notation for the Weingarten matrix [116].

[119] We are furthermore implicitly assuming that, while we take the limit $D \to \infty$, $\langle A^n \rangle$ remains finite. This indeed is the case for many choices of operators in a many-body physics setting; i.e., this is true if $\langle A \rangle$ is a few-body or local operator.

[120] James A. Mingo and Roland Speicher, *Free Probability and Random Matrices* (Springer, New York, 2017), Vol. 35.

[121] In the literature of free probability, the expectation value is usually noted by $\phi(\bullet)$ [79].

[122] Mixed means that the random variables appearing in the free cumulant do not all come from the same subalgebra.

[123] Formally, $K^{-1}$ is distinct from the duality transformation $^*$.
In fact, applying the duality twice, one would generate a
cyclic permutation, corresponding to a clockwise rotation
by two units in the notation of Fig. 3.

[124] Felix M. Haehl, R. Loganayagam, Prithvi Narayan, Amin A. Nizami, and Mukund Rangamani, *Thermal out-of-time-order correlators, KMS relations, and spectral functions*, J. High Energy Phys. 12 (2017) 154.

[125] Silvia Pappalardi, Felix Fritzsch, and Tomaž Prosen, *Full eigenstate thermalization via free cumulants in quantum lattice systems*, arXiv:2303.00713.

[126] F. H. L. Essler and A. J. J. M. de Klerk, *Statistics of matrix elements of local operators in integrable models*, Phys. Rev. X **14**, 031048 (2014).

[127] As mentioned above, since powers of ETH operators
are ETH, the condition (60) implies that all mixed free
cumulants of the order of smaller than or equal to 2k vanish
as well.

[128] Zhengyan Darius Shi, Shreya Vardhan, and Hong Liu, *Local dynamics and the structure of chaotic eigenstates*, Phys. Rev. B **108**, 224305 (2023).

[129] Daniel K. Mark, Federica Surace, Andreas Elben, Adam L. Shaw, Joonhee Choi, Gil Refael, Manuel Endres, and
Soonwon Choi, A maximum entropy principle in deep thermalization and in Hilbert-space ergodicity, Phys. Rev. X **14**, 041051 (2024).

[130] Here, $|i_1,...,k\rangle$ ($|\tilde{i}_1,...,k\rangle$) represents the basis of $\mathcal{H}_X$ ($\mathcal{H}_{\tilde{X}}$).

[131] Rather than being concerned with the average Rényi entropies of an ensemble of states, we consider the average of their exponentials, which directly probe the kth moment of the state; cf. Eq. (64).

[132] Don N. Page, Average entropy of a subsystem, Phys. Rev. Lett. **71**, 1291 (1993).

[133] Karol Zyczkowski and Hans-Jürgen Sommers, Induced measures in the space of mixed quantum states, J. Phys. A **34**, 7111 (2001).

[134] MuSeong Kim, Mi-Ra Hwang, Eylee Jung, and DaeKil Park, Average Rényi entropy of a subsystem in random pure state, Quantum Inf. Process. **23**, 37 (2024).

[135] Benoît Collins and Ion Nechita, Random quantum channels I: Graphical calculus and the Bell state phenomenon, Commun. Math. Phys. **297**, 345 (2010).

[136] Benoît Collins and Ion Nechita, Gaussianization and eigenvalue statistics for random quantum channels (III), Ann. Appl. Probab. **21**, 1136 (2011).

[137] Zi-Wen Liu, Seth Lloyd, Elton Zhu, and Huangjun Zhu, Entanglement, quantum randomness, and complexity beyond scrambling, J. High Energy Phys. 07 (2018) 041.

[138] Chaitanya Murthy and Mark Srednicki, Structure of chaotic eigenstates and their entanglement entropy, Phys. Rev. E **100**, 022131 (2019).

[139] Lev Vidmar and Marcos Rigol, Entanglement entropy of eigenstates of quantum chaotic Hamiltonians, Phys. Rev. Lett. **119**, 220603 (2017).

[140] Harshank Shrotriya and Wen Wei Ho, Nonlocality of deep thermalization, arXiv:2305.08437.

[141] Elliott H. Lieb and Derek W. Robinson, The finite group velocity of quantum spin systems, Commun. Math. Phys. **28**, 251 (1972).

[142] Joel L. Lebowitz and Oliver Penrose, Modern ergodic theory, Phys. Today **26**, No. 2, 23 (1973).

[143] Jun Li, Ruihua Fan, Hengyan Wang, Bingtian Ye, Bei Zeng, Hui Zhai, Xinhua Peng, and Jiangfeng Du,

Measuring out-of-time-order correlators on a nuclear magnetic resonance quantum simulator, Phys. Rev. X **7**, 031011 (2017).

[144] Martin Gärttner, Justin G. Bohnet, Arghavan Safavi-Naini, Michael L. Wall, John J. Bollinger, and Ana Maria Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nat. Phys. **13**, 781 (2017).

[145] K. A. Landsman, C. Figgatt, T. Schuster, N. M. Linke, B. Yoshida, N. Y. Yao, and C. Monroe, Verified quantum information scrambling, Nature (London) **567**, 61 (2019).

[146] Manoj K. Joshi, Andreas Elben, Benoît Vermersch, Tiff Brydges, Christine Maier, Peter Zoller, Rainer Blatt, and Christian F. Roos, Quantum information scrambling in a trapped-ion quantum simulator with tunable range interactions, Phys. Rev. Lett. **124**, 240505 (2020).

[147] M. S. Blok, V. V. Ramasesh, T. Schuster, K. O’Brien, J. M. Kreikebaum, D. Dahlen, A. Morvan, B. Yoshida, N. Y. Yao, and I. Siddiqi, Quantum information scrambling on a superconducting qutrit processor, Phys. Rev. X **11**, 021010 (2021).

[148] Benoît Collins, James A. Mingo, Piotr Śniady, and Roland Speicher, Second order freeness and fluctuations of random matrices. III: Higher order freeness and free cumulants, Doc. Math. **12**, 1 (2024).

[149] The Schur-Weyl duality states that an operator $A$ on $\mathcal{H}^{\otimes k}$ commutes with all operators $V^{\otimes k}$ with $V \in \mathcal{U}(\mathcal{H})$ if and only if it is a linear superposition of permutation operators.

[150] We assume the existence of the inverse $Q^{-1}$, valid for $k \le N$, which is the case under consideration.

[151] Alternatively, one might have tried to apply the result of Sec. IV by starting with the state in Eq. (65) and applying the formula for the Haar channel in Eqs. (42)–(44) with $A_1 = \cdots = A_k = |0\rangle\langle 0|$. Note, however, that the approach in Sec. IV was developed for operators $A_j$ such that $\text{Tr}(A_j) \sim O(D)$; hence, those results cannot be applied directly in this case, and it is necessary to proceed as described in the main text.