#  The statistical properties of eigenstates in chaotic many-body quantum
systems

Dominik Hahn<EMAIL_ADDRESS>Max Planck Institute for the Physics of Complex
Systems, Nöthnitzer Straße 38, 01187-Dresden, Germany  David J. Luitz
<EMAIL_ADDRESS>Physikalisches Institut, Universität Bonn, Nussallee
12, 53115 Bonn, Germany  J. T. Chalker<EMAIL_ADDRESS>Theoretical Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United
Kingdom

###### Abstract

We consider the statistical properties of eigenstates of the time-evolution
operator in chaotic many-body quantum systems. Our focus is on correlations
between eigenstates that are specific to spatially extended systems and that
characterise entanglement dynamics and operator spreading. In order to isolate
these aspects of dynamics from those arising as a result of local conservation
laws, we consider Floquet systems in which there are no conserved densities.
The correlations associated with scrambling of quantum information lie outside
the standard framework established by the eigenstate thermalisation hypothesis
(ETH). In particular, ETH provides a statistical description of matrix
elements of local operators between pairs of eigenstates, whereas the aspects
of dynamics we are concerned with arise from correlations amongst sets of four
or more eigenstates. We establish the simplest correlation function that
captures these correlations and discuss features of its behaviour that are
expected to be universal at long distances and low energies. We also propose a
maximum-entropy Ansatz for the joint distribution of a small number $n$ of
eigenstates. In the case $n=2$ this Ansatz reproduces ETH. For $n=4$ it
captures both the growth with time of entanglement between subsystems, as
characterised by the purity of the time-evolution operator, and also operator
spreading, as characterised by the behaviour of the out-of-time-order
correlator. We test these ideas by comparing results from Monte Carlo sampling
of our Ansatz with exact diagonalisation studies of Floquet quantum circuits.

##  I Introduction

Although textbook approaches to the thermodynamic equilibrium of quantum
systems rely on invoking a weak coupling to a heat bath, it was understood
over the course of the last four decades that this is not strictly necessary.
In a large class of many-body quantum systems, the interactions between its
constituents enable an isolated system to act as its own heat bath and to
reach a thermal equilibrium state at long times when starting from most
nonequilibrium initial states. While the late-time state remains a pure state,
it approaches a _typical state_ , representative of the Gibbs ensemble with
small fluctuations away from this state for large system size [1]. This
phenomenon emerges from the pseudo-random nature of physical observables in
the energy-eigenbasis, which was suggested as a criterion for quantum chaos
[1] and demonstrated numerically early on [2], with diagonal matrix elements
clustering around equilibrium expectation values [3]. Integrable systems can
evade this behaviour, but are not robust in the sense that very weak
perturbations suffice to recover thermalisation [4].

These observations were subsequently formalised and now constitute the
eigenstate thermalisation hypothesis (ETH) [1, 2, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13], which can be derived from the assumption that the behaviour of the
quantum system within a narrow energy window is essentially captured by random
matrix theory. This generically leads to Gaussian distributions of the matrix
elements of local operators in the energy eigenbasis [14, 8, 9, 10], with
possible exceptions in tails of the distribution for systems with slow
particle transport [15, 16, 17, 18]. In these Gaussian distributions, the mean
values of diagonal matrix elements are fixed to reproduce statistical
mechanical averages of observables, while off-diagonal matrix elements have
zero mean and an energy structure in their variance that governs the dynamics
of autocorrelation functions of the local operators [19, 12].

It was pointed out recently on general grounds that this picture cannot be
complete [20, 21]. If matrix elements of local operators in the eigenbasis
were independent Gaussian random variables, then their mean and variance would
determine not only autocorrelation functions, but also all higher-order
correlators. In particular, the implications for the out-of-time-order
correlators (OTOC) [22, 23] are in stark contrast to the known behaviour of
ergodic many-body quantum systems. This leads to the conclusion that there are
necessarily correlations between matrix elements, which contain information
characterised by higher order cumulants [20, 21, 24, 25].

More specifically, in quantum systems with local interactions, there are
strong bounds on how fast correlations can spread [26], and this limits for
instance the rate of growth of the entanglement entropy [27, 28]. This
behaviour is reflected in the typically linear growth of the operator
entanglement entropy of the time evolution operator of local systems [29, 30],
and is captured by light-cone structures of out-of-time-order correlators
[31]. Our focus here is on the resulting correlations between matrix elements
of observables, and related correlations between eigenstates of the time-
evolution operator.

A separate potential source of correlations between matrix elements is
provided by locally conserved densities. Such correlations were considered
early on [32], have been investigated via transport timescales [33] and
subsequently observed in larger systems [34, 35]. To isolate features arising
because of the dynamics of quantum information from those due to locally
conserved densities, we focus in the following on Floquet systems, for which
time-dependence of the Hamiltonian eliminates energy conservation, and in
which there are no other local conservation laws.

Recent work has presented a generalised formulation of ETH using Free
Probability theory and numerical tests for higher-order correlations between
matrix elements [36, 37]. That perspective considers matrix elements of local
operators as fundamental objects, revealing the frequency structure of the
higher-order free cumulants, particularly the fourth-order free cumulant,
which encodes the leading correlations of matrix elements beyond standard ETH.

An alternative perspective, adopted in Ref. [38], is to consider eigenstates
of the time-evolution operator, rather than matrix elements of observables, as
the relevant set of variables. In particular, the typical time evolution of
Renyi entropy in local systems implies non-trivial correlations between
eigenstates. Separately, in Ref. [39] a derivation is given of ETH via a study
of eigenstates in random Floquet quantum circuits using a field-theoretic
approach.

In the present paper, we discuss correlations in chaotic many-body quantum
systems that are associated with the dynamics of quantum information
scrambling and that can be expected to be universal in spatially extended
systems with local interactions. While ETH is formulated as a statement about
the statistical properties of matrix elements of operators, we find that it is
more transparent to consider instead eigenstates and correlators constructed
from them, without reference to particular operators. The first main
contribution of our work is to identify the leading-order eigenstate
correlator that contains information about this scrambling dynamics, and to
discuss its behaviour at long distances and low energy differences. The second
main contribution is to introduce an Ansatz for the joint probability
distribution of a small number $n$ of eigenstates of the time-evolution
operator. In the case of a pair of eigenstates, this reproduces the Gaussian
distribution for matrix elements of local observables that constitutes ETH.
Extending this Ansatz, we show that a constraint on the joint distribution of
four eigenstates is sufficient to capture the essential features of the OTOC
and of the operator entanglement entropy of the time-evolution operator. As a
third main contribution, we test these ideas by comparing results from exact
diagonalisation of the time-evolution operator with those from Monte Carlo
sampling of our Ansatz for the joint eigenstate distribution.

The ideas and results that we set out in this paper are consistent with, but
broadly complementary to, the recent discussion of a generalised version of
ETH, formulated to describe the average of products of arbitrary numbers of
matrix elements of observables using the language of Free Probability [36,
37]. In particular, our emphasis is different from that of [37] in two ways.
First, we focus on correlations at large distances and small energy
separations that we expect to be a universal consequence of the dynamics of
quantum information. Second, we find that it is advantageous to consider
correlations between eigenstates in place of matrix elements.

The remainder of this paper is organised as follows. In Sec. II we provide a
compact overview of our main results. In Sec. III we give details of our
calculations, including the microscopic models we use for numerical studies,
the determination of parameters in our Ansatz for the joint probability
distribution of eigenstates, and a summary of the numerical methods used. We
develop a treatment of our Ansatz based on a perturbative expansion in Sec. IV
and we present numerical results for additional models in Sec. V. We conclude
with a summary and outlook in Sec. VI.

##  II Overview

In this section we provide an overview of our results. We introduce the class
of model studied and the eigenstate correlators of interest in Sec. II.1. We
set out the relationship between these correlators and the OTOC in Sec. II.2,
and indicate generalisations in Sec. II.3. We review the sense in which the
original form of ETH fails to capture these correlators in Sec. II.4 and we
propose a representation for them in terms of the joint distribution function
for sets of four eigenstates in Sec. II.5. We present results for the
correlators from exact diagonalisation and from Monte Carlo sampling of this
distribution in Sec. II.6.

###  II.1 Models and Correlation Functions

We start by setting out some essential notation. We consider a one-dimensional
Floquet system consisting of $L$ sites, each with a local Hilbert space
dimension $q$, coupled by local interactions. We use $W$ to denote the Floquet
operator for the system, which is then a unitary $q^{L}\times q^{L}$ matrix
that generates evolution through one time period. Defining eigenstates
$\ket{a}$ and quasienergies $\theta_{a}$ satisfying
${W}|a\rangle=\mathrm{e}^{-\mathrm{i}\theta_{a}}|a\rangle$, the evolution
operator $W(t)$ for an integer number of time-steps $t$ has the spectral
decomposition

$W(t)=\sum_{a}\mathrm{e}^{-\mathrm{i}\theta_{a}t}|a\rangle\langle a|\,.$ (1)

The dynamics can be characterised in terms of correlators of local operators
$X_{\alpha},Y_{\beta}\ldots$. Here we use upper case letters $X,Y,\ldots$ as
labels for subsystems on which local operators act, with subscripts
$\alpha,\beta\ldots$ to distinguish different operators acting on a given
subsystem. Later we will use $\overline{X}$, $\overline{Y},\ldots$ to denote
the complements of these subsystems. The time evolution in the Heisenberg
picture is $X_{\alpha}(t)=W^{\dagger}(t)X_{\alpha}W(t)$, and for a Floquet
system it is natural to evaluate correlators using the infinite-temperature
density matrix. In addition, it is useful to consider an ensemble of
realisations of $W$ and to average physical quantities over the ensemble. We
indicate this average by $[\ldots]_{\rm av}$. An alternative average is over a
Haar distribution of vectors, which we indicate by $[\ldots]_{0}$.

The simplest correlator is the autocorrelation function of a single operator,
which has the spectral decomposition

$q^{-L}{\rm Tr}[X_{\alpha}(t)X_{\alpha}]=q^{-L}\sum_{ab}|\langle
a|X_{\alpha}|b\rangle|^{2}\mathrm{e}^{\mathrm{i}(\theta_{a}-\theta_{b})t}\,.$
(2)

In an ergodic system this is expected to decay on a timescale that is
microscopic in the sense that it is of order a few Floquet periods. Evidently,
its behaviour reflects statistical properties of pairs $|a\rangle,|b\rangle$
of eigenstates, which are therefore expected to show features as a function of
the quasienergy difference $\theta_{a}-\theta_{b}$ that vary on a scale only a
few times smaller than the spectral width $2\pi$ [21].

The OTOC has the definition and spectral decomposition

$\displaystyle q^{-L}{\rm
Tr}[X_{\alpha}(t)Y_{\beta}X_{\alpha}(t)Y_{\beta}]=q^{-L}$
$\displaystyle\sum_{abcd}\langle a|X_{\alpha}|b\rangle\langle
b|Y_{\beta}|c\rangle\times$ $\displaystyle\times\langle
c|X_{\alpha}|d\rangle\langle d|Y_{\beta}|a\rangle$ $\displaystyle
e^{i(\theta_{a}-\theta_{b}+\theta_{c}-\theta_{d})t}\,.$ (3)

If the subsystems $X$ and $Y$ are separated by a large distance $s$, the main
features of the OTOC appear on a large timescale. More specifically, the
support of the operator $X_{\alpha}(t)$ is expected [26] to grow with a
butterfly velocity $v_{\rm B}$; the OTOC is constant and non-zero if this
support is disjoint from that of $Y_{\beta}$, but falls towards zero when the
support of $X_{\alpha}(t)$ expands to contain that of $Y_{\beta}$. Clearly,
behaviour of the OTOC reflects statistical properties of sets of four
eigenstates, $|a\rangle,|b\rangle,|c\rangle$ and $|d\rangle$, which must show
features on a quasienergy scale $2\pi v_{\rm B}/s$ that is much smaller for
large $s$ than the spectral width.

Besides the OTOC, a second way to characterise quantum information dynamics is
via the spread of entanglement. Consider an initial state $|\psi\rangle$ with
low entanglement in the site basis. Although the corresponding density matrix
$\rho(t)=W(t)|\psi\rangle\langle\psi|W^{\dagger}(t)$ remains pure at all $t$,
for a subsystem $X$ that is much smaller than its complement $\overline{X}$,
the reduced density matrix $\rho_{X}(t)={\rm Tr}_{\overline{X}}\rho(t)$ is
expected to evolve towards an equilibrium one. This is probed at the simplest
level by considering the purity ${\rm Tr}_{X}[\rho_{X}(t)]^{2}$. Since the
definition of the purity involves two powers of $W(t)$ and two of
$W^{\dagger}(t)$, its behaviour, like that of the OTOC, reflects correlations
amongst sets of four eigenstates.

Both the OTOC and the purity require choices in their definitions – of the
operators denoted by $X_{\alpha}$ and $Y_{\beta}$ in Eq. (II.1) for the
former, and of the initial wavefunction $|\psi\rangle$ for the latter. This
arbitrariness can be eliminated by averaging the OTOC over two complete sets
of operators $\\{X_{\alpha}\\}$ and $\\{Y_{\beta}\\}$ with given supports $X$
and $Y$, and by considering the operator entanglement entropy of $W(t)$ [40,
29, 30] in place of the purity of $\rho_{X}(t)$. Both routes lead to an
identical correlator which is defined solely in terms of sets of four
eigenstates and the choice of $X$ and $Y$. We defer discussion of details and
present first an alternative argument that singles out the same correlator.

As a starting point, consider the Schmidt decomposition of an eigenstate in
terms of tensor products of orthonormal basis states $|i_{X}\rangle$ and
$|i_{\overline{X}}\rangle$ for subsystem $X$ and its complement
$\overline{X}$, which we write as

