††thanks: T.Y. and S.J.G. contributed equally to this work.††thanks: T.Y. and
S.J.G. contributed equally to this work.

# Operator dynamics in Floquet many-body systems

Takato Yoshimura  All Souls College, Oxford OX1 4AL, UK  Rudolf Peierls Centre
for Theoretical Physics, University of Oxford, 1 Keble Road, Oxford OX1 3NP,
UK  Samuel J. Garratt  Department of Physics, University of California,
Berkeley, California 94720, USA  J. T. Chalker  Rudolf Peierls Centre for
Theoretical Physics, University of Oxford, 1 Keble Road, Oxford OX1 3NP, UK

(August 28, 2024)

###### Abstract

We study operator dynamics in many-body quantum systems, focusing on generic
features of systems which are ergodic, spatially extended, and lack conserved
densities, as exemplified by spin chains with Floquet time evolution. To
characterise dynamics we examine, in solvable models and numerically, the
behaviour of operator autocorrelation functions, as a function of time and the
size of the operator support. The standard expectation is that operator
autocorrelation functions in such systems are maximum at time zero and decay,
over a few Floquet periods, to a fluctuating value that reduces to zero under
an average over an ensemble of statistically similar systems. Our central
result is that ensemble-averaged correlation functions also display a second
generic feature, which consists of a peak at a later time. In individual many-
body systems, this peak can also be revealed by averaging autocorrelation
functions over complete sets of operators supported within a finite spatial
region, thereby generating a partial spectral form factor. The duration of the
peak grows indefinitely with the size of the operator support, and its
amplitude shrinks, but both are essentially independent of system size
provided this is sufficiently large to contain the operator. In finite
systems, the averaged correlation functions also show a further feature at
still later times, which is a counterpart to the so-called ramp and plateau of
the spectral form factor; its amplitude in the autocorrelation function
decreases to zero with increasing system size. Both the later-time peak and
the ramp-and-plateau feature are specific to models with time-translation
symmetry, such as Floquet systems or models with a time-independent
Hamiltonian, and are absent in models with an evolution operator that is a
random function of time, such as the extensively-studied random unitary
circuits.

##  I Introduction

The essential phenomenology of dynamics in chaotic many-body quantum systems
describes the evolution of operators and of state vectors. Operators that are
initially simple evolve into increasingly complicated operators, and states
that initially have low entanglement develop increasing entanglement under
time evolution. These long-standing ideas [1, 2] have recently been
illustrated by calculations using random unitary circuits (RUCs) [3, 4, 5, 6],
in which the dynamics of a spin chain consists of a series of discrete steps
and at each step pairs of spins are coupled by randomly chosen gates.

In comparison with evolution under a time-independent Hamiltonian, RUCs embody
two related but distinct simplifications. Energy is eliminated as a conserved
quantity, and time-translation symmetry is broken by taking the gates at every
step to be statistically independent. The consequences of each of these
simplifications can be investigated separately, by constructing RUCs that have
an internal symmetry and an associated conserved density but no time-
translation symmetry [7, 8], and by studying random Floquet circuits (RFCs)
that have discrete time-translation symmetry but no conserved densities [9,
10]. Both alternatives show dynamical features that are absent in the simplest
RUCs. In the first case, the existence of a density that spreads diffusively
is imprinted on the correlation functions of observables to which it couples.
In the second case, time-translation symmetry gives rise to a distinctive ramp
and plateau in the spectral form factor (SFF) [11], obtained from the trace of
the evolution operator.

A natural question, and the subject of this paper, is what consequences time-
translation symmetry has for operator dynamics in chaotic systems. Strikingly,
we find that autocorrelation functions of observables — the simplest
quantities that characterise such dynamics — show a generic feature that we
believe has not previously been identified and that arises from this symmetry.
In contrast, operator spreading, as characterised by the out-of-time order
correlator (OTOC) [2], has exactly the same behaviour in the RFCs we study as
in RUCs examined previously. While the calculations presented in this paper
are restricted to Floquet models, we expect similar behaviour in Hamiltonian
systems for observables that do not couple to conserved densities.

Figure 1: Schematic behaviour of the average autocorrelation function
$\overline{C_{\alpha\alpha}(t)}$ ($P_{\alpha}\neq I$) for spatially local
Floquet many-body systems without conservation laws (taking $a<L/2$).

In order to summarise our results more precisely, we introduce some notation.
Consider a spin chain of $L$ sites with local Hilbert space dimension $q$ and
hence full Hilbert space dimension $D=q^{L}$. The ($D\times D$) unitary
evolution operator is defined for integer times $t$ as $W(t)\equiv W^{t}$,
with $W$ the Floquet operator. Now take a complete orthonormal basis of
$q^{2}$ Hermitian operators at each site, with the identity as one member of
the basis, and construct from them as direct products the $D^{2}$ observables
$P_{\alpha}$ satisfying $D^{-1}{\rm
Tr}[P_{\alpha}P_{\beta}]=\delta_{\alpha\beta}$. The correlation function of
the observables $P_{\alpha}$ and $P_{\beta}$ is

$C_{\alpha\beta}(t)=\langle P_{\alpha}(t)P_{\beta}\rangle,\quad
P_{\alpha}(t)=W^{\dagger}(t)P_{\alpha}W(t)$ (1)

with $\langle\,\cdots\rangle\equiv D^{-1}{\rm Tr}[\,\cdots]$, the infinite-
temperature average appropriate for a Floquet system. Autocorrelation
functions correspond to the case $\alpha=\beta$. We denote the average of the
correlation function over an ensemble of similar systems by
$\overline{C_{\alpha\beta}(t)}$.

It is useful to have a simple picture of the dynamics in mind. An
autocorrelation function takes its maximum value at time zero and is expected
to decay on a short timescale, typically set by the period in Floquet systems
and by the inverse of coupling energies in Hamiltonian models. After this
initial decay, it may exhibit irregular fluctuations, for example because of
spin precession in evolving local exchange fields. One expects such
fluctuations to be a specific fingerprint of the particular observable and
system. Averaging over an ensemble of similar systems, a natural possibility
is that fluctuations would be reduced to zero. We find this is not the case:
instead, a second peak in the autocorrelation function emerges from the
average. The duration of this peak grows indefinitely with the size of the
support of the operator concerned, and it is therefore a parametrically
distinct feature for ‘mesoscopic’ operators, which have large but finite
support.

This structure of the ensemble-averaged autocorrelation function
$\overline{C_{\alpha\alpha}(t)}$ is shown schematically in Fig. 1. It exhibits
three generic features as a function of time. The most prominent one is a peak
at $t=0$ with amplitude $\overline{C_{\alpha\alpha}(0)}=1$ fixed by our choice
of normalisation for the observable. In a system without locally conserved
densities (or for observables that do not couple to conserved densities) the
autocorrelation function relaxes from this initial value over timescale
$t_{\rm d}$, which we refer to as the decay time. Of the three features, we
believe that until recently only this first one has been recognised. By
contrast, our concern in this paper is with the two additional generic
features in the autocorrelation function which appear at longer timescales.

The second generic feature is an asymmetric peak which extends from a time of
order $t_{\rm d}$ up to a timescale that we call the Thouless time
$t_{\text{Th},\alpha}$ of the operator $P_{\alpha}$ (we comment below on this
terminology). For an operator having nontrivial support on $a$ sites, we will
see below that a natural scale for $\overline{C_{\alpha\alpha}(t)}$ is
$q^{-2a}$, and we define $t_{\text{Th},\alpha}$ as the time at which the tail
drops to a value of order $q^{-2a}$. For example, if the sites on which
$P_{\alpha}$ has non-trivial support are contiguous, we find that
$t_{\text{Th},\alpha}\sim\ln a$. Roughly, $t_{\text{Th},\alpha}$ can be
identified as the time beyond which the dynamics of $P_{\alpha}$ resembles
that with an unstructured random unitary matrix as Floquet operator.

At the latest times we find the third generic feature: a linear growth and
eventual plateau of the autocorrelation function, with an onset time $t_{\rm
ramp}\propto L-L_{A}$. The amplitude of this growth and the plateau value are
exponentially small in system size, and this feature is therefore completely
suppressed in the thermodynamic limit. This late-time feature in
autocorrelation functions is the counterpart of the ramp and plateau, familiar
in the behaviour of the SFF as a representation of random matrix correlations
between eigenvalues of the Floquet operator or Hamiltonian for a system [11].
Its behaviour has recently been analysed [12] in the context of the partial
spectral form factor (PSFF).

We arrive at the picture outlined in the preceding paragraphs using a
combination of analytical and numerical calculations in RFCs with short-range
interactions. Our analytical results are obtained for the random phase model
(RPM) introduced in Ref. [13], which is exactly solvable for $q\to\infty$,
while for numerical calculations we use the circuit studied in Ref. [14].

Both the asymmetric peak in autocorrelation functions, and the ramp-and-
plateau structure, are features visible only after a suitable average. One
possibility is to average autocorrelation functions over an ensemble of random
systems. An alternative is to average by summing the autocorrelation functions
over choices of the initial operator. This generates the PSFF [12] and for
sufficiently large subregions it reveals the asymmetric peak. The ramp and
plateau, on the other hand, are obscured by large temporal fluctuations in
individual systems, and are not visible without an explicit average over
different systems or over time. The absence in this sense of self-averaging
has long been recognised in the case of the SFF [15, 16].

Our analysis of dynamics is rooted in the deconstruction of the SFF as a sum
over autocorrelation functions of a complete set of observables [17, 18].
Under this correspondence the asymmetric peak in autocorrelation functions at
intermediate times is the source of deviations of the SFF from random matrix
behaviour, identified recently in ergodic, spatially extended many-body
systems [13, 19, 20]. These deviations are absent in dual-unitary circuits
[10, 21, 22], where many autocorrelation functions also vanish. We note that
the behaviour of OTOCs [23] has been studied in dual unitary circuits, and
more recently in systems where this restriction is weakly violated [24].

Our use of the term Thouless time deserves discussion, since it and the term
Thouless energy for its reciprocal have been employed in a variety of
contexts. As originally introduced in discussions of single-particle models
for disordered conductors [25], the terms have both a dynamical significance
(as the time for a particle to diffuse over a distance set by the system size)
and a spectral significance (as the energy scale below which eigenvalue
correlations match the predictions of random matrix theory [26]). In many-body
systems the terms have been used in both in a spectral sense [27, 13, 20, 28]
and a dynamical sense [29, 30] but the link between the two has been unclear.
A key result of the current work is to elucidate the relation between these
scales. We will identify a Thouless time in three different settings: the
dynamics of individual operators; the behaviour of the PSFF; and the behaviour
of the SFF. These are linked by the fact that the PSFF and SFF can be obtained
as sums over autocorrelation functions of operators. In all three settings the
Thouless time has a dynamical significance, while in the third it also has a
spectral significance.

The effects of time-translation symmetry can be understood in a unified way by
interpreting an autocorrelation function as a sum of amplitudes associated
with paths in the space of operators. In general, these amplitudes are real
but not necessarily positive. We identify a subset of paths with positive
amplitudes that make the dominant contribution to the averaged autocorrelation
functions. Interestingly, these paths have rather narrow trajectories in the
full space of operators. For example, in the RPM at large $q$, the boundaries
of operator paths form straight vertical lines in the space-time plane, as
illustrated in Fig. 2. This picture is to be contrasted with that for the
averages of squared autocorrelation functions, and also of OTOCs (both of
which are nonzero even in RUCs). In the first case the spatial widths of the
dominant operator paths grow as $t^{1/2}$ [31], whereas in the second they
grow ballistically [4, 5]. In systems with conservation laws, we note that the
operator paths relevant to dynamics have recently been studied in Refs. [32,
33].

(a)  (b)

Figure 2: Operator paths contributing to ensemble-averaged autocorrelation
functions $\overline{C_{\alpha\alpha}(t)}$. Coloured regions represent regions
of space where the operator acts non-trivially (i.e. not as the identity). (a)
Examples of the dominant paths for various choices of
$\overline{C_{\alpha\alpha}(t)}$ distinguished by colour. (b) Examples of
paths that are sub-dominant for the same set of
$\overline{C_{\alpha\alpha}(t)}$.

The asymmetric peak is striking even for operators supported on a few sites.
This feature is therefore straightforward to observe in experiments involving
only conventional dynamics and a number of samples that is independent of $L$.
We should compare this situation with that for OTOCs and the SFF. For example,
one method for determining OTOCs involves an explicit time-reversal of many-
body dynamics [34], while randomised measurement approaches for extracting
either the OTOC [35] or the SFF [12] suffer from a sample complexities that
scale exponentially with the system size.

This paper is organised as follows. In Sec. II we provide an overview of our
results, and in Sec. III we present calculations of averaged autocorrelation
functions and PSFFs. Following this, in Sec. IV, we calculate the OTOC in the
RPM. In Sec. V we provide an interpretation of our results in the language of
operator paths. Finally, in Sec. VI, we discuss statistical fluctuations
relevant to dynamics. We summarise and discuss our results in Sec. VII.

##  II Overview

This section is organised as follows. In Sec. II.1 we introduce our notation
and define the basic physical quantities of interest. Following this, in Sec.
II.2, we describe the models we use for analytical as well as numerical
calculations. In Secs. II.3 and II.4 we then outline the behaviour of
autocorrelation functions and PSFFs, respectively, before discussing OTOCs in
Sec. II.5.

###  II.1 Definitions

In this work we study chaotic many-body dynamics in interacting spin chains.
Our focus will be on the effects of locality in systems with fixed evolution
operators, where there is time-translation symmetry in the dynamics. Our
minimal models for this kind of dynamics are RFCs without conservation laws.
As an orthonormal basis for the Hilbert space of many-body quantum states we
choose the tensor products $|a\rangle=|a_{1}\ldots a_{L}\rangle$, where
$a_{x}=0,\cdots,q-1$, and $a=0,\cdots,D-1$. Since our focus is primarily on
the dynamics of operators, it is convenient to choose a Hermitian operator
basis $P_{\alpha}$; here the label $\alpha=0,\ldots,(q^{2L}-1)$. Explicitly,
$P_{\alpha}=P_{\alpha_{1}}\ldots P_{\alpha_{L}}$ where $P_{\alpha_{x}}$, with
$\alpha_{x}=0,\ldots,(q^{2}-1)$, acts on site $x$. For example, for $q=2$, we
can choose $P_{\alpha_{x}}$ to be Pauli matrices. We reserve the label
$\alpha_{x}=0$ to denote the $q\times q$ identity at site $x$, and $\alpha=0$
to denote the $D\times D$ identity. Note that in our convention
$q^{-1}\text{Tr}[P_{\alpha_{x}}P_{\beta_{x}}]=\delta_{\alpha_{x}\beta_{x}}$,
where the trace is over states at $x$, and
$D^{-1}\text{Tr}[P_{\alpha}P_{\beta}]=\delta_{\alpha\beta}$, where the trace
is over the full Hilbert space of states. This condition implies that all
$P_{\alpha\neq 0}$ are traceless.

We will make multiple uses of the completeness relation for our operator
basis. In terms of the matrix elements $\langle a|P_{\alpha}|b\rangle$ it has
the form

$\displaystyle D^{-1}\sum_{\alpha=0}^{D^{2}-1}\langle
a|P_{\alpha}|b\rangle\langle c|P_{\alpha}|d\rangle=\delta_{ad}\delta_{bc}.$
(2)

Our probes of the dynamics are (infinite-temperature) correlation functions;
see Eq. (1). Note that $C_{\alpha\beta}(t)$ is real number, and that the
normalisation condition $D^{-1}\text{Tr}P_{\alpha}^{2}(t)=1$ implies
$\sum_{\beta\in\mathcal{P}}C^{2}_{\alpha\beta}(t)=1$ for $\alpha\not=0$, where
$\mathcal{P}$ is the set of non-identity $P_{\alpha}$. This means that
$C_{\alpha\beta}(t)$ is an orthogonal matrix. Using Eq. (2) the
autocorrelation function $C_{\alpha\beta}(t)$ can be recast as

$\displaystyle C_{\alpha\beta}(t)$
$\displaystyle=\sum_{\alpha(1),\cdots,\alpha(t-1)\in\mathcal{P}}C_{\alpha\alpha(1)}\cdots
C_{\alpha(t-1)\beta},$ (3)

where in the correlation function for a single time step we omit the time
argument, $C_{\alpha\beta}\equiv C_{\alpha\beta}(1)$. The expansion given in
Eq. (3) can be thought of as a sum of the amplitudes of paths
$(\alpha,\alpha(1),\ldots,\beta)$. The normalisation condition
$D^{-1}\text{Tr}P_{\alpha}^{2}(t)=1$ can now be rewritten as

$\displaystyle\sum_{\alpha(1)\cdots\alpha(t)\in\mathcal{P}}C^{2}_{\alpha\alpha(1)}\ldots
C^{2}_{\alpha(t-1)\alpha(t)}=1,$ (4)

and so we can interpret the squared amplitudes $C^{2}_{\alpha\alpha(1)}\ldots
C^{2}_{\alpha(t-1)\alpha(t)}$ as probabilities for the various operator paths.

Our focus in Sec. III will be on autocorrelation functions
$C_{\alpha\alpha}(t)$, the PSFF

