ij\rangle}\hat{G}_{i}\hat{G}_{j}\hat{U}_{ij}(\theta,\phi)$ applied to non-
overlapping pairs of nearest neighbour qubits $\langle i,j\rangle$. For each
pair a cycle consists of two single qubit gates $\hat{G}_{i},\hat{G}_{j}$ and
an entangling two-qubit gate,

$\displaystyle\hat{U}_{i,j}(\theta,\phi)=e^{-\frac{i}{2}\theta\left(X_{i}X_{j}+Y_{i}Y_{j}\right)-\frac{i\phi}{2}Z_{i}Z_{j}},$
(S7)

parameterized by fixed angles $\theta,\phi$. The random instances are
generated by drawing single qubit gates from a specific finite set
$\left\\{\hat{G}_{i}\right\\}$. We consider two sets corresponding to
generators of single qubit Clifford group $\\{\sqrt{X^{\pm 1}},\sqrt{Y^{\pm
1}}\\}$ and the set which in conjunction with any entangling two-qubit gate
generates a universal set $\\{\sqrt{X^{\pm 1}},\sqrt{Y^{\pm 1}},\sqrt{W^{\pm
1}},\sqrt{V^{\pm 1}}\\}$ introduced in Sec. I.

Average OTOC can be calculated by first considering the dynamics of the pair
of butterfly operators $\hat{O}(t)=U^{\dagger}\hat{O}U$ evolved under circuit
$U$. We introduce the average of the pair of the operators acting in the
direct product of two Hilbert spaces of the two replicas of the same circuit:

$\displaystyle\hat{\cal
O}^{(2)}(t)=\overline{\hat{O}(t)\otimes\hat{O}(t)}\equiv\overline{\hat{O}(t)^{\otimes
2}}$ (S8)

where averaging over ensemble of the circuits is denoted by
$\overline{(...)}$, and the rightmost equation defines the short-hand
notation.

Analogously, one can introduce the higher order averages

$\displaystyle\hat{\cal
O}^{(4)}(t)=\overline{\hat{O}(t)\otimes\hat{O}(t)\otimes\hat{O}(t)\otimes\hat{O}(t)}\equiv\overline{\hat{O}(t)^{\otimes
4}},$ (S9)

as a starting point to analyze the circuit to circuit fluctuations of $C(t)$.

Then the value of OTOC is obtained by taking the matrix element of
$C(t)=\bra{\psi}\overline{\hat{M}\hat{O}(t)\hat{M}\hat{O}(t)}|\psi\rangle$
with the initial state that in our experimental setup is chosen to be
$|\psi\rangle=\bigotimes_{i=1}^{n}|+\rangle_{i}$ where $|+\rangle_{i}$ is the
symmetric superposition of the computational basis states. The average and the
second moment of $C(t)$ are obtained as a straightforward convolution of the
indices in Eqs. (S8) and (S9) respectively. Operator $\hat{\cal O}^{(4)}$ can
be also used to study the effect of the initial state to the OTOC (clearly its
average is not sensitive to the initial conditions).

###  IV.1 Symmetric single qubit gate set

We first consider average over uniformly random single qubits gate such that
for any Pauli matrix,

$\displaystyle\overline{\hat{\alpha}_{i}}\equiv\overline{\hat{G}^{\dagger}\hat{\alpha}_{i}\hat{G}}=0,$
(S10)

We will analyze the specific discrete gate sets used in the experiment in
Section IV.4. In this section, we use the Latin indices to label a qubit, and
the Greek ones to denote the corresponding Pauli matrices
$\hat{\alpha}_{i}=\\{\hat{X}_{i},\hat{Y}_{i},\hat{Z}_{i}\\}$

For a pair of operators as in Eq. (S8) there exist an analog of scalar product
– a spherically symmetric combination that does not vanish after averaging,

$\displaystyle\overline{\hat{\alpha}_{i}\otimes\hat{\beta}_{j}}=\hat{G}^{\dagger}\hat{\alpha}_{i}\hat{G}\otimes\hat{G}^{\dagger}\hat{\beta}_{j}\hat{G}=\delta_{\alpha\beta}\delta_{ij}\mathcal{B}_{i},$
(S11)

where we introduced “bond” notation,

$\displaystyle\mathcal{B}_{i}\equiv\frac{1}{3}\hat{\alpha}_{i}\otimes\hat{\alpha}_{i}\equiv\frac{1}{3}\left(\hat{X}_{i}^{\otimes
2}+Y_{i}^{\otimes 2}+Z_{i}^{\otimes 2}\right),$ (S12)

and the summation over the repeated Greek indices is always implied.

Analogously one can find non-vanishing averages of the Pauli matrices acting
in four replicas space of the same qubit (we will omit the qubit label). There
are in total six bilinear combinations
$\displaystyle\overline{\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\openone}\otimes\hat{\openone}}=\delta_{\alpha\beta}\mathcal{B}^{(12)},$
$\displaystyle\mathcal{B}^{(12)}=\frac{1}{3}\hat{\alpha}\otimes\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\openone},$
(S13a)
$\displaystyle\overline{\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\beta}\otimes\hat{\openone}}=\delta_{\alpha\beta}\mathcal{B}^{(13)},$
$\displaystyle\mathcal{B}^{(12)}=\frac{1}{3}\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\openone},$
$\displaystyle\vdots$
$\displaystyle\overline{\hat{\openone}\otimes\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\beta}}=\delta_{\alpha\beta}\mathcal{B}^{(34)},$
$\displaystyle\mathcal{B}^{(34)}=\frac{1}{3}\hat{\openone}\otimes\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\alpha}.$

Four cubic invariants are possible because there is no inversion operation for
the spin:

$\displaystyle\overline{\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma}}=\epsilon^{\alpha\beta\gamma}\mathcal{C}^{(1)},$
$\displaystyle\quad\mathcal{C}^{(1)}\equiv\frac{\epsilon^{\alpha\beta\gamma}}{6}\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma},$
(S13b)
$\displaystyle\overline{\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\beta}\otimes\hat{\gamma}}=\epsilon^{\alpha\beta\gamma}\mathcal{C}^{(2)},$
$\displaystyle\quad\mathcal{C}^{(2)}\equiv\frac{\epsilon^{\alpha\beta\gamma}}{6}\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\beta}\otimes\hat{\gamma},$
$\displaystyle\vdots$
$\displaystyle\overline{\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma}\otimes\hat{\openone}}=\epsilon^{\alpha\beta\gamma}\mathcal{C}^{(4)},$
$\displaystyle\quad\mathcal{C}^{(4)}\equiv\frac{\epsilon^{\alpha\beta\gamma}}{6}\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma}\otimes\hat{\openone}.$

where $\epsilon_{\alpha\beta\gamma}$ is the three-dimensional Levi-Civita
symbol. Because they would change sign under the inversion operation, the
cubic invariants can appear only in pairs on neighboring qubits.

Finally, the three quartic invariants are

$\displaystyle\overline{\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma}\otimes\hat{\delta}}=\mathcal{D}^{(2)}\Upsilon_{\alpha\beta\gamma\delta}+\mathcal{D}^{(3)}\Upsilon_{\alpha\gamma\beta\delta}+\mathcal{D}^{(4)}\Upsilon_{\alpha\gamma\delta\beta},$
(S13c)
$\displaystyle\Upsilon_{\alpha\beta\gamma\delta}=\delta_{\alpha\beta}\delta_{\gamma\delta}-\frac{1}{4}\delta_{\alpha\gamma}\delta_{\beta\delta}-\frac{1}{4}\delta_{\alpha\delta}\delta_{\beta\gamma},$
$\displaystyle\mathcal{D}^{(2)}=\frac{2}{15}\hat{\alpha}\otimes\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\beta},$
$\displaystyle\mathcal{D}^{(3)}=\frac{2}{15}\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\alpha}\otimes\hat{\beta},$
$\displaystyle\mathcal{D}^{(4)}=\frac{2}{15}\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\beta}\otimes\hat{\alpha}.$

This simple form implies the spherical symmetry of the single qubit averaging.
For a lower symmetry (e.g. all the rotations of a cube) other quartic and
cubic invariants are possible (like $X^{\otimes 4}+Y^{\otimes 4}+Z^{\otimes
4}$ or
$|\epsilon_{\alpha\beta\gamma}|\hat{\openone}\otimes\alpha\otimes\beta\otimes\gamma$)
but we will ignore them for the sake of simplicity.

###  IV.2 Efficient Population Dynamics for the Averaged OTOC

We expand average evolution of the pair of butterfly operators (S8) as,

$\displaystyle\hat{\cal
O}^{(2)}(t)=\sum_{\\{v_{i}\\}}P_{\\{v_{i}\\}}\bigotimes_{i=1}^{n}\left(v_{i}\hat{\mathcal{B}}_{i}+u_{i}\hat{\openone}_{i}^{\otimes
2}\right).$ (S14)

where normalization condition reads,

$u_{i}+v_{i}=1,$ (S15)

the variable $v_{i}=\\{0,1\\},i=1,2,...,n$ indicates whether or not Pauli
matrices occupy qubit $i$, and $P_{\\{v_{i}\\}}$ are formfactors. In other
words, each $i$th is characterized either by
$\rho^{(0)}_{i}\equiv\openone_{i}^{\otimes 2}$, ”vacuum”, or the “bond”
$\mathcal{B}_{i}$. All other terms in the right hand side of Eq. (S14) vanish
upon averaging over single qubit gates see Sec. IV.1.

Application of a two-qubit Sycamore gate to a pair of qubits $\\{i,j\\}$ is
then described by $2^{2}\times 2^{2}$ matrix in the space $\\{v_{i},v_{j}\\}$,

$\displaystyle\Omega=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&1-a-b&a&b\\\
0&a&1-a-b&b\\\ 0&\frac{b}{3}&\frac{b}{3}&\left(1-\frac{2}{3}b\right)\\\
\end{array}\right),$ (S20) $\displaystyle
a=\frac{1}{3}\left(2\sin^{2}\theta+\sin^{4}\theta\right),$ $\displaystyle
b=\frac{1}{3}\left(\frac{1}{2}\sin^{2}2\theta+2\left(\sin^{2}\theta+\cos^{2}\theta\right)\right).$

where the matrix $\Omega$ acts from the right on four dimensional row vector
with the basis $(00),(01),(10),(11)$ The standard $\sqrt{\text{iSWAP}}$ gate
corresponds to $\theta=\frac{\pi}{4}$, $a=\frac{1}{12},b=\frac{1}{2}$, and
iSWAP $\theta=\frac{\pi}{2}$ with $a=\frac{1}{3},b=\frac{2}{3}$.

Each time when the two-qubit gate is applied the formfactor $P$ is updated
according to the rules

$\displaystyle
P_{v_{1}\text{...}v_{i}v_{j}\text{...}v_{n}}(t+1)=\sum_{v_{i}^{\prime}v_{j}^{\prime}}P_{v_{1}\text{...}v_{i}^{\prime}v_{j}^{\prime}\text{...}v_{n}}(t)\Omega_{v_{i}^{\prime}v_{j}^{\prime},v_{i}v_{j}}.$
(S21)

Equations (S20)–(S21) are obtained by an application of a two-qubit gate (S7)
with $\phi=0$ to a pair $(i,j)$ of factors in Eq. (S14) and averaging the
result using Eqs. (S10) – (S11).

Some additional constraints can be extracted from the exact condition

$\left[\hat{O}(t)\right]^{2}=\bigotimes_{i=1}^{n}\hat{\openone}_{n}.$ (S22)

Convoluting operators in Eq. (S14), using $\hat{\openone}_{i}^{\otimes
2}\to\hat{\openone}_{i}$, $\hat{\mathcal{B}}_{i}\to\hat{\openone}_{i}$ and the
normalization condition (S15), we obtain the requirement

$\sum_{\\{v_{i}\\}}P_{\\{v_{i}\\}}(t)=1.$ (S23)