$|a\rangle=\sum_{i_{X}i_{\overline{X}}}[C_{X}(a)]_{i_{X}i_{\overline{X}}}|i_{X}\rangle\otimes|i_{\overline{X}}\rangle\,.$
(4)

Here $C_{X}(a)$ is the matrix version of the eigenstate $\ket{a}$, separating
states on subsystem $X$ into row indices and states on $\overline{X}$ into
column indices. It hence has dimensions $q^{L(X)}\times q^{L(\overline{X})}$,
where $L(X)=|X|$, the number of sites in $X$, and similarly for
$L(\overline{X})$.

The problem of constructing correlators from sets of eigenstates is equivalent
to one of forming scalars from sets of matrices $C_{X}(a)$. This can be done
by taking the trace of products of an even number of terms, in which the
matrices alternate with their Hermitian conjugates. The matrices within such a
trace must all refer to a given choice of subsystem $X$ but may refer to
multiple eigenstates $|a\rangle$, $|b\rangle\ldots$. At lowest order this
recipe simply yields the quantity

${\rm Tr}[C_{X}(a)C^{\dagger}_{X}(b)]=\delta_{ab},$ (5)

which has a value fixed by orthonomality of the eigenstates. At next order it
gives

$M_{X}(abcd)={\rm
Tr}[C_{X}(a)C^{\dagger}_{X}(b)C_{X}(c)C_{X}^{\dagger}(d)]\,.$ (6)

Such quantities can be represented diagrammatically as shown in Fig. 1

Figure 1: Ingredients for a diagrammatic representation: from left to right
$C_{X}(a)$, $C_{X}(b)^{\dagger}$, $X_{\alpha}$, $\langle
a|X_{\alpha}|b\rangle$, and $M_{X}(abcd)$.

We now invoke two guiding ideas. One follows from the fact that we want to
characterise dynamics in space and time, which suggests that we should
consider more than one way of subdividing the system into subsystems, with
different alternatives labelled $X$, $Y$, $\ldots$. The other follows from the
fact that the overall phases of individual eigenstates can be chosen
arbitrarily, but physical quantities should be invariant under the
transformation $|a\rangle\to e^{i\phi_{a}}|a\rangle$. To eliminate the phases
$\phi_{a}$, each $C_{X}(a)$ appearing in a correlator must be accompanied by a
Hermitian conjugate $C_{Y}^{\dagger}(a)$, referring to the same eigenstate but
possibly with a different division into subsystems. Employing both these
ideas, we are led to the main quantity of interest in this paper,
$M^{\phantom{*}}_{X}(abcd)M^{*}_{Y}(abcd)$ and its appropriately normalised
ensemble average

$\displaystyle F_{4}(X,Y,\theta)=$ $\displaystyle
q^{-L(X,Y)}\Big{[}\sum_{abcd}M^{\phantom{*}}_{X}(abcd)M_{Y}^{*}(abcd)\times$
(7)
$\displaystyle\times\delta(\theta-\theta_{a}+\theta_{b}-\theta_{c}+\theta_{d})\Big{]}_{\rm
av}\,.$

Here and in the following, the argument of the $\delta$-function on
quasienergy differences is taken modulo $2\pi$. The length $L(X,Y)$ appearing
in the normalisation is defined by

$\displaystyle L(X,Y)$ $\displaystyle=$ $\displaystyle
2L-|\overline{X}\setminus\overline{Y}|-|X\setminus Y|$ (8) $\displaystyle=$
$\displaystyle 2L-|\overline{Y}\setminus\overline{X}|-|Y\setminus
X|=L(Y,X)\,,$

where notation of the form $|A\setminus B|$ indicates the number of sites that
are in subsystem $A$ but not in $B$. The choice of subscript on
$F_{4}(X,Y,\theta)$ indicates that this quantity characterises correlations
within sets of four eigenstates.

Somewhat surprisingly, at this order in powers of the eigenstates, Eq. (7) is
the unique outcome of interest from the approach we have sketched for
constructing correlators. To see this, consider potential alternatives. Any
such alternative should involve two factors, each consisting of a trace over
products of Schmidt matrices $C_{X}(a)$, since each trace carries a subsystem
label and we are interested in correlations between a pair of subsystems.
Moreover, at this order the two factors together involve four such matrices
and four Hermitian conjugates. If these matrices are equally distributed
between the two traces, then alternatives to Eq. (7) must all be generated by
replacing $M_{Y}^{*}(abcd)$ with a similar factor that preserves invariance
under changes of eigenstate phases. These can all be obtained from Eq. (7)
using the equalities
$M_{Y}(abcd)=M_{Y}(cdab)=M_{\overline{Y}}(adcb)=M^{*}_{Y}(badc)$ and so are
equal to $F_{4}(X,Y,\theta)$ or $F_{4}(X,\overline{Y},\theta)$. Finally, one
might regroup matrices under the trace, so that one trace involves six
matrices and the other trace involves only two. Then, however, the value of
the second trace is fixed via Eq. (5) and is independent of subsystem label,
eliminating the spatial dependence of interest.

Amongst correlators involving only a single $X$, the lowest order quantity
that is independent of eigenstate phases is $M_{X}(abba)$ (equal to
$M_{\overline{X}}(aabb)$). From this we define

$\displaystyle F_{2}(X,\theta)$ $\displaystyle=$ $\displaystyle
q^{-(L+L(X))}\times$ (9) $\displaystyle\times$
$\displaystyle\Big{[}\sum_{ab}M_{X}(abba)\delta(\theta-\theta_{a}+\theta_{b})\Big{]}_{\rm
av}\,.$

The subscript on $F_{2}(X,\theta)$ indicates that this correlator
characterises correlations between pairs of eigenstates.

The two correlators $F_{2}(X,\theta)$ and $F_{4}(X,Y,\theta)$ are the central
quantities of interest in the following, together with their counterparts in
the time domain, defined by

$f_{2}(X,t)=\int_{-\pi}^{\pi}{\rm d}\theta\,F_{2}(X,\theta)e^{i\theta t}$ (10)

and

$f_{4}(X,Y,t)=\int_{-\pi}^{\pi}{\rm d}\theta\,F_{4}(X,Y,\theta)e^{i\theta
t}\,.$ (11)

The initial values $f_{2}(X,0)=f_{4}(X,Y,0)=1$ follow from completeness of the
set of eigenstates and the choices of normalisation in Eqns. (7) and (9) (see
Eqns. (28) and (29) for a discussion). The late-time limits are also system
independent. For $f_{4}(X,Y,t)$ this limit comes (assuming no degeneracies)
from terms in (7) with pairwise equal labels $a{=}b,c{=}d$ or $a{=}d,b{=}c$.
For $f_{2}(X,t)$ it comes from terms in (9) with $a{=}b$. All these terms can
be simplified by noting that $C_{X}(a)C_{X}^{\dagger}(a)\equiv{\rm
Tr}_{\overline{X}}|a\rangle\langle a|$ is the reduced density matrix on
subsystem $X$ formed from the eigenstate $|a\rangle$. If $L(X)\ll L$, one
expects from ETH that to an excellent approximation
$C_{X}(a)C^{\dagger}_{X}(a)=q^{-L(X)}\openone_{X}$, where $\openone_{X}$ is
the identity on $X$. This implies that
$\lim_{|t|\to\infty}f_{2}(X,t)=q^{-2L(X)}$ and that
$\lim_{t\to\infty}f_{4}(X,Y,t)=q^{2L-L(X,Y)-L(X)-L(Y)}$ for $L(X),L(Y)\ll L$.

###  II.2 Relation to autocorrelation functions of observables and OTOC

As we now discuss, these correlators are related respectively to the
autocorrelation function [Eq. (2)] and the OTOC [Eq. (II.1)] via averages over
the operators appearing in the latter two quantities. We begin by stating a
key relation between $M_{X}(abcd)$ and an operator average of matrix elements.
Given a subsystem $X$, choose a complete basis of $q^{2L(X)}$ Hermitian
operators $\\{X_{\alpha}\\}$ that act on the subsystem and obey the
orthonormality condition

$q^{-L(X)}{\rm Tr}_{X}[X_{\alpha}X_{\beta}]=\delta_{\alpha\beta}\,.$ (12)

Using the resolution of the identity in the vector space of operators, one
finds

$q^{-L(X)}\sum_{\alpha}\langle a|X_{\alpha}|b\rangle\langle
c|X_{\alpha}|d\rangle=M_{X}(abcd)\,.$ (13)

The operator resolution of the identity and the relation given in Eq. (13) are
represented diagrammatically in Fig. 2.

Figure 2: Diagrammatic representation of: (top) the operator resolution of the
identity; (bottom) Eq. (13).

We first apply this to the simpler case of the autocorrelation function. The
autocorrelation function averaged over all choices of operator (and over the
ensemble of systems) is

$q^{-2L(X)}\sum_{\alpha}\Big{[}q^{-L}{\rm
Tr}[X_{\alpha}(t)X_{\alpha}]\Big{]}_{\rm av}=f_{2}(X,t)\,.$ (14)

The special case $X_{\alpha}=\openone_{X}$ is the only contribution to this
average that survives at late times, giving a value of
$\lim_{|t|\to\infty}f_{2}(X,t)$ consistent with the discussion above.

Similar arguments apply to the OTOC. In this case we choose two complete sets
of operators. Operators in one set act on the subsystem labelled $X$ and are
denoted by $X_{\alpha}$. Operators in the other set act on the subsystem
labelled $Y$ and are denoted by $Y_{\beta}$. Then

$\displaystyle q^{-(L(X)+L(Y))}$ $\displaystyle\sum_{\alpha\beta}\langle
a|X_{\alpha}|b\rangle\langle b|Y_{\beta}|c\rangle\langle
c|X_{\alpha}|d\rangle\langle d|Y_{\beta}|a\rangle$ (15)
$\displaystyle=M_{X}(abcd)M_{Y}(bcda)\,.$

This can be rewritten using $M_{Y}(bcda)=M^{*}_{\overline{Y}}(abcd)$, and so
the average of the OTOC over both sets of operators is

$\displaystyle q^{-2(L(X)+L(Y))}$
$\displaystyle\sum_{\alpha\beta}\Big{[}q^{-L}{\rm
Tr}[X_{\alpha}(t)Y_{\beta}X_{\alpha}(t)Y_{\beta}]\Big{]}_{\rm av}$ (16)
$\displaystyle=q^{-S(X,Y)}f_{4}(X,\overline{Y},t)$

with
$S(X,Y)=L(X)+L(Y)+L-L(X,\overline{Y})=|X\setminus\overline{Y}|+|Y\setminus\overline{X}|$.
Contributions to this average from the special cases $X_{\alpha}=\openone_{X}$
and/or $Y_{\beta}=\openone_{Y}$ survive at long times and are responsible for
the limiting value given above.

The correlator $f_{4}(X,Y,t)$ also arises from a discussion of the operator
entanglement entropy of the evolution operator. This quantity stems from
considering the operator $W(t)$ as a state on a doubled Hilbert space, with
components given by the matrix elements
$[W(t)]_{i_{X}i_{\overline{X}},j_{Y}j_{\overline{Y}}}$. The corresponding
reduced operator density matrix, obtained by tracing out the degrees of
freedom in the subsystems $\overline{X}$ and $\overline{Y}$, is

$\displaystyle[\rho(X,Y,W(t))$
$\displaystyle]_{i_{X}j_{Y},l_{X}m_{Y}}=\sum_{ab}[C_{X}(a)C_{X}^{\dagger}(b)]_{i_{X}l_{X}}\times$
$\displaystyle\times$
$\displaystyle[C_{Y}(b)C_{Y}^{\dagger}(a)]_{m_{Y}j_{Y}}e^{i(\theta_{b}-\theta_{a})t}\,.$
(17)

The ensemble-averaged operator purity arising from this reduced density matrix
is simply

$\big{[}{\rm Tr}[\rho(X,Y,W(t))^{2}]\big{]}_{\rm
av}=q^{L(X,Y)}f_{4}(X,Y,t)\,.$ (18)

The proportionality between $f_{4}(X,Y,t)$ and the operator purity of the
evolution operator implies a straightforward link to the idea of an
entanglement membrane, which has been proposed as a coarse-grained description
of entanglement dynamics in chaotic many-body quantum systems [41, 42]. For
the one-dimensional models we are considering, the entanglement membrane is a
curve in space-time, and to discuss the link to operator purity we build on
the exposition of Ref. [42]. In outline, coarse-grained features of
entanglement dynamics are determined by the line tension ${\cal E}(v)$ of this
membrane, which is a function of a velocity $v$. In our notation, $v=s/t$,
where $s$ is the distance between the ends of subsystems $X$ and $Y$, defined
in Fig. 3(d). For a fixed choice of $X$ and $Y$ with $s$ large, the operator
purity of the evolution operator is proportional to the line tension, and so
traces out the function ${\cal E}(v)$ as $t$ varies. Hence the correlators
$f_{4}(X,Y,t)$ and $F_{4}(X,Y,\theta)$ can be seen as representations of the
line tension ${\cal E}(v)$.

###  II.3 Multi-time and multi-quasienergy correlators

An obvious generalisation [21] of the OTOC [Eq. (II.1)] is to introduce three
time arguments, by considering the quantity $q^{-L}{\rm
Tr}[X_{\alpha}(t+t_{2})Y_{\beta}(t_{1})X_{\alpha}(t)Y_{\beta}]$.
Correspondingly, in the quasienergy domain we have a generalisation of the
correlator $F_{4}(X,Y,\theta)$, defined by