$\displaystyle K_{A}(t)=1+\sum_{\alpha\in\mathcal{P}_{A}}C_{\alpha\alpha}(t),$
(5)

where $\mathcal{P}_{A}$ is the set of $P_{\alpha}$ with non-trivial support in
the region $A$ (and trivial support in its complement $\overline{A}$), and the
SFF

$\displaystyle K(t)=1+\sum_{\alpha\in\mathcal{P}}C_{\alpha\alpha}(t).$ (6)

These quantities have alternative expressions in terms of the Floquet
operator, as given in Eq. (28) in the case of the PSFF, and as $K(t)=|{\rm
Tr}W(t)|^{2}$ in the case of the SFF. The link between these expressions and
Eqns. (5) and (6) is provided by the operator resolution of the identity, Eq.
(2).

In Sec. IV we study the OTOC,

$\displaystyle\mathcal{C}(x,t)$
$\displaystyle=-\frac{1}{2}\langle[\mathcal{O}(0,t),\mathcal{O}(x,0)]^{2}\rangle$
(7) $\displaystyle=1-\sum_{\alpha,\beta\in\mathcal{P}}\langle
P_{\alpha}\mathcal{O}(x,0)P_{\beta}\mathcal{O}(x,0)\rangle
C_{\mathcal{O}\alpha}(t)C_{\mathcal{O}\beta}(t),$

where $\mathcal{O}(x,0)$ is a local operator supported at site $x$, and
$C_{\mathcal{O}\alpha}(t):=\langle\mathcal{O}(0,t)P_{\alpha}\rangle$.

From the perspective highlighted in this work, the crucial difference between
correlation functions and OTOCs is that the former involve sums over
amplitudes of paths whereas, in both RUCs and RFCs, the dominant contribution
to the latter is a sum over probabilities. Meanwhile, a difference between
RUCs and RFCs with Haar-random local unitary operations is that, in the
former, ensemble averages of operator path amplitudes necessarily vanish.
Their sums, from autocorrelations functions to the SFF, are therefore trivial.
This is not the case in RFCs.

###  II.2 Models

We perform analytical calculations for the RPM [13], a minimal model for
Floquet dynamics without conservation laws which is tractable in the limit of
large local Hilbert space dimension $q$. It has the Floquet operator
$W=W_{2}W_{1}$, with

$\displaystyle W_{1}=\bigotimes_{x=1}^{L}U_{x}$ (8)

a tensor product of single-site ($q\times q$) Haar-random unitary operators,
and $W_{2}$ a product of diagonal two-site unitary operators which couple
neighbouring sites,

$\displaystyle W_{2}|a\rangle=\exp\left({\rm
i}\sum_{x}\varphi_{x;a_{x},a_{x+1}}\right)|a\rangle.$ (9)

Each $\varphi_{x;a_{x},a_{x+1}}$ is an independent Gaussian random variable
with mean zero and variance $\varepsilon$. The structure of the Floquet
operator is illustrated in Fig. 3. Previously it has been shown that, in this
model and in the limit of large $q$, the average SFF behaves as in RMT,
$\overline{K(t)}\simeq t$, beyond the spectral Thouless time
$t_{\text{Th}}\sim\log L/\varepsilon$. Due to the local Haar-random
operations, the relaxation time $t_{\mathrm{d}}=1$, but this is non-generic;
we explore a generalisation of the RPM which is not locally Haar random in
Appendix A.

Figure 3: Floquet operator in the RPM. The blue squares represent the local
unitary operations $U_{x}$ [Eq. (8)], and the red ovals represent the
interactions in Eq. (9).

For numerical calculation we do not use the RPM because, for $q=2$, it is
known [13] to be in or close to a many-body localised regime for all
$\varepsilon$, while setting $q>2$ limits the system sizes that are
accessible. Instead, we study the kicked random field Heisenberg model of
Refs. [14] in its ergodic phase. This model involves single-site Haar-random
unitary operations as above, but the coupling between sites is described by an
exponentiated SWAP operator, which for $q=2$ is equivalent to a Heisenberg
interaction. We will refer to this as the random SWAP model, or RSM.

The Floquet operator $W$ for this model can be expressed as $W=W_{2}W_{1}$
where $W_{1}=\bigotimes_{x\text{ odd}}V_{x,x+1}$ and
$W_{2}=\bigotimes_{x\text{ even}}V_{x,x+1}$; with periodic boundary
conditions, which necessitates $L$ even, there are $L/2$ gates in each layer,
and we identify $V_{L,L+1}=V_{L,1}$. With open boundary conditions and $L$
even there are $L/2$ gates in $W_{1}$ and $L/2-1$ in $W_{2}$, while with $L$
odd there are $(L-1)/2$ gates in each layer. The operators $U_{x,x+1}$ take
the form

$\displaystyle V_{x,x+1}=\exp\big{(}iJ\,\text{SWAP}\big{)}U_{x}\otimes
U^{\prime}_{x+1},$ (10)

where SWAP is the two-site swap operator, and each of $U_{x}$ and
$U^{\prime}_{x+1}$ are independent $2\times 2$ Haar-random unitary matrices,
which can be viewed as describing the precession of qubits in random fields.
For all numerical calculations we use the $q=2$ RSM with $J=0.2$. The boundary
conditions are periodic except in Fig. 9, where they are open.

###  II.3 Autocorrelation functions

Here we outline our results on the behaviour of autocorrelation functions.
First we sketch the behaviour in systems whose Floquet operators are $D\times
D$ Haar-random unitary matrices — we will refer to these as ‘random matrix
Floquet’ (RMF) systems — and then we offer a qualitative picture of the
behaviour we should expect in a system with local interactions. Following this
we describe our results in the RPM at large $q$.

####  II.3.1 Qualitative approach

In RUC systems with either global or local Haar-random unitary operations, and
without time-translation invariance, ensemble averages of autocorrelation
functions vanish. In a RMF system, which by construction has no spatial
structure, for any operator $P_{\alpha}\in\mathcal{P}$ the autocorrelation
function behaves as

$\overline{C_{\alpha\alpha}(t)}=\frac{\overline{K(t)}-1}{D^{2}-1},$ (11)

while the SFF $\overline{K(t)}$ is given by

$\overline{K(t)}=\begin{dcases*}D^{2}&$t=0$\\\ t&$1\leq t<D$\\\ D&$D\leq
t$.\end{dcases*}$ (12)

That is, the average autocorrelation function drops from unity to zero over
the first time step, and then increases linearly (at a rate that is
exponentially small in the number of degrees of freedom) before an eventual
plateau beyond $t=D$.

As a step towards understanding the behaviour in a system with spatial
structure, first consider a system of uncoupled sites. That is, take the
evolution operator for the $L$-spin system to be a tensor product of $q\times
q$ Haar-random unitary operators (note that this corresponds to the
$\varepsilon=0$ limit of the RPM). The evolution operator for each spin is
fixed in time, and so the dynamics corresponds simply to the precession of
spins in local fields. For single-site operators we can apply the above
results for RMF systems, replacing the global Hilbert space dimension $D$ with
the local one, giving

$\displaystyle\overline{C_{\alpha\alpha}(t)}=\begin{cases}1&t=0\\\
(t-1)/(q^{2}-1)&1\leq t\leq q\\\ (q-1)/(q^{2}-1)&q\leq t.\end{cases}$ (13)

For $q=2$ the result is particularly simple:
$\overline{C_{\alpha\alpha}(0)}=1$, $\overline{C_{\alpha\alpha}(1)}=0$ and
$\overline{C_{\alpha\alpha}(t\geq 2)}=1/3$.

Now we can ask how this behaviour is modified by intersite coupling that is
weak but chosen sufficiently large that the coupled system is in the ergodic
rather than the many-body localised phase (a condition that is easy to satisfy
if $q\gg 1$ [13]). The behaviour should be approximately described by Eq. (13)
at early times and by the RMF result in Eq. (11) at late times. Physically, we
can understand this interpolation as corresponding to the dephasing of spins’
precession by interactions with their neighbors. For this reason, we should
expect that averaged autocorrelation functions first increase with time as in
Eq. (13), but that the late-time plateau which occurs there for $t\geq q$ is
replaced by gradual decay of the autocorrelation function to a value that is
exponentially small in $L$, as in Eq. (11).

As we will demonstrate, both in the RPM at large $q$ and numerically at $q=2$,
the ideas above provide an qualitative characterisation of the dynamics of
local observables. However, it is not immediately clear how they should be
extended to describe the behaviour of observables with support on multiple
sites, and how the autocorrelation function will depend on the spatial
structure of the observable. This can be understood through a study of the
RPM.

####  II.3.2 RPM at large $q$

In the RPM we can calculate ensemble-averaged autocorrelation functions
analytically in the limit of large $q$. Our basis operators $P_{\alpha}$ are
tensor products of identity and non-identity operators; it is useful to view
each $P_{\alpha}$ as consisting of ‘clusters’ of non-identity operators which
are separated by identities. We denote the number of clusters in $P_{\alpha}$
by $n$, use $m=1,\ldots,n$ to label the clusters, and let $a_{m}$ be the
length of the cluster $m$. The total number of sites on which the operator
$P_{\alpha}$ differs from the identity is then $a=\sum_{m=1}^{n}a_{m}$; we
call this the weight of the operator. As we will show in Sec. III, in the
limit of large $q$ the ensemble-averaged autocorrelation functions depend only
on $a$ and $n$ and are given by

$\overline{C_{\alpha\alpha}(t)}=q^{-2a}e^{-2n\varepsilon
t}(t-1)^{n}\left[1+(t-2)e^{-\varepsilon t}\right]^{a-n}.$ (14)

Clearly $\overline{C_{\alpha\alpha}(t)}=0$ at the decay time
$t_{\mathrm{d}}=1$, but from Eq. (14) we see that this is followed by a peak.
In the simplest case of $n=1$, i.e. an operator $P_{\alpha}$ with support on
$a$ contiguous sites, the late-time tail of this peak behaves as

$\displaystyle\overline{C_{\alpha\alpha}(t)}=$
$\displaystyle\,q^{-2a}e^{-2\varepsilon t}(t-1)$ (15)
$\displaystyle\times[1+e^{-\varepsilon t}(a-1)(t-2)+\ldots].$

From this expression we can identify the time scale at which
$\overline{C_{\alpha\alpha}(t)}$ falls to a value of order $q^{-2a}$ as
$t_{\text{Th},\alpha}$, where

$t_{\text{Th},\alpha}=\frac{\ln a}{\varepsilon},$ (16)

and we have assumed that $a\gg n$.

More generally, the late-time behaviour of $\overline{C_{\alpha\alpha}(t)}$
can be understood by viewing each cluster (of size $a_{m}$) as having random
matrix dynamics. Approximating the evolution of a cluster of size $a_{m}$ by
its evolution in a $q^{a_{m}}$-dimensional RMF system, its contribution to
autocorrelation functions is $q^{-2a_{m}}(t-1)e^{-2\varepsilon t}$ in the
large-$q$ limit, where factors $e^{-\varepsilon t}$ at each boundary capture
dephasing induced by coupling to the remainder of the system. Combining the
contributions from each cluster multiplicatively, we recover the late time
asymptotics of autocorrelation functions in the RPM.

Note that the value $t_{\mathrm{d}}=1$ for the relaxation time is a
consequence of Haar-random local unitary rotations included in the Floquet
operator, and is not generic. In Appendix A we modify the RPM so that the
single-site operations are not Haar random, giving a value for
$t_{\mathrm{d}}$ that varies with model parameters.

We can gain more intuition about the late-time behaviour of autocorrelation
functions at $q\to\infty$ by interpreting them as sums over amplitudes of
operator paths using Eq. (3). In Sec. V we will show that in the RPM at large
$q$ the leading-order contribution to ensemble-averaged autocorrelation
functions comes from operator paths that, despite the ergodicity of the model,
are localised in the sense that all operators in the path have clusters of the
same length and in the same spatial locations.

###  II.4 The PSFF in the RPM

Summing over autocorrelation functions of all operators supported within a
region $A$, we arrive at the PSFF via Eq. (5). We will restrict ourselves to a
contiguous regions $A$ of length $L_{A}$, and we are primarily interested in
the regime $L_{A}\ll L$. In the limit of large $q$, the average PSFF is

$\overline{K_{A}(t)}=\frac{1}{t}\left[\lambda_{0}^{L_{A}+1}+(t-1)\lambda^{L_{A}+1}\right]$
(17)

where $\lambda_{0}=1+(t-1)e^{-\varepsilon t}$ and $\lambda=1-e^{-\varepsilon
t}$. Like autocorrelation functions, the PSFF exhibits an asymmetric peak
after the decay time $t_{\rm d}$, and the height of the peak grows
exponentially with $L_{A}$. For $L_{A}\ll L$ we then have $K_{A}(t)\simeq 1$
beyond a time scale which we denote by $t_{\text{Th},A}$. The Thouless time
$t_{\text{Th},A}$ for the region $A$ is given by

$t_{\text{Th},A}=\frac{\ln L_{A}}{\varepsilon}.$ (18)

Although this time scale coincides with $t_{\text{Th},\alpha}$ for operators
having $n=1$ and $a=L_{A}$, the late-time behaviour of the PSFF has
contributions from all $n=1$ operators with $1\leq a\leq L_{A}$. For example
the contribution from single-site operators is $L_{A}e^{-2\varepsilon t}$,
where the factor of $L_{A}$ comes from summing over all possible locations of
the operator within $A$. In the case $L_{A}=L$, i.e. $K_{A}(t)=K(t)$, we
recover $t_{\text{Th},A}=t_{\text{Th}}=\varepsilon^{-1}\ln L$ from Ref. [13].

The behaviour above should be contrasted with that in RMF systems without
spatial structure. In that setting one finds a shift-ramp-plateau structure
[12]

$\overline{K_{A}(t)}=\frac{D_{A}^{2}-1}{D^{2}-1}\overline{K(t)}+\frac{D^{2}-D_{A}^{2}}{D^{2}-1},$
(19)

where we recall that the SFF behaves as in Eq. (12). Note that Eq. (19) can
also be obtained by summing over averaged autocorrelation functions of the
form in Eq. (11). So far we have provided only analytic results, but in Secs.
III.4 and V.3, we show numerically that the $q=2$ RSM has similar dynamics to
the RPM at large $q$.

###  II.5 The OTOC in the RPM

The notion of operator spreading provides a different perspective on dynamics,
and here the central quantity is the OTOC in Eq. (7). The basic phenomenology,
supported by calculations in RUCs [4, 5], is that operators grow
ballistically, with the left- and right-hand operator ‘fronts’ broadening
diffusively in time, and so

$\overline{\mathcal{C}(x,t)}\simeq\Phi\left(\frac{v_{B}t-|x|}{\sqrt{2\mathcal{D}t}}\right).$
(20)

Here $\Phi(u)$ is the error function, with $\Phi(u\to-\infty)\to 0$ and
$\Phi(u\to\infty)\to 1$, and the quantity $v_{B}$ is known as the butterfly
velocity. In RUCs with Haar-random two-site gates, $v_{B}$ and $\mathcal{D}$
are functions of $q$. A pathology of the large-$q$ limit in these models [4,
5] (and in the RFC with Haar-distributed two-site gates [9]) is that $v_{B}$
approaches the lightcone velocity $v_{L}$ [36] from below, and the constant
$\mathcal{D}$ goes to zero. One also finds $v_{B}=v_{L}$ in dual unitary
circuits [23], although the generic behaviour $v_{B}<v_{L}$ has recently been
recovered in a perturbative expansion around the dual unitary limit [24].

An attractive feature of the RPM is that diffusive broadening of the operator
front survives even at $q\to\infty$, and $v_{B}<v_{L}$ (see also Ref. [37]).
To be more precise, again on a coarse-grained space-time scale, the OTOC in
the RPM at large $q$ has exactly the same functional form as in Eq. (20) where
the butterfly velocity $v_{\mathrm{B}}$ and the diffusion constant
$\mathcal{D}$ are in this case

$v_{\mathrm{B}}=1-e^{-2\epsilon},\quad\mathcal{D}=\frac{1}{2}e^{-2\epsilon}(1-e^{-2\epsilon}).$
(21)

Note that the structure above shows that the operator paths contributing to
OTOCs are qualitatively different to those contributing to autocorrelation
functions. In particular, while operator paths contributing to autocorrelation
functions do not change their spatial structure in time, those contributing to
OTOCs grow ballistically; we elaborate on this difference in Sec. V.

##  III Autocorrelation functions

Our main source of insight into the behaviour of autocorrelation functions
comes from the RPM. In this model the effects of locality on the SFF are
evident even for $q\to\infty$ [13], and through Eq. (6) we can infer that the
phenomena encountered in that setting should also manifest themselves in
dynamics. In this section we calculate autocorrelation functions, finding
features that are a consequence of locality and are most dramatic for local
operators. First, in Sec. III.1, we summarise the diagrammatic techniques we
use for evaluating Haar averages in autocorrelation functions and the PSFF. We
then calculate averaged autocorrelation functions in the RPM at large $q$ in
Sec. III.2, generalising these results to a model with a parameter-dependent
relaxation time in Appendix A. As noted above, PSFFs can be expressed as
partial sums over autocorrelation functions, and we calculate the averages of
these sums in Sec. III.3. Guided by these calculations, in Sec. III.4 we turn
to numerics on qubits chains ($q=2$), and show that the behaviour of averaged
autocorrelation functions in these chains is very similar to that at large $q$
in the RPM. In Sec. III.5, we discuss an approximate description of the
behaviour of the PSFF of the RPM at finite $q$, including a treatment of
systems of finite size. Finally, in Sec. III.6, we connect our discussion of
autocorrelation functions to the behaviour of the SFF, and in this way provide
an intuitive perspective on the deviations of spectral statistics from RMT
which are known to arise from locality.