Preserving this condition in the update rule (S21) requires the elements in
each row of the matrix $\Omega$ from Eq. (S20) to add up to one. Moreover, all
the elements of $\Omega$ are non-negative and therefore not only is
$P_{\\{v_{i}\\}}(t)$ normalized but it is also non-negative.

Therefore, rules (S21) and (S23) defines the Markov process and variables
$v_{i}=0,1$ correspond to the classical population of each qubit. As the non-
populated states $v_{i}=v_{j}=0$ do not evolve and the reproducing
$(01)\to(11)$ and destruction $(11)\to(01)$ processes are allowed the problem
(S21) is nothing but the classical population dynamics. The formfactors
$P\left(v_{1}\text{...}v_{i}^{\prime}v_{j}^{\prime}\text{...}v_{n},t\right)$
are interpreted as a distribution function over an $n$ bit register
$\left\\{v_{1}\text{...}v_{n}\right\\}$, subject to the Markov process defined
by the update, Eq. (S20).

Direct solution of Eq. (S21) would require $2^{n}$ real numbers. However,
unlike the original unitary evolution, the classical population dynamics
involves only positive numbers constrained by normalization (S23). Such
dynamics is very efficiently stimulated using a classical Monte Carlo type
algorithm.

The butterfly operator is the starting point of the Markov process Eq. (S21).
At long times the distribution $P_{\\{v_{i}\\}}$ converges to the stationary
state $\bigotimes_{i=1}^{n}\rho_{i}^{(erg)}$, where
$\rho_{i}^{(erg)}=\frac{1}{4}\left(\openone_{i}^{\otimes
2}+3\mathcal{B}_{i}\right)$, which corresponds to ”vacuum” occurring on each
site $i$ with probability $p\left(\openone_{i}^{\otimes 2}\right)=\frac{1}{4}$
and ”bond” occurring with probability $p\left(\mathcal{B}_{i}\right)=3/4$.
OTOC in this limit takes the value corresponding to the random matrix
statistics Yoshida and Kitaev (2017). Intermediate dynamics of OTOC between
these two limits is fully described by the Markov process Eq. (S21). It has a
form of shock wave spreading from the initial butterfly operator to cover the
whole system. Note that the choice of the two qubit gate parameter $\theta$
has a dramatic effect on the intermediate OTOC dynamics. As discussed in the
main text $\theta=\pi/2$ corresponds to the probability of butterfly operator
to spread equal one, and therefore spreading with maximum velocity equal to
the light cone velocity and saturating the Lieb-Robinson bound. At any other
value of $\theta$ probability to spread is less then one which results in
diffusive broadening of the front, and the center of the front propagates with
a butterfly velocity that is smaller than the light cone velocity.

###  IV.3 Sign Problem in the Population Dynamics for OTOC Fluctuations

As we already mentioned, the analytic calculation of the OTOC for an
individual circuit is impossible. One can expect, however, that it is possible
to express the variance of the OTOC in terms of some products of classical
propagators similarly to the analysis of the mesoscopic fluctuations in the
disordered metals. The purpose of this subsection is to show that even such a
modest task can not be efficiently undertaken.

The starting point is the expansion in terms of the single qubit rotations
invariants (S13). Similarly to Eq. (S14), we write for $\hat{\cal O}^{(4)}$ of
Eq. (S9)

$\begin{split}\hat{\cal
O}^{(4)}(t)=\sum_{\\{{\mathbf{V}}_{i}\\}}P_{\\{{\mathbf{V}}_{i}\\}}\bigotimes_{i=1}^{n}{\mathbf{V}}_{i}\cdot\hat{\bf\cal{Q}}_{i},\end{split}$
(S24)

where $\hat{\bf\cal{Q}}_{i}$ is the vector with $14=1+6+4+3$ operator
components given by invariants of Eq. (S9):

$\begin{split}&\hat{\bf\cal{Q}}_{i}=\left[\hat{\openone}^{\otimes
4}_{i},\mathbf{b},\mathbf{c},\mathbf{d}\right].\\\
&\mathbf{b}=\left(\hat{\mathcal{B}}^{(12)}_{i},...,\hat{\mathcal{B}}^{(34)}_{i}\right),\\\
&\mathbf{c}=\left(\hat{\mathcal{C}}^{(1)}_{i},...,\hat{\mathcal{C}}^{(4)}_{i}\right),\\\
&\mathbf{d}=\left(\hat{\mathcal{D}}^{(2)}_{i},\hat{\mathcal{D}}^{(3)}_{i},\hat{\mathcal{D}}^{(4)}_{i}\right),\end{split}$
(S25)

and ${\mathbf{V}}_{i}$ are the $14$ basis unit vectors so that $13$ components
equal to zero and the remaining component is $1$. In other words, each site
can be in $14$ possible states.

A straightforward generalization of the evolution equation (S21) reads

$\displaystyle
P_{{\mathbf{V}}_{1}\cdot{\mathbf{V}}_{i}{\mathbf{V}}_{j}\cdot{\mathbf{V}}_{n}}(t+1)=\sum_{{\mathbf{V}}_{i}^{\prime}{\mathbf{V}}_{j}^{\prime}}P_{{\mathbf{V}}_{1}\cdot{\mathbf{V}}_{i}^{\prime}{\mathbf{V}}_{j}^{\prime}\cdot{\mathbf{V}}_{n}}(t){\mathbf{\Omega}}_{{\mathbf{V}}_{i}^{\prime}{\mathbf{V}}_{j}^{\prime},{\mathbf{V}}_{i}{\mathbf{V}}_{j}},$
(S26)

where ${\mathbf{\Omega}}$ is now $14^{2}\times 14^{2}$ whose explicit form is
known but not quite important for the further consideration.

Equation (S26) involves $14^{n}$ real numbers. The only feasible path to the
solution would be an efficient Monte Carlo sampling. Naively, one can hope to
map the problem to the multicolored population dynamics. However, it is not
possible as we explain below.