$\displaystyle F_{4}(X,Y;\theta,\theta_{1},\theta_{2})$ $\displaystyle=$
$\displaystyle
q^{-L(X,Y)}\Big{[}\sum_{abcd}M^{\phantom{*}}_{X}(abcd)M_{Y}^{*}(abcd)\times$
(19) $\displaystyle\times$
$\displaystyle\delta(\theta-\theta_{a}+\theta_{b}-\theta_{c}+\theta_{d})\times$
$\displaystyle\times$
$\displaystyle\delta(\theta_{2}-\theta_{a}+\theta_{b})\times\delta(\theta_{1}-\theta_{b}+\theta_{c})\Big{]}_{\rm
av}$

with the Fourier transform

$\displaystyle f_{4}(X,Y;t,t_{1},t_{2})$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\\!\\!\\!\\!\\!{\rm
d}\theta\int_{-\pi}^{\pi}\\!\\!\\!\\!\\!{\rm
d}\theta_{1}\int_{-\pi}^{\pi}\\!\\!\\!\\!\\!{\rm
d}\theta_{2}\,\,F_{4}(X,Y;\theta,\theta_{1},\theta_{2})\times$ (20)
$\displaystyle\times$ $\displaystyle e^{i(\theta
t+\theta_{1}t_{1}+\theta_{2}t_{2})}$

which is related to the generalised OTOC by

$\displaystyle f_{4}(X,\overline{Y};t,t_{1},t_{2})$ $\displaystyle=$
$\displaystyle q^{S(X,Y)-2L(X)-2L(Y)}\times$
$\displaystyle\times\sum_{\alpha\beta}q^{-L}$ $\displaystyle{\rm Tr}$
$\displaystyle[X_{\alpha}(t+t_{2})Y_{\beta}(t_{1})X_{\alpha}(t)Y_{\beta}]\,.$
(21)

The single-quasienergy correlator can be recovered from the multi-quasienergy
version using

$F_{4}(X,Y,\theta)=\int_{-\pi}^{\pi}{\rm d}\theta_{1}\int_{-\pi}^{\pi}{\rm
d}\theta_{2}\,\,F_{4}(X,Y;\theta,\theta_{1},\theta_{2})\,.$ (22)

For large separation $s$ between the ends of $X$ and $\overline{Y}$ these
correlation functions factorise as follows (see Ref. [21]). Consider behaviour
in the time domain. Let $v_{\rm B}$ be the butterfly velocity and $D$ the
diffusion constant for the spread of an operator front. Then we expect