###  III.1 Diagrammatics and Haar averages

In the RPM, the single-site unitary operators $U_{x}$ are independent for
different sites, so Haar averages can be performed independently for each
site. To evaluate single-site Haar averages, it is convenient to express the
averages diagrammatically. In this section, we outline the diagramatic scheme
of Brouwer and Beenakker [38]. In brief, an average over a product of matrix
elements of a unitary and its conjugate can be expressed as a sum over
‘pairings’ of their indices, and a subset of pairings contributes to physical
quantities at large $q$.

To illustrate the idea, consider first a simple example
$q^{-1}\overline{\text{Tr}[UUU^{\dagger}U^{\dagger}]}$, where $U$ is a
$q\times q$ Haar-distributed unitary matrix. This quantity is trivial because,
without taking a Haar average, the unitarity of $U$ implies that it is simply
unity. It is however instructive to reproduce this result by explicitly Haar
averaging the unitaries rather than resorting to the unitarity condition
$UU^{\dagger}=1$. The calculation highlights the importance of pairings that
are subleading and hence are often neglected at large $q$.

We represent matrix elements $U_{aa^{\prime}}$ of a unitary $U$ by a pair of
black and white dots, with the first index ($a$) black and the second
($a^{\prime}$) white, connected by a dotted line. $U^{*}_{bb^{\prime}}$ is
also represented in the same way, but to distinguish it from
$U_{aa^{\prime}}$, an asterisk is shown alongside the dotted line. With these
rules, Haar averaging induces pairings of indices used for unitaries $U$ and
their conjugates $U^{*}$, and there are four possible pairings as depicted in
Fig. 4. These pairings can be classified into two categories: Gaussian
pairings and non-Gaussian pairings. In Fig. 4, the first two pairings are
referred to as Gaussian: the first and second indices of a unitary $U$ are
paired with those of a unique $U^{*}$. This is not the case in the last two
pairings in Fig. 4, where at least one unitary $U$ is paired with two
different conjugates $U^{*}$. These pairings are thus called non-Gaussian
pairings.

Figure 4: The average of
$q^{-1}\overline{\text{Tr}[UUU^{\dagger}U^{\dagger}]}$ is a sum of two
Gaussian (first and second) and two non-Gaussian (third and fourth) pairings.

It is known that averages of unitaries, such as those in Fig. 4, can be
expressed in terms of the Weingarten function [38]. The pairings on the right
side of the equality in Fig. 4 have respective values

$\frac{q^{2}}{q^{2}-1},\quad\frac{1}{q^{2}-1},\quad-\frac{1}{q^{2}-1},\quad-\frac{1}{q^{2}-1}.$
(22)

The four pairings sum to one; although the first Gaussian pairing is the
leading contribution at large $q$, taking into account the others is necessary
for the result to be compatible with unitarity at finite $q$. Below we will
see that non-Gaussian pairings are important for the calculation of the OTOC
in Sec. IV even in the large-$q$ limit. In the current section, however, we
will only have to consider Gaussian pairings. For this reason, it will be
convenient to simplify our notation further, and below we will represent both
the first and second indices of a unitary together as a single black dot.

###  III.2 Autocorrelation functions in the RPM

Having the diagrammatics involved in Haar averaging in mind, we move to
evaluate autocorrelation functions $\overline{\langle
P_{\alpha}(t)P_{\alpha}\rangle}$ at large $q$. Since unitaries are acting on
single sites in the RPM, we can take an independent Haar average at each site
and study which pairings are relevant in the large-$q$ limit. As we alluded to
earlier, contributions from non-Gaussian diagrams are suppressed at large $q$.
Moreover, of the $t!$ Gaussian pairings, only the $t$ cyclic pairings will
contribute.

To start, suppose the operator $P_{\alpha}$ acts on the site $x$ nontrivially,
in the sense that $P_{\alpha_{x}}\not=I_{q}$ (the identity). Figs. 5a-5c then
show the leading pairings when $t=3$. Note that every unitary $U_{x}$ and
$U^{*}_{x}$ is now represented by a single black dot, and for brevity, we have
also removed the asterisk that was previously used for indicating conjugation.
If we label the unitaries $U_{x}$ and $U^{*}_{x}$ on the left and right edge
of the diagrams as
$[U_{x}]_{a(0)a(1)},[U_{x}]_{a(1)a(2)},\cdots,[U_{x}]_{a(t-1)a(t)}$ and
$[U^{*}_{x}]_{b(0)b(1)},[U^{*}_{x}]_{b(1)b(2)},\cdots,[U^{*}_{x}]_{b(t-1)b(t)}$
from bottom to top, these pairings can be conveniently parameterised by
$s_{x}=0,\cdots,t-1$, corresponding to indices paired cyclically as
$a(r)=b(r+s_{x})$.

(a)  (b)  (c)  (d)

Figure 5: Leading local pairings in autocorrelation functions at $t=3$. When
$P_{\alpha}$ acts on the site nontrivially, we have three pairings: (a) $s=0$,
(b) $s=1$, and (c) $s=2$. Note that the first pairing has a vanishing
contribution as mentioned in the main text. (d) The only dominant local
pairing on the site where the operator acts on trivially, $s=0$. Black dots
represent on-site unitaries $U_{x}$ and $U^{*}_{x}$.

Importantly, when $P_{\alpha_{x}}\not=I_{q}$, the contribution from the
$s_{x}=0$ pairing vanishes because our basis operators are traceless.
Therefore $s_{x}\neq 0$ pairings become the leading pairings at site $x$. With
such pairings, the contribution from the choice of operator is fixed by our
normalisation condition $q^{-1}{\rm Tr}[P_{\alpha_{x}}^{2}]=1$. On the other
hand, if $P_{\alpha_{x}}=I_{q}$, the opposite happens: $s_{x}=0$ is the only
pairing (with the weight $q$) compatible with the boundary condition, and as
such $s_{x}\neq 0$ pairings are all subleading at large $q$. As a result,
averaged autocorrelation functions of basis operators can be expressed as sums
over ‘pairing configurations’ $s_{1},\ldots,s_{L}$, with $s_{x}=0$ at any site
$x$ where the operator of interest acts as the identity, and
$s_{x}=1,\ldots,(t-1)$ otherwise.

In other words, the averaged autocorrelation function then has the form of a
partition function for a one-dimensional system of $L$ sites, whose on-site
degrees of freedom are pairings $s_{x}$. Averaging over the phases $\varphi$
on each bond [see Eq. (9)] we find a factor of $e^{-\varepsilon t}$ wherever
$s_{x}\neq s_{x+1}$, and a factor of unity otherwise. To connect with the
calculation of the SFF in Ref. [13] we define a $t\times t$ transfer matrix
which acts on a space in which the $t$ cyclic pairings correspond to
orthonormal basis vectors. The matrix elements of $T$ in this basis are

$\displaystyle
T_{ss^{\prime}}=\delta_{ss^{\prime}}+(1-\delta_{ss^{\prime}})e^{-\varepsilon
t},$ (23)

where for brevity we have dropped the site label $x$ from the pairings
$s,s^{\prime}$. This object is related to the SFF via
$\overline{K(t)}=\text{Tr}T^{L}$ in a system with periodic boundary conditions
[13]. A key concept which emerges from the structure of the transfer matrix is
the notion of a domain wall, which is simply a configuration of unequal
pairings on two adjacent sites. Domain walls do not exist in RMF systems, and
they provide a way to understand the effects of locality in dynamics and on
many-body spectra [13, 20].

To calculate autocorrelation functions, in addition to the transfer matrix we
require the on-site matrices

$\displaystyle D_{0}=|0\rangle\langle 0|,\quad
D_{\mathrm{c}}=I_{t}-|0\rangle\langle 0|=\sum_{s=1}^{t-1}|s\rangle\langle s|,$
(24)

where $|s\rangle\langle s|$ is the projector onto the pairing $s$, and $I_{t}$
is the identity matrix of size $t$. If the operator acts nontrivially on the
site we will insert a matrix $D_{\mathrm{c}}$ into the product of transfer
matrices, whereas if it acts trivially we insert $D_{0}$. Suppose now that the
operator $P_{\alpha}$ is made from $n$ clusters of non-identity operators and
each cluster contains $a_{m}$ on-site operators ($m=1,\cdots,n$). The
autocorrelation function of the operator $P_{\alpha}$ at large $q$ is then
given by

$\displaystyle\overline{C_{\alpha\alpha}(t)}$
$\displaystyle=q^{-2a}\prod_{m=1}^{n}\left[(TD_{\mathrm{c}})^{a_{i}}T\right]_{00}$
(25) $\displaystyle=q^{-2a}e^{-2n\varepsilon
t}(t-1)^{n}\left[1+(t-2)e^{-\varepsilon t}\right]^{a-n},$

provided we do not simultaneously have $n=1$ and $a=L$. In that case
$\overline{C_{\alpha\alpha}(t)}=q^{-2L}\overline{K(t)}$ at leading order in
$q$.

Eq. (25) is one of the main technical results of this paper. For
$t>t_{\mathrm{d}}$ (here $t_{\mathrm{d}}=1$) correlation functions increase
due to the above factors of $t$, and at later times this behaviour crosses
over into a decrease controlled by the factors $e^{-\varepsilon t}$ which come
from domain walls. The timescale that separates these two behaviours is
$2+\varepsilon^{-1}$, at which $(t-2)e^{-\varepsilon t}$ is maximized. Before
this time, Eq. (25) can be approximated by

$\overline{C_{\alpha\alpha}(t)}\simeq
q^{-2a}(t-1)^{n}(t-2)^{a-n}e^{-(a+n)\varepsilon t}.$ (26)

At late times, on the other hand, we have

$\displaystyle\overline{C_{\alpha\alpha}(t)}\simeq q^{-2a}$ $\displaystyle
e^{-2n\varepsilon t}(t-1)^{n}$
$\displaystyle\times\left[1+(a-n)(t-2)e^{-\varepsilon t}+\ldots\right],$ (27)

where the ellipsis denotes terms decaying faster than $e^{-\varepsilon t}$.
While the early-time behaviour of autocorrelation functions depends on the
internal structure of the operator, the late-time behaviour is much simpler.
Equation (III.2) shows that autocorrelation functions approach
$q^{-2a}e^{-2n\varepsilon t}(t-1)^{n}$; up to an overall $a$-dependent
prefactor, in this regime the autocorrelation function becomes insensitive to
the sizes of individual clusters.

The approach to this asymptotic value is controlled by the timescale
$t_{\mathrm{Th,\alpha}}$, which can be inferred from the late-time expansion
in Eq. (III.2). Over the timescale $t_{\mathrm{Th,\alpha}}$ the average
autocorrelation function drops to a value of order $q^{-2a}$, which implies
$t_{\mathrm{Th,\alpha}}=\varepsilon^{-1}\ln(a-n)$. As mentioned in Sec. II.3,
$t_{\mathrm{Th,\alpha}}$ bears a concrete physical meaning: beyond this
timescale the structure within each cluster becomes unimportant, and a cluster
of size $a_{m}$ behaves as though it evolves as a RMF system of dimension
$q^{a_{m}}$ that is subjected to dephasing induced by coupling at the boundary
to the rest of the system.

###  III.3 PSFF in the RPM

In this section, we provide a derivation of the exact PSFF in the RPM at
$q\to\infty$ reported in Sec. II.4. For the exact computation, it is
convenient to express the PSFF as

$K_{A}(t)=q^{-(L-L_{A})}\mathrm{Tr}_{\bar{A}}\left[(\mathrm{Tr}_{A}W(t)^{\dagger})(\mathrm{Tr}_{A}W(t))\right].$
(28)

Averaging this object, we find that in the region $A$ the leading pairings are
cyclic with $s=0,\ldots,(t-1)$; see Figs. 6a-6d. In region $\bar{A}$, on the
other hand, the leading pairing is $s=0$. The average PSFF is a sum over
configurations of such pairings, and hence can be expressed in terms of the
transfer matrix $T$ introduced above,

$\overline{K_{A}(t)}=[T^{L_{A}+1}]_{00}=\frac{1}{t}\left[\lambda_{0}^{L_{A}+1}+(t-1)\lambda^{L_{A}+1}\right]$
(29)

where $\lambda_{0}=1+(t-1)e^{-\varepsilon t}$ and $\lambda=1-e^{-\varepsilon
t}$ are eigenvalues of $T$, with degeneracies $1$ and $t-1$ respectively, and
$L_{A}<L$. Note that, unlike in autocorrelation functions, here we do not have
a $q$-dependent prefactor. The PSFF approaches unity at late times, with the
asymptotic behaviour

$\overline{K_{A}(t)}=1+\frac{t-1}{2}e^{-2\varepsilon(t-t_{\mathrm{Th},A})}+\cdots.$
(30)

The Thouless time $t_{\mathrm{Th},A}$ for the region $A$ is
$t_{\mathrm{Th},A}=\varepsilon^{-1}\ln L_{A}$. The infinite-time value is
controlled by the identity operator which does not evolve in time, hence
$\lim_{t\to\infty}\overline{K_{A}(t)}=q^{-L}{\rm Tr}I_{q^{L}}=1$.

(a)  (b)  (c)  (d)

Figure 6: Leading local pairings in the PSFF at $t=3$. On $A$ we have again
three cyclic pairings: (a) $s=0$, (b) $s=1$, and (c) $s=2$. (d) The only
leading local pairing $s=0$ on $\bar{A}$.

In fact, the decomposition of the PSFF as a sum over autocorrelation functions
in Eq. (5), along with our result in Eq.(25), provides us a more detailed
understanding of its behaviour. At first $\overline{K_{A}(t)}$ increases due
to the factor $t-1$, but this is eventually cut off by the domain wall cost
$e^{-\varepsilon t}$. The timescale associated with this corresponds to the
maximum of $(t-1)e^{-\varepsilon t}$, and we can infer which operators
contribute to the build-up of this peak from Eq. (25). First, recall that the
relevant operator strings are supported only in region $A$. From Eq. (26), it
is clear that the contributions from operator strings are larger when $a$ is
large and $n$ is small. Note, however, that we also need to take into account
the combinatorial factor associated with the number of ways of distributing
clusters when $n>1$. Since for a given $n=k$ (with $k=1,\cdots,\lfloor
n/2\rfloor-1$) the largest possible $a$ is $a=L_{A}+1-k$, the total
contribution from the operator strings with $a=L_{A}+1-k$ and $n=k$ is

$\displaystyle\begin{pmatrix}L_{A}-2\\\
k-1\end{pmatrix}(t-2)^{L_{A}+1-k}e^{-(L_{A}+1)\varepsilon t}$
$\displaystyle\quad\approx\frac{1}{k!}\times\left(\frac{L_{A}}{t-2}\right)^{k-1}\left[(t-2)e^{-\varepsilon
t}\right]^{L_{A}+1},$ (31)

where the factor of $q^{-2a}$ in individual autocorrelation functions has, at
leading order in $q$, cancelled with the number $q^{2a}-1$ of contributing
operators. Assuming that $t\ll L_{A}$, it is readily seen that the expression
is maximized when $k=L_{A}/(t-2)$. Thus the operator strings that contribute
to the onset of the second peak in the PSFF vary as the dynamics proceeds, and
at a given time $t$ they are specified by

$a=1+\left\lfloor L_{A}\frac{t-3}{t-2}\right\rfloor,\quad
n=\left\lfloor\frac{L_{A}}{t-2}\right\rfloor.$ (32)

Similarly, the late-time behaviour Eq. (30) can be recovered by recalling that
the operator strings consisting of a single cluster decay in the slowest way.
These operators have $a$ from $1$ to $L_{A}$, and summing over all such
operators we find factors of $L_{A}-a+1$ corresponding to their possible
locations within $A$. Combining these with the factor $q^{2a}$, the net
contribution of these operator strings to $\overline{K_{A}(t)}$ is

$\sum_{a=1}^{L_{A}}(L_{A}-a+1)e^{-2\varepsilon
t}(t-1)=\begin{pmatrix}L_{A}+1\\\ 2\end{pmatrix}e^{-2\varepsilon t}(t-1),$
(33)

which matches Eq. (30).

###  III.4 Numerics

Figure 7: Averaged autocorrelation functions $\overline{C_{\alpha\alpha}(t)}$
for operators consisting of a single cluster ($n=1$). The various supports $a$
are shown on the legend (with the main panel showing $a\leq 5$ only), and we
use system sizes $L=8$ (dotted), $L=10$ (dashed) and $L=12$ (solid). The inset
shows the $L=8$ data with the vertical axis scaled by $D^{2}-1=2^{2L}-1$ in
order to highlight the ramp. With this scaling, we would find a ramp $t$
within RMT; for clarity we show this behavior shifted relative to the data as
$t-10$ (dotted black line).