Reliable Monte Carlo sampling requires (1) normalizable and non-negative
weights, and (2) absence of the correlations in contribution of different
configurations. Let us demonstrate that both conditions do not hold for the
evolution of formfactors in Eq. (S24). Once again, we use exact Eq. (S22).
Convoluting third and fourth replica in $\hat{\cal O}^{(4)}(t)$ we obtain

$\hat{\cal O}^{(4)}(t)\to\hat{\cal
O}^{(2)}(t)\otimes\hat{\mathbf{\openone}},\quad\hat{\mathbf{\openone}}\equiv\left[\bigotimes_{i=1}^{n}\hat{\openone}_{i}\right].$
(S27)

Convolutions of the four operators (S13) yield

$\displaystyle\mathcal{B}^{(12)}$
$\displaystyle\to\mathcal{B}^{(3)}\equiv\frac{1}{3}\hat{\alpha}\otimes\hat{\alpha}\otimes\hat{\openone},$
(S28) $\displaystyle\hat{\openone}^{\otimes 4},\mathcal{B}^{(34)}$
$\displaystyle\to\hat{\openone}^{\otimes 3},$
$\displaystyle\mathcal{B}^{(13)},\ \mathcal{B}^{(14)}$
$\displaystyle\to\mathcal{B}^{(2)}\equiv\frac{1}{3}\hat{\alpha}\otimes\hat{\openone}\otimes\hat{\alpha},$
$\displaystyle\mathcal{B}^{(23)},\ \mathcal{B}^{(24)}$
$\displaystyle\to\mathcal{B}^{(1)}\equiv\frac{1}{3}\hat{\openone}\otimes\hat{\alpha}\otimes\hat{\alpha},$
$\displaystyle\mathcal{C}^{(3)},\mathcal{C}^{(4)}$
$\displaystyle\to\mathcal{C}\equiv\frac{\epsilon^{\alpha\beta\gamma}}{6}\hat{\alpha}\otimes\hat{\beta}\otimes\hat{\gamma},$
$\displaystyle\mathcal{C}^{(1)}$ $\displaystyle\to
i{B}^{(1)},\qquad\mathcal{C}^{(2)}\to i{B}^{(2)},$
$\displaystyle\mathcal{D}^{(2)}$
$\displaystyle\to\frac{6}{5}\mathcal{B}^{(3)},$
$\displaystyle\mathcal{D}^{(3)}$
$\displaystyle\to\frac{2}{5}\mathcal{B}^{(3)}+\frac{4i}{5}\mathcal{C},$
$\displaystyle\mathcal{D}^{(4)}$
$\displaystyle\to\frac{2}{5}\mathcal{B}^{(3)}-\frac{4i}{5}\mathcal{C}.$

Let us convolute both sides of Eq. (S24) using the rules (S27) – (S28). We
find

$\begin{split}\hat{\cal
O}^{(2)}&(t)\otimes\hat{\mathbf{\openone}}=\sum_{\\{{\mathbf{V}}_{i}\\}}P_{\\{{\mathbf{V}}_{i}\\}}\bigotimes_{i=1}^{n}\left(q_{i}^{(1)}+q_{i}^{(2)}+q_{i}^{(3)}+q_{i}^{(4)}\right),\\\
q_{i}^{(1)}&\equiv V_{i}^{(1)}\hat{\openone}^{\otimes
3}_{i}+V_{i}^{(2)}\mathcal{B}^{(3)}_{i},\\\ q_{i}^{(2)}&\equiv
V_{i}^{(7)}\hat{\openone}^{\otimes
3}_{i}+\frac{2}{5}\left(3V_{i}^{(12)}+V_{i}^{(13)}+V_{i}^{(14)}\right)\mathcal{B}^{(3)}_{i},\\\
q_{i}^{(3)}&\equiv\left(V_{i}^{(3)}+V_{i}^{(4)}+iV_{i}^{(8)}\right)\mathcal{B}^{(1)}_{i}\\\
+&\left(V_{i}^{(5)}+V_{i}^{(6)}+iV_{i}^{(9)}\right)\mathcal{B}^{(2)}_{i},\\\
q_{i}^{(4)}&\equiv\left(V_{i}^{(10)}+V_{i}^{(11)}+\frac{4i}{5}V_{i}^{(13)}-\frac{4i}{5}V_{i}^{(14)}\right)\mathcal{C}_{i}.\end{split}$

Here, the superscript in vector $V_{i}^{(\cdot)}$ enumerates components
according to Eq. (S25). The total result is real as the imaginary terms are
always generated in pairs.

Configurations with $V_{i}^{(1)}=u_{i}$, $V_{i}^{(2)}=v_{i}$,
$V_{i}^{(3,..,14)}=0$, exactly reproduce averaged result (S14). It means that
all the other terms exactly cancel each other. It is possible only if the
corresponding formfactors can be of both signs. Therefore, any finite
unconstrained sampling leads to an arbitrary result (sign problem). Moreover,
the constraints are non-local. Consider, e.g. cancellation of contribution
proportional to $\mathcal{B}^{(1)}_{i}\mathcal{B}^{(1)}_{j}$. Cancellation
occurs only if the formfactors different by replacement
$V_{i}^{(8)}V_{j}^{(8)}\to V_{i}^{(a)}V_{j}^{(b)}$, $a,b=\\{3,4\\}$ must be
kept precisely the same. The requirement quickly becomes intractable with the
increasing of the number of non-trivial matrices involved into cancellation.

To summarize, the general requirements on the evolution of the $\hat{\cal
O}^{(4)}(t)$ lead to the sign and locality problems. Those two facts render a
brute force classical Monte Carlo algorithm impossible. At present, we are not
aware of any algorithm enabling to circumvent those obstacles.

###  IV.4 Population Dynamics for iSWAP Gate Sets Implemented in the Main Text

Circuits with $\theta=\pi/2$ were used to realize Clifford as well as a
universal ensemble. It is therefore instructive to study dynamics of the
specific gate sets used in the experiment in more detail. Consider conjugation
of a pair of Pauli operators $\hat{\alpha}_{i}\hat{\beta}_{j}$ by the iSWAP
gate $\hat{S}$. It maps the pair onto another pair of Pauli operators
$\hat{\gamma}_{i}\hat{\delta}_{j}$ according to the following rules,