$\displaystyle q^{-L}$ $\displaystyle{\rm
Tr}[X_{\alpha}(t+t_{2})\overline{Y}_{\beta}(t_{1})X_{\alpha}(t)\overline{Y}_{\beta}]$
$\displaystyle\approx$ $\displaystyle\left\\{\begin{array}[]{lll}q^{-L}{\rm
Tr}[X_{\alpha}(t+t_{2})X_{\alpha}(t)]\times&&\\\ \quad\times q^{-L}{\rm
Tr}[\overline{Y}_{\beta}(t_{1})\overline{Y}_{\beta}]&&s-v_{\rm
B}|t|\gg\sqrt{Dt}\\\ 0&&v_{\rm B}|t|-s\gg\sqrt{Dt}\end{array}\right.$

This motivates the approximation

$f_{4}(X,{Y};t,t_{1},t_{2})\approxeq
f_{4}(X,{Y},t)f_{2}(X,t_{2})f_{2}({Y},t_{1})\,,$ (26)

which in the quasienergy domain is

$F_{4}(X,{Y};\theta,\theta_{1},\theta_{2})\approxeq
F_{4}(X,{Y},\theta)F_{2}(X,\theta_{2})F_{2}({Y},\theta_{1})\,.$ (27)

For the models we study in this paper, $F_{2}(X,\theta_{2})$ is only weakly
dependent on $\theta$ and so the multi-quasienergy correlator carries only
limited extra information compared to the single-quasienergy version. For this
reason we leave study of $F_{4}(X,{Y};\theta,\theta_{1},\theta_{2})$ for
future work.

###  II.4 Existence of correlations beyond ETH

Our objective in the remainder of this work is to find a form for the joint
distribution function (JDF) of a small number of eigenstates that reproduces
these correlations. We do this using a maximum entropy Ansatz with a final
form that we build up by considering first individual vectors, then pairs of
vectors, and finally sets of four vectors.

To place our approach in context, it is useful to recall (following Refs. [20]
and [21]) the limitations of ETH in its standard formulation when applied to
the OTOC. As a starting point, consider the spectral decomposition of the OTOC
in terms of operator matrix elements, as displayed in Eq. (II.1). ETH asserts
that matrix elements of the form $\langle a_{1}|X_{\alpha}|a_{2}\rangle$ and
$\langle a_{3}|Y_{\beta}|a_{4}\rangle$ appearing in this expression are
Gaussian random variables, and are independent apart from the constraint
implied by Hermiticity of the operators $X_{\alpha}$ and $Y_{\beta}$. The mean
values of off-diagonal matrix elements are automatically zero, and those of
diagonal matrix elements are zero for traceless operators in the Floquet
setting of interest. Finally, the variance of these matrix elements is set by
the Hilbert space size and is ${\cal O}(q^{-L})$.

Applying these ideas to Eq. (II.1), the OTOC is given by $q^{-L}$ times a sum
of $q^{4L}$ random ${\cal O}(q^{-2L})$ terms. Of these, only the $q^{L}$ terms
with $a{=}b{=}c{=}d$ are expected from ETH to have a non-zero average. This
would imply an average value for the OTOC of ${\cal O}(q^{-2L})$. Treating the
remaining terms as independent random variables, one expects ${\cal
O}(q^{-L})$ fluctuations around this average. In contrast, the true value is
${\cal O}(1)$ at short times. To resolve this discrepancy it is necessary that
a product of four matrix elements of the form $\langle
a|X_{\alpha}|b\rangle\langle b|Y_{\beta}|c\rangle\langle
c|X_{\alpha}|d\rangle\langle d|Y_{\beta}|a\rangle$ has a non-zero average that
is ${\cal O}(q^{-3L})$ in addition to the ${\cal O}(q^{-2L})$ fluctuations
captured by the standard version of ETH [21]. These additional correlations
are the central concern in this paper and in the generalisation of ETH
discussed in Refs. [36, 37].

A simple demonstration that such correlations must be present, regardless of
details of the dynamics, is provided by a sum rule related to the value of the
OTOC at $t=0$. From the left-hand side of Eq. (II.1), assuming for simplicity
that the subsystems $X$ and $Y$ do not overlap and using the operator
normalisation of Eq. (12), we have

$q^{-L}{\rm
Tr}[X_{\alpha}(t)Y_{\beta}X_{\alpha}(t)Y_{\beta}]\Big{|}_{t=0}=1\,.$ (28)

Using this in Eq. (16) with Eq. (11) we have

$\int_{-\pi}^{\pi}{\rm d}\theta\,F_{4}(X,\overline{Y},\theta)=1\,.$ (29)

This sum rule for $F_{4}(X,\overline{Y},\theta)$ is automatically satisfied if
eigenstates are Haar distributed vectors, and in that case
$F_{4}(X,\overline{Y},\theta)$ is independent of $\theta$. The eigenstate
correlations that we are concerned with generate a dependence of
$F_{4}(X,\overline{Y},\theta)$ on $\theta$ but do not alter the fact that,
with the normalisation of Eq. (7), it has an order of magnitude that is
independent of the Hilbert space dimension $q^{L}$.

###  II.5 Describing correlations beyond ETH

Some constraints on the eigenstate JDF are implied by ETH, which we now
consider. ETH specifies statistical properties of both diagonal and off-
diagonal matrix elements of local observables between eigenstates, and we
discuss the two classes of matrix elements separately.

For a system with a time-independent Hamiltonian, a key part of ETH is that
diagonal matrix elements of observables vary smoothly with energy, taking
average values compatible with a thermal ensemble at the same energy density,
and with fluctuations of a characteristic size that vanishes rapidly as the
thermodynamic limit is approached. By contrast, for Floquet systems,
statistical properties of diagonal matrix elements of observables are
independent of quasienergy. We capture this property of diagonal matrix
elements in a Floquet system by taking individual eigenstates to be
isotropically distributed vectors in the Hilbert space for the model. We
outline in Sec. VI the alternative choice required to model the energy
dependence of diagonal matrix elements in Hamiltonian systems within our
approach. We denote the isotropic (Haar) distribution for one, two or four
orthonormal vectors by $P_{1}^{(0)}(a)$, $P_{2}^{(0)}(a,b)$ and
$P_{4}^{(0)}(a,b,c,d)$ respectively, and set out to modify these distributions
in a way that introduces the correlations of interest.

Statistical properties of off-diagonal matrix elements determine the approach
to equilibrium and the autocorrelation functions of observables. ETH applied
to Floquet systems asserts that these matrix elements are independent Gaussian
random variables with a variance that depends only on quasienergy separation.
A central idea in our work is that strict independence is incompatible with
the correlations implied by the dynamics of quantum information. Instead there
are correlations (albeit weak) and these are better handled by considering
distributions for eigenstates rather than matrix elements. We take the joint
distribution of a pair of eigenstates to have the Maximum Entropy form

$P_{2}(a_{1},a_{2})=Z_{2}^{-1}P_{2}^{(0)}(a_{1},a_{2})e^{-S_{2}(a_{1},a_{2})}$
(30)

with $Z_{2}$ a normalisation constant and

$\displaystyle
S_{2}(a,b)=\sum_{X}G_{2}(X,\theta_{a}-\theta_{b})M_{X}(abba)\,,$ (31)

where the coefficients $G_{2}(X,\theta)$ act as a Lagrange multipliers and
should be chosen to reproduce the behaviour of $F_{2}(X,\theta)$ as determined
for a particular system.

Extending this pattern, we take the joint distribution of four eigenstates to
have the form

$\displaystyle P_{4}(a_{1},a_{2},a_{3},a_{4})$
$\displaystyle=Z_{4}^{-1}P_{4}^{(0)}(a_{1},a_{2},a_{3},a_{4})\times$ (32)
$\displaystyle\times
e^{-\sum_{j<k}S_{2}(a_{j},a_{k})-S_{4}(a_{1},a_{2},a_{3},a_{4})}$

with $Z_{4}$ a normalisation constant and

$\displaystyle S_{4}(a,b,c,d)$ $\displaystyle=$
$\displaystyle\sum_{XY}G_{4}(X,Y,\theta_{a}-\theta_{b}+\theta_{c}-\theta_{d})\times$
(33) $\displaystyle\times M_{X}(abcd)M^{*}_{Y}(abcd)\,.$

Here the Lagrange multipliers $G_{4}(X,Y,\theta)$ should be chosen to
reproduce the behaviour of $F_{4}(X,Y,\theta)$.

Two further ingredients are required. One is to establish a practical method
for deducing the values of the Lagrange multipliers from information on the
correlators. The other is to test the approach by sampling $P_{2}(a,b)$ or
$P_{4}(a,b,c,d)$ and comparing the results with correlators calculated from
exact diagonalisation (ED) of $W$. In this work we focus on the correlator
$F_{4}(X,Y,\theta)$ since it contains the long-distance, low-energy
information related to the dynamics of quantum information. Moreover, in the
models we study, the lower-order correlator $f_{2}(X,t)$ decays rapidly in
time. This implies that $F_{2}(X,\theta)$ is approximately independent of
quasi-energy and so we simply set $G_{2}(X,\theta)$ to zero in our initial
treatment. We return to consideration of a $\theta$-dependent
$F_{2}(X,\theta)$ and non-zero $G_{2}(X,\theta)$ immediately after our
discussion of $F_{4}(X,Y,\theta)$.

The determination of the Lagrange multipliers from ED data for eigenstate
correlator can be seen simply as a fitting problem, but since this involves a
high-dimensional parameter space, alternative approaches are desirable.
Fortunately, as we describe in Sec. III.2, we have been able to find direct
and straightforward methods to derive $G_{4}(X,Y,\theta)$ from
$F_{4}(X,Y,\theta)$ and $G_{2}(X,Y,\theta)$ from $F_{2}(X,Y,\theta)$.

###  II.6 Results

Figure 3: Overview of main results, calculated for the brickwork circuit model
defined in Sec. III.1 with $L=12$, $q=2$ and open boundary conditions:
behaviour from ED in (a) and (b); comparison of ED and MC in (c); and
geometries of subsystems $X$ and $Y$ in (d). (a) The correlator $f_{4}(X,Y,t)$
[Eq. (9)] vs $t$, obtained using ED. Recall [Eq. (16)] that this correlator is
proportional to the OTOC $q^{-L}{\rm
Tr}[X_{\alpha}(t)\overline{Y}_{\beta}X_{\alpha}(t)\overline{Y}_{\beta}]$
averaged over operators $X_{\alpha}$ and $\overline{Y}_{\beta}$ with support
on subsystems $X$ and $\overline{Y}$ respectively. Decay of the correlator
reflects operator spreading, and the onset time for decay increases with $s$.
(b) The correlator $F_{4}(X,Y,\theta)$ [Fourier transform of $f_{4}(X,Y,t)$:
see Eq. (7)] vs quasienergy difference $\theta$, obtained using ED
[contributions to Eq. (7) in which any of the state labels $a$, $b$, $c$ and
$d$ are equal have been omitted: they are atypical and carry vanishing weight
in the thermodynamic limit]. It has a peak centred on $\theta=0$ which grows
narrower and higher with increasing $s$, reflecting the short-time plateau in
(a); black dashed line: behaviour when the Floquet operator $W$ is modelled
using a $q^{L}\times q^{L}$ Haar unitary, showing for this structureless case
that $F_{4}(X,Y,\theta)$ is non-zero but $\theta$-independent. (c) Comparison
of ED results with MC results from the Ansatz for the JDF of Eq. (32) fitted
to behaviour in the geometries of (d), showing excellent agreement between
$F^{\rm MC}_{4}(X,Y,\theta)$ (open circles from MC) and $F_{4}(X,Y,\theta)$
(lines from ED) vs $\theta$ for various $s$. (d) Illustration of two ways of
dividing the 12-site system with open boundary conditions into subsystems by
means of a single spatial cut. In one case the subsystems are labelled $X$ and
$\overline{X}$; in the other the labels are $Y$ and $\overline{Y}$. The
distance between the spatial cuts in the two cases is denoted by $s$.

We implement and test these ideas using an open chain with the Floquet
operator defined by a brickwork circuit (see Sec. III.1). Our fastest method
(see Sec. III.2.1) for determining the Lagrange multipliers
$G_{4}(X,Y,\theta)$ is effective if the probability distribution of
$M_{X}(abcd)$ is well-approximated by a Gaussian, which is the case for the
model studied provided the subsystem sizes $L(X)$ and $L(\overline{X})$ are
not too small. In order to satisfy this requirement, and in order to limit the
total number of Lagrange multipliers under consideration, we include
$G_{4}(X,Y,\theta)$ for all $L-5$ subsystems $X$ with $L(X)>2$ that can be
obtained from the full system by means of a single spatial cut, and similarly
for $Y$. Taking account of symmetry under interchange of $X$ and $Y$, this
gives $(L-5)(L-4)/2$ (i.e. 28 for the example studied of $L=12$) independent
Lagrange multipliers. Using as input the values of $F_{4}(X,Y,\theta)$ for
these subsystem choices from ED, we arrive at a final form for the JDF.

We use Monte Carlo (MC) sampling of this distribution to compute the
correlator [denoted by $F^{\rm MC}_{4}(X,Y,\theta)$] with two objectives.
First, for the simplest choices of $X$ involving single spatial cuts (and
similarly for $Y$), comparison with ED is a test of our Ansatz for the JDF and
of our procedure to determine the Lagrange multipliers. Second, it is
interesting to see whether this input alone is sufficient to capture
correlations more generally. To probe this we compare ED and MC results for
the correlator, making choices of $X$ (and also $Y$) that are defined by more
than one spatial cut. This is a test of the extent to which the proposed JDF
captures long-distance, low energy correlations in general. In particular,
taking $X$ and $Y$ each to consist of a small number of sites acting as the
support for a local observable, we test the implications of the JDF for the
(operator-averaged) OTOC.

Figure 4: Test of JDF fitted to behaviour in the geometries of Fig. 3(d) but
applied to geometries of Fig. 4(b). (a) Comparison of data from MC (open
circles) and ED (solid lines). (b) Partition used for (a), in which the
12-site system is divided by two spatial cuts into a two-site subsystem $X$
and its complement $\overline{X}$, or a two-site subsystem $\overline{Y}$ and
its complement $Y$. Calculations are for the brickwork circuit model defined
in Sec. III.1 with $L=12$, $q=2$ and open boundary conditions.

Some of our principal results are shown in Fig. 3 and discussed in the figure
caption. The main conclusions are as follows. (i) As expected from its
relation to the OTOC, the correlator $f_{4}(X,Y,t)$ is time-independent at
short times and falls off at a timescale that is long if the spatial
separation $s$ between subsystems $X$ and $\overline{Y}$ is large. All aspects
of this behaviour are apparent in Fig. 3(a). (ii) In turn, this implies
structure in $F_{4}(X,Y,\theta)$ at small quasienergies $\theta$, as is
visible in Fig. 3(b); the width in quasienergy of this structure decreases
with increasing $s$. (iii) Monte Carlo sampling of the JDF of Eq. (32), with
Lagrange multipliers determined as described in Sec. III.2, generates results
for $F_{4}(X,Y,\theta)$ that are in excellent agreement with those from ED, as
demonstrated in Fig. 3(c).

Further important results are shown in Fig. 4. Here we examine how well the
JDF constructed using correlators for single-cut subsystems can capture
correlators for two-cut subsystems. It is apparent from Fig. 4(a) that MC
sampling of the JDF generates a moderately good representation of
$F_{4}(X,Y,\theta)$ for the two-site choices of $X$ and $Y$ shown in Fig.
4(b). Equivalently, the JDF determined using information from the geometries
of Fig. 3(d) reproduces the main features of the OTOC, as a function of time
and spatial separation, for two operators, each supported on a pair of sites
in the geometry of Fig. 4(b).

As a final indication of the effectiveness of our approach, we return to the
behaviour of $F_{2}(X,\theta)$, which characterises the correlations that are
incorporated by ETH. For simplicity we consider only the joint distribution of
pair of eigenstates, in this way treating $F_{2}(X,\theta)$ separately from
$F_{4}(X,Y,\theta)$. By determining the Lagrange multiplier $G_{2}(X,\theta)$
from ED data as described in Sec. III.2 we generate and sample from this joint
distribution, with results that are shown in Fig. 5. As is evident, our MC
data are in excellent agreement with ED results for all values of quasienergy
difference $\theta$ and all subsystem choices $X$.

The fact that deviations are small from the value
$F_{2}(X,\theta)=(2\pi)^{-1}\approxeq 0.159$ for a pair of Haar-distributed
orthogonal unit vectors is justification for our omission of $F_{2}(X,\theta)$
in our discussion of $F_{4}(X,Y,\theta)$. A more complete treatment would
require the simultaneous inclusion of both $G_{2}(X,\theta)$ and
$G_{4}(X,Y,\theta)$, following Eq. (32). Some consequences have been discussed
previously in Ref. [21] and we do not consider them further here.

The overall aim of MC sampling from our Maximum Entropy Ansatz for the JDF of
a small number of eigenstates is to test whether, with a suitable choice of
Lagrange multipliers, the JDF reproduces correlations measured from ED. This
test of our approach is a crucial one, and we believe Fig. 3(c) offers
excellent evidence that the JDF can indeed reproduce the required correlations
for the model and parameter range investigated there. Further discussion of
the determination of Lagrange multipliers, including a treatment of other
models, is presented in Sec. V and the Appendix.

$X$$k$

Figure 5:  Numerical results for $F_{2}(X,\theta)$ in a system with $L=12$ and
$q=2$ as a function of the left cut position $k$: ED results (solid lines) vs.
Monte-Carlo (circles). [Delta-function contributions at $\theta=0$ have been
omitted; they are responsible for ensuring that the sum rule
$\int_{-\pi}^{\pi}{\rm d}\theta\,F_{2}(X,\theta)=1$ is satisfied for all $X$.]
The almost perfect agreement between both sets of data indicates that the
correlations between pairs of eigenstates, as captured by ETH, may also be
represented accurately within our approach.

##  III Models, Lagrange multipliers and numerical methods

In Sec. III.1 we describe the microscopic models we use to generate the
numerical results shown in this paper. In Sec. III.2 we set out efficient
methods to determine the Lagrange multipliers $G_{2}(X,\theta)$ and
$G_{4}(X,Y,\theta)$ that appear in the JDF for eigenstates. In Sec. III.3 we
give details of the numerical methods used in this paper.

###  III.1 Models

We use the brick-wall circuit depicted in Fig. 6 as a simple model of a
periodically driven many-body quantum system with local interactions in one
dimension. Such Floquet models have been studied extensively in past work: see
for example [43, 21, 44].

$w_{0}$$w_{2}$$w_{4}$$w_{6}$$w_{8}$$w_{1}$$w_{3}$$w_{5}$$w_{7}$$w_{9}$ Figure
6: Unitary time evolution operator $W$ of the Floquet circuit written as a
tensor network.

This circuit is defined in terms of two-site gates
$w_{i}\in\mathbb{C}^{q^{2}\times q^{2}}$ and the driving period is decomposed
into two parts: in the first half of the period couplings are active only on
even bonds of the system, while in the second half the couplings are active
only on odd bonds. Thus, the time evolution operator for the first half period
is $W_{1}=w_{0}\otimes w_{2}\otimes w_{4}\dots$ and for the second half period
is $W_{2}=\openone\otimes w_{1}\otimes w_{3}\otimes\dots$. The evolution
operator over one full period is $W=W_{2}W_{1}$, and for $t$ periods we write
$W(t)\equiv W^{t}$.

In order to define an ensemble of systems, a natural choice would be to draw
each unitary matrix $w_{i}$ independently from the Haar distribution. We find,
however, that in this case determination of the Lagrange multipliers
$G_{4}(X,Y,\theta)$ is complicated by effects that we attribute to
realisations containing weak links $i$ on which the gate $w_{i}$ is close to
the identity (especially in small systems or with small subsystems). To avoid
such weak links, we draw the $w_{i}$ from a truncated version of the Haar
distribution in which all gates with an operator purity above a cutoff are
discarded. With operator purity defined as in Eq (18) (so that the two-site
identity operator has a purity of $q^{4}$) we take the cutoff to be $0.3\times
q^{4}$ for local Hilbert space dimension $q=2$. The consequences of changing
or omitting this cutoff are discussed in Sec. V.1 and Appendix A.2, and
results for $q=3$ are given in Sec. V.2.

###  III.2 Determining Lagrange multipliers in the JDF

We now discuss the problem of determining the Lagrange multipliers that appear
in Eq. (31) and Eq. (33), and that define the JDFs, Eq. (30) and Eq. (32), of
a small number of eigenstates. We present a method specific to
$G_{4}(X,Y,\theta)$ in Sec. III.2.1 and one specific to $G_{2}(X,\theta)$ in
Sec. III.2.2. These are both single-shot methods that make explicit use of
information about the probability distributions of $M_{X}(abcd)$ and
$M_{X}(abba)$ respectively, and for this reason they are particularly
efficient. In Sec. III.2.3 we outline a third, more generally applicable
iterative approach, that is agnostic to the probability distributions
involved.

####  III.2.1 Determination of $G_{4}(X,Y,\theta)$

As noted above, pairwise correlations between eigenstates are very weak in the
Floquet model we consider, in the sense that the correlator $F_{2}(X,\theta)$
is only weakly dependent on $\theta$, taking values close to that for a Haar-
distributed pair of orthogonal vectors. As a simplifying approximation, we
therefore set to zero the Lagrange multipliers $G_{2}(X,\theta)$ [see Eq.
(31)] that control these pairwise correlations. We make further choices
concerning the set of subsystems $X$ and $Y$ for which Lagrange multipliers
$G_{4}(X,Y,\theta)$ are included in the JDF. Without restrictions there are
$2^{L}-1$ distinct subsystems: to reduce this number we include only Lagrange
multipliers for connected subsystems – those that can be obtained from the
full system (which has open boundary conditions) by means of a single cut –
and we omit them for all subsystems obtained using multiple cuts.

Our objective is to find the values of the Lagrange multipliers
$G_{4}(X,Y,\theta)$ for which the JDF reproduces the eigenstate correlator
$F_{4}(X,Y,\theta)$ (known from ED) as accurately as possible. We are able to
simplify this task and avoid attacking directly a high-dimensional fitting
problem if the quantities $M_{X}(abcd)$ [Eq. (6)] are Gaussian distributed.
Our motivation for treating a model with a truncated Haar distribution of gate
unitaries is that in this system $M_{X}(abcd)$ is well-approximated by a
Gaussian.

Evidence for this is presented in Fig. 7. Here we consider only the magnitude
$|M_{X}(abcd)|$ since the phase of $M_{X}(abcd)$ is dependent on the phase
convention used for the eigenstates, and we compare the distribution of
$|M_{X}(abcd)|$ with one in which the real and imaginary parts of
$M_{X}(abcd)$ are assumed to be uncorrelated Gaussian variables with equal
variances and zero means. An analysis of the dependence of these distributions
on the system size $L$ and subsystem size $L(X)$ (see Appendix A) suggests
that deviations from a Gaussian vanish in the limit of large $L$ and $L(X)$.
Since deviations are significant for small subsystem size, we omit Lagrange
multipliers $G_{4}(X,Y,\theta)$ for subsystems $X$ with $L(X)$ or
$L(\overline{X})\leq 2$. This leaves $(L-5)^{2}$ Lagrange multipliers, or
$(L-5)(L-4)/2$ independent quantities after taking account of symmetry
relations.

Figure 7: Probability distribution of $q^{L}|M^{\phantom{*}}_{X}(abcd)|$ in
the Floquet model of Sec. III.1 (coloured data) compared with fitted Gaussian
distributions (black dashed lines). Data are for $L=12$, $q=2$ and a subsystem
$X$ consisting of the four sites closest to the end of an open system, and are
shown for the four indicated values of the relative phase
$\theta=\theta_{a}-\theta_{b}+\theta_{c}-\theta_{d}$. Results for other $q$
and $X$ are shown in Appendix A.

Similarly, a Gaussian distribution for $M_{X}(abcd)$ also arises from the Haar
distribution for four orthogonal vectors [$P_{4}^{(0)}(a,b,c,d)$ in the
notation of (32)] in the large $q^{L}$, $q^{L(X)}$ limit. In this case, and in
contrast to the Floquet model, the covariance is independent of the
quasienergy difference $\theta$ and has the value

$[M_{X}(abcd)M^{*}_{Y}(abcd)]_{0}=q^{L(X,Y)-4L}\,.$ (34)

In order to represent these Gaussian distributions in a compact way, it is
convenient to introduce notation in which the $L-1$ single-cut subsystems $X$
used to define our Lagrange multipliers are labels for basis states in a
$(L-1)$-dimensional vector space. Then the values of $M_{X}(abcd)$ for
different $X$ are components of an $(L-1)$-component column vector
$\mathsf{M}$, so that $[\mathsf{M}]_{X}=M_{X}(abcd)$, while the Lagrange
multipliers are entries in the $(L-1)\times(L-1)$ matrix $\mathsf{G_{4}}$ with
rows and columns labelled by $X$ and $Y$ respectively, so that
$[\mathsf{G_{4}}]_{X,Y}=G_{4}(X,Y,\theta)$. Similarly, the eigenstate
correlators $F_{4}(X,Y,\theta)$ are elements of a matrix $\mathsf{F_{4}}$.
Then Eq. (33) can be rewritten in the compact form

$S_{4}=\mathsf{M}^{\rm T}\mathsf{G_{4}}\mathsf{M}^{*}\,.$ (35)

The fact that the distributions $P_{4}(a,b,c,d)$ and $P^{(0)}_{4}(a,b,c,d)$
are [at large $L$ and $L(X)$] Gaussian for $\mathsf{M}$ suggests that a
convenient coordinate system consists of the components of $\mathsf{M}$
together with additional variables $\Omega$ that we do not specify explicitly.
We indicate this by writing $P_{4}({\mathsf{M}},\Omega)$ and
$P^{(0)}_{4}({\mathsf{M}},\Omega)$. The result given in Eq. (34) for the
covariance within the Haar distribution implies that

$\displaystyle\int{\rm d}\Omega P^{(0)}_{4}({\mathsf{M}},\Omega)$
$\displaystyle=$ $\displaystyle[{\cal Z}_{4}^{(0)}]^{-1}e^{-\mathsf{M}^{\rm
T}\mathsf{G_{4}^{(0)}}\mathsf{M}^{*}}$ (36)
$\displaystyle\mbox{with}\quad[(\mathsf{G_{4}^{(0)}})^{-1}]_{X,Y}$
$\displaystyle=$
$\displaystyle\big{[}[\mathsf{M}]_{X}[\mathsf{M}^{*}]_{Y}\big{]}_{\rm
0}=q^{L(X,Y)-4L}\,.$

and ${\cal Z}_{4}^{(0)}=\pi^{L-1}/\det{\mathsf{G_{4}^{(0)}}}$. This in turn
implies that

$\int{\rm d}\Omega P_{4}({\mathsf{M}},\Omega)=[{\cal
Z}_{4}]^{-1}e^{-\mathsf{M}^{\rm
T}[\mathsf{G_{4}}+\mathsf{G_{4}^{(0)}]}\mathsf{M}^{*}}$ (37)

and hence that

$\big{[}[\mathsf{M}]_{X}[\mathsf{M}^{*}]_{Y}\big{]}_{\rm
av}=\left[(\mathsf{G_{4}}+\mathsf{G^{0}_{4}})^{-1}\right]_{X,Y}.$ (38)

In addition, we have from Eq. 7

$[\mathsf{F_{4}}]_{X,Y}=(2\pi)^{-1}{q^{4L-L(X,Y)}}\big{[}[\mathsf{M}]_{X}[\mathsf{M}^{*}]_{Y}\big{]}_{\rm
av}\,.$ (39)

Eq. (38) allows the determination of the Lagrange multipliers in Eq. (32) in
terms of the matrix $\mathsf{F_{4}}$, which is obtained using ED. This method
was used to generate Fig. 3(c).

####  III.2.2 Determination of $G_{2}(X,\theta)$

Next we describe the method we use to determine the Lagrange multipliers
$G_{2}(X,\theta)$, treating explicitly the case of $L-1$ subsystems $X$
generated by single cuts. A different approach is required to that for
$G_{4}(X,Y,\theta)$ because the probability distribution of $M_{X}(abba)$ is
quite different to that of $M_{X}(abcd)$ for $a\not=b\not=c\not=d$. Indeed,
while (as discussed) $M_{X}(abcd)$ has a complex Gaussian distribution with
zero mean, $M_{X}(abba)$ is from its definition real and non-negative.

To understand the distribution of $M_{X}(abba)$ it is useful to start from Eq.
(13), which specialises here to

$M_{X}(abba)=q^{-L(X)}\sum_{\alpha}|\langle a|X_{\alpha}|b\rangle|^{2}\,.$
(40)

With $a\not=b$ we expect from ETH that $\langle a|X_{\alpha}|b\rangle$ for
each $\alpha$ is an independent complex Gaussian random variable. From this we
can conclude, first, that the two quantities $M_{X}(abba)$ and
$M_{X^{\prime}}(abba)$ are correlated if the sets $\\{X_{\alpha}\\}$ and
$\\{X^{\prime}_{\alpha}\\}$ have operators in common, and second, that
statistically independent quantities can be constructed using a transformation
that organises the operators into disjoint sets.

In detail, this transformation is defined recursively by considering the $L-1$
subsystems in order of increasing size. Again we introduce an
$(L-1)$-dimensional vector space, and define the vector $\mathsf{M}$ to have
components $[\mathsf{M}]_{X}=M_{X}(abba)$. Similarly, we introduce the
notation $\mathsf{T}$ for our target vector with statistically independent
components $[\mathsf{T}]_{X}$, and $\mathsf{G}$ with
$[\mathsf{G}]_{X}=G_{2}(X,\theta)$ for the vector of Lagrange multipliers. To
write the transformation we abuse notation and substitute in place of the
component label $X$ the value $\ell=L(X)$. Then using Eq. (40) we write

$\displaystyle\begin{split}\mathsf{T}_{1}&=\mathsf{M}_{1}\\\
\mbox{and}\quad\mathsf{T}_{\ell}&=\mathsf{M}_{\ell}-q^{-1}\mathsf{M}_{\ell-1}\end{split}$
(41)

for $\ell=2$ to $L-1$. This can be recast in the matrix form

$\mathsf{T}=\mathsf{V}\mathsf{M}$ (42)

where

$\mathsf{V}=\begin{pmatrix}1&0&0&\dots&0\\\ -q^{-1}&1&0&\dots&0\\\
\vdots&&\ddots&\\\ 0&0&\dots&-q^{-1}&1\end{pmatrix}\,.$ (43)

The effect of this transformation is that Eq. (40) is replaced by
$\mathsf{T}_{\ell}=q^{-\ell}\sum^{\prime}_{\alpha}|\langle
a|X_{\alpha}|b\rangle|^{2}\,,$ where the sum runs over the subset of
$n_{\ell}\equiv(q^{2}-1)q^{2(\ell-1)}$ operators $X_{\alpha}$ that act non-
trivially at the rightmost site in $X$ [see illustration in Fig. 3(d)]. Since
$\langle a|X_{\alpha}|b\rangle$ is Gaussian with mean zero, the variable
$s_{\alpha}\equiv|\langle a|X_{\alpha}|b\rangle|^{2}$ has the distribution
$p_{\alpha}(s_{\alpha})=\sigma_{\alpha}e^{-\sigma_{\alpha}s_{\alpha}}$, where
$\sigma_{\alpha}\equiv\big{[}|\langle a|X_{\alpha}|b\rangle|^{2}]_{\rm
av}\big{]}^{-1}$. We make the approximation that $\sigma_{\alpha}$ takes the
same value $\sigma_{\ell}$ for all $X_{\alpha}$ that contribute to a given
$\mathsf{T}_{\ell}$. This framework applies not only to the true eigenstate
distribution under consideration, but also to vectors with a Haar
distribution, and in the latter case we denote the value of $\sigma_{\alpha}$
by $\sigma^{(0)}_{\ell}$. Then