As shown earlier in this section, in the limit of large $q$ autocorrelation
functions develop considerable structure as a consequence of time-translation
symmetry. It is natural to ask which of these features carry over to finite
$q$, and in particular to the dynamics of chains of qubits, and so here we
discuss the behavior of autocorrelation functions in the RSM with $q=2$. Our
numerical methods are described in Appendix B.

First, in Fig. 7, we show $\overline{C_{\alpha\alpha}(t)}$ for operator
strings consisting of single clusters ($n=1$) of various lengths $a$. As
predicted by our large-$q$ calculations in the RPM, there is a significant
peak in $\overline{C_{\alpha\alpha}(t)}$ at early times, and the height of
this peak is independent of $L$. At finite $q$ the autocorrelation function
additionally has a ramp $\sim 2^{-2L}(t-1)$ at late times, and we emphasize
this feature in the inset for small $L$. For small $a$ we see a significant
vertical shift of the ramp relative to the RMT prediction, although even for
moderate large $L$ this feature is negligible compared with the early-time
peak.

In Fig. 8, we show the average PSFF. Perhaps the most notable feature, which
corroborates the large-$q$ result in Eq. (29), is the growth of the PSFF with
$L_{A}$. In the RPM, at any finite $t$ and at sufficiently large $L_{A}$, the
fact that $\lambda_{0}(t)\geq\lambda(t)$ implies that the growth with $L_{A}$
is exponential. This behaviour can be seen to persist at finite $q$ by
expressing the averaged PSFF in terms of the transfer matrix generating the
averaged SFF [20, 14] (see also Refs. [39, 10, 19] for studies of this
transfer matrix in related models) and the influence matrix introduced in Ref.
[40]. At finite $t$ the leading eigenvalue $\lambda_{0}(t)$ of the transfer
matrix is unique so, at large $L_{A}$,
$\overline{K_{A}(t)}\sim\lambda^{L_{A}}_{0}(t)$ in a system with periodic
boundary conditions.

While the data shown in Figs. 7 and 8 are for systems with periodic boundary
conditions, verifying the exponential growth of $\overline{K_{A}(t)}$ with
$L_{A}$ numerically is simplest in a system with open boundary conditions.
This is because, with periodic boundary conditions, $\overline{K_{A}(t)}$ has
contributions from all $t$ of the leading eigenvalues of the transfer matrix,
and large values of $L_{A}$ are necessary to resolve the gaps between them
even at moderate $t$. With open boundary conditions, on the other hand, taking
$A$ to be e.g. the left-hand $L_{A}$ sites of the system, we only have a
contribution from the leading eigenvalue (see the discussion of the effects of
boundary conditions on the SFF in Ref. [20]). For example, in the RPM at large
$q$ we have the equality $\overline{K_{A}(t)}=\lambda^{L_{A}}_{0}(t)$. At
finite $q$, Refs. [20, 14] and Ref. [40] imply that
$\overline{K_{A}(t)}\simeq\mathcal{B}(t)\lambda^{L_{A}}_{0}(t)$ for large $L$
and $L_{A}$. Here $\mathcal{B}(t)$ is controlled by the overlap of leading
eigenvectors of transfer and influence matrices, and so is independent of $L$
and $L_{A}$, while the ellipses denotes contributions from subleading
eigenvalues. This relation implies that for large $L_{A}$ we should observe
the scaling $\log\overline{K_{A}(t)}\propto L_{A}\log\lambda_{0}(t)$, and we
demonstrate this numerically in Fig. 9. At large $q$ we have
$\mathcal{B}(t)=1$, and although at $q=2$ there are small variations of
$\mathcal{B}(t)$ with $t$, the changes they induce in $\overline{K_{A}(t)}$
are much less than dramatic than those induced by the evolution of
$\lambda_{0}(t)$, which was studied in this model in Ref. [14]: at $t=1$ we
have $\lambda_{0}(1)=1$, with the eigenvalue then passing through a broad
maximum extending from $t=2$ to $t=10$, and then gradually decreasing to unity
at late times.

Figure 8: Averaged PSFFs $\overline{K_{A}(t)}$ for various $L_{A}$ (legend).
Here $L=12$.

(a)(b)

Figure 9: Scaling of the PSFF with open boundary conditions, and with
subregions $A$ and $\bar{A}$ each involving a system boundary. (a)
$\overline{K_{A}(t)}$ versus $t$ for various $L_{A}$ (with the same color
indexing as in Fig. 8). (b) Exponential growth of $\overline{K_{A}(t)}$ with
increasing $L_{A}$ for various $t$ (see legend); the rate of exponential
growth is set by the leading eigenvalue of the transfer matrix which also
generates the average SFF [see Ref. [14] for the analysis of the transfer
matrix for this model]. Here $L=14$ and the vertical scale of the figure is
logarithmic.

###  III.5 PSFF at finite $q$

So far we have focused exclusively on the large $q$ limit in the RPM. In the
present section we attempt to go beyond this limit and capture in a
quantitative way effects for $q$ large but finite. We find that finite $q$ has
important consequences in finite-size systems for both autocorrelation
functions and the PSFF as it restores their late-time ramps.

Accounting for the full ramifications of finite $q$ in a systematic way is a
rather demanding task and we leave it for future studies. Rather, here we
propose a possible mechanism by which the the ramp is induced in the PSFF. To
this end, it is illuminating to recall how the diagonal approximation [20],
which we define shortly, gives rise to the ramp at large but finite $D$ in RMF
systems.

####  III.5.1 Diagonal approximation in RMF systems

We start from the expression Eq. (28) with $W$ is a $D\times D$ Haar random
matrix. To perform the average it will be convenient to make the various index
summations explicit; here indices $a(r)=0,\ldots,(D_{A}-1)$ with
$r=0,\ldots,(t-1)$ label basis states in region $A$, while
$\bar{a}_{r}=0,\ldots,(D/D_{A}-1)$ label basis states its complement; we use
$b$ rather than $a$ for indices in the product of $W^{*}$. Using this notation
we have

$\displaystyle K_{A}(t)=\frac{D_{A}}{D}\sum_{\\{a\\},\\{b\\}\in
A}\sum_{\\{\bar{a}\\},\\{\bar{b}\\}\in\bar{A}}\delta_{\bar{a}(0)\bar{b}(0)}\delta_{\bar{a}(t)\bar{b}(t)}$
$\displaystyle\quad\times
W_{a(0)\bar{a}(0)}^{a(1)\bar{a}(1)}W_{a(1)\bar{a}(1)}^{a(2)\bar{a}(2)}\cdots
W_{a(t-1)\bar{a}(t-1)}^{a(0)\bar{a}(t)}$
$\displaystyle\quad\times[W^{*}]_{b(0)\bar{b}(0)}^{b(0)\bar{b}(1)}[W^{*}]_{b(1)\bar{b}(1)}^{b(2)\bar{b}(2)}\cdots[W^{*}]_{b(t-1)\bar{b}(t-1)}^{b(0)\bar{b}(t)}$
(34)

where $W_{ab}^{cd}=\langle a,b|W|c,d\rangle$. The factor
$\delta_{\bar{a}(0)\bar{b}(0)}\delta_{\bar{a}(t)\bar{b}(t)}$ implements the
contraction of $\bar{A}$ indices between $W(t)$ and $W^{*}(t)$ at the initial
and final times. Haar averaging then induces pairings among indices, and as
above we are interested in the cyclic pairings, i.e. those parameterized by
$(a(r),\bar{a}(r))=(b(r+s),\bar{b}(r+s))$ for all $r$ and fixed
$s=0,\cdots,t-1$; it can be shown that in this setting the cyclic pairings are
dominant for $D_{A},D/D_{A}\gg 1$. Note that, since RMF systems have no
spatial structure the pairing $s$ is here uniform throughout the system. By
contrast, in the RPM at large $q$ we found $s=0$ in region $\bar{A}$ while
$s=0,\ldots,(t-1)$ in region $A$.

In the computation of the SFF, the restriction to cyclic pairings without
spatial structure is known as the diagonal approximation. The sum over the $t$
cyclic pairing generates the characteristic ramp $\overline{K(t)}=t$ in RMF
systems, with each pairing contributing unity. In the computation of the PSFF,
on the other hand, the $s=0$ and $s\neq 0$ pairings are inequivalent. First,
the contribution from $s=0$ pairing in the summation in Eq. (III.5.1) simply
gives

$\frac{1}{D^{t}}\delta_{\bar{a}(t)\bar{b}(t)}\prod_{r=0}^{t-1}\delta_{a(r)b(r)}\delta_{\bar{a}(r)\bar{b}(r)}.$
(35)

Summing over $\\{b\\}$ and $\\{\bar{b}\\}$, we are left with $t$ and $t+1$
free indices on $A$ and $\bar{A}$. Hence after the sum over the remaining
indices, we get $D/D_{A}$. When pairings are $s\neq 0$, the indices used in
imposing the boundary condition on $\bar{A}$, i.e.
$\bar{a}(0),\bar{b}(0),\bar{a}(t),\bar{b}(t)$ also appear in the pairing
conditions $\delta_{\bar{a}(r)\bar{a}(r+s)}$, leaving only $t-1$ free indices
on $\bar{A}$. This means that the net contribution upon taking the sum is
$(t-1)D_{A}/D$. Combining both cases we obtain the diagonal approximation to
the PSFF

$\overline{K_{A}(t)}\simeq\frac{D^{2}_{A}}{D^{2}}(t-1)+1.$ (36)

Comparing this with Eq. (19) it can be verified that the diagonal
approximation is appropriate in this setting for $1\ll t\ll D$ and $1\ll
D_{A}\ll D$.

####  III.5.2 Local diagonal approximation in the RPM

Motivated by the above result, we make a similar approximation in computing
the finite $q$ PSFF in the RPM, which we call the local diagonal
approximation. Namely, we include the contribution from cyclic Gaussian
pairings on $\bar{A}$ with $s\neq 0$, as in the diagonal approximation for RMF
systems. Of course, at finite $q$ there are contributions from many other
possible pairings both on $A$ and $\bar{A}$ which we simply neglect. The
approximation generates a ramp at late times.

We first notice that on-site diagrams on $\bar{A}$ with $s\neq 0$ local
pairing (imagine replacing the $s=0$ pairing in Fig. 6d with that in Fig. 6b
and Fig. 6c) carry the factor $q^{-1}$. Further, accounting for the extra
$q^{-1}$ from the normalisation of trace, each $s\neq 0$ pairing carries a
factor of $q^{-2}$ whereas the $s=0$ pairing carries a factor of unity.
Writing the average PSFF as a sum over all cyclic pairings, including all
values of $s$ at each site, we find

$\overline{K_{A}(t)}=\mathrm{Tr}(T^{L_{A}}\hat{T}^{L-L_{A}}),$ (37)

where $\hat{T}=Td$ and the on-site diagonal matrix $d$ has entries
$d_{ss^{\prime}}=\delta_{ss^{\prime}}(\delta_{s0}+(1-\delta_{s0})q^{-2})$
where $s,s^{\prime}=0,\cdots,t-1$. It is immediate to see that this gives the
desired RMT behaviour at $t\to\infty$. At infinite time, the transfer matrix
becomes the identity matrix, hence

$\lim_{t\to\infty}\overline{K_{A}(t)}=\mathrm{Tr}(d^{L-L_{A}})=1+\frac{t-1}{q^{2(L-L_{A})}}.$
(38)

The result in Eq. (38) corresponds to uniform pairing throughout the system.

We must now ask when the behaviour in Eq. (38) sets in; this defines the time
scale $t_{\mathrm{ramp}}$. In fact, that can be inferred readily from an
argument based on domain walls. As we already know, at large $q$, the decay of
the PSFF is governed by two domain walls supported on $A$ where each domain
wall comes with a statistical cost $e^{-\varepsilon t}$ (see Fig. 10b, where
each domain wall is indicated by a red cross). This is still the case at
finite $q$ within the local diagonal approximation, but now domain walls can
be supported anywhere in the system.

The next-to-leading order corrections within the local diagonal approximation
correspond to two different kinds of domain wall configurations. The first
corresponds to a single-site $s\neq 0$ domain within $\bar{A}$, with a $s=0$
domain covering the rest of the system (see Fig. 10c). The number of such
pairing configurations is $(L-L_{A})(t-1)$. Thus, their total contribution to
the average PSFF is

$(L-L_{A})(t-1)\frac{e^{-2\varepsilon t}}{q^{2}}.$ (39)

In the second kind of configuration, there is a $s\neq 0$ domain having one of
its walls in $A$ and one in $\bar{A}$. Since each site having $s\neq 0$ in
$\bar{A}$ costs a factor $q^{-2}$, the leading contribution of this kind
corresponds to a domain that has only a single site in $\bar{A}$ [see Fig.
10d]. Summing over the $L_{A}$ possible domain wall locations we find a
contribution

$2L_{A}(t-1)\frac{e^{-2\varepsilon t}}{q^{2}},$ (40)

(a)  (b)  (c)

(d)

Figure 10: Leading domain wall configurations to the late-time behaviour of
the PSFF at large but finite $q$ with periodic boundary conditions. A red
cross indicates the location of a domain wall and the shaded blue region
corresponds to region $A$. (a), (b), and (c) and (d) give the second, third,
and fourth terms in Eq. (41), respectively.

and for the PSFF we then have

$\displaystyle\overline{K_{A}(t)}\simeq 1+\frac{t-1}{q^{2(L-L_{A})}}$ (41)
$\displaystyle\quad+\begin{pmatrix}L_{A}+1\\\
2\end{pmatrix}(t-1)e^{-2\varepsilon t}+(L+L_{A})(t-1)\frac{e^{-2\varepsilon
t}}{q^{2}}.$

at late times within the local diagonal approximation. The first and second
terms come from configurations without domain walls, and correspond to $s=0$
and $s\neq 0$, respectively. The third term, which is leading order in $q$
comes from domain walls in or on the boundary of $A$. The fourth term, which
we now see is always smaller than the third, corresponds to domain wall
configurations of the kind described above.

From Eq. (41) we can identify $t_{\mathrm{ramp}}$. In brief, the domain wall
contribution to the PSFF, which arises even at leading order in $q$, decays as
$e^{-2\epsilon t}$, and eventually this becomes smaller than the ramp expected
from RMT, which here has prefactor $q^{2(L_{A}-L)}$. Therefore

$t_{\mathrm{ramp}}\sim\frac{\ln q}{\varepsilon}(L-L_{A}).$ (42)

Although this timescale is proportional to the distance $L-L_{A}$, it should
not be interpreted in terms of a velocity since it is simply the time at which
two distinct contributions to the PSFF exchange dominance.

###  III.6 Spectral form factor

In closing this section, we offer a brief overview of the relation between the
spectral statistics of many-body Floquet systems and the dynamics of
operators. Expressing the SFF as a sum over all autocorrelation functions, we
can understand the ramp as the regime where all operators contribute equally.
That is, each of the $D^{2}-1$ non-identity operators contributes
$(D^{2}-1)^{-1}(t-1)$, with the identity contributing unity. At early times,
on the other hand, there is a peak in $\overline{K(t)}$, and we can now
understand this peak from the perspective of operator dynamics. At late times,
we have seen that the dynamics of an operator string that is dense over a
given set of $a$ sites can be understood by modelling the dynamics with (i) a
$q^{a}$-dimensional RMF system and (ii) dephasing through interactions with
neighbouring regions. This picture suggests
$\overline{C_{\alpha\alpha}(t)}\sim q^{-2a}(t-1)e^{-2\varepsilon t}$. The
number of such operators is $\sim q^{2a}$, and therefore each $a$-site region
contributes $\sim(t-1)e^{-2\varepsilon t}$ to the SFF (where we neglect
faster-decaying contributions). Summing over all of these processes, we
recover the SFF of the RPM
$\overline{K(t)}=1+\frac{1}{2}L(L-1)e^{-2\varepsilon t}+\ldots$.

The time scale at which RMT behavior of the SFF sets in, which we call the
spectral Thouless time, is given by $t_{\text{Th}}=\varepsilon^{-1}\log L$.
Let us compare this with the dynamical Thouless times $t_{\text{Th},\alpha}$
associated with individual operators $P_{\alpha}$. For an operator consisting
of $n$ clusters, autocorrelation functions decay as $e^{-2n\varepsilon t}$ at
late times; from the above expressions we can infer that the late-time
behaviour of the SFF is dominated by operators with $n=1$. For large operator
support $a$, the dynamical Thouless time
$t_{\text{Th},\alpha}\simeq\varepsilon^{-1}\log a$, and so naively the
spectral Thouless time appears to be the dynamical Thouless time for $a=L$,
i.e. for single cluster operators with support over the entire system. This is
not quite correct, however, because for any $a<L$ the SFF picks up an entropic
factor from the sum over spatial locations of operators, and this also leads
to a time scale $\varepsilon^{-1}\log L$. The implication is the spectral
Thouless time emerges from the sum over autocorelation functions of all single
cluster operators. This picture is complementary to the discussions of orbit
pairing in Refs. [13, 20].

The situation above should be compared with that in Floquet systems having
scalar charge conservation. Working within the diagonal approximation, Refs.
[41, 42, 43] found $t_{\text{RMT}}\sim L^{2}$, and so there is a very clear
connection between time scales appearing the spectrum and in the (diffusive)
dynamics of a local charge densities.

##  IV Out-of-time-order correlators