$\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|}\hline\cr\hat{\alpha}_{i}\hat{\beta}_{j}&\hat{X}_{i}\hat{\openone}_{j}&\hat{Y}_{i}\hat{\openone}_{j}&\hat{Z}_{i}\hat{\openone}_{j}&\hat{Z}_{i}\hat{X}_{j}&\hat{Z}_{i}\hat{Y}_{j}&\hat{\alpha}_{i}\hat{\alpha}_{j}&\hat{X}_{i}\hat{Y}_{j}\\\
\hline\cr
S^{\dagger}\hat{\alpha}_{i}\hat{\beta}_{j}S&-\hat{Z}_{i}\hat{Y}_{j}&\hat{Z}_{i}\hat{X}_{j}&\hat{\openone}_{i}\hat{Z}_{j}&-\hat{Y}_{i}\hat{\openone}_{j}&\hat{X}_{i}\hat{\openone}_{j}&\hat{\alpha}_{i}\hat{\alpha}_{j}&\hat{Y}_{i}\hat{X}_{j}\\\
\hline\cr\end{array}$

####  IV.4.1 Clifford Gate Set

Clifford gate set used to obtain the data in the main text is drawn form the
single qubit gate set
$\\{\hat{G}\\}=\\{\sqrt{X},\sqrt{X^{-1}},\sqrt{Y},\sqrt{Y^{-1}}\\}$. Averaging
over this gates set of a symmetric pair of Pauli operators
$\alpha_{i}\otimes\alpha_{i}$ reads,

$\displaystyle\overline{\hat{G}^{\dagger}\hat{\alpha}\hat{G}\otimes\hat{G}^{\dagger}\hat{\alpha}\hat{G}}=\Xi_{\alpha},$
(S29)

where we introduce the basis
$\Xi_{\alpha}=\\{\openone,\mathbb{X},\mathbb{Y},\mathbb{Z}\\}$,

$\displaystyle\mathbb{X}=\frac{1}{2}(\hat{Z}^{\otimes 2}+\hat{Y}^{\otimes
2}),$ $\displaystyle\mathbb{Y}=\frac{1}{2}(\hat{X}^{\otimes
2}+\hat{Z}^{\otimes 2}),$
$\displaystyle\mathbb{Z}=\frac{1}{2}(\hat{X}^{\otimes 2}+\hat{Y}^{\otimes
2}).$

Single qubit gate average can be described in terms of the invariants
$\Xi_{\alpha}\to M^{(S)}_{\alpha\beta}\Xi_{\beta}$,

$\displaystyle M^{(S)}=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&0&\tfrac{1}{2}&\tfrac{1}{2}\\\ 0&\tfrac{1}{2}&0&\tfrac{1}{2}\\\
0&\tfrac{1}{2}&\tfrac{1}{2}&0\\\ \end{array}\right).$ (S34)

In this new basis $\\{\Xi_{\gamma}\Xi_{\delta}\\}$ instead of a $4\times 4$
matrix $\Omega$ in Eq. (S20) effect of iSWAP gate is described by $16\times
16$ matrix constructed straightforwardly from the rules Eqs. (S34) and the
rules for iSWAP.

####  IV.4.2 Universal Gate Set

Figure S19: Extended data for Fig. 2 in the main text: Average OTOC for four
different locations of butterfly operator, $X$. Solid lines correspond to
experimental data and dashed lines to the population dynamics simulation with
Eq. (S39).

The universal gates set used in the main text consists of eight choices of
single qubit gates $\\{\hat{G}\\}=\\{\sqrt{X^{\pm 1}},\sqrt{Y^{\pm
1}},\sqrt{W^{\pm 1}},\sqrt{V^{\pm 1}}\\}$. For the butterfly operator we
choose $\hat{O}_{i}=\hat{X}_{i}$ the dynamics is described in the reduced
subspace spanned by the basis $\Lambda=\\{\openone^{\otimes
2},\hat{X}^{\otimes 2},\hat{Y}^{\otimes 2},\hat{Z}^{\otimes 2}\\}$. Average of
a pair of Pauli operators over this ensemble reads $\Lambda_{\alpha}\to
M^{(U)}_{\alpha\beta}\Lambda_{\beta}$,

$\displaystyle M^{(U)}=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\tfrac{3}{8}&\tfrac{1}{8}&\tfrac{1}{2}\\\
0&\tfrac{1}{8}&\tfrac{3}{8}&\tfrac{1}{2}\\\ 0&\tfrac{1}{2}&\tfrac{1}{2}&0\\\
\end{array}\right),$ (S39)

Together with the rules for iSWAP gate it is straightforward to generate
$16\times 16$ matrix defining the Markov population dynamics process.
Comparison of this noise-free population dynamics prediction for several
different circuits with different locations of the butterfly operator $X$
throughout the 53 qubit chip implemented experimentally is shown in Fig. S19.

##  V Efficient Population Dynamics for Noisy Circuits

###  V.1 Inversion Error

The circuit to measure OTOC require inversion of gates
$\hat{U}_{i,j}(\theta,\phi)$ which is not perfect. An important source of
error is non-invertible phase $\phi$ such that instead of
$\hat{U}_{i,j}^{\dagger}(\theta,\phi)$ the gate $\hat{U}_{i,j}(-\theta,\phi)$
is implemented. This inversion error can be included in the Markov population
dynamics process in the following way,

$\displaystyle\Omega=\left(\begin{array}[]{cccc}\cos^{2}\phi&0&0&\sin^{2}\phi\\\
0&1-\tilde{a}-b&\tilde{a}&b\\\ 0&\tilde{a}&1-\tilde{a}-b&b\\\
\tfrac{\sin^{2}\phi}{9}&\frac{b}{3}&\frac{b}{3}&\tfrac{8+\cos^{2}\phi}{9}-\frac{2}{3}b\\\
\end{array}\right),$
$\displaystyle\tilde{a}=\frac{1}{3}\left(\cos^{2}\phi\sin^{4}\theta+\sin^{2}\phi\cos^{4}\theta\right).$