$[\mathsf{T}_{\ell}]_{\rm
av}=\frac{n_{\ell}}{\sigma_{\ell}}\quad\mbox{and}\quad[\mathsf{T}_{\ell}]_{\rm
0}=\frac{n_{\ell}}{\sigma^{(0)}_{\ell}}\,.$ (44)

With these ingredients in hand, Eq. (31) can be written in the form

$S_{2}=\mathsf{G}^{\rm T}\mathsf{V}^{-1}\mathsf{T}\quad\mbox{so
that}\quad\sigma_{\ell}=\sigma^{(0)}_{\ell}+[\mathsf{G}^{\rm
T}\mathsf{V}^{-1}]_{\ell}\,.$ (45)

Substituting Eq. (44) into Eq. (45) and rearranging, we obtain

$\mathsf{G}_{\ell}=\sum_{\ell^{\prime}}\left\\{\frac{n_{\ell^{\prime}}}{[(\mathsf{V}\mathsf{M})_{\ell^{\prime}}]_{\rm
av}}-\frac{n_{\ell^{\prime}}}{[(\mathsf{V}\mathsf{M})_{\ell^{\prime}}]_{\rm
0}}\right\\}\mathsf{V}_{\ell^{\prime}\ell}\,.$ (46)

We employ Eq. (46) to determine the Lagrange multiplier $G_{2}(X,\theta)$ in
terms of $[\mathsf{M}]_{\rm av}$ obtained from ED and the Haar average

$[\mathsf{M}_{\ell}]_{\rm 0}=q^{-L}(q^{\ell}-q^{-\ell})\,.$ (47)

This approach is used to obtain the data shown in Fig. 5.

Note [from the definition of $\mathsf{T}_{\ell}$ below Eq. (43) and the
discussion following Eq. (44)] that the probability distribution of
$\mathsf{T}_{\ell}$ is consistent with a Gaussian distribution for the matrix
elements $\langle a|X_{\alpha}|b\rangle$, as expected from ETH, with a
variance controlled by the Lagrange multipliers $G_{2}(X,\theta)$.

####  III.2.3 Iterative method for determining Lagrange multipliers

We next describe a straightforward method for determining the Lagrange
multipliers without making use of information about the probability
distributions of $M_{X}(abcd)$ and $M_{X}(abba)$. We treat the case of
$G_{4}(X,Y,\theta)$; the necessary modifications for $G_{2}(X,\theta)$ are
obvious.

Our starting point is a perturbative expansion of Eq. (32) to first order in
$G(X,Y,\theta)$, which yields

$\displaystyle\begin{split}[M_{X}M_{Y}^{*}]_{\rm av}=&[M_{X}M_{Y}^{*}]_{0}\\\
&-\sum_{X^{\prime}Y^{\prime}}[\mathsf{K}_{4}]_{XY,X^{\prime}Y^{\prime}}G(X^{\prime},Y^{\prime},\theta)\,,\end{split}$
(48)

where $[\mathsf{K}_{4}]_{XY,X^{\prime}Y^{\prime}}$ is the connected correlator

$\displaystyle[\mathsf{K}_{4}]_{XY,X^{\prime}Y^{\prime}}$
$\displaystyle=[M_{X}M_{Y}^{*}M_{X^{\prime}}M_{Y^{\prime}}^{*}]_{0}$
$\displaystyle\quad-[M_{X}M_{Y}^{*}]_{0}[M_{X^{\prime}}M_{Y^{\prime}}^{*}]_{0}\,.$
(49)

Introducing the abbreviations

$\displaystyle\begin{split}[\mathsf{V}_{4}]_{XY}&=[M_{X}M_{Y}^{*}]_{0}-[M_{X}M_{Y}^{*}]_{\rm
av}\\\ \mbox{and}\quad[\mathsf{G}_{4}]_{XY}&=G_{4}(X,Y,\theta),\end{split}$
(50)

this gives at first order in perturbation theory

$\mathsf{G}_{4}=\mathsf{K}_{4}^{-1}\mathsf{V}_{4}.$ (51)

To go beyond first order perturbation theory, we define an iterative procedure
based on Eq. (51). Let $[M_{X}M^{*}_{Y}]^{(n)}_{\rm MC}$ denote the MC result
obtained with the $n$th approximant $G_{4}^{(n)}(X,Y,\theta)$ as Lagrange
multiplier, and let $[M_{X}M^{*}_{Y}]_{\rm ED}$ be the value from ED. Then
iterate for $n=1,2\ldots$

$\displaystyle[\mathsf{V}_{4}^{(n+1)}]_{XY}=[M_{X}M^{*}_{Y}]^{(n)}_{\rm
MC}-[M_{X}M^{*}_{Y}]_{\rm ED}$ (52)

with

$\displaystyle[\mathsf{V}_{4}^{(1)}]_{XY}=[M_{X}M^{*}_{Y}]_{0}-[M_{X}M^{*}_{Y}]_{\rm
ED}\,.$ (53)

and

$\displaystyle\mathsf{G}^{(n)}_{4}=\mathsf{G}^{(n-1)}_{4}+\mathsf{K}_{4}^{-1}\mathsf{V}^{(n)}_{4}\,.$
(54)

with

$\displaystyle\mathsf{G}^{(0)}_{4}=0\,.$ (55)

A possible refinement is to replace Eq. (54) with

$\displaystyle\mathsf{G}^{(n)}_{4}=\mathsf{G}^{(n-1)}_{4}+\alpha\mathsf{K}_{4}^{-1}\mathsf{V}^{(n)}_{4}\,.$
(56)

where $0<\alpha\leq 1$ is a real parameter. Small $\alpha$ reduces the risk of
overshooting the solution at the expense of a slower convergence rate.

This method is used to produce the data shown in Fig. 14.

###  III.3 Numerical methods

In this section, we give details of the ED, MC and ensemble averaging
procedures used to obtain the data presented in this paper.

####  III.3.1 Exact diagonalization

In order to study eigenstate correlations, we use ED of the Floquet operator
$W$ to compute exact eigenvectors and hence averages of
$M^{\phantom{*}}_{X}(abcd)M_{Y}^{*}(abcd)$ for all choices of our selected
subsystems $X$ and $Y$. We obtain phase-resolved averages by dividing the
phase interval $[-\pi,\pi]$ into $64$ bins. Each realisation of $W$ with $L$
spins generates approximately $q^{4L}/4!$ different quadruples of eigenvectors
and from these we randomly choose up to $10^{6}$ quadruples. We average over
between 150 (for $L=12$) and 1000 (for $L=10$ and $L=8$) realisations of $W$
to obtain the ED data in Fig. 3, and over 25000 realisations for the results
shown in Fig. 4. The symmetry relation $F_{4}(X,Y,\theta)=F_{4}(X,Y,-\theta)$
allows us to restrict calculations to $\theta\geq 0$. Similarly, in the case
of $M_{X}(abba)$ we take $10^{6}$ tuples out of the approximately $q^{2L}/2$
possibilities for each realisation of $W$ and $L=12$. We average over 1000
realisations for $L=12$ to obtain the data in Fig. 5.