In this section we calculate the average OTOC for the RPM. In RUCs the
standard calculation of the OTOC [4, 5] involves mapping the dynamics of
operators to a kind of classical stochastic dynamics. There one finds that
operators grows ballistically, and that the operator fronts broaden
diffusively as $\sim(\mathcal{D}t)^{1/2}$. Here we recover this behaviour,
including the broadening of the front, in the large-$q$ limit of a Floquet
model [c.f. Ref. [9]].

In generic Floquet systems it is not possible to directly map the calculation
of the OTOC to a property of a stochastic process, simply because unitary
operations cannot be averaged independently at different time steps, and so
here we will instead evaluate the OTOC using a spatial transfer matrix.

###  IV.1 Diagrammatics and Haar averages

Compared to the computation of averaged autocorrelation functions and the PSFF
at large $q$, that of the OTOC is more complicated. There are two reasons for
this. First, the fact that the OTOC involves two copies of both $W(t)$ and its
conjugate means that the transfer matrix has a dimension that is twice as
large as that for the PSFF. Second, the pairings of unitary operators which
contribute at leading order in $q$ can be non-Gaussian [38], whereas for the
PSFF are all contributing pairings are Gaussian at this order.

In Fig. 11a we illustrate the basic structure of the OTOC, as defined in Eq.
(7), and in Fig. 11b-11d we show the contractions of indices relevant to the
different sites.

(a)  (b)  (c)  (d)

Figure 11: (a): Diagrammatic representation of $\mathcal{C}(x,t)$. (b, c, d):
Single-site diagram on $0$, $x$, and other sites.

The strategy we follow is the same as before. We first evaluate averages over
one-site unitary operators $U_{x}$ and enumerate all relevant leading-order
pairings. This will allow us to express the ensemble-averaged OTOC as a
product of transfer matrices, where these objects can be viewed as acting on
the space of pairings. However, since non-Gaussian pairings now contribute at
leading order, it will be useful to represent diagrams in a slightly different
way.

Let us start with a site $x^{\prime}\neq 0,x$, which contains no local
operators [see Fig. 11d]. The on-site diagram consists of two Floquet
operators $W$ and their conjugates $W^{*}$, with indices contracted as in the
figure. An equivalent way of representing this is to draw a square and assign
a Floquet operator (or its conjugate $W^{*}$) to each edge. To illustrate
this, consider first $t=1$: in Fig. 12a we show the contractions of unitary
operators using both “U-shape” and “square-shape” diagrams, where indices
$a(0),b(0),c(0),d(0)$ label the boundaries of the edges in both diagrams.
Comparing the location of the labels in both diagrams, one can also infer the
correspondence between edges.

At $t=1$, there are four pairings: two are Gaussian (see Figs. 12a and 12b),
and two are non-Gaussian (see Figs. 12c and 12d). From Ref. [38] the values of
these pairings can be readily obtained for Figs. 12a and 12b as
$q^{3}/(q^{2}-1)$, for 12c as $-q^{3}/(q^{2}-1)$, and for 12d
as$-q/(q^{2}-1)$. At large $q$ these become $q,-q$, and $-q^{-2}$,
respectively. This indicates that, at large $q$, we have contributions from
the Gaussian diagrams 12a and 12b, as well as non-Gaussian diagrams 12c which
have no crossings in the square representation. Non-Gaussian diagrams with
crossings, such as 12d, are subleading. Note that, if we keep track of only
the Gaussian diagrams and non-crossing non-Gaussian diagrams, we can again
avoid distinguishing first and second indices of unitaries $U_{x}$ provided we
use square diagrams. For this reason, as in the previous section, we will just
use black dots to indicate the unitaries $U_{x}$. As we will see shortly, the
inclusion of non-Gaussian diagrams turn out to be crucial for the OTOC to
behave in the way we expect.

(a)  (b)  (c)  (d)

Figure 12: (a,b): Gaussian pairings carrying the value $q$ and (c, d): non-
Gaussian pairings carrying the values $-q$ and $-q^{-2}$, respectively, at
$t=1$.

The analysis above can be readily carried over to generic $t>1$, in which case
we have $t$ unitaries on each edge. This means that on each edge there are
$t-1$ bonds that connect unitaries. On the left and right edges of the square
diagrams (which each correspond to $W(t)$) we assign indices
$a(1),a(2),\cdots,a(t-1)$ and $c(1),c(2),\cdots,c(t-1)$ to these bonds,
whereas on the upper and lower edges (which each correspond to $W^{*}(t)$) we
assign indices $b(t-1),b(t-2),\cdots,a(1)$ and $d(t-1),d(t-2),\cdots,d(1)$;
our ordering convention is shown in Fig. 13a. The two ends of each edge also
carry indices, which we denote e.g. $a(0)$ and $a(t)$, but the trace structure
imposes conditions $a(0)=d(0),a(t)=b(t),b(0)=c(0),c(t)=d(t)$, leaving four
free indices out of eight edge-attached indices.

With the above convention the relevant Gaussian diagrams, which we denote by
$g_{p}$ with $p=1,\cdots,t-1$, correspond to indices paired as

$\displaystyle a(\tau)=b(\tau),\quad c(\tau)=d(\tau),$ $\displaystyle\quad
p\leq\tau\leq t-1$ (43) $\displaystyle a(\tau)=d(\tau),\quad b(\tau)=c(\tau),$
$\displaystyle\quad 1\leq\tau\leq p.$ (44)

This parametrisation, however, does not cover two diagrams that are given by
$a(\tau)=b(\tau),\,c(\tau)=d(\tau)$ and $a(\tau)=d(\tau),\,b(\tau)=c(\tau)$
for $1\leq\tau\leq t-1$, which we denote by $g_{t}$ and $g_{0}$, respectively.
We therefore have in total $t+1$ Gaussian diagrams $\\{g_{p}\\}_{p=0}^{t}$,
and in Figs. 13a-13d we depict $g_{0},\cdots,g_{3}$ at $t=3$. Evaluating the
various diagrams, we find that each carries an overall weight $q$.

The above labelling convention is also useful for enumerating the contributing
non-Gaussian diagrams. Those which generate contributions that are leading
order in $q$ are obtained by stacking two consecutive Gaussian diagrams (e.g.
$p=p_{1},p_{1}+1$) on top of each other; see for example the stacking of Figs.
13a and 13b to generate the non-Gaussian diagram Fig. 13e. The non-Gaussian
diagrams, which we denote by $n_{p}$ for $0\leq p\leq t-1$, correspond to the
pairings

$\displaystyle a(\tau)=b(\tau),\quad c(\tau)=d(\tau),\quad$ $\displaystyle
p<\tau\leq t$ (45) $\displaystyle a(\tau)=d(\tau),\quad b(\tau)=c(\tau),\quad$
$\displaystyle 0\leq\tau\leq p.$ (46)

The non-Gaussian pairings $n_{1},n_{2}$, and $n_{3}$ are depicted in Figs.
13e, 13e, and 13g, and we find that each $n_{p}$ carries a weight $-q$.

(a)  (b)

(c)  (d)

(e)  (f)

(g)

Figure 13: (a)-(d): Gaussian pairings $g_{0},\cdots,g_{3}$ and (e)-(g): the
non-Gaussian pairings $n_{0},n_{1}$, and $n_{2}$ obtained as superpositions of
the Gaussian diagrams at $t=3$.

Haar averaging for each site on $y\neq 0,x$ therefore yields $2t+1$ leading
pairings where $t+1$ of them are Gaussian carrying the factor $q$ and $t$ of
them are non-Gaussian with the factor $-q$. We organise these pairings by
introducing a new set of labels $\\{h_{p}\\}_{p=0}^{2t}$ defined by

$h_{p}:=\begin{cases}g_{p}&p=0,1,\cdots,t\\\
n_{p}&p=t+1,t+2,\cdots,2t\end{cases}$ (47)

so that the first $t+1$ pairings are Gaussian whereas the latter $t$ pairings
are non-Gaussian. For our transfer matrix calculation it will be convenient to
associate the $2t+1$ pairings with an orthonormal basis $|p\rangle$ of vectors
labelled by the index $p=0,\ldots,2t$.

We now move on to evaluate the contributions from nearest-neighbour
interactions in the model. First, if the sites $y$ and $y+1$ have local
pairings corresponding to the diagram Fig. 13a and 13g at $t=3$, the non-
vanishing contribution to the phase is

$\displaystyle\Phi$
$\displaystyle=\varphi_{y;a_{y}(1),a_{y+1}(1)}+\varphi_{y;a_{y}(2),a_{y+1}(2)}-\varphi_{y;b_{y}(1),b_{y+1}(1)}$
$\displaystyle\quad-\varphi_{y;b_{y}(2),b_{y+1}(2)}+\varphi_{y;c_{y}(1),c_{y+1}(1)}+\varphi_{y;c_{y}(2),c_{y+1}(2)}$
$\displaystyle\quad-\varphi_{y;d_{y}(1),d_{y+1}(1)}-\varphi_{y;d_{y}(2),d_{y+1}(2)}.$
(48)

Upon averaging over the phases, we find $\langle e^{{\rm
i}\Phi}\rangle=\rho^{2}$ where $\rho=e^{-2\varepsilon}$. The pattern can be
readily inferred for arbitrary $t$, from which we obtain the
$(2t+1)\times(2t+1)$ transfer matrix

$S=\begin{pmatrix}S_{1}&S_{2}\\\ S_{2}^{\mathrm{T}}&S_{3}\end{pmatrix},$ (49)

where $S_{1}$, $S_{2}$, and $S_{3}$ are $(t+1)\times(t+1)$, $(t+1)\times t$,
and $t\times t$ matrices, respectively. Their matrix elements are given by
$[S_{1}]_{pp^{\prime}}=\delta_{pp^{\prime}}+\rho^{|p-p^{\prime}|-1}(1-\delta_{pp^{\prime}})$,
$[S_{2}]_{pp^{\prime}}=\delta_{pp^{\prime}}+\rho^{|p-p^{\prime}|-1}(1-\delta_{pp^{\prime}}+\Theta(p^{\prime}-p)(\rho-1))$,
and
$[S_{3}]_{pp^{\prime}}=\delta_{pp^{\prime}}+\rho^{|p-p^{\prime}|}(1-\delta_{pp^{\prime}})$,
where indices $p,p^{\prime}$ run from $0$ to $t$ and $\Theta(p)$ is the step
function with $\Theta(0)=0$. For instance, the transfer matrix for $t=3$ reads

$S=\begin{pmatrix}1&1&\rho&\rho^{2}&1&\rho&\rho^{2}\\\ 1&1&1&\rho&1&1&\rho\\\
\rho&1&1&1&\rho&1&1\\\ \rho^{2}&\rho&1&1&\rho^{2}&\rho&1\\\
1&1&\rho&\rho^{2}&1&\rho&\rho^{2}\\\ \rho&1&1&\rho&\rho&1&\rho\\\
\rho^{2}&\rho&1&1&\rho^{2}&\rho&1\end{pmatrix}.$ (50)

The transfer matrix as defined above does not account for the different signs
of pairings (non-Gaussian pairings can make negative contributions), so we
will also introduce an on-site $(2t+1)\times(2t+1)$ matrix $qE$ that encodes
this information

$E=\begin{pmatrix}I_{t+1}&0\\\ 0&-I_{t}\end{pmatrix}.$ (51)

Finally, we turn to Haar averages over unitary operators on sites $0$ and $x$.
In this case the trace structures, shown in Fig. 11a and 11c, respectively,
admit only one leading pairing because $\mathcal{O}_{0}$ and
$\mathcal{O}_{x}=0$ are traceless. More precisely, on site $0$ only the
diagonal pairing $a(\tau)=d(\tau),\,b(\tau)=c(\tau)$ contributes, whereas on
site $x$ the anti-diagonal pairing $a(\tau)=b(\tau),\,c(\tau)=d(\tau)$
contributes; each of these pairings carries a weight $q$. Due to the different
index contractions at sites $0$ and $x$ it will be convenient to introduce new
on-site matrices for the transfer matrix calculation. These are
$E_{0}=q|0\rangle\langle 0|$ and $E_{t}=q|t\rangle\langle t|$, respectively.

(a)  (b)

Figure 14: The only leading pairing on (a): site $0$ and (b): site $x$ at
$t=3$.

###  IV.2 Computation of the OTOC at large $q$

Now we are ready to compute the OTOC at large $q$. Note first that a factor
$q^{-L}$ in the definition of the OTOC cancels the $L$ factors of $q$
corresponding to on-site pairings. This leads us to the expressions

$\overline{\mathcal{C}(x,t)}=\begin{dcases*}1-\hat{S}^{L-x}_{t0}\hat{S}^{x}_{0t}&$x>0$\\\
1-\hat{S}^{L-|x|}_{0t}\hat{S}^{|x|}_{t0}&$x<0$\end{dcases*},$ (52)

where $\hat{S}=ES$, $\hat{S}_{pp^{\prime}}=\langle
p|\hat{S}|p^{\prime}\rangle$, $E|0\rangle=|0\rangle$, and
$E|t\rangle=|t\rangle$.

The transfer matrix $\hat{S}$ has a number of interesting properties. First,
it is readily seen that the matrix always has an eigenvalue 1, as well as $2t$
degenerate eigenvalues equal to $0$. Since there are only two linearly
independent eigenvectors, $v_{1}=|t\rangle+|2t\rangle$ and
$v_{2}=|0\rangle+|t+1\rangle$, associated with the eigenvalue 0 (i.e. the
geometric multiplicity of the eigenvalue $0$ is 2), the matrix $\hat{S}$ is
not diagonalisable. However, it becomes diagonalisable if we take higher
enough powers:

$\hat{S}^{k+t}=\hat{S}^{t}=P,$ (53)

where $k\geq 0$ and the matrix elements of $P$ are given by
$P_{pp^{\prime}}=1-2\Theta(p^{\prime}-t)$. Eq. (53) implies that when the
power $p$ reaches $t$, $S^{p}$ becomes a rank-one projection matrix $P$, which
is diagonalisable. This is in accordance with what we expect: up to $t=|x|$
the operator front of $\mathcal{O}(0,t)$ will not have reached the site $x$,
hence, for large enough $L$, the OTOC is zero. Note also that, in the
thermodynamic limit, $\lim_{L\to\infty}\hat{S}^{L-|x|}_{t0}=P_{t0}=1$. The
question now is how to access the matrix elements of $\hat{S}^{\ell}$ when
$\ell<t$.

While obtaining matrix elements is in general a formidable task, some patterns
in their structure are apparent already for small $p$ and $t$. This motivates
an ansatz $\hat{S}^{\ell}_{pt}=F(\ell,t-p)$, where

$F(\ell,k)=1-(1-\rho)^{\ell}\sum_{i=0}^{k-\ell-1}\begin{pmatrix}\ell+i-1\\\
i\end{pmatrix}\rho^{i},$ (54)

or equivalently

$F(\ell,k)=1-\sum_{i=0}^{k-\ell-1}\begin{pmatrix}k-1\\\
i\end{pmatrix}\rho^{i}(1-\rho)^{k-i-1},$ (55)

which we now show is correct by induction. First consider $p=0$: using the
above matrix elements of $\hat{S}$, the matrix element $\hat{S}_{0t}^{\ell+1}$
is

$\hat{S}_{0t}^{\ell+1}=F(\ell,t)+\sum_{i=0}^{t-1}\left[F(\ell,t-i-1)-F(\ell,t-i)\right]\rho^{i}.$
(56)

Plugging Eq. (54) into the square bracket, it is readily seen that

$\displaystyle
F(\ell,t)+\sum_{i=0}^{t-1}\left[F(\ell,t-i-1)-F(\ell,t-i)\right]\rho^{i}$
$\displaystyle=F(\ell,t)+\begin{pmatrix}t-1\\\
t-\ell-1\end{pmatrix}\rho^{t-\ell-1}(1-\rho)^{\ell}.$ (57)

Using Eq. (55), the right hand side precisely coincides with $F(\ell+1,t)$,
establishing that $\hat{S}_{0t}^{\ell+1}=F(\ell+1,t)$. The induction for $p>0$
can be carried out in the same way, and we can similarly show that
$\hat{S}^{\ell}_{t0}=\hat{S}^{\ell}_{t0}$. From these results we obtain the
exact OTOC at $q\to\infty$. Outside of the light cone (i.e. for $|x|>t$) we
have $\overline{\mathcal{C}(x,t)}=0$, whereas for $|x|\leq t$ we have

$\overline{\mathcal{C}(x,t)}=\sum_{i=0}^{t-|x|-1}\begin{pmatrix}t-1\\\
i\end{pmatrix}\rho^{i}(1-\rho)^{t-i-1}.$ (58)

Note that the same function appeared in the computation of the OTOC in RUCs
[4]. For large $x$ and $t$ with $|x|<t$ we have

$\displaystyle\overline{\mathcal{C}(x,t)}$
$\displaystyle\simeq\frac{1}{\sqrt{4\pi\mathcal{D}t}}\int_{0}^{t-|x|}{\rm
d}y\,\exp[-\frac{(y-t\rho)^{2}}{4\mathcal{D}t}]$
$\displaystyle\simeq\Phi\left(\frac{v_{B}t-|x|}{\sqrt{2\mathcal{D}t}}\right)$
(59)