###  V.2 Generic Error Model

Figure S20: Average fidelity of OTOC, i.e. the circuit with no butterfly
applied shown as red line. OTOC fidelity from Monte Carlo simulation of noisy
population dynamics Sec. V. The Pauli error rate $r_{p}=0.013$ is determined
by fitting the numerics to the experimental data. This value is within $10\%$
of the gate fidelity determined in a separate experiment benchmarking two-
qubit gates, see Sec. II. This is the single parameter of the model which is
used to reproduce the rest of the OTOC data for 51 different locations of
butterfly operator on the chip. Black dashed line shows a naive estimate of
circuit fidelity via product of fidelities of gates within the light cone.
Note that naive fidelity decays much faster with time than OTOC fidelity.

For an ensemble that satisfies Eq. (S10) the effect of relaxation and
dephasing in the quantum processor can be captured by a one and two qubit
depolarizing channel noise model,

$\displaystyle\rho\to(1-p_{1})\rho+\frac{p_{1}}{3}\sum_{\hat{\alpha}=\hat{X},\hat{Y},\hat{Z}}\hat{\alpha}\rho\hat{\alpha},$
(S40)
$\displaystyle\rho\to(1-p_{2})\rho+\frac{p_{2}}{15}\sum_{\hat{\alpha},\hat{\beta}}\hat{\alpha}\hat{\beta}\rho\hat{\alpha}\hat{\beta}.$
(S41)

In this case the Markov process Eq. (S20) can be modified in a straight
forward way. For a single qubit depolarizing channel the Markov process Eq.
(S20) is supplemented by the exponential decay rate as follows,

$\displaystyle\openone_{i}^{\otimes 2}\rightarrow\openone_{i}^{\otimes 2},~{}$
(S42) $\displaystyle\mathcal{B}_{i}\rightarrow
e^{-4p_{1}/3}\mathcal{B}_{i},~{}$ (S43)

for each bond per gate cycle. Two-qubit depolarizing channel noise is
accounted by supplementing the Markov process by the following,

$\displaystyle\openone_{i}^{\otimes 2}\openone_{j}^{\otimes
2}\rightarrow\openone_{i}^{\otimes 2}\openone_{j}^{\otimes 2},~{}$ (S44)
$\displaystyle\mathcal{B}_{i}\openone_{j}^{\otimes 2}\rightarrow
e^{-\frac{16}{15}p_{2}}\mathcal{B}_{i}\openone_{j}^{\otimes 2},~{}$ (S45)
$\displaystyle\mathcal{B}_{i}\mathcal{B}_{j}\rightarrow
e^{-\frac{16}{15}p_{2}}\mathcal{B}_{i}\mathcal{B}_{j}.~{}$ (S46)

Note that the effect of noise on the Markov population dynamics cannot be
described by the global depolarizing channel model that is often conjectured
for ergodic circuits. Instead the time dependence of OTOC will demonstrate
characteristic time dependence that can be used to verify the experimental
data.

Figure S21: OTOC circuits (with no butterfly applied) for universal iSWAP gate
set. Solid red line shows experimental data and dashed blue line is the
theoretical fit with a single parameter, Pauli error rate. The extracted Pauli
error rate is $r_{p}=0.012$ within $20\%$ of the value obtained via two-qubit
gate benchmarking.

The described procedure allows us to verify the experimental results for
average fidelity of OTOC circuits, the circuit with no butterfly applied, by
direct comparison to the noisy population dynamics, see Fig. S20 and Fig. S21.
We use the two-qubit Pauli error as a fitting parameter. Best fit corresponds
to errors within $10\%$ of the average error of two qubit cycle measured
independently. We use the extracted error to predict values of OTOC for every
position of the butterfly with no additional fitting. Comparison of the values
of normalized OTOC predicted in this way with experimental data is shown in
Fig. S22. This procedure introduces substantial noise-dependent bias into the
observed OTOC values, which is illustrated in Fig. S23.

Figure S22:  Extended data for Fig. 2 of the main text: Normalized OTOC for
four different locations of butterfly operator $Z$. Solid lines show
experimental data, and dashed lines correspond to the noisy population
dynamics.

###  V.3 Error Mitigation for Average OTOC

We estimate circuit fidelity from the circuit shown in, Fig. 1 A, B, without
the butterfly operator. This fidelity estimate is used for the error
mitigation procedure aplied to the data in Fig. 2. In the error-free circuit
absence of the butterfly operator means that $U$ and $U^{\dagger}$ cancel each
other exactly resulting in $C_{z0}(t)=1$. In practice, inversion is imperfect
as detailed in Sec. V.1. Both errors in the unitary parameters and the effect
of noise are reflected in the time dependence of $C_{z0}(t)<1$ which serves as
circuit fidelity, an analog of Loschmidt echo for local operator $\hat{M}$.

Average over circuits reads,

$\displaystyle\overline{C}_{0z}=\mathcal{F}_{\openone^{\otimes
2}_{1}}+\mathcal{F}_{\mathcal{B}_{1}},$ (S47)
$\displaystyle\overline{C}_{zz}=\mathcal{F}_{\openone^{\otimes
2}_{1}}-\frac{1}{3}\mathcal{F}_{\mathcal{B}_{1}}.$ (S48)

where the probabilities of vacuum $\mathcal{F}_{\openone^{\otimes 2}_{1}}$ and
bond $\mathcal{F}_{\mathcal{B}_{1}}$ at the measurement site are described by
the population dynamics Eqs. (S20) with respective decay rates Eqs. (S42-S46).
Note that the decay rate grows with the extent of spreading of the butterfly
operator, resulting in the decay of fidelity that is not a simple exponent.
Moreover, in general $\mathcal{F}_{\openone^{\otimes
2}_{1}},\mathcal{F}_{\mathcal{B}_{1}}$ is not a simple probability no error
occurred, as one could naively expect, see Fig. S20 for comparison.