From these ED results we compute the Lagrange multipliers $G_{4}(X,Y,\theta)$
and $G_{2}(X,\theta)$ using the procedures described in Sec. III.2. We find
that a further symmetrisation of the data, using the spatial symmetry under
the interchange of $X$, $Y$ with $\overline{X}$, $\overline{Y}$, improves
stability.

####  III.3.2 Monte-Carlo sampling

In order to determine the eigenstate correlator $F_{4}(X,Y,\theta)$ from the
JDF [Eq. (32)] for four eigenstates $|a\rangle$, $|b\rangle$, $|c\rangle$ and
$|d\rangle$ we use Monte Carlo sampling with $e^{-S_{4}(a,b,c,d)}$ [Eq. (33)]
as the weighting term. Similarly, to determine the correlator
$F_{2}(X,\theta)$ from the JDF [Eq. (30)] for two eigenstates $|a\rangle$ and
$|b\rangle$ we use Monte Carlo sampling with $e^{-S_{2}(a,b)}$ [Eq. (31)] as
the weighting term.

We follow the Metropolis-Hastings algorithm to obtain a Markov chain of vector
quadruples $(a,b,c,d)$ distributed according to $P_{4}$. To generate the next
quadruple $(a^{\prime},b^{\prime},c^{\prime},d^{\prime})$, we use a random
unitary rotation of the vectors in the previous sample of the form
$V=\exp(i\epsilon A)$, where $A$ is a random Hermitian matrix drawn from the
Gaussian unitary ensemble with unit variance. We set $\epsilon=0.1$ for $L=12$
and $q=2$, and $\epsilon=0.8$ for $L=8$ and $q=2$. Since $V$ is unitary,
orthonormal vectors retain this property after the transformation.

This choice of update rule with tuning parameter $\epsilon$ allows us to
perform effective importance sampling, since we only propose relatively small
changes to the sample. The new sample is then accepted in the Markov chain
with probability

$P_{\text{accept}}=\min\left(1,\mathrm{e}^{-S_{4}(a,b,c,d)+S_{4}(a^{\prime},b^{\prime},c^{\prime},d^{\prime})}\right).$
(57)

For each $\theta$ we perform 2000 Monte-Carlo runs in parallel with up to
$10^{6}$ samples per run. The results are obtained by averaging over all runs.
In the case of small $\theta$, we find in rare cases (one out of 1000 runs)
instabilities towards local maxima of the weighting function during the
sampling process. This is visible by tracking $S_{4}(a,b,c,d)$ or the
acceptance rate. This problem can be circumvented by decreasing the size of
the update steps, but only at the expense of longer autocorrelation times. As
a compromise, we discard runs where $S_{4}$ falls below a threshold value of
$-50$.

Figure 8: Probability distribution of $S_{4}(a,b,c,d)$ [Eq. (33)] for
quasienergy differences $\theta$ as indicated: comparison between ED results
(black dashed lines) and MC results (solid coloured lines). Parameters as in
Fig. 3.

To provide an overall test of our form [Eq. (32)] for the JDF, we show in Fig.
8 a comparison of the distributions of $S_{4}(a,b,c,d)$ obtained respectively
from ED and from MC sampling of the JDF. The excellent agreement between the
two distributions over a range of values for $\theta$ is evidence of the
internal consistency of our approach.

##  IV Perturbative treatment of Lagrange multipliers

An obvious approach to calculations based on the joint eigenstates
distributions of Eqns. (30) and (32) is a perturbative expansion in powers of
the Lagrange multipliers $G_{2}(X,\theta)$ and $G_{4}(X,Y,\theta)$. In this
section we set out a general framework for such an expansion and apply it in
several ways. While the expansion does not generate fundamentally new results,
it provides a useful perspective that is complementary to the one set out in
Sec. II and Sec. III.

We use the expansion to provide an alternative justification of the fitting
procedure for $G_{4}(X,Y,\theta)$ to the one described in Sec. III.2.1. We
also use it to consider Eq. (32) without the simplification employed in Sec.
II.6 of setting $G_{2}(X,\theta)=0$. We show that for large $q^{L}$ and
$q^{L(X)}$ that if $F_{2}(X,\theta)$ is calculated from the JDF for four
eigenstates rather than two, the influence of $G_{4}(X,Y,\theta)$ on the
result is small. This means that the predictions of ETH for matrix elements
between pairs of eigenstates are only weakly affected by the correlations
between sets of four eigenstates that are introduced with our Ansatz for the
JDF. Finally, we show that the effect of the Lagrange multipliers
$G_{2}(X,\theta)$ and $G_{4}(X,Y,\theta)$ on the normalisation of eigenstates
drawn from the joint distributions is small in large systems.

An exact perturbative expansion requires the evaluation of averages over a
Haar distribution of vectors, denoted by $P_{1}^{0)}(a)$, $P_{2}^{(0)}(a,b)$
and $P_{4}^{(0)}(a,b,c,d)$ in Sec. II.5. While the formalism required for this
is well developed (see e.g. [45]), it is quite cumbersome. Moreover, we
require results only at leading order for $q^{L}$ large. These can be obtained
by substituting in place of the Haar distribution, one in which each vector
component is an independent Gaussian random variable. Specifically, consider a
computational basis $\\{\ket{k_{c}}\\}$ and denote the overlap of the
eigenstate $\ket{a}$ with the basis state $\ket{k_{c}}$ by
$a(k)=\braket{k_{c}}{a}$. Define $S_{1}(a)$ by

$\displaystyle S_{1}(a)=q^{L}\sum_{k}|a(k)|^{2}\,.$ (58)

We replace $P_{4}^{(0)}(a_{1},a_{2},a_{3},a_{4})$ by

$P^{(\rm
G)}_{4}(a_{1},a_{2},a_{3},a_{4})=\left(\frac{q^{L}}{\pi}\right)^{4q^{L}}e^{-\sum_{i}S_{1}(a_{i})}\,.$
(59)

Although different vectors drawn from this distribution are not in general
exactly orthonormal, orthonormality is recovered in the limit
$q^{L}\rightarrow\infty$. We denote averages with respect to this Gaussian
distribution by $[\ldots]_{\rm G}$.

###  IV.1 Diagrammatic notation

It is useful to employ diagrammatic notation for these Gaussian averages. As
in Fig. 1, eigenstates are represented by circles and index contractions
following from the definition of $M_{X}(abcd)$ are indicated by solid lines
carrying arrows that run from states $|a\rangle$ towards conjugate states
$\langle b|$; these lines carry labels to indicate the subsystem within which
the contraction is done. The combination of circles and full lines, with
labels for states and subsystems, is fixed by the choice of quantity we
average, and diagrams are generated by making all possible Wick pairing of
circles. These pairings are represented by dashed lines. The contribution of a
diagram is a product of two factors. One factor stems from Eq. (59) and
consists of $q^{-L}$ for every dashed line. The other factor arises from sums
over the Hilbert space at each site. To evaluate this factor we form closed
loops in the diagram consisting alternately of dashed lines and full lines
traversed in the direction of the arrows. These full lines must carry the
label that corresponds to a subsystem in which the site lies. For each site we
have a factor of $q$ for every such closed loop.

For example, the evaluation of $[M_{X}(abba)]_{\rm G}$ is shown in Fig. 9.

Figure 9: Diagrammatic representation and evaluation of $[M_{X}(abba)]_{\rm
G}$: (a) full diagram; (b) decomposition into loops for a site in subsystem
$X$; (c) decomposition into loops for a site in subsystem $\overline{X}$.

It yields

$[M_{X}(abba)]_{\rm G}=q^{-2L}q^{2L(X)}q^{L(\overline{X})}=q^{L(X)-L}\,.$ (60)

###  IV.2 Relating $F_{4}(X,Y,\theta)$ and $G_{4}(X,Y,\theta)$

Figure 10: Leading connected contributions to
$\big{[}[\mathsf{M}]_{X}[\mathsf{M}]_{Y}[\mathsf{M}^{\rm
T}\mathsf{G}_{4}\mathsf{M}^{*}]^{n}\big{]}_{\rm G}$ for: (i) $n=0$, (ii) $n=1$
and (iii) $n=2$. The pair of squares in (i) and the left-most pairs of squares
in (ii) and (iii) represent $[\mathsf{M}]_{X}[\mathsf{M}]_{Y}$; other pairs of
squares represent $\mathsf{M}^{\rm T}\mathsf{G}_{4}\mathsf{M}^{*}$. The dashed
lines in (ii) and (iii) that leave the diagrams on the left are joined to the
dashed lines that leave on the right.

Our first objective is to obtain a relationship between $F_{4}(X,Y,\theta)$
and $G_{4}(X,Y,\theta)$ by expanding $e^{-S_{4}(a,b,c,d)}$ in a power series,
averaging with respect to the Gaussian distribution, and then resumming the
terms at each order in the expansion that are leading for $q^{L}$ and
$q^{L(X)}$ large. By this means we will recover Eq. (38).

Using the notation of Eq. (35) we require the connected contributions to
$\big{[}[\mathsf{M}]_{X}[\mathsf{M}]_{Y}[\mathsf{M}^{\rm
T}\mathsf{G}_{4}\mathsf{M}^{*}]^{n}\big{]}_{\rm G}$ at each order $n$ in
perturbation theory that are leading for large system and subsystem sizes.
These come from the Wick contractions that generate the largest number of
loops in a decomposition of the type illustrated in Fig. 9. These contractions
are illustrated for $n=0,1$ and $2$ in Fig. 10, establishing an obvious
pattern for general $n$. Retaining only these terms we have

$\displaystyle\big{[}[\mathsf{M}]_{X}[\mathsf{M}^{*}]_{Y}\big{]}_{\rm av}$
$\displaystyle\approx$
$\displaystyle\sum_{n=0}^{\infty}(-1)^{n}[(\mathsf{G}^{(0)}_{4})^{-1}(\mathsf{G}_{4}(\mathsf{G}^{(0)}_{4})^{-1})^{n}]_{X,Y}$
(61) $\displaystyle=$
$\displaystyle\big{[}(\mathsf{G}^{(0)}_{4})^{-1}[\openone+\mathsf{G}_{4}(\mathsf{G}^{(0)}_{4})^{-1}]^{-1}\big{]}_{X,Y}$
$\displaystyle=$
$\displaystyle\big{[}(\mathsf{G}_{4}+\mathsf{G}^{(0)}_{4})^{-1}\big{]}_{X,Y}$

where $\mathsf{G}^{(0)}_{4}$ is as defined following Eq. (36). [Note that for
these terms the factor of $(n!)^{-1}$ arising from the power series expansion
of $e^{-S_{4}(a,b,c,d)}$ is cancelled by a combinatorial factor arising in the
pairing of terms in Fig. 10.] Comparing Eqns. (38) and (61) we see that this
diagrammatic resummation provides an alternative derivation of the main
results of Sec. III.2.1.

###  IV.3 Cross-correlations

In our MC studies of the eigenstate JDF [Eq. (32)] we have made two
simplifications: one is to omit $G_{4}(X,Y,\theta)$ when studying
$F_{2}(X,\theta)$, and the other is to omit $G_{2}(X,\theta)$ when studying
$F_{4}(X,Y,\theta)$. The magnitudes of the resulting errors can be assessed
using perturbation theory, as we now discuss.

We begin by considering the effect of $G_{4}(X,Y,\theta)$ on the value of
$F_{2}(X,\theta)$, or equivalently on the value of $[M_{X}(abba)]_{\rm av}$,
as follows. To first order in perturbation theory in $G_{4}(X,Y,\theta)$ we
have