where $\Phi(x):=\int_{-\infty}^{x}{\rm d}y\,e^{-y^{2}/2}/(\sqrt{2\pi})$ with
the butterfly velocity $v_{B}=1-\rho$ and the diffusion constant
$\mathcal{D}=\frac{1}{2}\rho(1-\rho)$ as highlighted above in Eq. (21). As
mentioned above, even though the form of the OTOC is quite similar to that in
finite-$q$ RUCs with Haar-random two-site gates, sending $q$ to infinity in
that setting causes $\mathcal{D}\to 0$ and the butterfly velocity $v_{B}\to
1$, the geometrical lightcone velocity. This is markedly different from the
case in the RPM where, even in this limit, the OTOC retains the
characteristics expected for chaotic quantum many-body systems in one
dimension.

##  V Operator paths

A useful perspective on our results in the previous sections comes from
expressing quantities in terms of sums over paths in the space of operators.
In RUCs with Haar-random gates and without time-translation symmetry, the
structure of the paths contributing to correlation functions was studied in
Ref. [31]. In that setting ensemble-averaged correlation functions vanish, but
rich structures arise in their magnitudes; for squared autocorrelation
functions it was shown that $\overline{C^{2}_{\alpha\alpha}(t)}$ is dominated
by operator paths whose widths grow as $\sim t^{1/2}$. These paths should be
contrasted with those contributing to the OTOC, whose widths grow
ballistically [4, 5]. Separately, we note that Refs. [32, 33] have studied the
operator paths relevant for the description of hydrodynamics in RUCs with
conserved quantities; see also Refs. [7, 8] for a discussion of their OTOCs. A
related idea, which involves analysing operator paths in a Krylov space rather
than in real space, has been introduced to characterise Hamiltonian systems
[44].

In Sec. V.1 we study the operator paths contributing to autocorrelation
functions in random Floquet circuits. Our main result is to show that time-
translation symmetry in the dynamics causes a specific subset of these paths
to have real non-negative amplitudes, and therefore to survive an ensemble
average. The paths we identify have a spatially localised support that does
not change with time. Strikingly, contributions from these localised paths
alone generate the full ensemble-averaged autocorrelation function
$\overline{C_{\alpha\alpha}(t)}$ of the RPM at large $q$. Following this, in
Sec. V.2 we show that, in the RPM at large-$q$, the operator paths
contributing to the ensemble-averaged OTOC are independent of time-translation
symmetry, in the sense that they are the same as in RUCs. In Sec. V.3 we show
numerically that these pictures of the paths contributing to the
autocorrelation function and to the OTOC hold also for qubit chains.

To streamline our discussion we will focus on systems of qubits, although the
following discussion can be generalised. For $q=2$ we can use a tensor product
basis of Pauli strings $P_{\alpha}=P_{\alpha_{1}}\cdots P_{\alpha_{L}}$, where
each of $P_{\alpha_{x}}$ is a Pauli operator or the identity. The operators
$P_{\alpha}$ are unitary, Hermitian, and have simple commutation and
anticommutation relations $P_{\alpha}P_{\beta}=\pm P_{\beta}P_{\alpha}$. Where
we discuss the limit of large $q$, we can for simplicity restrict ourselves to
the case where $q$ is a power of $2$, since then each site $x=1,\ldots,L$ can
viewed as a set of $\log_{2}q$ qubits, and we can continue to use an operator
basis of Pauli strings.

###  V.1 Amplitudes

Here we ask which kinds of operator paths give rise to structure in ensemble-
averaged autocorrelation functions. First we discuss the case where the
Floquet operator is a single $D\times D$ Haar-random unitary matrix, and hence
the only nontrivial feature is the ramp. Since the SFF and PSFF are sums over
autocorrelation functions, this discussion will also provide an explanation
for the ramps appearing in those quantities. Following this we study the role
of locality via the RPM at large-$q$. In this setting the operator paths
contributing to autocorrelation functions of local operators will turn out to
be themselves localised.

Our discussion is based on Eq. (3). It will be useful to expand the
correlators $C_{\alpha_{r}\alpha_{r+1}}$ for individual time steps as

$\displaystyle C_{\alpha(r)\alpha(r+1)}$
$\displaystyle=D^{-1}\sum_{\begin{subarray}{c}a(r),a^{\prime}(r)\\\
b(r),b^{\prime}(r)\end{subarray}}W^{*}_{b(r)b^{\prime}(r)}W_{a(r)a^{\prime}(r)}$
$\displaystyle\times[P_{\alpha(r)}]_{b(r)a(r)}[P_{\alpha(r+1)}]_{a^{\prime}(r)b^{\prime}(r)},$
(60)

and we now turn to the Haar random problem.

####  V.1.1 Paths in RMF systems

First consider the amplitude $C_{\alpha(0)\alpha(1)}\ldots
C_{\alpha(t-1)\alpha(0)}$ for a single operator contributing to the
autocorrelation function. Inserting Eq. (60) into this expression, we see that
some terms in the sum have non-negative contributions for all Floquet
operators $W$. These are the ones where $a(\sigma(r))=b(r)$ and
$a^{\prime}(\sigma(r))=b^{\prime}(r)$ for $r=0,\ldots,(t-1)$, with $\sigma$ a
permutation of $t$ elements. There are $t!$ such permutations, and for each of
them we find that the corresponding contribution to the amplitude involves a
factor $\prod_{r=0}^{t-1}|W_{a(r)a^{\prime}(r)}|^{2}$.

In an unstructured RMF system we will now see that $t-1$ of the $t$ cyclic
permutations dominate for $1<t\ll D$. It is helpful to isolate the cyclic
permutations, $\sigma(r)=r+s$, with $s=1\ldots(t-1)$ and addition of indices
defined modulo $t$. The cyclic permutation with $s=0$ can be verified to give
a vanishing contribution as a consequence of the fact that all operators in
the path are traceless. We will additionally write

$\displaystyle|W_{a(r)a^{\prime}(r)}|^{2}=D^{-1}+\ldots,$ (61)

where the ellipsis denotes a fluctuation with average zero that is of the same
order as $D^{-1}$. Making the restriction that $t$ is prime for simplicity (we
will relax this below), and summing over states $a(r),a^{\prime}(r)$, then
gives

$\displaystyle C_{\alpha(0)\alpha(1)}\ldots C_{\alpha(t-1)\alpha(0)}$ (62)
$\displaystyle=D^{-2t}\sum_{s=1}^{t-1}|\mathrm{Tr}(P_{\alpha(0)}P_{\alpha(s)}\cdots
P_{\alpha([t-1]s)})|^{2}+...$

In this expression the ellipsis denotes (i) contributions from non-cyclic
pairings of indices and (ii) fluctuations inherited from Eq. (61). The
appearance of a single squared trace is a consequence of the fact that $s$ and
$t$ have no common factors. It is clear that the contribution to the operator
path amplitude displayed on the right-hand side of Eq. (62) is non-negative
for all unitary evolution operators $W$. If we average $W$ over the Haar
distribution of $D\times D$ unitary matrices, we find that only the displayed
term survives at large $D$. For simplicity, we will restrict to such an
average in the remainder of this subsection.

(a)  (b)

(c)

Figure 15: Cyclic pairings in $C_{\alpha_{0}\alpha_{1}}\cdots
C_{\alpha_{t-1}\alpha_{0}}$ at $t=3$.

Observe now that $|\mathrm{Tr}(P_{\alpha(0)}P_{\alpha(s)}\cdots
P_{\alpha([t-1]s)})|^{2}=D^{2}$ only if

$\displaystyle P_{\alpha(0)}P_{\alpha(s)}\cdots P_{\alpha([t-1]s)}=\pm I$ (63)

and is otherwise zero. We will refer to as a selection rule on operator paths.
Note that, if a given operator path satisfies this condition for one value of
$s$, it also satisfies it for all other values of $s$, simply because
$P_{\alpha}P_{\beta}=\pm P_{\beta}P_{\alpha}$ for our operator basis of Pauli
strings. Additionally, if a given operator path has nonzero amplitude, within
this model of Haar-random matrices any reordering of that path also has
nonzero amplitude.

The implication is that operator paths satisfying Eq. (63), which make up a
fraction of just $D^{-2}$ of all possible operator paths, have averaged
amplitude

$\displaystyle\overline{C_{\alpha(0)\alpha(1)}\ldots
C_{\alpha(t-1)\alpha(0)}}=D^{-2(t-1)}(t-1),$ (64)

while all others have zero averaged amplitude. So far we have shown that there
is a ramp at the level of individual operator path amplitudes, and that this
ramp can be understood as arising from a restricted set of paths.

The analysis is similar for non-prime $t$. For cyclic pairings of state
indices $b(r)=a(r+s)$ and $b^{\prime}(r)=a^{\prime}(r+s)$ we found a
contribution proportional to $|\mathrm{Tr}(P_{\alpha(0)}P_{\alpha(s)}\cdots
P_{\alpha([t-1]s)})|^{2}$ whenever $s$ and $t$ share no common factors [Eq.
(62)]. If $s$ and $t$ share a common factor, then instead of a single squared
trace we find a product of them. For example, summing over $a(r)$ and $b(r)$
for $t=6$ and $s=2$ we find the product of
$|\text{Tr}[P_{\alpha(0)}P_{\alpha(2)}P_{\alpha(4)}]|^{2}$ and
$|\text{Tr}[P_{\alpha(1)}P_{\alpha(3)}P_{\alpha(5)}]|^{2}$. This object is
zero for all but a fraction $D^{-4}$ of paths, but when it is nonzero we find
a contribution to the amplitude $D^{-2(t-2)}$, c.f. Eq. (64). More generally,
if we denote by $m$ and $n$ the smallest integers satisfying $ms=nt$, we find
a contribution

$\displaystyle\prod_{r=0}^{t/m-1}|\text{Tr}[P_{\alpha(r)}P_{\alpha(r+s)}\ldots
P_{\alpha(r+(m-1)s)}]|^{2},$ (65)

and summing over $s$ we find symmetry factors analogous to the factor $t-1$ in
Eq. (64). For example if for $t=6$ the $s=2$ contribution to the amplitude is
nonzero, the property $P_{\alpha}P_{\beta}=\pm P_{\beta}P_{\alpha}$ implies
that the $s=4$ contribution is also nonzero.

By summing over non-negative contributions to operator amplitudes, such as on
the right-hand side of Eq. (64), we will now recover

$\displaystyle\overline{C_{\alpha(0)\alpha(0)}(t)}=D^{-2}(t-1)$ (66)

which is the leading contribution to the autocorrelation function for $t\ll
D$. The implication is that the late-time ramp in autocorrelation functions
comes from operator paths satisfying the above selection rules.

To see how Eq. (66) arises, let us sum the right-hand side of Eq. (62) over
the $D^{2(t-1)}$ paths which begin and end at $P_{\alpha(0)}$. Out of these,
$D^{2(t-2)}$ satisfy Eq. (63), and so for prime $t$ we find an overall
contribution $D^{-2}(t-1)$ from Eq. (64). For $t$ that is not prime the result
is the same at large $D$, but it is simpler to consider each cyclic
permutation $s=1\ldots(t-1)$ independently. For $s$ and $t$ sharing no common
factors the $D^{2(t-2)}$ paths satisfying Eq. (63) make an overall
contribution $D^{-2}$ as above, whereas if there exist integers $n<s$ and
$m<t$ for which $ms=nt$ we find contributions of the form Eq. (65). The number
of operator paths for which this expression is nonzero is $D^{2t(1-1/m)-2}$,
but for each path we find a contribution $D^{2t(1/m-1)}$, and therefore each
of the $t-1$ allowed values of $s$ gives an overall contribution $D^{-2}$.
Their sum is the term displayed on the right-hand side of Eq. (66).

To recover the ramps in the PSFF and SFF we need only sum over autocorrelation
functions. Writing
$K_{A}(t)=1+\sum_{\alpha\in\cal{P}_{A}}C_{\alpha\alpha}(t)$, we immediately
find $\overline{K_{A}(t)}=1+(D_{A}/D)^{2}(t-1)$ and $\overline{K(t)}=t$ for
$1\leq t\ll D$. All of these ramps therefore have their origin in operator
paths satisfying the selection rule in Eq. (63), and its generalization which
is implicit in Eq. (65).

####  V.1.2 Paths in the RPM

We can adapt aspects of this analysis to the RPM in the $q\to\infty$ limit,
although first some general comments are in order. Most importantly, the
operator paths that can have nonzero amplitudes are severely restricted by
locality. For example, in a system with local interactions, and over a single
time step, an operator cannot change the distance between its left- and right-
hand end-points by more than a number of sites of order unity. Subject to this
constraint, we must ask which kinds of paths necessarily make real, positive
contributions to autcorrelations functions (for all but a measure zero set of
disorder realisations).

Remarkably, in the RPM at $q\to\infty$, the paths which dominate
autocorrelation functions are ‘localised’. By this we mean that the operators
$P_{\alpha(r)}$ in the path all have the same support. In fact, this can be
seen already from our calculations of averaged autocorrelation functions:
$\overline{C_{\alpha\alpha}(t)}$ is independent of the system size at large
$q$, and depends only on the support of the operator $P_{\alpha}$. The
mechanism by which this arises is quite general, and later we will show that a
related feature is present in the dynamics of qubit chains.

It is simplest to consider single-site operators, and so here we will restrict
ourselves to this example. For a closed operator path
$P_{\alpha(0)},\ldots,P_{\alpha(t-1)}$ in which all $P_{\alpha(r)}$ are
nontrivial on just a single site, we can write $C_{\alpha(r+1)\alpha(r)}$ in
Eq. (60) as

$\displaystyle
C_{\alpha(r+1)\alpha(r)}=q^{-3}[P_{\alpha(r+1)}]_{b(r)a(r)}[P_{\alpha(r)}]_{a^{\prime}(r)b^{\prime}(r)}$
(67) $\displaystyle\times
V_{a(r)a^{\prime}(r)}V^{*}_{b(r)b^{\prime}(r)}\sum_{c,c^{\prime}=0}^{q-1}e^{i[\varphi_{ca(r)}-\varphi_{cb(r)}+\varphi^{\prime}_{a(r)c^{\prime}}-\varphi_{b(r)c^{\prime}}]}$

where $V$ is a Haar-random unitary acting at the same site as the single-site
operators $P_{\alpha(r)}$, and the diagonal matrices $\varphi$ and $\varphi$,
whose diagonal entries we denote $\varphi_{ca(r)}$ and
$\varphi^{\prime}_{a(r)c^{\prime}}$, describe interactions between this site
and its left- and right-hand neighbours, respectively. Here
$a(r),a^{\prime}(r),b(r),b^{\prime}(r)=0\ldots(q-1)$ are state indices for the
single site at which the operators $P_{\alpha(r)}$ act, while $c$ and
$c^{\prime}$ are state indices for this site’s neighbours. The overall
prefactor $q^{-3}$ has arisen from contracting the circuit away from the
three-site region centered on our operator path.

At large $q$, observe that we have $a(r)\neq b(r)$ and $b_{r}\neq b_{r}^{*}$
for the vast majority of terms. Summing over $c$ and $c^{\prime}$, Eq. (67)
then reduces to

$\displaystyle C_{\alpha_{r+1}\alpha_{r}}=q^{-1}e^{-2\varepsilon}$
$\displaystyle[P_{\alpha(r+1)}]_{b(r)a(r)}[P_{\alpha(r)}]_{a^{\prime}(r)b^{\prime}(r)}$
$\displaystyle\times V_{a(r)a^{\prime}(r)}V^{*}_{b(r)b^{\prime}(r)}+\ldots$
(68)

Restricting to localised operator paths has at this stage left us with a
single-site problem. The factors $e^{-2\epsilon}$ at each time step describe
the dephasing induced by interactions with neighbouring sites. We can now
follow essentially the same line of reasoning as in the previous subsection.

That is, we first insert Eq. (68) into the expression Eq. (3) and ask about
contributions from terms where the single-site indices are paired cyclically
as $b(r)=a(r+s)$ and $b^{\prime}(r)=a^{\prime}(r+s)$. For these terms we find
contributions $\prod_{r=0}^{t-1}|V_{a(r)a^{\prime}(r)}|^{2}$ and, mirroring
Eq. (61), we will separate out non-negative parts

$\displaystyle|V_{a(r)a^{\prime}(r)}|^{2}=q^{-1}+\ldots.$ (69)

Summing over state indices then generates a non-negative contribution to the
operator path. For prime $t$ this is

$\displaystyle C_{\alpha(0)\alpha_{1}}\ldots C_{\alpha(t-1)\alpha(0)}$ (70)
$\displaystyle=q^{-2t}e^{-2\varepsilon
t}\sum_{s=1}^{t-1}|\mathrm{Tr}(P_{\alpha(0)}P_{\alpha_{s}}\cdots
P_{\alpha([t-1]s)})|^{2}+...$

We therefore have a single-site analogue of the selection rule in Eq. (63),
and the contributions displayed in Eq. (62) will turn out to dominate the
averaged autocorrelation functions of local operators at large $q$.

If $P_{\alpha(0)}\ldots P_{\alpha(t-1)}=\pm I$, we find that the expression
displayed above is the dominant contribution to the averaged amplitude