The ratio of $\overline{C}_{zz}/\overline{C}_{0z}$ is compared to normalized
data in Fig. S22. Note that this ratio does not correspond to the noise free
OTOC, Fig. S23. This is because the probability of vacuum at measure site
$\mathcal{F}_{\openone^{\otimes 2}_{1}}$ decays slower with respect to its
noise free value $p(\openone^{\otimes 2}_{1})$ than the probability of the
bond $\mathcal{F}_{\mathcal{B}_{1}}$. The latter corresponds to the weight of
the operators which span the distance between measure and butterfly qubit
which are more susceptible to noise. As a result normalized OTOC from noisy
circuits overestimates the value of noise-free OTOC. For individual circuits
such a simple description is no longer valid, as demonstrated by the
dependence of fidelity estimate $C_{z0}(t)$ on the circuit instance, see Sec.
VI. In many cases circuit dependent corrections are relatively small and the
normalization procedure still works for individual circuits as well.
Nonetheless, for our OTOC fluctuations data we develop a more efficient error
mitigation procedure described in the following section.

Figure S23: Comparison of noisy OTOC obtained by normalization procedure,
dashed lines, to noise-free population dynamics for the same circuit, solid
lines.

###  V.4 Theory of Error Mitigation for Individual Circuits

We use a more precise error mitigation procedure for individual OTOC
measurements presented in Figs. 3, 4 of the main text. These circuits contain
iSWAP entangling gate and OTOC value can be expanded in terms of contributions
from Clifford circuits. In presence of non-Clifford gates the butterfly
operator can be conveniently expanded into Pauli strings $B_{i}$,

$\displaystyle O(t)=\sum
w_{\alpha_{1}...\alpha_{n}}B_{\alpha_{1}...\alpha_{n}},$ $\displaystyle
B_{\alpha_{1}...\alpha_{n}}=\hat{\alpha}_{1}\otimes...\otimes\hat{\alpha}_{n}$

For the initial state used in the experimental protocol
$|\psi\rangle=\bigotimes|+\rangle_{i}$, the value of OTOC is expanded as,

$\displaystyle C=\sum
w_{\alpha\alpha_{2}...\alpha_{n}}w_{\beta\alpha_{2}...\alpha_{n}}\bra{+_{1}}\sigma^{z}_{1}\sigma^{\alpha}_{1}\sigma^{z}_{1}\sigma^{\beta}_{1}|+\rangle_{1}.$
(S49)

The real part of the OTOC corresponds to $\alpha=\beta$,

$\displaystyle\mathrm{Re}C=\sum
w_{\alpha_{1}...\alpha_{n}}^{2}\kappa_{\alpha_{1}},$ (S50)

where $\kappa_{\alpha}=\\{1,-1,-1,1\\}$ for $\alpha=\\{0,x,y,z\\}$.

For Clifford circuits ideal value of OTOC can be calculated efficiently. It
can then be used to calculate circuit fidelity by comparing data with the
expected value. The fidelity of the OTOC in presence of non-Clifford gates is
then calculated by sampling a subset of OTOCs for Cliffords that appear in its
expansion, see Eq. (S49). The fidelity is calculated by averaging over
fidelities of individual Clifford contributions. Averaging reduces the circuit
specific effect of noise and gives a more accurate estimate of fidelity.

##  VI Numerical Simulation of Error-Limiting Mechanisms for OTOC

In this section, we provide numerical simulation results aimed at identifying
the potential error sources that limit our experimental accuracy in resolving
OTOCs. As demonstrated in Section I, shot noise from finite statistical
sampling is unlikely the dominant mechanism. The remaining known error
channels are: 1. Incoherent, depolarizing noise in the quantum circuits, which
can arise from qubit dephasing or relaxation. 2. Coherent errors in the
quantum gates, e.g. the remnant conditional-phase $\phi$ demonstrated in
Section II. We study each of these errors below.

###  VI.1 Coherent and Incoherent Contributions

Figure S24: Simulation of incoherent and coherent OTOC errors. (A) OTOCs, $C$,
of 47 quantum circuit instances simulated in three different ways: with $\phi$
error, where a conditional-phase of $0.136$ rad is added to every two-qubit
gate; with $r_{\text{p}}$ error, where a Pauli error of 0.015 is added to each
two-qubit gate; ideal, where no error is introduced. The qubit configuration
used here is a 1D chain of 10 qubits wherein $Q_{1}$ and $Q_{\text{b}}$ reside
at opposite ends of the chain. The two-qubit gate is iSWAP and
$N_{\text{D}}=12$ non-Clifford gates are used for each instance. The number of
circuit cycles is 34 ($N_{\text{S}}=250$). (B) The scaling of OTOC error (i.e.
RMS deviation from the ideal OTOC values) against $\phi$ (red) and
$r_{\text{p}}$ (blue). Dashed line shows the size of the OTOC signal (i.e. RMS
value of the ideal OTOCs).

We first describe how OTOC error from depolarizing noise may be simulated. We
consider a depolarizing channel model parameterized by an error probability
$p$,

$\mathcal{E}(\rho)=(1-p)\rho+p\frac{\mathbf{I}}{d},$ (S51)

where $d=2^{n}$ is the dimension, and $\mathbf{I}$ the identity operator. The
Kraus operators of this map are all Pauli strings of length $n$, where each
non-trivial string has weight $\frac{p}{d^{2}}$. We consider a model where
after each two qubit gate (these typically dominate the loss in fidelity
compared to single-qubit gates), a two-qubit depolarizing channel is applied.

In Fig. S24A, we show a number of instance-dependent OTOC values $C$ simulated
using full density matrix calculations and an experimentally measured Pauli
error rate $r_{\text{p}}$ of 0.015 ($r_{\text{p}}=15p/16$ in Eq. (S51), due to
the trivial Pauli string in the Kraus map from the identity matrix). Here we
have used a smaller number of qubits due to the high cost of density-matrix
simulation. To match with experiment, we have adjusted the number of circuit
cycles to yield a total number of iSWAP gates close to those shown in Fig. 4
of the main text. The same normalization protocol as the experiment is also
adopted, such that $\braket{\hat{\sigma}_{\text{y}}}$ is simulated with and
without the butterfly operator and their ratio is recorded as $C$. The
results, compared to the ideal OTOC values also plotted in Fig. S24B, show
little deviation.