$\displaystyle[M_{X}(abba)]_{\rm av}=$ $\displaystyle[M_{X}(abba)]_{\rm G}$
$\displaystyle-\sum_{X^{\prime}Y^{\prime}}G_{4}(X^{\prime},Y^{\prime},$
$\displaystyle\theta)[M_{X}(abba)M_{X^{\prime}}(abcd)M^{*}_{Y^{\prime}}(abcd)]_{G,c}\,,$
(62) $\displaystyle+{\cal O}([G_{4}(X^{\prime},Y^{\prime},\theta]^{2})\,,$

where $[\ldots]_{\rm G,c}$ denotes the connected average, defined by

$\displaystyle[M_{X}(abba)M_{X^{\prime}}(abcd)M_{Y^{\prime}}^{*}(abcd)]_{\rm
G,c}=$
$\displaystyle{[M_{X}(abba)M_{X^{\prime}}(abcd)M_{Y^{\prime}}^{*}(abcd)]_{\rm
G}}$ $\displaystyle-{[M_{X}(abba)]_{\rm
G}[M_{X^{\prime}}(abcd)M_{Y^{\prime}}^{*}(abcd)]_{\rm G}}$ $\displaystyle\,.$
(63)

The diagrams that contribute to this connected average are shown in Fig. 11.
Evaluating these diagrams for the representative case
$X=X^{\prime}=Y^{\prime}$, we obtain

$\displaystyle[M_{X}(abba)$ $\displaystyle M_{X}(abcd)M_{X}^{*}(abcd)]_{\rm
G,c}$ $\displaystyle=q^{-3L}$
$\displaystyle(q^{-L(X)}+2q^{-L(\overline{X})})\,.$ (64)

We compare the zeroth-order term, which is

$[M_{X}(abba)]_{\rm G}=q^{L(X)-L}$ (65)

from Eq. (60), with the first-order term, using Eq.(IV.3) and
$G_{4}(X,X,\theta)\sim q^{2L}$ from Eq. (36), which gives

$\displaystyle\sum_{X^{\prime}Y^{\prime}}G_{4}(X^{\prime},Y^{\prime}$
$\displaystyle,\theta)[M_{X}(abba)M_{X^{\prime}}(abcd)M^{*}_{Y^{\prime}}(abcd)]_{G,c}\,,$
$\displaystyle\sim$ $\displaystyle q^{-L}(q^{-L(X)}+2q^{-L(\overline{X})})\,.$
(66)

Hence we see that the effect of $G_{4}(X^{\prime},Y^{\prime},\theta)$ on
$F_{2}(X,\theta)$ is small provided that $q^{L(X)}$ and $q^{L(\overline{X})}$
are large.

We have not systematically investigated higher order terms in Eq. (IV.3), but
we expect them individually to be small: note that although
$G_{4}(X^{\prime},Y^{\prime},\theta)\sim q^{2L}$, each such contribution is
accompanied by a factor of $q^{-4L}$ from extra dashed lines, as well as
diagram-dependent factors from sums over the Hilbert space at each site.

A similar study of the effect of $G_{2}(X,\theta)$ at first order in
perturbation theory on $F_{4}(X,Y,\theta)$, or equivalently on
$[M_{X^{\prime}}(abcd)M_{Y^{\prime}}^{*}(abcd)]_{\rm av}$ reaches a different
conclusion: the influence is not small in $q^{L}$ or $q^{L(X)}$, but only in
powers of $G_{2}(X,\theta)$. This is not unexpected: the functional form for a
correlator involving four eigenstates and depending on three quasienergy
differences has been discussed previously [46] and in Sec. II.3; it involves
factors of both $F_{2}(X,\theta)$ and $F_{4}(X,Y,\theta)$ (see Eq. (27)). We
do not pursue this further since we have chosen here to study models in which
$F_{2}(X,\theta)$ lies close to its value for Haar-distributed pairs of
eigenstates, and $G_{2}(X,\theta)$ is therefore small.

Figure 11: Contributions to the connected average defined in Eq. (IV.3)

###  IV.4 Renormalisation of propagators and vertices

An alternative perspective is provided by considering the perturbative
renormalisation of the propagators and vertices appearing in the JDF. We begin
with the former: the bare propagator in the theory, represented using dashed
lines in the figures, is generated by Eqns. (58) and (59), and carries a
factor of $q^{-L}$. More generally, it acquires a self-energy $\Sigma$ and the
factor becomes $(q^{L}+\Sigma)^{-1}$. Our aim is to evaluate $\Sigma$ at
leading order in $G_{2}(X,\theta)$ and $G_{4}(X,Y,\theta)$ and compare it to
the bare inverse propagator $q^{L}$: if it is small under this comparison,
then the effects of the vertices on the normalisation of eigenstates can be
neglected. The diagrams contributing to $\Sigma$ at this order are shown in
Fig. 12. They are diagonal matrices in the site basis and the magnitudes of
their largest entries are respectively

$\Sigma^{\rm(i)}_{kk}=q^{-L}\sum_{X}G_{2}(X,\theta)q^{L(X)}\sim q^{2L(X)}$
(67)

[where we have used Eqns. (46) and (47) to estimate $G_{2}(X,\theta)\sim
q^{L+L(X)}$] and

$\Sigma^{\rm(ii)}_{kk}=q^{-4L}\sum_{X,Y}q^{L(X,Y)}G_{4}(X,Y,\theta)\sim{\cal
O}(1)\,,$ (68)

[where we have used an estimate of $G_{4}(X,Y,\theta)$ given above Eq.
(IV.3)]. From this we see that the renormalisation of the propagator is indeed
small if $q^{L}$ is large and if $q^{2L(X)}\ll q^{L}$. In the opposite regime
($q^{2L(X)}>q^{L}$) we believe that $G_{2}(X,\theta)$ is small since
$F_{2}(X,\theta)$ approaches the value it takes for Haar-distributed
eigenstates, and that $\Sigma^{\rm(i)}_{kk}\ll q^{L}$ notwithstanding the
estimate of Eq. 67.

Figure 12: Contributions to the self-energy $\Sigma$ at first order: (i) from
$G_{2}(X,\theta)$, and (ii) from $G_{4}(X,Y,\theta)$.

A similar discussion can be developed of renormalised vertices generated by
combining contributions from $S_{2}(a,b)$ and $S_{4}(a,b,c,d)$ in (32) and
forming sufficient internal contractions to generate a new effective
contribution to either $S_{2}(a,b)$ or $S_{4}(a,b,c,d)$. For example, one (out
of three possible terms) contributing in this way to $S_{2}(a,b)$ at first
order in both $G_{2}(X,\theta)$ and $G_{4}(X,Y,\theta)$ is shown in Fig. 13.
We do not pursue this further because this is exactly the same phenomenon as
has been discussed from a different perspective in Sec. IV.3. Note that it is
not necessary to consider renormalisation of $S_{2}(a,b)$ in powers of
$G_{2}(X,\theta)$ alone, or of $S_{4}(a,b,c,d)$ in powers of
$G_{4}(X,Y,\theta)$ alone, since these effects (which are not generally small)
are covered in full by the approaches described in Sec. IV.2 and in Sec.
III.2.

Figure 13: One of three contributions to the renormalisation at first order in
both $G_{2}(X,\theta)$ and in $G_{4}(X,Y,\theta)$ of the vertex that appears
in $S_{2}(a,b)$.

##  V Results for further models

In this section we provide numerical results for brickwork models additional
to the one treated in Sec. II. In Sec. V.1 we omit the cutoff in the operator
purity of gates introduced in Sec. III.1 and present results for gates drawn
from a Haar-distribution. In Sec. V.2 we give results with local Hilbert space
dimension $q=3$ and Haar-distributed gates.

###  V.1 Results for $q=2$ with Haar gates

Figure 14: $F^{\rm MC}_{4}(X,Y,\theta)$ (open circles from MC) and
$F_{4}(X,Y,\theta)$ (lines from ED) vs $\theta$ for various $s$, for a
brickwork model with Haar-distributed gates, $q=2$ and $L=8$.

As discussed in Sec. III, the determination of the Lagrange multipliers
$G_{4}(X,Y,\theta)$ in the brickwork model with Haar-distributed gates is
complicated by the presence of weak links. Specifically, we find that the
simplifying assumption used in Sec. III.2.1, namely that the probability
distribution of $M_{X}(abcd)$ is approximately Gaussian, does not hold for
Haar-distributed gates. Instead, this distribution exhibits long tails, as we
demonstrate in Appendix A.2.

To determine $G_{4}(X,Y,\theta)$ under these circumstances, we use the
iterative fitting procedure introduced in Sec. III.2.3, which is not
predicated on a particular form for the probability distribution of
$M_{X}(abcd)$. We take all subsystems $X$ and $Y$ that can be obtained using a
single cut from a system with open boundary conditions. In contrast to Sec.
II, where the requirement of an approximately Gaussian distribution for
$M_{X}(abcd)$ led to the restriction $2<L(X)<L-2$, here we include all
subsystem sizes $1\leq L(X)\leq L-1$. A disadvantage of the iterative fitting
procedure is that it requires multiple MC evaluations of $F_{4}^{\rm
MC}(X,Y,\theta)$, which is slow if the total number of degrees of freedom
$q^{L}$ involved is large; this restricts us to $L=8$ with $q=2$. In this case
estimates $F^{\rm MC}_{4}(X,Y,\theta)$ are obtained using between 60 and 300
iterations and a step size $\alpha=0.2$. The MC results shown in Fig. 14
display very good agreement with ED data.

###  V.2 Results for $q=3$ brickwork model

Figure 15: $F^{\rm MC}_{4}(X,Y,\theta)$ (open circles from MC) and
$F_{4}(X,Y,\theta)$ (lines from ED) vs $\theta$ for various $s$, for a
brickwork model with Haar-distributed gates, $q=3$ and $L=8$.  Figure 16: Test
of JDF fitted to behaviour in the geometries of Fig. 15 (b) but applied to
geometries of Fig. 16 (b), for the brickwork model with $q=3$ and $L=8$. (a)
Comparison of data from MC (open circles) and ED (solid lines). (b) Partition
used for (a), in which the 8-site system is divided by two spatial cuts into a
two-site subsystem $X$ and its complement $\overline{X}$, or a two-site
subsystem $\overline{Y}$ and its complement $Y$.

In our study of the brickwork model with $q=3$ and Haar-distributed gates we
use system size $L=8$ so that ED calculations are straightforward. The effect
of weak links becomes less pronounced with increasing bond dimension $q$ and
as a consequence the single-shot approach to obtain $G_{4}(X,Y,\theta)$ of
Sec. III.2.1 becomes more accurate. Conversely, for a given system size the
iterative method of Sec. III.2.3 is more difficult to apply with increasing
$q$, because a large number of samples is required to obtain accurate
estimates for $F^{\rm MC}_{4}(X,Y,\theta)$ when the Hilbert space dimension
$q^{L}$ of the system is large (see the discussion of Sec. II.4 and Ref.
[21]). We therefore determine Lagrange multipliers for this model using the
single-shot approach, taking subsystems $X$ obtained from a system with open
boundary conditions by making a single cut, and with $L(X)\geq 2$,
$L(\overline{X})\geq 2$. The MC results shown in Fig. 15 display excellent
agreement with ED data.

Finally, we test how well the JDF fitted to the geometries of Fig. 14(b) can
reproduce the behaviour of the OTOC in the geometries of Fig. 16(b). The
results shown in Fig. 16(a) display excellent agreement between ED and MC
data.

##  VI Summary and outlook

In this work we have analysed the interplay between the statistical properties
of the time evolution operator for chaotic many-body quantum systems, and the
quantum information dynamics that these systems display. The eigenstate
thermalisation hypothesis provides an accurate description of some key
aspects, in terms of the probability distribution of matrix elements of
observables between eigenstates of the time-evolution operator. In its
original form, however, it does not capture the consequences of a finite speed
for quantum information spreading. To remedy this, we advocate a change of
viewpoint. In place of the matrix elements considered within ETH, we construct
correlators from eigenstates. We show that the simplest such correlator
capturing signatures of quantum information dynamics at long times and large
distances is unique and involves a set of four eigenstates. We also propose an
Ansatz for the joint probability distribution function of small numbers of
eigenstates, in which the values of correlators are controlled by Lagrange
multipliers. We support these general ideas using numerical studies of Floquet
quantum circuits, showing firstly that correlators have the expected features,
and secondly that it is possible to choose values for the Lagrange multipliers
so that our Ansatz reproduces these features accurately. We believe our
approach is complementary to recent work [36, 37] that generalises ETH using
the language of Free Probability theory.

An advantage of viewing this problem in terms of correlations between
eigenstates rather than matrix elements is that even the simplest assumption,
of a Haar distribution for eigenstates, yields correlations with the correct
order of magnitude. In this instance, however, the correlations are
independent of eigenvalues. By contrast, a microscopic model for local quantum
dynamics yields a characteristic dependence of correlators on differences in
quasienergies, as we have summarised in Fig. 3. We have shown that our Ansatz
for the joint distribution of eigenstates, together with a suitable choice for
the Lagrange multipliers, captures this dependence on eigenvalues differences.

The quantity that characterises correlations between sets of four eigenstates,
denoted $F_{4}(X,Y,\theta)$ in this paper, is a function not only of the
eigenphase difference $\theta$, but also of subsystem choices $X$ and $Y$.
This constitutes a large set of possibilities if there are no restrictions on
the way the subsystems are selected. We have focussed on the simple set of
subsystems that can be obtained from an overall system with open (rather than
periodic) boundary conditions by making a single spatial cut, choosing
Lagrange multipliers in our Ansatz for the JDF so that $F_{4}(X,Y,\theta)$ is
reproduced accurately for these subsystems. In addition we have shown that the
eigenstate correlations imposed in this way are sufficient to reproduce to a
good approximation the correlators for some other choices of $X$ and $Y$. In
particular, we have demonstrated that the behaviour of the OTOC for operators
supported on a few sites of the system (requiring choices of $X$ and $Y$ each
involving two cuts) is well approximated by our JDF.

Several obvious directions remain open for future work. Perhaps most
importantly, while our discussion in this paper has been restricted to Floquet
systems, it would be desirable to extend the approach to systems with a time-
independent Hamiltonian. In our context, the significance of such an extension
is that diagonal matrix elements of local observables in the basis of
Hamiltonian eigenstates are generically functions of energy, a feature absent
from Floquet systems. Indeed, ETH is formulated in part to describe this
energy dependence. A natural modification of our JDF to impose such an energy
dependence is to supplement Eqns. (30) and (32) with an extra factor
$e^{-\sum_{k}S_{1}(a_{k})}$ chosen to bias selected eigenstates towards a pre-
specified energy shell. An obvious choice is to take $S_{1}(a)=\beta\langle
a|H|a\rangle$, where $H$ is the Hamiltonian and $\beta$ is a Lagrange
multiplier.

It is worth emphasising that the correlations induced respectively by
$S_{1}(a)$, $S_{2}(a,b)$ and $S_{4}(a,b,c,d)$ are associated with widely
separated energy scales. Taking the characteristic strength of local
interactions as the unit of energy, $S_{1}(a)$ selects vectors from a broad
energy window, which has a width that increases with system size. In turn,
$S_{2}(a,b)$ generates correlations between pairs of vectors within an energy
window with width of order unity. Finally $S_{4}(a,b,c,d)$ generates
correlations between groups of four vectors with energy or quasienergy
differences lying in a narrow energy window, whose width decreases
indefinitely with increasing separation between the subsystems $X$ and $Y$ in
$F_{4}(X,Y,\theta)$.

A further direction for future work is to examine correlations between sets of
$n$ eigenstates with $n>4$. This opens up many possibilities, since
generalisations of the correlators $F_{2}(X,\theta)$ and $F_{4}(X,Y,\theta)$
may involve multiple subsystems in their definition. Restricting to a pair of
subsystems $X$ and $Y$, these higher-order correlators are required, for
example, to describe the higher-order Rényi entropies of the operator
entanglement for the time-evolution operator. Beyond this, one can ask whether
there are physical phenomena exposed only by higher-order correlators.

## Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft through the
cluster of excellence ML4Q (EXC2004, project-id 390534769) and by the UK
Engineering and Physical Sciences Research Council through Grants EP/N01930X/1
and EP/S020527/1. We also acknowledge support from the QuantERA II Programme,
which has received funding from the European Union’s Horizon 2020 research
innovation programme (GA 101017733), and from the DFG through the project
DQUANT (project-id 499347025), and from the National Science Foundation under
Grant No. NSF PHY-1748958. We thank S. Parameswaran for useful comments.

##  Appendix A Distribution of $M_{X}(abcd)$

$X$$k$

Figure 17: Probability distribution of $|M_{X}(abcd)|$ for different
subsystems $X$ with $L(X)=k$ and $L=12$, $q=2$, and fitted complex Gaussian
distributions (dashed black) as guide for the eye: (a) $\theta=0.1$ and (b)
$\theta=1.7$. Brickwork model with gates sampled from the distribution defined
in III.1 and a cutoff $0.3\times q^{4}$ for the operator purity. With
decreasing partition size $L(X)$, the distribution of $M_{X}(abcd)$ exhibits
non-Gaussian tails.

$X$$k$

Figure 18: Probability distribution of $|M_{X}(abcd)|$ at partition $k=2$ and
$L=12$, $q=2$, and a complex Gaussian distribution (dashed black) as guide for
the eye: (a) $\theta=0.1$ and (b) $\theta=1.7$. Brickwork model with unitary
gates sampled from the distribution defined in III.1 and values of the cutoff
$c$ as indicated. With increasing cutoff, the distributions show non-Gaussian
tails for small relative phases $\theta$.

$x$$y$$X$$\overline{X}$$\overline{Y}$$Y$

Figure 19: Relative deviation between $F^{MC}_{4}(X,Y,\theta)$ and
$F_{4}(X,Y,\theta)$ for $q=2$, and $L=8$ (left column), $L=10$ (center column)
and $L=12$ (right column), and for $\theta=0.1$ (top row) and $\theta=2.03$
(bottom row).

In this appendix we examine the probability distribution of $M_{X}(abcd)$,
providing further details beyond the information given in Fig. 7. This
information is important because one of the two methods that we use for
determining the Lagrange multipliers $G_{4}(X,Y,\theta)$ relies on this
distribution having a Gaussian form. In Appendix A.1 we show how the form of
the distribution varies with the value of $L(X)$. In Appendix A.2 we
investigate the effect on the distribution of the cutoff in the operator
entanglement purity of two-qubit gates introduced in Sec. III.1.

###  A.1 Effect of the subsystem size $L(X)$

The distribution of $M_{X}(abcd)$ for different partitions $X$ is shown in
Fig. 17 (a) and (b). With increasing partition size $L(X)$, the distribution
approaches a Gaussian. This supports our observation that the method we use to
determine $G_{4}(X,Y,\theta$ becomes more accurate with increasing $L(X)$.

###  A.2 Effect of the cutoff in the operator entanglement purity

In Sec. III we presented results for a model (see Sec. III.1) that is well
adapted to our procedure for finding $G_{4}(X,Y,\theta$). This model is
designed to have a probability distribution for $M_{X}(abcd)$ that is close to
Gaussian, and has gates drawn from a truncated version of the Haar
distribution, with a cutoff on the operator entanglement purity of $c\times
q^{4}$.

In Fig. 18 we examine the effect of the value of this cutoff on the
probability distribution of $M_{X}(abcd)$, considering the range from $c=0.3$
(the value used in Sec. III) to $c=1$ (an unrestricted Haar distribution). The
distribution shows non-Gaussian tails for small relative phase $\theta$ and a
Haar-random Floquet model. The tails are suppressed with decreasing cutoff in
the operator entanglement purity. We attribute this effect to weak links $i$
on which the gate $w_{i}$ is close to the identity (especially in small
systems or with small subsystems). The effect of such weakly entangling gates
on the dynamics of quantum information was studied in more detail in [47].

##  Appendix B Accuracy of our approach with increasing system size

Here we provide complete information on deviations between
$F^{MC}_{4}(X,Y,\theta)$ and $F_{4}(X,Y,\theta)$ with Lagrange Multipliers
$G_{4}(X,Y,\theta)$ determined using the single-shot procedure described in
Sec. III.2.1. We compare the results for system sizes $L=8,\,10,\,12$ and all
possible partitions defined by one cut that have $L(X),L(\overline{X})>2$.

The relative deviations are shown in Fig. 19. In all cases the relative
deviation is less than $10\%$, showing the accuracy of our approach.
Furthermore, the accuracy improves with increasing system size. We see the
largest deviations for small partition sizes $L(X)$ and $L(Y)$. This is
expected, as we discuss in Sec. IV: Correction terms in the perturbation
theory and contributions from $G_{2}(X,\theta)$ are most significant when
$L(X)$ and $L(\overline{X})$ are small.

## References

  * Peres [1984] A. Peres, Ergodicity and mixing in quantum theory. I, Phys. Rev. A 30, 504 (1984).
  * Feingold  _et al._ [1984] M. Feingold, N. Moiseyev, and A. Peres, Ergodicity and mixing in quantum theory. II, Phys. Rev. A 30, 509 (1984).
  * Feingold and Peres [1986] M. Feingold and A. Peres, Distribution of matrix elements of chaotic systems, Phys. Rev. A 34, 591 (1986).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Beugeling  _et al._ [2014] W. Beugeling, R. Moessner, and M. Haque, Finite-size scaling of eigenstate thermalization, Phys. Rev. E 89, 042112 (2014).
  * Beugeling  _et al._ [2015] W. Beugeling, R. Moessner, and M. Haque, Off-diagonal matrix elements of local operators in many-body quantum systems, Phys. Rev. E 91, 012144 (2015).
  * Mondaini and Rigol [2017] R. Mondaini and M. Rigol, Eigenstate thermalization in the two-dimensional transverse field Ising model. II. Off-diagonal matrix elements of observables, Phys. Rev. E 96, 012157 (2017).
  * Khaymovich  _et al._ [2019] I. M. Khaymovich, M. Haque, and P. A. McClarty, Eigenstate Thermalization, Random Matrix Theory, and Behemoths, Phys. Rev. Lett. 122, 070601 (2019).
  * Alba [2015] V. Alba, Eigenstate thermalization hypothesis and integrability in quantum spin chains, Phys. Rev. B 91, 155123 (2015).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From Quantum Chaos and Eigenstate Thermalization to Statistical Mechanics and Thermodynamics, Advances in Physics 65, 239 (2016).
  * Borgonovi  _et al._ [2016] F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky, Quantum chaos and thermalization in isolated systems of interacting particles, Physics Reports 626, 1 (2016).
  * Prosen [1999] T. Prosen, Ergodic properties of a generic nonintegrable quantum many-body system in the thermodynamic limit, Phys. Rev. E 60, 3949 (1999).
  * Luitz [2016] D. J. Luitz, Long tail distributions near the many-body localization transition, Phys. Rev. B 93, 134201 (2016).
  * Luitz and Bar Lev [2016] D. J. Luitz and Y. Bar Lev, Anomalous thermalization in ergodic systems, Phys. Rev. Lett. 117, 170404 (2016).
  * Serbyn  _et al._ [2017] M. Serbyn, Z. Papić, and D. A. Abanin, Thouless energy and multifractality across the many-body localization transition, Phys. Rev. B 96, 104201 (2017).
  * Roy  _et al._ [2018] S. Roy, Y. B. Lev, and D. J. Luitz, Anomalous thermalization and transport in disordered interacting Floquet systems, Phys. Rev. B 98, 060201 (2018).
  * Khatami  _et al._ [2013] E. Khatami, G. Pupillo, M. Srednicki, and M. Rigol, Fluctuation-Dissipation Theorem in an Isolated System of Quantum Dipolar Bosons after a Quench, Phys. Rev. Lett. 111, 050403 (2013).
  * Foini and Kurchan [2019] L. Foini and J. Kurchan, Eigenstate thermalization hypothesis and out of time order correlators, Phys. Rev. E 99, 042139 (2019).
  * Chan  _et al._ [2019a] A. Chan, A. De Luca, and J. Chalker, Eigenstate Correlations, Thermalization, and the Butterfly Effect, Phys. Rev. Lett. 122, 220601 (2019a).
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Quasiclassical Method in the Theory of Superconductivity, Soviet Physics JETP, Vol. 28, p.1200 28, 1200 (1969).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, Journal of High Energy Physics 2014, 67 (2014).
  * Brenes  _et al._ [2021] M. Brenes, S. Pappalardi, M. T. Mitchison, J. Goold, and A. Silva, Out-of-time-order correlations and the fine structure of eigenstate thermalization, Phys. Rev. E 104, 034120 (2021).
  * Jafferis  _et al._ [2023] D. L. Jafferis, D. K. Kolchmeyer, B. Mukhametzhanov, and J. Sonner, Matrix models for eigenstate thermalization, arXiv:2209.02130 (2023).
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, Commun.Math. Phys. 28, 251 (1972).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, J. Stat. Mech. 2005, P04010 (2005).
  * Kim and Huse [2013] H. Kim and D. A. Huse, Ballistic Spreading of Entanglement in a Diffusive Nonintegrable System, Phys. Rev. Lett. 111, 127205 (2013).
  * Zhou and D. J. Luitz [2017] T. Zhou and D. J. Luitz, Operator entanglement entropy of the time evolution operator in chaotic systems, Phys. Rev. B 95, 094206 (2017).
  * Dubail [2017] J. Dubail, Entanglement scaling of operators: a conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1d, Journal of Physics A: Mathematical and Theoretical 50, 234001 (2017).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, J. High Energ. Phys. 2016 (8), 106.
  * Prosen [1994] T. Prosen, Statistical Properties of Matrix Elements in a Hamilton System Between Integrability and Chaos, Annals of Physics 235, 115 (1994).
  * Dymarsky [2022] A. Dymarsky, Bound on Eigenstate Thermalization from Transport, Phys. Rev. Lett. 128, 190601 (2022).
  * Richter  _et al._ [2020] J. Richter, A. Dymarsky, R. Steinigeweg, and J. Gemmer, Eigenstate thermalization hypothesis beyond standard indicators: Emergence of random-matrix behavior at small frequencies, Phys. Rev. E 102, 042127 (2020).
  * Wang  _et al._ [2022] J. Wang, M. H. Lamann, J. Richter, R. Steinigeweg, A. Dymarsky, and J. Gemmer, Eigenstate Thermalization Hypothesis and Its Deviations from Random-Matrix Theory beyond the Thermalization Time, Phys. Rev. Lett. 128, 180601 (2022).
  * Pappalardi  _et al._ [2022] S. Pappalardi, L. Foini, and J. Kurchan, Eigenstate Thermalization Hypothesis and Free Probability, Phys. Rev. Lett. 129, 170603 (2022).
  * Pappalardi  _et al._ [2023] S. Pappalardi, F. Fritzsch, and T. Prosen, General Eigenstate Thermalization via Free Cumulants in Quantum Lattice Systems, arXiv:2303.00713 (2023).
  * Shi  _et al._ [2023] Z. D. Shi, S. Vardhan, and H. Liu, Local dynamics and the structure of chaotic eigenstates, arXiv:2306.08032 (2023).
  * Liao and Galitski [2022] Y. Liao and V. Galitski, Field theory approach to eigenstate thermalization in random quantum circuits, arXiv:2210.06480 (2022).
  * Prosen and Pižorn [2007] T. Prosen and I. Pižorn, Operator space entanglement entropy in a transverse Ising chain, Physical Review A 76, 032316 (2007).
  * Nahum  _et al._ [2017] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Quantum entanglement growth under random unitary dynamics, Phys. Rev. X 7, 031016 (2017).
  * Jonay  _et al._ [2018] C. Jonay, D. A. Huse, and A. Nahum, Coarse-grained dynamics of operator and state entanglement, arXiv:1803.00089 (2018).
  * Chan  _et al._ [2018] A. Chan, A. De Luca, and J. T. Chalker, Solution of a minimal model for many-body quantum chaos, Phys. Rev. X 8, 041019 (2018).
  * Garratt and Chalker [2021] S. J. Garratt and J. T. Chalker, Local pairing of feynman histories in many-body floquet models, Phys. Rev. X 11, 021051 (2021).
  * Brouwer and Beenakker [1996] P. W. Brouwer and C. W. J. Beenakker, Diagrammatic method of integration over the unitary group, with applications to quantum transport in mesoscopic systems, Journal of Mathematical Physics 37, 4904 (1996).
  * Chan  _et al._ [2019b] A. Chan, A. De Luca, and J. T. Chalker, Eigenstate correlations, thermalization, and the butterfly effect, Phys. Rev. Lett. 122, 220601 (2019b).
  * Nahum  _et al._ [2018] A. Nahum, J. Ruhman, and D. A. Huse, Dynamics of entanglement and transport in one-dimensional systems with quenched randomness, Phys. Rev. B 98, 035118 (2018).