$\displaystyle\overline{C_{\alpha(0)\alpha(1)}\ldots
C_{\alpha(t-1)\alpha(0)}}=q^{-2(t-1)}e^{-2\varepsilon t}(t-1)$ (71)

at large $q$. This can be generalized to non-prime $t$ along identical lines
to the previous section. Again, we find a ramp at the level of an individual
path’s amplitude, although in contrast to the previous subsection we have here
restricted to an extremely small subset of the possible operator paths. The
locality of the interactions compensates for this, and causes these paths to
have relatively large amplitudes. Summing over paths (for prime or non-prime
$t$) we pick up an overall factor $q^{2(t-2)}$ which leads us to

$\displaystyle\overline{C_{\alpha(0)\alpha(0)}(t)}=q^{-2}e^{-2\varepsilon
t}(t-1).$ (72)

Remarkably, considering only the localised paths only has given us the
dominant contribution to the averaged autocorrelation function of the RPM.

The above approach can be generalized to the autocorrelation functions of more
complicated operators, and one then finds that all operators $P_{\alpha(r)}$
in the contributing paths $P_{\alpha(0)},\ldots,P_{\alpha(t-1)}$ have the same
support. It is important to note, however, that our large-$q$ analysis does
not expose the late-time ramp in $C_{\alpha(0)\alpha(0)}(t)$, although an
approximate description of this feature can be formulated along similar lines
to Sec. III.5.

Equation (72) is the central result of this section: the operator paths which
contribute to ensemble-averaged autocorrelation functions of local operators
in the RPM are themselves localised. This feature arises from the combination
of (i) local interactions and (ii) the correlations between amplitudes of
different steps of the operator path which are imposed by time-translation
symmetry.

###  V.2 Squared amplitudes

Here we discuss the operator paths contributing to squared correlation
functions, and to the OTOC. First note that in this limit the OTOC is a
weighted sum over squared correlation functions; starting from Eq. (7) it can
be verified that in the RPM at large $q$,

$\displaystyle\overline{\mathcal{C}(x,t)}=1-\sum_{\alpha\in\mathcal{P}}\langle
P_{\alpha}\mathcal{O}(x,0)P_{\alpha}\mathcal{O}(x,0)\rangle\overline{C^{2}_{\mathcal{O}\alpha}(t)}.$
(73)

We can therefore restrict our attention to averages of squared correlation
functions.

The key observation is that, even with time-translation symmetry, in the
large-$q$ limit this does not play a role in the average. To see this let us
first write

$\displaystyle
C_{\alpha\beta}(t)=\sum_{\gamma\in\mathcal{P}}C_{\alpha\gamma}(t^{\prime})C_{\gamma\beta}(t-t^{\prime}),$
(74)

for $t^{\prime}\neq 0,t$, and then square this expression. The square involves
a sum over two sets of operators $\gamma,\gamma^{\prime}\in\mathcal{P}$, and
so we have both diagonal $\gamma=\gamma^{\prime}$ and off-diagonal
$\gamma\neq\gamma^{\prime}$ contributions. The sum over diagonal contributions
factorises at large $q$,

$\displaystyle\sum_{\gamma\in\mathcal{P}}\overline{[C_{\alpha\gamma}(t^{\prime})]^{2}[C_{\gamma\beta}(t-t^{\prime})]^{2}}$
(75)
$\displaystyle=\sum_{\gamma\in\mathcal{P}}\overline{[C_{\alpha\gamma}(t^{\prime})]^{2}}\,\,\overline{[C_{\gamma\beta}(t-t^{\prime})]^{2}}.$

It can then be shown that the off-diagonal contribution
$\gamma\neq\gamma^{\prime}$ is subleading relative to the above diagonal
contribution above. Therefore,

$\displaystyle\overline{[C_{\alpha\beta}(t)]^{2}}=\sum_{\gamma\in\mathcal{P}}\overline{[C_{\alpha\gamma}(t^{\prime})]^{2}}\,\,\overline{[C_{\gamma\beta}(t-t^{\prime})]^{2}}.$
(76)

This decomposition can be repeated until we are left with a product of average
squared correlation functions for individual time steps,

$\overline{C^{2}_{\alpha\beta}(t)}=\sum_{\alpha(1),\cdots,\alpha(t-1)\in\mathcal{P}}\overline{C^{2}_{\alpha\alpha(1)}}\,\overline{C^{2}_{\alpha(1)\alpha(2)}}\cdots\overline{C^{2}_{\alpha(t-1)\beta}}.$
(77)

The fact that the averages on the right-hand side of Eq. (77) can be performed
independently for each time step shows that averages of squared
autocorrelation functions, as well as of the OTOC, are insensitive to time-
translation symmetry at large $q$. We can therefore understand their behaviour
from Ref. [31] and Refs. [4, 5], respectively. For example, the average
squared autocorrelation function $\overline{C^{2}_{\alpha\alpha}(t)}$ is
dominated by operator paths whose extent grows as $t^{1/2}$ [31], to be
contrasted with the localised operator paths contributing to
$\overline{C_{\alpha\alpha}(t)}$ which we have discussed above.

###  V.3 Numerical results

In this section we have so far shown that, at large $q$ in the RPM, the kinds
of operator paths contributing to $\overline{C_{\alpha\alpha}(t)}$ and
$\overline{C^{2}_{\alpha\alpha}(t)}$ are very different. We now demonstrate
numerically that this difference carries over to the behavior of a Floquet
system with $q=2$.

The first object we consider is a deconstructed sum over autocorrelation
functions, and the second is a sum over squared correlation functions. To
motivate the first object, let us write an autocorrelation function as

$\displaystyle
C_{\alpha\alpha}(t)=\sum_{\beta\in\mathcal{P}}C_{\alpha\beta}(\lceil
t/2\rceil)C_{\beta\alpha}(\lfloor t/2\rfloor),$ (78)

which is a sum over all paths from $P_{\alpha}$ to $P_{\beta}$ over the first
half of the interval $t$, and over all paths back from $P_{\beta}$ to
$P_{\alpha}$ in the second half. For odd $t$ we split the intervals as $\lceil
t/2\rceil=(t+1)/2$ and $\lfloor t/2\rfloor=(t-1)/2$. Let us define the
contribution to $C_{\alpha\alpha}(t)$ from operator paths that, halfway
through the evolution, have support within a region $B$,

$\displaystyle
C^{B}_{\alpha\alpha}(t)=\sum_{\beta\in\mathcal{P}_{B}}C_{\alpha\beta}(\lceil
t/2\rceil)C_{\beta\alpha}(\lfloor t/2\rfloor).$ (79)

If we now define a sequence of regions $B_{0},B_{1},\ldots$, with $B_{j}$
strictly contained within $B_{j+1}$ we can write the contribution to
$C_{\alpha\alpha}(t)$ from operator paths that have support outside of $B_{j}$
but not outside of $B_{j+1}$ as
$C^{B_{j+1}}_{\alpha\alpha}(t)-C^{B_{j}}_{\alpha\alpha}(t)$. To make contact
with our discussions of the PSFF, it will be convenient to sum over
autocorrelation functions for sets of operators $P_{\alpha}$ that are
nontrivial within a region $A$. In analogy with $C^{B}_{\alpha\alpha}(t)$ we
define

$\displaystyle
K^{B}_{A}(t)=\sum_{\alpha\in\mathcal{P}_{A}}\sum_{\beta\in\mathcal{P}_{B}}C_{\alpha\beta}(\lceil
t/2\rceil)C_{\beta\alpha}(\lfloor t/2\rfloor),$ (80)

which is the sum over all amplitudes of operators paths that start in $A$,
have support entirely within $B$ at the halfway point, and then return to $A$.
Having defined a sequence of regions $B_{j}$ as above, with $B_{0}$ a region
of zero sites for which $K_{A}^{B_{0}}(t)=0$, and with the largest region
corresponding to entire system, we can write the PSFF as

$\displaystyle K_{A}(t)=\sum_{j=0}[K_{A}^{B_{j}}(t)-K_{A}^{B_{j-1}}(t)].$ (81)

At large $q$ in the RPM we have $K_{A}(t)=K_{A}^{A}(t)$, i.e. all operator
paths contributing to $K_{A}(t)$ are confined to the region $A$.

The second object we consider measures operator weight via squared correlation
functions. To construct this weight we write

$\displaystyle P_{\alpha}(\lceil
t/2\rceil)=\sum_{\beta\in\mathcal{P}}C_{\alpha\beta}(\lceil
t/2\rceil)P_{\beta},$ (82)

and observe that $D^{-1}\text{Tr}P^{2}_{\alpha}(\lceil t/2\rceil)=1$, which
corresponds to the normalization condition on the sum over squared correlation
functions $\sum_{\beta}C_{\alpha\beta}^{2}(\lceil t/2\rceil)=1$. We can then
ask how much of this weight comes from operators that have support only within
$B$, i.e. $\sum_{\beta\in\mathcal{P}_{B}}C_{\alpha\beta}^{2}(\lceil
t/2\rceil)$. To compare with $K^{B}_{A}(t)$, we will sum this object over all
operators within the region $A$,

$\displaystyle
N^{B}_{A}(t)=\sum_{\alpha\in\mathcal{P}_{A}}\sum_{\beta\in\mathcal{P}_{B}}C_{\alpha\beta}^{2}(\lceil
t/2\rceil).$ (83)

This is the total weight of the operator strings which evolve from $A$ to $B$
over the interval $\lceil t/2\rceil$, as opposed to the total amplitude of
operator strings which follow this route and ultimately return to $A$ at time
$t$. In analogy with Eq. (81) we can write the total weight of operators which
start in $A$ as

$\displaystyle N_{A}(t)=\sum_{j=0}[N_{A}^{B_{j}}(t)-N_{A}^{B_{j-1}}(t)].$ (84)

Since we have a normalisation condition $N_{A}(t)=D_{A}^{2}-1$, it will be
convenient to study the objects
$[N_{A}^{B_{j}}(t)-N_{A}^{B_{j-1}}(t)]/(D_{A}^{2}-1)$ below. The sum of these
objects over regions $B_{j}$ is then unity. Our analysis of the RPM at large
$q$ (along with earlier studies of the OTOC in RUCs) shows that the value of
$B$ for which $N^{B}_{A}(t)$ is maximized should grow ballistically with time.

Figure 16: Contributions $\overline{K^{B}_{A}(t)}$ to the average PSFF from
operator paths of varying support. We consider regions $B_{j}$ of
$L_{B_{j}}=2j-1$ sites for $j=1,2,3,4$ centred on the region $A$ of $L_{A}$
sites (see text in centres of panels). The data shown is
$K^{B_{j}}_{A}(t)-K^{B_{j-1}}_{A}(t)$, corresponding to operator paths which
at the halfway point have support outside of $B_{j-1}$ but within $B_{j}$.
Here $L=12$. Figure 17: Average contributions
$(D_{A}^{2}-1)^{-1}\overline{N^{B}_{A}(t)}$ to the operator weight for
operators which are initially within $A$ and at time $\lceil t/2\rceil$ are in
$B$. As in Fig. 16 we calculate differences between the contributions from
different regions $B_{j}$, with $B_{j}$ strictly within $B_{j+1}$. Here
$L=12$.

In Fig. 16 we calculate $K^{B_{j}}_{A}(t)-K^{B_{j-1}}_{A}(t)$ numerically. For
each of the values of $L_{A}$ shown it is clear that, at early times, the
largest contribution to $K_{A}(t)$ is from operators which have support $\sim
L_{A}$ at time $\lceil t/2\rceil$,, which is the behavior anticipated from the
RPM at large $q$. This behavior is to be contrasted with
$N^{B_{j}}_{A}(t)-N^{B_{j-1}}(t)$: in Fig. 17 we show that, by this measure,
the operators $P_{\alpha}$ develop weight on ever longer Pauli strings as time
progresses. At late times the difference between the structure of operators
contributing to autocorrelation functions, and to squared autocorrelation
functions, is particularly striking: only in the latter case is the operator
weight is in the longest operator strings (see $L_{B_{j}}=12$, $L_{B_{j-1}}=9$
data in Fig. 17).

##  VI Individual systems

Above we have studied ensemble-averaged dynamics, and so it is natural to ask
how much of the above phenomenology carries over the behaviour of individual
Floquet many-body systems. In the following we show first that sample-to-
sample fluctuations of autocorrelation functions obscure their second peak. We
then show that this structure is instead revealed by considering averages over
sets of autocorrelation functions such as the PSFF. In particular, we will
show that the sample-to-sample fluctuations of $K_{A}(t)$ are parametrically
smaller than its late-time tail for $1\ll L_{A}\ll L$.

First consider the ensemble variance
$\overline{C^{2}_{\alpha\alpha}(t)}-\overline{C_{\alpha\alpha}(t)}^{2}$ of the
autocorrelation function in the large-$q$ limit of the RPM. For $a$-site
single-cluster operators, it can be verified that at late times
$\overline{C^{2}_{\alpha\alpha}(t)}\simeq q^{-2a}(t-1)e^{-4\varepsilon t}$,
while $\overline{C_{\alpha\alpha}(t)}\simeq q^{-2a}(t-1)e^{-2\varepsilon t}$.
The sample-to-sample fluctuations are therefore of order $q^{-a}$, and so are
large compared with the ensemble average at large $q$. One might consider a
time average, but since for single-site operators the autocorrelation function
has nontrivial structure only over a time interval $1\leq t\lesssim t_{\rm
Th,\alpha}$, for an individual Floquet system there are too few time steps to
suppress the fluctuations. Although at finite $q$ the ratio of the variance to
the squared mean is of course finite, a mild ensemble average is nevertheless
required to observe the second peak in autocorrelation functions; in Fig. 18a
we illustrate the effects of averaging over $10$ random Floquet systems.

(a)(b)

Figure 18: (a) Autocorrelation functions of operators consisting of a single
cluster ($n=1$) for various $a$ (legend); the coloured regions are centred on
the average $\overline{C_{\alpha\alpha}(t)}$ and have vertical width equal to
$10^{-1/2}$ times the root-mean-square fluctuation of $C_{\alpha\alpha}(t)$;
they represent the effect of averaging over just $10$ random Floquet circuits.
(b) Suppressing fluctuations through an average over operators in individual
Floquet systems. The coloured regions are centred on $\overline{K_{A}(t)}$ and
have vertical width equal to the root-mean-square fluctuation of $K_{A}(t)$.
The legend shows $L_{A}$. Here $L=10$.

Instead we will sum $C_{\alpha\alpha}(t)$ over operators $P_{\alpha}$
supported within a region $A$, i.e. we consider the PSFF. To be concrete, we
will calculate the second moment $\overline{K_{A}^{2}(t)}$ in the RPM at large
$q$, and as above we use a transfer matrix technique. In calculating on-site
Haar averages, denote by $a(r)$ and $b(r)$ indices appearing in the first
factor of $K_{A}(t)$, as in Sec. III, and by $c(r)$ and $d(r)$ indices in the
second factor. Note that $a(r)$ and $c(r)$ appear in products of Floquet
operators $W^{t}$, while $b(r)$ and $d(r)$ appear in the conjugates. The index
$r=0,\ldots,t$, and for sites in $A$ we have the trace structure $a(0)=a(t)$,
whereas for sites in $\bar{A}$ we have $a(0)=b(0)$, $a(t)=b(t)$, $c(0)=d(0)$,
and $c(t)=d(t)$. We only need to account for cyclic Gaussian pairings, but
there are now two sets of such pairings in region $A$. The first set, which we
denote $p=0$, corresponds to $a(r)=b(r+s)$ and $c(r)=d(r+\tilde{s})$ for
$s,\tilde{s}=0,\ldots,(t-1)$. The second, which we denote $p=1$, corresponds
to $a(r)=c(r+s)$, $b(r)=d(r+s^{\prime})$. In region $\bar{A}$, the leading
order pairings are simply $a(r)=b(r)$ and $c(r)=d(r)$.

We identify each pairing $(p,s,\tilde{s})$ with one element of a set of
$2t^{2}$ orthonormal basis vectors, and define a transfer matrix $T^{(2)}$
with elements

$\displaystyle T^{(2)}_{ps\tilde{s},p^{\prime}s^{\prime}\tilde{s}^{\prime}}$
$\displaystyle=\delta_{pp^{\prime}}T_{ss^{\prime}}T_{\tilde{s}\tilde{s}^{\prime}}+[1-\delta_{pp^{\prime}}]e^{-2\varepsilon
t}.$ (85)

Here $T$ is the transfer matrix which generates the average SFF [13], as well
as the PSFF above. We note that $T^{(2)}$ also generates the second moment of
the SFF [20, 45]. In terms of $T^{(2)}$ we have the ensemble variance

$\displaystyle\overline{K_{A}^{2}(t)}-\overline{K_{A}(t)}^{2}=[(T^{(2)})^{L_{A}+1}]_{000,000}-[T^{L_{A}+1}]^{2}_{00},$
(86)

where the indices $p,s,\tilde{s}=0,0,0$ correspond to the boundary condition
in region $\bar{A}$. At late times the variance goes to zero, and it decays to
this value as

$\displaystyle\overline{K_{A}^{2}(t)}-\overline{K_{A}(t)}^{2}=L_{A}(L_{A}+1)t^{2}e^{-4\varepsilon
t}+\ldots,$ (87)