Next, we simulate OTOCs of the same circuits with $r_{\text{p}}=0$ but an
experimentally measured conditional-phase of $\phi=0.136$ rad on each iSWAP
gate and its inverse. The results are also plotted in Fig. S24A and seen to
deviate much more from the ideal values. To quantify these observations, we
have plotted the OTOC error as a function of both $r_{\text{p}}$ and $\phi$ in
Fig. S24B. Here we see that at the experimental limit ($r_{\text{p}}=0.015$
and $\phi=0.136$ rad), the OTOC error is dominated by the contribution from
$\phi$ (where it is $\sim$0.03) rather than $r_{\text{p}}$ (where it is
$\sim$0.006). Interestingly, the SNR is about 0.9 for $\phi=0.136$ rad, which
is close to the value measured in Fig. 4 of the main text for
$N_{\text{S}}\approx 250$.

Figure S24B also provides preliminary indication of how the OTOC accuracy for
our experiments may be improved as we decrease both $\phi$ and $r_{\text{p}}$.
We see that the OTOC error is linearly proportional to $r_{\text{p}}$ whereas
its scaling is steeper against $\phi$. In particular, reducing $\phi$ by a
factor of 2 leads to an OTOC error that is 3 times lower. Although these
results may depend on the number of qubits and the structure of circuits,
their similarity to experimental data nevertheless provides hints that
reducing the conditional-phases in our iSWAP gates can potentially lead to
large improvements in the OTOC accuracy. This can, for example, be achieved
through concatenated pulses demonstrated in Ref. Foxen  _et al._ (2020).

###  VI.2 Perturbative Expansion of OTOC Error

Figure S25: Relative error $\epsilon$ of the 0’th and 1’st order error
approximation (from Eq. (LABEL:eq:single-error-approx)) calculated by a pure
state simulation, comparing to the exact value from a density matrix
calculation, as a function of depolarizing probability $p$ (Eq.(S51)). Here we
compute $C_{\text{zx}}$, varying the number of iSWAP ($N_{s}$) and non-
Clifford gates ($N_{D}$), with 11 qubits on a chain. The dots are the median
over 48 circuits with lines of best fit given by the legend. The shaded region
is the middle 50% of the data. The scaling of the error is close to
$\epsilon\propto p^{2}$ as expected for the one-error approximation in all
three cases.

The exact density-matrix calculation used in the preceding section is costly
to implement and becomes quickly intractable as the number of qubits increase.
One can more systematically estimate the instance specific noise contribution
from a perturbation theory expansion of the quantum map Eq. (S51), with an
entirely pure state calculation. Let us adopt the notation
$\sigma_{m_{1},m_{2}}^{(i,j)}=\sigma_{m_{1}}^{(i)}\sigma_{m_{2}}^{(j)}$, where
$\sigma_{m}^{(i)}$ is the $m$’th Pauli operator ($m\in\\{0,x,y,z\\}$) applied
on qubit $i$. We will also call $C[\sigma_{m_{1},m_{2}}^{(i,j)}(d)]$ the OTOC
value with the additional ‘error gate’ $\sigma_{m_{1},m_{2}}^{(i,j)}$ inserted
at layer $d$ in the circuit. Then, to accuracy $O(p^{2})$ one has

$\begin{split}&C_{p}=(1-p)^{n_{2}}C_{\mathrm{ideal}}\\\
&+\frac{p}{15}(1-p)^{n_{2}-1}\sum_{\begin{subarray}{c}(m_{1},m_{2})\neq(0,0),\\\
\langle
i,j\rangle,d\end{subarray}}C[\sigma_{m_{1},m_{2}}^{(i,j)}(d)],\end{split}$
(S52)

where $C_{\mathrm{ideal}}$ is the OTOC value in the limit of no noise, $C_{p}$
the first order approximation of $C_{\mathrm{ideal}}$ in parameter $p$, and
$n_{2}$ the total number of two-qubit gates in the circuit, occurring over
pairs defined by $\langle i,j\rangle$ (i.e. Eq. (LABEL:eq:single-error-approx)
is a sum over all error terms in the circuit). For convenience, we have also
redefined $p\rightarrow 15p/16$ from Eq. (S51) due to the trivial contribution
from the all zero Pauli string.

Eq. (LABEL:eq:single-error-approx) can be used to separate out the ideal OTOC
value of an individual circuit, from the noise contribution. This noise
contribution can be computed using $15\times n_{2}$ circuit simulations
(inserting each Pauli pair at all two-qubit gate locations in the circuit),
though in some cases symmetries can be used to reduce this. For example, for
the normalization curve $C_{\text{z0}}$, only around half of this is required,
since for each error in the reverse circuit $U^{\dagger}$, there is an
equivalent one in $U$. Of course one can continue this expansion to arbitrary
order, however the number of terms quickly becomes infeasible to compute (the
$k$’th order contains $15\times{n_{2}\choose k}$ terms).

In Fig. S25 we show the error scaling (comparing to an exact density matrix
computation) for the zero’th (i.e. only the first term in Eq.
(LABEL:eq:single-error-approx)) and first order approximation for
$C_{\text{zx}}$, giving scaling in the error as expected; the error $\epsilon$
is computed as
$|C_{p}(\text{exact})-C_{p}(\text{approx})|/|C_{p}(\text{exact})|$, where
$C_{p}(\text{approx})$ is from Eq. (LABEL:eq:single-error-approx), and
$C_{p}(\text{exact})$ the value from the density matrix calculation with noise
rate $p$. Moreover, the accuracy of the approximation remains fairly
consistent, for a fixed noise level ($p$) for the three different circuit
sizes (using a similar number of iSWAP gates as in the main text).

Here we have outlined a protocol of simulating the instance specific error
contribution to the OTOC. This can be used to more accurately separate the
instance specific OTOC value, systematic errors, and contributions from gate
noise. Of course, the overhead in the simulation is a factor of $15n_{2}$,
which can itself become challenging in deep enough circuits, although
statistical sampling may be feasible in some cases.