which corresponds to a domain with $p\neq 0$ and any
$s,\tilde{s}=0,\ldots,(t-1)$. For $L_{A}\gg 1$, the scale of the fluctuations
around the average are therefore of order $L_{A}te^{-2\varepsilon t}$.

We should compare the fluctuations with the average PSFF itself, and at late
times we have shown that $\overline{K_{A}(t)}\simeq
1+L_{A}(L_{A}+1)(t-1)e^{-2\varepsilon t}/2$. For large $L_{A}$, the sample-to-
sample fluctuations in $K_{A}(t)$ are therefore parametrically smaller than
its late-time tail, with the signal-to-noise ratio proportional to $L_{A}$ for
$L_{A}\gg 1$. The implication is that structure in the PSFF arising from
locality can be observed in individual systems.

In Fig. 18b  we verify this numerically for the $q=2$ RSM. There we show the
root-mean-square fluctuations of $K_{A}(t)$ compared with
$\overline{K_{A}(t)}$, and observe a clear separation between the different
values of $L_{A}$ displayed.

##  VII Discussion

There are fundamental differences between the dynamics of operators as viewed
through correlation functions and OTOCs. This is clear even in RUCs having
gates which are Haar-random in both space and time: there the average of the
OTOC can be understood through the hydrodynamic picture of operator spreading
[4, 5], while averages of correlation functions vanish.

In this work we have studied the effects of time-translation invariance (i.e.
the effect of having a fixed time-evolution operator) on these objects. To
simplify the analysis, we focused on Floquet quantum circuits having (i) Haar-
random local unitary operations and (ii) no conservation laws. A key result
has been to show that, in locally-interacting spin chains, time-translation
symmetry of the evolution operator leads to observable and nontrivial features
in autocorrelation functions. In contrast, we demonstrated that the structure
of the OTOC is insensitive to time-translation invariance, at least at large
$q$ in the RPM.

This difference can be understood by representing correlation functions and
OTOCs in terms of sums over paths in the space of operators. While correlation
functions can be expressed as sums over amplitudes of operator paths [see Eq.
(3)], OTOCs involve sums over (non-negative) ‘probabilities’ [see Eq. (4)]. In
RUCs the amplitudes of operator paths have zero mean, and so averaged
correlation functions vanish. Introducing time-translation symmetry generates
correlations between the amplitudes of different segments of operator paths,
and these give rise to nonzero averages of autocorrelation functions. In the
local systems we have considered, this perspective has also revealed a
striking difference between the morphologies of operator paths which
contribute to ensemble-averaged autocorrelation functions relative to those
dominating the OTOC. In particular, we find that the operator paths
contributing to autocorrelation functions are much narrower in space; see
Figs. 2a and 2b. We note that the morphologies of operators paths contributing
to correlation functions were also recently discussed in Ref. [31].

A second result has been to show how the spectral statistics of many-body
Floquet systems can be understood in terms of the dynamics of operators. The
SFF can be expressed as a sum over autocorrelation functions of all operators
[18], and a local analogue (the PSFF) has been defined as the restriction of
this sum to operators with support only within a subregion [12]. In brief, at
late times autocorrelation functions match their behavior in RMF systems, and
their sum reproduces the famous ramp in the SFF. At the earliest times
$t<t_{\mathrm{d}}$ autocorrelation functions relax, and strikingly this
relaxation is followed by an increase for $t>t_{\mathrm{d}}$, which can be
understood as arising from local unitary scrambling in the region of the
operator’s support (a residue of a ‘local ramp’). This local unitary
scrambling competes with dephasing induced by interactions with neighbouring
regions, and the result of this competition is a finite-time peak in the
autocorrelation function. The decay following the peak ultimately gives way to
late-time ramp and plateau, whose amplitudes are exponentially small in the
system size.

These findings were supported by analytical as well as numerical calculations.
Our analytical investigations were centered on the RPM at large $q$ [13], as
well as a spinful variant of it which we have introduced (see Appendix A). In
these models we calculated OTOCs, autocorrelation functions, and PSFFs. We
also provided a heuristic argument to understand the finite-$q$ physics of
PSFFs, although evidence for the generality of our analytic results in
describing finite-$q$ systems can already be seen in Ref. [20], which studied
spectral statistics numerically at $q=2$. Here we have demonstrated that
calculations in the large-$q$ RPM also provide an accurate depiction of
autocorrelation functions at $q=2$ (focusing on the ergodic phase of the model
in Ref. [14]).

Some aspects of autocorrelation functions remain to be understood. The local
diagonal approximation has provided an approximate description of the
emergence of the ramp $(t-1)/(q^{2L}-1)$ expected from RMT. However, our
numerical calculations indicate that the gradient of the ramp may depend on
the structure of the operator, and this feature is absent in RMT.

It is worth commenting on the difference between our results and those which
would be found in ergodic dual unitary circuits (having brickwork geometry).
In that setting autocorrelation functions of site-local operators are exactly
zero for finite times $t>0$ in the limit of large $L$ [46, 23, 47], while in
this regime the SFF matches the RMT result [10, 22]. Representing the SFF as a
sum over autocorrelation functions, these two properties are consistent with
the ‘perfect dephaser’ property of dual unitary circuits [40]; the effective
dephasing time for the autocorrelation function of a local operator is zero
(i.e. it is as if the effective domain wall tension $\varepsilon\to\infty$).
However, autocorrelation functions of multi-site operators in dual-unitary
circuits do not necessarily vanish, and in fact some of them must be non-zero
in order to generate the anticipated behaviour of the SFF from their sum. It
would be interesting to understand how this arises.

Moving away from the ergodic phase of Floquet circuits, it is helpful to
consider first the extreme limit of decoupled qubits evolving under local
Haar-random fields (i.e. $q=2$ and $\varepsilon=0$ in the RPM, or $J=0$ in
RSM). There it is straightforward to show that
$\overline{C_{\alpha\alpha}(t\geq 2)}=(1/3)^{a}$, with $a$ is the number of
non-identity Pauli operators in the string $P_{\alpha}$, and
$\overline{K(t\geq 2)}=2^{L}$, corresponding to Poisson spectral statistics.
Weakly-coupled random spin chains are expected to be many-body localised [48,
49], and in that setting corrections to the above behaviour of autocorrelation
functions can be understood in the energy domain through the theory of many-
body resonances [50, 51, 52, 53, 54]; the connection to spectral statistics
has been discussed in Ref. [52]. Moving to the time domain, it would be
interesting to understand how that picture manifests itself in the structure
of the operator paths contributing to correlation functions.

Discussions of operator paths are relevant to questions about the complexity
of simulating quantum dynamics on classical computers (and hence also for
discussions of quantum advantage [55, 56]). For example, in a one-dimensional
system of qubits, a generic operator having support on $\ell$ contiguous sites
can be represented using matrix-product operator bond dimension $\chi\leq
2^{\ell}$. When studying an infinite-temperature autocorrelation function of a
local operator, it is therefore important to understand how rapidly it grows
in time. One might expect that the operator paths relevant to autocorrelation
functions have $\ell\simeq 2v_{B}t$ at their widest, although the results of
Ref. [31] suggest that the dominant paths may be confined to a parametrically
smaller region of space with $\ell\sim t^{1/2}$. Our work suggests that the
contributions to evolving operators which survive the ensemble average are
confined to a region with $\ell\sim t^{0}$.

In closing, we note that the operator paths relevant to the simulation of
diffusion have recently been investigated in [33, 32]. Building on these
ideas, as well as previous work on the spectral statistics of Floquet systems
with conservation laws [41, 42, 43], it would be interesting to understand the
connection between dynamics and spectra in Hamiltonian systems. This would
provide an alternative to ideas based on interference effects [57, 13, 20,
58].

##  VIII Acknowledgements

We thank Remy Dubertrand, Michele Fava, Max McGinley, and Adam Nahum for
useful discussions. T.Y. acknowledges hospitality at University of Tokyo where
part of the manuscript was written. S.J.G. was supported by the Gordon and
Betty Moore Foundation and by the U.S. Department of Energy, Office of
Science, Office of High Energy Physics, under QuantISED Award DE-SC0019380.
J.T.C. was supported in part by EPSRC Grant EP/S020527/1.

##  Appendix A Spinful RPM

The RPM as defined in Sec. II.2 has the property that
$\overline{C_{\alpha\alpha}(t=1)}=0$ for all operators $P_{\alpha}$. This
feature is shared with other models in which the ensemble of Floquet operators
is invariant under independent unitary rotations at each site, but is not
expected to be generic. For that reason it is interesting to construct a
solvable model without such an invariance. We do so here by introducing into
the RPM defined previously an extra degree of freedom, which has dynamics that
does not consist simply of Haar-distributed rotations. We call the model with
this extra degree of freedom the spinful RPM.

We define the spinful RPM by extending the local Hilbert space to have
dimension $2q$, identifying this space as the tensor product of the original
$q$-dimensional space (the ‘colour’ degrees of freedom) $\mathbb{C}^{q}$ and
the $2$-dimensional space of a ‘spin’ degree of freedom $\mathbb{C}^{2}$. The
single-site unitary operations (now $2q\times 2q$ matrices) mix the two spin
states on a timescale set by a parameter $\theta$ (when $\theta=0$ they reduce
to two independent copies of $q\times q$ Haar-distributed unitary matrices),
and this mixing introduces new features in autocorrelation functions and the
OTOC. To be more precise, the on-site unitaries consisting of the first layer
of the Floquet operator $W_{1}=\bigotimes_{x=1}^{L}U_{x}$ are promoted to
$2q\times 2q$ matrices

$U_{x}=(R\otimes I_{q})\begin{pmatrix}V_{1,x}&0\\\ 0&V_{2,x}\end{pmatrix}$
(88)

with

$R=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}$ (89)

where $\theta\in[0,\pi/4)$ and $V_{1,x},V_{2,x}$ are independent $q\times q$
matrices that are Haar random. On the other hand, even in the spinful RPM the
second year $W_{2}$ still acts on the computational basis
$\mathbb{C}^{q}\otimes\mathbb{C}^{2}$ diagonally.

###  A.1 PSFF of the spinful RPM

The PSFF for the spinful RPM can be computed in the same way as shown in Sec.
III.1 for the spinless case. To compute Haar averages diagrammatically, let us
first focus on a single site $x$ in $A$. The on-site diagram consists of two
orbits, one a product of $U_{x}$ and the other a product of $U_{x}^{*}$. We
need to evaluate

$\displaystyle\Big{<}[U_{x}]_{\mathtt{s}(0),a(0)}^{\mathtt{s}(1),a(1)}[U_{x}]_{\mathtt{s}(1),a(1)}^{\mathtt{s}(2),a(2)}\cdots[U_{x}]_{\mathtt{s}(t-1),a(t-1)}^{\mathtt{s}(0),a(0)}$
$\displaystyle\times[U^{*}_{x}]_{\mathtt{r}(0),b(0)}^{\mathtt{r}(1),b(1)}[U^{*}_{x}]_{\mathtt{r}(1),b(1)}^{\mathtt{r}(2),b(2)}\cdots[U^{*}_{x}]_{\mathtt{r}(t-1),b(t-1)}^{\mathtt{r}(0),b(0)}\Big{>},$
(90)

where indices $\mathtt{s}(u),\mathtt{r}(u)=1,2$ are associated with the spin
degree of freedom, while $a(u),b(u)=1,\cdots,q$ are for the colour degree of
freedom for $u=0,\cdots,t-1$. It turns out that the leading pairings upon Haar
averaging are again $t$ cyclic pairings, which can be parameterised by
$\underline{a}(u)=\underline{b}(u+s)$ for fixed $u$ and $s=0,\cdots,t-1$,
where $\underline{a}=(\mathtt{s},a)$ (see Fig. 6a-6c in the case of $t=3$).
Note that the parameterization for the spin degree of freedom
$\mathtt{s}(u)=\mathtt{r}(u+s)$ comes from the fact that $V_{1}$ and $V_{2}$
are independent unitaries. This means that the phase factor due to having
different local pairings on neighboring sites is still the same in this case.
Finally, we also have to account for the contribution from the spin degree of
freedom, which happens to be decoupled from the rest. Namely, for each site
$x$ on $A$ we have

$\displaystyle\prod_{u=1}^{t}\sum_{\mathtt{s}_{x}(u)=1,2}R_{\mathtt{s}_{x}(u),\mathtt{s}_{x}(u+1)}R_{\mathtt{s}_{x}(u+s_{x}),\mathtt{s}_{x}(u+1+s_{x})}$
$\displaystyle=\mathrm{Tr}\,\tilde{R}^{t}=1+\cos^{t}2\theta,$ (91)

where $\tilde{R}_{\mathtt{s}\mathtt{r}}=(R_{\mathtt{s}\mathtt{r}})^{2}$ for
any $s_{x}$. Thus $K_{A}(t)$ can be recast into a product of two partition
functions as explained in Sec. II.4 and we get

$K_{A}(t)=K^{\mathrm{s}}_{A}(t)K^{\mathrm{c}}_{A}(t),$ (92)

where $K^{\mathrm{s}}_{A}(t)=(1+\cos^{t}2\theta)^{L_{A}}$ and

$K^{\mathrm{c}}_{A}(t)=\frac{1}{t}\left[\lambda_{0}^{L_{A}+1}+(t-1)\lambda^{L_{A}+1}\right]$
(93)

with $\lambda_{0}=1+(t-1)e^{-\varepsilon t}$ and $\lambda=1-e^{-\varepsilon
t}$.

Two comments are in order. First, the PSFF no longer drops to $1$ at
$t_{\mathrm{d}}=1$ but rather takes a nontrivial value
$\overline{K_{A}(1)}=(1+\cos 2\theta)^{L_{A}}$. We thus observe that the
additional spin degree freedom in the model introduces nontrivial
$\overline{K_{A}(t_{\mathrm{d}}=1)}\neq 1$. Another observation is that, in
the spiful RPM, there are two sources that influence the Thouless time
$t_{\mathrm{Th},A}$ of the PSFF. The first one is the domain-walls that are
encoded in $K^{\mathrm{c}}_{A}(t)$ and another one is the extra on-site spin
structure accounted for by $K^{\mathrm{s}}_{A}(t)$. Since
$K^{\mathrm{s}}_{A}(t)$ decays over the timescale $-\ln L_{A}/\log\cos
2\theta$, which again grows with $\ln L_{A}$, the Thouless time of the PSFF of
the spinful RPM is given by

$t_{\mathrm{Th},A}=\max\\{\ln L_{A}/\varepsilon,-\ln L_{A}/\ln\cos 2\theta\\}$
(94)

###  A.2 Autocorrelation functions of the spinful RPM

As in Sec. II.3, let us start with accounting for the spin degree of freedom
by parameterising an on-site operator on site $x$ as
$\tilde{P}_{\mu_{x},\alpha_{x}}=\sigma_{\mu_{x}}\otimes P_{\alpha_{x}}$ where
$\sigma_{\mu}$ are Pauli matrices
($\sigma_{0}=I_{2},\sigma_{1}=X,\sigma_{2}=Y,\sigma_{3}=Z$) and normalised as
$\mathrm{Tr}(\sigma_{\mu}\sigma_{\nu})=2\delta_{\mu\nu}$. The on-site non-
identity operators $\tilde{P}_{\mu_{x},\alpha_{x}}$ can be classified into two
classes, which we call class I and class II, according to whether they act
nontrivially on the colour degree of freedom (class I) or not (class II). It
turns out that the only on-site operator in class II that produces a
nontrivial contribution is $Z\otimes I_{q}$, as $X$ and $Y$ are nondiagonal.
Let us consider the operator string $\tilde{P}_{\mu,\alpha}$ that consists of
$n$ clusters of class I operators, each of which contains $a_{m}$ class I on-
site operators for $m=1,\cdots,n$. Every pair of neighboring clusters is
connected by a string of on-site operators that are either class II or the
identity. We denote the total number of class II operators in the operator
string $\tilde{P}_{\mu,\alpha}$ by $b$. This means that the total number of
single-site identity operators is $L-b-a$ where $a=\sum_{m=1}^{n}a_{m}$, and
the weight of $\tilde{P}_{\mu,\alpha}$ is $a+b$. Further, let us denote the
total number of class I operators with $\sigma_{\mu}$ in the entire operator
string by $\mathtt{a}_{\mu}$. See Fig. 19 for an example of an operator
string.

Figure 19: An operator string where $L=11$, $n=2$, $b=2$, and $a=5$ (with
$\mathtt{a}_{0}=\mathtt{a}_{1}=\mathtt{a}_{3}=1$ and $\mathtt{a}_{2}=2$). The
letter in each on-site operator corresponds to the spin component of the
operator $\sigma_{\mu}$, and those in red colour act nontrivially on the
colour degree of freedom. The weight of the operator string is $a+b=7$.

Next we compute the Haar averages of single-site diagrams at $q\to\infty$.
There are three types of single-site diagrams: those with class I operators,
class II operators, and diagrams without local operators. Let us start with
the first type. We can repeat the same argument as in the RPM without the spin
degree of freedom, and observe that the leading pairings are again given by
the cyclic ones with the factor $q^{-1}$ except the one associated with the
identity permutation. The extra spin space, however, generates a factor
$f_{\mu}(\theta)$, which turns out to be different depending on which Pauli
matrix $\sigma_{\mu}$ is used in the on-site operator. A simple calculation
shows that the factor is