aainstitutetext: Joint Quantum Institute, Department of Physics, University of
Maryland, College Park, MD 20742, USAbbinstitutetext: Condensed Matter Theory
Center, Department of Physics, University of Maryland, College Park, MD 20742,
USA

#  Effective field theory of random quantum circuits

Yunxiang Liao a Victor Galitski<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantum circuits have been widely used as a platform to simulate generic
quantum many-body systems. In particular, random quantum circuits provide a
means to probe universal features of many-body quantum chaos and ergodicity.
Some such features have already been experimentally demonstrated in the noisy
intermediate-scale quantum (NISQ) devices. On the theory side, properties of
random quantum circuits have been studied on a case-by-case basis and for
certain specific systems, a hallmark of quantum chaos - universal Wigner-Dyson
level statistics - has been derived. This work develops an effective field
theory for a large class of random quantum circuits. The theory has the form
of a replica sigma model and is similar to the low-energy approach to
diffusion in disordered systems. The method is used to explicitly derive
universal random matrix behavior of a large family of random circuits. In
particular, we rederive Wigner-Dyson spectral statistics of the brickwork
circuit model by Chan, De Luca, and Chalker [Phys. Rev. X 8, 041019 (2018)]
and show within the same calculation that its various permutations and higher-
dimensional generalizations preserve the universal level statistics. Finally,
we use the replica sigma model framework to rederive the Weingarten calculus,
which is a method to evaluate integrals of polynomials of matrix elements with
respect to the Haar measure over compact groups and has many applications in
the studies of quantum circuits. The effective field theory, derived here,
provides both a method to quantitatively characterize quantum dynamics of
random Floquet systems (e.g., calculating operator and entanglement spreading)
and also a path to understanding the general fundamental mechanism behind
quantum chaos and thermalization in these systems.

##  1 Introduction

Recent years have seen a surge of interest in random quantum circuits, which
can be used to simulate various properties of interacting many-body quantum
systems, including universal ergodic dynamics. Random quantum circuits QC-0 ;
QC-1 ; QC-2 ; QC-3 ; QC-4 ; QC-5 ; QC-6 ; QC-7 ; QC-8 ; QC-9 ; QC-10 ; QC-11
consist of qubits (or qudits) evolving under successive applications of
unitary quantum gates drawn randomly from ensembles of unitaries. They capture
several general properties of many-body quantum chaotic systems. In
particular, universal random matrix theory (RMT) statistics of the quasi-
energy spectra has been derived for some time-periodic (Floquet) circuits
Chalker-1 ; Chalker-2 ; Chalker-3 ; Chalker-4 ; Chalker-5 ; Chalker-6 ;
Prosen-SFF3 ; Prosen-SFF2 ; Prosen-GSFF ; Huse2021 ; chan2021 . Quantum
systems with underlying classical chaotic dynamics are conjectured to exhibit
level statistics identical to that of a suitably chosen random matrix ensemble
BGS , and the RMT statistics has been used as one of the indicators for
quantum chaos. For decades, numerous efforts have been made to understand
theoretically this Bohigas-Giannoni-Schmit (BGS) quantum chaos conjecture
Prosen-SFF1 ; Prosen-Fermion ; Prosen-Boson ; Prosen-KIsing ; BerrySFF ;
Sieber_2001 ; Muller_2005 ; SYK-Altland ; SSS ; PRB ; PRR ; PRL . Random
quantum circuits, due to their fine-tuned structure, allow exact analysis of
the spectral statistics, which sheds light on the underlying mechanism
responsible for the emergence of RMT structure. In addition to the RMT
spectral statistics, random quantum circuits also exhibit other fundamental
properties of many-body quantum chaotic systems, such as the decay of
correlation functions of local observables Prosen-DUQC-corr ; Prosen-DUQC-dyn
, ballistic spreading of the local operators Nahum-2018 ; Sondhi ;
Keyserlingk-2 ; Huse-2018 ; HJ-2018 ; Chalker-1 ; Bertini , ballistic growth
of the entanglement Nahum-2017 ; Nahum-2018-2 ; Sondhi ; Prosen-Spin-EE ;
Prosen-DUQC-dyn ; Prosen-QDUC-EE-I ; Prosen-QDUC-EE-II ; Chalker-1 ;
Gopalakrishnan ; Keyserlingk-2019 ; Sondhi ; Chalker-1  and Gaussian
distribution of the matrix elements of observables in the energy eigenbasis
(as expected from the eigenstate thermalization hypothesis srednicki1994chaos
; srednicki1999 ; Deutsch ) Chalker-ETH ; Prosen-ETH . Experimentally, random
quantum circuits can be simulated in the noisy intermediate-scale quantum
(NISQ) devices Preskill  built with superconducting qudits SC-1 ; SC-2 ,
trapped ions ion-1 ; ion-2 , and Rydberg atoms Rydberg , and some of these
generic features for quantum chaotic systems have been observed google-
scrambling . Due to their high controllability, random quantum circuits
provide a useful tool to explore fundamental principles underlying chaos and
thermalization.

At a high level, the approaches to study the conventional interacting many-
body systems can be broadly categorized as either (i) The microscopic
approach, where the specific Hamiltonian of a particular system is studied in
the ab initio or bottom up way and where all microscopic details enter the
calculation or (ii) The low-energy (usually field-theory-based) approach where
only the universal low-energy features enter the theory (e.g., gapless
collective modes). In the context of the quantum circuits used to simulate
many-body systems, it was mostly the former-type of approach that has been
used to date. In particular, theoretical studies of random quantum circuits
relied heavily on fine-tuned features of specific models. Here we construct a
field-theory approach, where universal features of ergodic dynamics of a large
family of random quantum circuits can be studied in a unified way.

This field-theory method generalizes a supersymmetric sigma model developed by
Zirnbauer Zirnbauer_1996 ; Zirnbauer_1998 ; Zirnbauer_1999 , Altland QKR  and
others making use of a generalized Hubbard-Stratonovich transformation - the
color-flavor transformation Zirnbauer_1996 ; Zirnbauer_1998 ; Zirnbauer_1999 ;
Zirnbauer_2021 . This model is similar to the sigma model for disordered
electron systems Wegner-1979 ; Efetov , but is formulated for systems modeled
by an evolution operator or a scattering matrix, rather than a Hamiltonian. In
particular, it has been used to investigate the connection between the
appearance of universal RMT statistics and the underlying classical chaotic
dynamics (i.e., the BGS quantum chaos conjecture BGS ) in generic quantum
chaotic maps Zirnbauer_1998 ; Zirnbauer_1999 , localization in the quantum
kicked rotors QKR ; Tian_2010 ; Tian_2011 ; Tian_2012 , eigenenergy and
eigenfunction statistics of quantum graphs graph-1 ; graph-2 ; graph-eigenfun
; graph-eigenfun-2 , and the plateau transition in quantum Hall systems
Zirnbauer-QH ; Zirnbauer-QH2 . See Refs. Altland-rev ; Haake  for a review of
this model and a more complete list of references. For quantum circuits, we
find it convenient to reformulate this sigma model in the replica formalism
(in contrast to supersymmetry), and use it to derive an ensemble-averaged
effective field theory.

The spectral statistics of the Floquet operator of time-periodic systems can
be extracted from this effective field theory. In contrast to the
aforementioned model-specific methods, the field-theoretic approach is
applicable to a wide class of Floquet systems, and is particularly useful for
random quantum circuits composed of local random gates. To demonstrate the
usefulness of this method, we apply it to study the spectral statistics of a
family of Floquet random quantum circuits, some of which are shown in Figs. 1
and 2. These Floquet circuits are composed of random unitary matrices which
are drawn randomly and independently from the Circular Unitary Ensemble (CUE)
and are applied to all pairs of neighboring qudits at various time substeps
during one period. Among these circuits, the brickwork circuit depicted in
Fig. 1(a) was studied earlier by Chan, De Luca, and Chalker in Ref. Chalker-1
(see also a related work Ref. Chalker-4 ), and the RMT spectral statistics
were derived in the limit of a large on-site Hilbert space dimension
$q\rightarrow\infty$, using a generalized diagrammatic approach initially
constructed in Ref. Beenakker  for the Weingarten calculus Samuel ; Weingarten
; collins2003 ; collins2006 ; collins2021 ; Weingarten-rev ; Zee . Applying
the field theoretical approach, we rederive this result and also show that the
universal RMT statistics is preserved under an arbitrary reordering of the
two-qudit gates, for both periodic and open boundary conditions. In
particular, we show that the effective field theory describing the Floquet
quantum circuit of this type is identical to that of the CUE ensemble.
Moreover, we consider higher-dimensional generalizations of this family of the
Floquet circuits, and find that their quasi-energy spectra all exhibit the RMT
statistics (see Fig. 2 for some of the 2D examples).

Figure 1:  Floquet operators of 1D Floquet random quantum circuits with
different orderings of the two-qudit unitary gates: (a) brickwork circuit
Chalker-1 , (c) staircase circuit, and (b) circuit obtained from rearranging
the local gates (blue boxes) in (a) or (b). The horizontal direction
represents the space coordinates, and each qudit is indicated by a black dot.
The vertical direction shows the discrete time evolution within one period,
and different layers represent different substeps. Each blue box represents an
independent CUE matrix $w^{(n,n+1)}$ acting on the Hilbert space of two
neighboring sites $n$ and $n+1$. The two half blue boxes at the boundaries
constitute a random CUE matrix $w^{(L,1)}$ (identical matrix) acting on the
sites $L$ and $1$, for periodic (open) boundary condition, with $L$ being the
total number of qudits. For any of these Floquet quantum circuits, which
differ only by the ordering of local gates and are subject to either periodic
or open boundary condition, we prove that the statistical properties of the
quasi-energy spectra are identical to that of the CUE ensemble in the limit of
large onsite Hilbert space dimension $q\rightarrow\infty$. Figure 2:  Floquet
operators of 2D Floquet random quantum circuits with different orderings of
the two-qudit gates. Within one period, two-qudit unitaries (orange boxes)
drawn independently from the CUE ensemble are applied to all pairs of
neighboring qudits (black dots) in a 2D lattice at various time substeps
(layers). Each qudit is coupled to all of its neighbors at different substeps.
The half box at the boundary together with its neighbor on the opposite side
combine to give a random unitary (identical matrix) acting on the
corresponding pair of qudits, for periodic (open) boundary condition. We prove
that, as in the 1D case (Fig. 1), the quasi-energy spectrum for this type of
Floquet random quantum circuits with arbitrary ordering of the quantum gates
exhibits universal statistical properties described by the CUE ensemble, in
the limit of large onsite Hilbert space dimension $q\rightarrow\infty$, for
both periodic and open boundary conditions. This statement extends to
arbitrary dimensions (see appendix B).

The field theoretical approach is not only useful in the study of the spectral
statistics of Floquet random quantum circuits, but can also be applied to the
analysis of more generic properties of a wide class of random quantum circuits
(which are not necessarily time periodic). To show this, within the same sigma
model framework, we rederive the known results for Weingarten calculus, which
has extensive applications in studies of quantum circuits (see for example
Refs. Nahum-2018 ; Sondhi ; Nahum-2017 ; HJ-2018 ; Huse ; Ludwig-1 ; Chalker-1
). When considering a quantum circuit where all or part of the quantum gates
are random unitaries drawn from some ensembles with Haar probability, one
usually encounters polynomial functions of the matrix elements integrated over
the relevant group with Haar measure (or equivalently averaged over the
relevant ensemble). This type of Haar integrals can be computed using the
Weingarten calculus Samuel ; Weingarten ; collins2003 ; collins2006 ;
collins2021 ; Weingarten-rev ; Beenakker . In the present paper, we rederive
the results for the Haar integrals of products of the matrix elements of a
unitary matrix over the unitary group (i.e., moments of CUE random matrices)
using the field theoretical approach. The derivation can be generalized to
other compact Lie groups or compact symmetric spaces collins2006 ; collins2009
; Wg-COE ; Wg-SymSpace .

The rest of the paper is organized as follows: In section 2, we introduce a
standard diagnostic of the statistics of quasi-energies of time-periodic
systems and briefly review the spectral statistical properties of the CUE
ensemble. In section 3, we present an effective field theory which is
formulated as a replica sigma model and can be used to investigate the
spectral statistics of a wide class of Floquet quantum systems. This effective
field theory is applied in section 4 to study a family of Floquet random
quantum circuits related to the brickwork circuit in Ref. Chalker-1  by a
reordering of the Haar random unitary gates (which couple the neighboring
qudits) and by the higher-dimensional generalization. We show that the
statistical properties of the quasi-energy spectra of the Floquet random
quantum circuits of this type are universally described by the CUE ensemble,
for both the periodic and open boundary conditions. Using the same sigma model
method, we rederive the results for the Weingarten calculus for the unitary
group in section 5. Finally, in section 6 we conclude with a discussion of
directions for future studies. The appendices are devoted to technical
details. In appendices A and B, we provide the calculation of the second- and
fourth- order moments of the Floquet operator for the family of Floquet random
quantum circuits studied in the current paper, in 1D and in higher dimensions,
respectively. We prove that, for arbitrary ordering of the two-qudit unitary
gates, and in any dimension, these moments are identical to that of the CUE
ensemble. In appendix C, we study a noninteracting Floquet model and examine
the higher oder fluctuations in the effective field theory. We show that the
quartic order fluctuations of this noninteracting model give rise to a larger
contribution to the level correlation function, compared with their chaotic
counterparts, and are therefore no longer negligible. In appendix D, we
rederive the asymptotic behavior of the Weingarten function for a unitary
group of dimension $q$ in the large $q\rightarrow\infty$ limit. Appendix E
contains a derivation of the recursion relation for the Weingarten function.

##  2 Diagnostic of spectral statistics of time-periodic systems

For time-periodic (Floquet) systems, the statistics of the quasi-energy
spectra serve as a diagnostic of quantum chaos. The quasi-energies
$\left\\{\theta_{i}\right\\}$ are the eigen-phases of the time evolution
operator over one period $U$ ( i.e., the Floquet operator):
$U\left|n\right\rangle=e^{i\theta_{n}}\left|n\right\rangle$. The quasi-energy
density can be expressed as

$\displaystyle\rho(\phi)=\sum_{i=1}^{N}\delta_{2\pi}(\phi-\theta_{i}),$ (1)

where $N$ is the Hilbert space dimension, and $\delta_{2\pi}(\phi)$ represents
the $2\pi$-periodic delta function defined as
$\delta_{2\pi}(\phi)=\sum_{n=-\infty}^{\infty}e^{-in\phi}/2\pi$. For brevity,
the subscript $2\pi$ will be omitted from now on.

For an ensemble of random Floquet systems, the two-point correlation function
of the quasi-energy density is usually introduced as:

$\displaystyle\begin{aligned}
R_{2}(\phi_{1},\phi_{2})=&\left\langle\rho(\phi_{1})\rho(\phi_{2})\right\rangle,\end{aligned}$
(2)

where the angular bracket denotes the ensemble averaging. The correlation
function, $R_{2}(\phi_{1},\phi_{2})$, measures the probability of finding two
quasi-energies at $\phi_{1}$ and $\phi_{2}$, and is one of the widely used
probes of the spectral statistics. In the present paper, in addition to the
ensemble averaging, we also perform the averaging over the entire quasi-energy
spectrum:

$\displaystyle\begin{aligned}
\bar{R}_{2}(\Delta\phi)=&\int_{0}^{2\pi}\frac{d\phi_{0}}{2\pi}R_{2}(\phi_{0}+\Delta\phi/2,\phi_{0}-\Delta\phi/2).\end{aligned}$
(3)

This averaging Altland-rev  is similar to the energy integration employed in
Refs. Andreev-Ballistic ; Andreev-Ballistic-2  to derive a ballistic sigma
model for individual Hamiltonian systems, and allows us to extract information
about the correlation function of an arbitrary pair of quasi-energy levels
with the separation $\Delta\phi$, irrespective of their positions in the
quasi-energy spectrum on the Floquet circle (see also, Ref. Zoller ).

The Fourier transform of the two-level correlation function is known as the
spectral form factor:

$\displaystyle
K(t)=\int_{0}^{2\pi}d\phi_{1}\int_{0}^{2\pi}d\phi_{2}R_{2}(\phi_{1},\phi_{2})e^{-i(\phi_{1}-\phi_{2})t}=2\pi\int_{0}^{2\pi}d(\Delta\phi)\,\bar{R}_{2}(\Delta\phi)e^{-i\Delta\phi\,t},$
(4)

which is equivalent to

$\displaystyle
K(t)=\left\langle\left|\operatorname{Tr}\left(U^{t}\right)\right|^{2}\right\rangle=\left\langle\sum_{i,j=1}^{N}e^{-it(\theta_{i}-\theta_{j})}\right\rangle.$
(5)

From this definition, one can immediately see that $K(t)=K(-t)$ and
$K(t=0)=N^{2}$. Moreover, assuming no degeneracy in the quasi-energy spectrum,
in the large time limit $t\rightarrow\infty$, the off-diagonal ($i\neq j$)
terms in the summation in Eq. (5) vanish upon ensemble averaging due to the
random phases, and the SFF acquires the value of $K(t\rightarrow\infty)=N$
(plateau).

Chaotic Floquet systems without the time-reversal invariance are expected to
exhibit the same spectral statistics properties as that of the Circular
Unitary Ensemble (CUE) Dyson-I ; Haake ; Mehta , which is an ensemble of
unitary matrices with Haar probability measure. For the CUE ensemble, the
spectrum is statistically homogeneous, and the mean level density is given by
$\bar{\rho}=N/2\pi$. The two-level correlation function
$R_{2}(\phi_{1},\phi_{2})$ depends only on the level separation
$\Delta\phi=\phi_{1}-\phi_{2}$ and therefore is equal to its energy-averaged
value Mehta :

$\displaystyle\begin{aligned}
R_{2}(\phi_{1},\phi_{2})=\bar{R}_{2}(\Delta\phi)=-\frac{1}{4\pi^{2}}\dfrac{\sin^{2}({N\Delta\phi/2})}{\sin^{2}\left(\Delta\phi/2\right)}+\frac{N^{2}}{4\pi^{2}}+\frac{N}{2\pi}\delta(\Delta\phi).\end{aligned}$
(6)

At nonzero $\Delta\phi$, the first term in the equation above corresponds to
the connected part of the two-level correlation function defined as

$\displaystyle R^{\rm
con}_{2}(\phi_{1},\phi_{2})=\left\langle\rho(\phi_{1})\rho(\phi_{2})\right\rangle-\left\langle\rho(\phi_{1})\right\rangle\left\langle\rho(\phi_{2})\right\rangle.$
(7)

In the large $N$ limit, after rescaling the quasi-energy
$\varepsilon=\bar{\rho}\Delta\phi$ and keeping it finite, the connected two-
level correlation function assumes the form identical to that of the Gaussian
Unitary Ensemble (GUE) Mehta ; Haake :

$\displaystyle\begin{aligned} \frac{1}{\bar{\rho}^{2}}R_{2}^{\rm
con}(\Delta\phi)=-\left(\dfrac{\sin(\pi\varepsilon)}{\pi\varepsilon}\right)^{2}.\end{aligned}$
(8)

From Eq. (6), one finds that the SFF of the CUE ensemble acquires the form

$\displaystyle\begin{aligned}
K(t)=\min(|t|,N)+N^{2}\delta_{t,0}.\end{aligned}$ (9)

It exhibits a linear ramp until a plateau sets in at $t=N$. This linear ramp
reflects the repulsion between the quasi-energies and is expected to be a
universal feature of the quantum chaotic systems with broken time-reversal
symmetry.

##  3 Replica sigma model for generic Floquet systems

In this section, we reformulate the sigma model initially constructed in Refs.
Zirnbauer_1996 ; Zirnbauer_1998 ; Zirnbauer_1999 ; QKR  using the replica
trick instead of the supersymmetric method, and obtain an ensemble averaged
effective field theory from which one can extract information about the
statistical properties of the quasi-energy spectra for generic time-periodic
quantum systems.

###  3.1 Generating function for level correlation function

Consider now an ensemble of Floquet systems whose Floquet operators $U$ are
from an arbitrary ensemble of unitary matrices. Note that, for many-body
systems, $U$ represents the many-body Floquet operator that acts in the many-
body Hilbert space. The statistical properties of the quasi-energies of this
ensemble can be obtained from the following generating function using the
replica trick

$\displaystyle\begin{aligned}
\mathcal{Z}^{(R)}(\alpha,\beta)=&\left\langle\int_{0}^{2\pi}\frac{d\phi}{2\pi}\left[\det\left(1-\alpha
e^{i\phi}U\right)\det\left(1-\beta
e^{-i\phi}U^{\dagger}\right)\right]^{R}\right\rangle.\end{aligned}$ (10)

Here $\alpha$ and $\beta$ are two complex numbers, and $R$ represents the
replica number which will be set to zero (the replica limit) at the end of the
calculation. Taking derivatives of the generating function
$\mathcal{Z}^{(R)}(\alpha,\beta)$ with respect to the complex variables
$\alpha$ and $\beta$, multiplying the result by $\alpha\beta/R^{2}$ and then
taking the replica limit $R\rightarrow 0$, we find

$\displaystyle\begin{aligned}
\mathcal{C}(\alpha,\beta)\equiv&\lim_{R\rightarrow
0}\frac{\alpha\beta}{R^{2}}\dfrac{\partial^{2}\mathcal{Z}^{(R)}(\alpha,\beta)}{\partial\alpha\partial\beta}=\left\langle\int_{0}^{2\pi}\frac{d\phi}{2\pi}\operatorname{Tr}\left(\dfrac{\alpha
e^{i\phi}U}{1-\alpha e^{i\phi}U}\right)\operatorname{Tr}\left(\dfrac{\beta
e^{-i\phi}U^{\dagger}}{1-\beta e^{-i\phi}U^{\dagger}}\right)\right\rangle\\\
=&\left\langle\int_{0}^{2\pi}\frac{d\phi}{2\pi}\sum_{n,m=1}^{\infty}\operatorname{Tr}\left(\alpha
e^{i\phi}U\right)^{n}\operatorname{Tr}\left(\beta
e^{-i\phi}U^{\dagger}\right)^{m}\right\rangle\\\
=&\sum_{n=1}^{\infty}K(n)(\alpha\beta)^{n}.\end{aligned}$ (11)

In this equation, the powers of $U$ and $U^{\dagger}$ (i.e., $n$ and $m$) need
to be identical to have a nonvanishing contribution to
$\mathcal{C}(\alpha,\beta)$ after the integration over $\phi$.

The two-level correlation function $\bar{R}_{2}(\Delta\phi)$ can be obtained
from $\mathcal{C}(\alpha,\beta)$, which is a weighted summation of the SFF
$K(n)$ (Eq. (5)) at discrete time $n$, by setting
$\alpha=\beta=e^{i\Delta\phi/2}$,

$\displaystyle\begin{aligned}
\bar{R}_{2}(\Delta\phi)=\frac{1}{(2\pi)^{2}}\left[2\operatorname{Re}\mathcal{C}(\alpha=\beta=e^{i\Delta\phi/2})+N^{2}\right].\end{aligned}$
(12)

Here we have used the fact that $K(-n)=K(n)$ and $K(0)=N^{2}$, with $N$ being
the dimension of the Hilbert space. We note that the higher-order correlation
functions of the quasi-energy density can be evaluated in an analogous manner
making use of a similar generating function Zirnbauer_1996 .

###  3.2 Replica sigma model for level correlation function

The generating function defined in Eq. (10) can be expressed as a Grassmann
path integral:

$\displaystyle\begin{aligned}
\mathcal{Z}^{(R)}(\alpha,\beta)=&\left\langle\int_{0}^{2\pi}\frac{d\phi}{2\pi}\int\mathcal{D}(\bar{\psi},\psi)\exp\left\\{-\bar{\psi}^{+,u}_{i}\left(\delta_{ij}-\alpha
e^{i\phi}U_{ij}\right)\psi^{+,u}_{j}\right.\right.\\\
&\left.\left.\qquad-\bar{\psi}^{-,u}_{i}\left(\delta_{ij}-\beta
e^{-i\phi}U^{\dagger}_{ij}\right)\psi^{-,u}_{j}\right\\}\right\rangle.\end{aligned}$
(13)

The Grassmann $\psi^{s,u}_{i}$ carries three different indices: $u=1,2,...,R$
is the replica index, $i=1,2,...,N$ labels the Hilbert space, and $s=+/-$
distinguishes the contributions from the forward ($U$) and backward
($U^{\dagger}$) evolution operators. Throughout the paper, we employ the
convention that repeated indices imply the summation.

Using the color flavor transformation Zirnbauer_1996 ; Zirnbauer_1998 ;
Zirnbauer_1999 ; Zirnbauer_2021 , the integration over the center phase $\phi$
in Eq. (13) can be converted into an integration over $RN\times RN$ complex
matrix field $Z$: 111Here we call matrix $Z$ a ‘field’, but note that it does
not have time dependence and carries indices in the Hilbert and replica spaces
only.

$\displaystyle\begin{aligned}
&\mathcal{Z}^{(R)}(\alpha,\beta)=c_{0}\left\langle\int\mathcal{D}(\bar{\psi},\psi)\exp\left\\{-\bar{\psi}^{+,u}_{i}\psi^{+,u}_{i}-\bar{\psi}^{-,u}_{i}\psi^{-,u}_{i}\right\\}\right.\\\
&\times\left.\int\mathcal{D}(Z,Z^{\dagger})\det(1+Z^{\dagger}Z)^{-2RN-1}\exp\left\\{\alpha\bar{\psi}^{+,u}_{i}Z^{uv}_{ij}U^{\dagger}_{jj^{\prime}}\psi^{-,v}_{j^{\prime}}-\beta\bar{\psi}^{-,u}_{i}(Z^{\dagger})^{uv}_{ij}U_{jj^{\prime}}\psi^{+,v}_{j^{\prime}}\right\\}\right\rangle\\\
=&c_{0}\left\langle\int\mathcal{D}(Z,Z^{\dagger})\det(1+Z^{\dagger}Z)^{-2RN-1}\right.\\\
&\qquad\times\left.\int\mathcal{D}(\bar{\psi},\psi)\exp\left\\{-\begin{bmatrix}\bar{\psi}^{+}&\bar{\psi}^{-}\end{bmatrix}\begin{bmatrix}1&-\alpha
ZU^{\dagger}\\\ \beta Z^{\dagger}U&1\end{bmatrix}\begin{bmatrix}\psi^{+}\\\
\psi^{-}\end{bmatrix}\right\\}\right\rangle.\end{aligned}$ (14)

Here $Z^{uv}_{ij}$ is a complex matrix which carries indices in both the
replica space (labeled by $u,v$) and the Hilbert space (labeled by $i,j$).
$c_{0}$ is an indefinite unessential normalization factor, which can be
determined in the $R\rightarrow 0$ limit.

Performing the Gaussian integration over $\psi$, we arrive at a sigma model
representation of the generating function

$\displaystyle\begin{aligned}
&\mathcal{Z}^{(R)}(\alpha,\beta)=c_{0}\left\langle\int\mathcal{D}(Z,Z^{\dagger})e^{-S[Z^{\dagger},Z]}\right\rangle,\\\
&S[Z^{\dagger},Z]=(2RN+1)\operatorname{Tr}\ln\left(1+Z^{\dagger}Z\right)-\operatorname{Tr}\ln\left(1+\alpha\beta
ZU^{\dagger}Z^{\dagger}U\right).\end{aligned}$ (15)

It is sometimes convenient to make the transformation

$\displaystyle\begin{aligned} &Q=T\Lambda T^{-1},\qquad
T=\begin{bmatrix}1&-Z\\\
Z^{\dagger}&1\end{bmatrix},\qquad\Lambda=\begin{bmatrix}1&0\\\
0&-1\end{bmatrix},\end{aligned}$ (16)

after which the action becomes

$\displaystyle\begin{aligned}
&S[Q]=-2RN\operatorname{Tr}\ln\left(\frac{Q\Lambda+1}{2}\right)-\operatorname{Tr}\ln\left[\frac{1}{2}\left(1-\begin{bmatrix}\alpha
U&0\\\ 0&\beta
U^{\dagger}\end{bmatrix}\right)Q\Lambda+\frac{1}{2}\left(1+\begin{bmatrix}\alpha
U&0\\\ 0&\beta U^{\dagger}\end{bmatrix}\right)\right].\end{aligned}$ (17)

Similar to the sigma model for disordered systems Wegner-1979 ; Efetov ;
Kamenev-book , the matrix field $Q$ stays on the manifold with the constraints
$\operatorname{Tr}Q=0$ and $Q^{2}=1$. Eq. (15) (or equivalently Eq. (17)) is a
replica version of the supersymmetric sigma model derived earlier to study the
spectral statistics of the circular ensembles and quantum chaotic maps
Zirnbauer_1996 ; Zirnbauer_1998 ; Zirnbauer_1999 ; Altland-rev .

###  3.3 Ensemble averaged effective theory

Starting from Eq. (15), we then perform the ensemble averaging and derive an
effective field theory for the matrix field $Z$ for any ensemble of Floquet
systems. Note that Eq. (15) can be rewritten as

$\displaystyle\begin{aligned}
\mathcal{Z}^{(R)}(\alpha,\beta)=\,&c_{0}\int\mathcal{D}(Z,Z^{\dagger})e^{-S_{\rm
eff}[Z^{\dagger},Z]},\\\ S_{\rm
eff}[Z^{\dagger},Z]\equiv\,&(2RN+1)\operatorname{Tr}\ln\left(1+Z^{\dagger}Z\right)-\ln\left\langle\exp\left[\operatorname{Tr}\ln\left(1+\alpha\beta
ZU^{\dagger}Z^{\dagger}U\right)\right]\right\rangle.\end{aligned}$ (18)

Until now, no approximations have been made and the expression above is exact.
However, it is difficult to perform the ensemble averaging in the second term
in the action $S_{\rm eff}$, especially for many-body systems whose Floquet
operator has a complicated structure in the many-body Hilbert space. To
proceed, we expand $S_{\rm eff}$ around the saddle point $Z=0$ in powers of
fluctuation $Z$ and carry out the ensemble averaging term by term. Up to
quartic order in $Z$, the action $S_{\rm eff}$ is given by

$\displaystyle S_{\rm eff}[Z^{\dagger},Z]=S_{\rm
eff}^{(2)}[Z^{\dagger},Z]+S_{\rm eff}^{(4)}[Z^{\dagger},Z]+O(Z^{6}),$ (19a)
$\displaystyle\begin{aligned} S_{\rm
eff}^{(2)}[Z^{\dagger},Z]=&\sum_{i_{1},...,i_{4}}\operatorname{tr}\left(Z^{\dagger}_{i_{1}i_{2}}Z_{i_{3}i_{4}}\right)\left[(2RN+1)\delta_{i_{2},i_{3}}\delta_{i_{4},i_{1}}-\alpha\beta\left\langle
U_{i_{2}i_{3}}U^{\dagger}_{i_{4}i_{1}}\right\rangle\right],\end{aligned}$
(19b) $\displaystyle\begin{aligned} S_{\rm
eff}^{(4)}[Z^{\dagger},Z]=&-\frac{(\alpha\beta)^{2}}{2}\sum_{i_{1},...,i_{8}}\operatorname{tr}\left(Z^{\dagger}_{i_{1}i_{2}}Z_{i_{3}i_{4}}\right)\operatorname{tr}\left(Z^{\dagger}_{i_{5}i_{6}}Z_{i_{7}i_{8}}\right)\\\
&\times\left[\left\langle
U_{i_{2}i_{3}}U^{\dagger}_{i_{4}i_{1}}U_{i_{6}i_{7}}U^{\dagger}_{i_{8}i_{5}}\right\rangle-\left\langle
U_{i_{2}i_{3}}U^{\dagger}_{i_{4}i_{1}}\right\rangle\left\langle
U_{i_{6}i_{7}}U^{\dagger}_{i_{8}i_{5}}\right\rangle\right]\\\
&+\frac{1}{2}\sum_{i_{1},...,i_{8}}\operatorname{tr}\left(Z^{\dagger}_{i_{1}i_{2}}Z_{i_{3}i_{4}}Z^{\dagger}_{i_{5}i_{6}}Z_{i_{7}i_{8}}\right)\\\
&\times\left[-(2RN+1)\delta_{i_{2},i_{3}}\delta_{i_{4},i_{5}}\delta_{i_{6},i_{7}}\delta_{i_{8},i_{1}}+(\alpha\beta)^{2}\left\langle
U_{i_{2}i_{3}}U^{\dagger}_{i_{4}i_{5}}U_{i_{6}i_{7}}U^{\dagger}_{i_{8}i_{1}}\right\rangle\right].\end{aligned}$
(19c)

In the present paper, we use ‘$\operatorname{tr}$’ to denote the trace
operation that acts on the replica space only and ‘$\operatorname{Tr}$’ to
trace over both the replica and Hilbert spaces. We note that the last term in
$S_{\rm eff}^{(4)}$ (Eq. (19c)) gives rise to a contribution of higher order
in the replica number $R$, compared with the remaining term in $S_{\rm
eff}^{(4)}$, and is not important in the replica limit $R\rightarrow 0$.

The specific forms of the moments of the Floquet operator in the effective
theory are model dependent. Once they are known, one can insert them into the
expression for the effective action Eq. (19) to obtain an effective field
theory which encodes the information about the statistical properties of the
quasi-energy spectrum. This field theoretical approach is therefore applicable
to a wide class of Floquet systems, including many-body and single-particle
systems. One advantage of this method, compared with other model specific
ones, is that it can be used to investigate why or when universal statistical
behaviors emerge. For Floquet random quantum circuits consisting of
independent local unitary gates, the moments of the many-body Floquet operator
are given by products of moments of the local unitaries, and therefore are
usually not difficult to evaluate. This makes the current field theoretical
approach especially useful for random quantum circuits with local gates.

##  4 Application to Floquet random quantum circuits

###  4.1 Floquet random quantum circuits

As an example, we apply the effective field theory (Eq. (19)) derived in the
previous section to study the spectral statistics of a family of Floquet
quantum circuits composed of random local unitary gates. Figs. 1 and 2 show
the Floquet operators of some of these Floquet random quantum circuits,
including the brickwork circuit (Fig. 1(a)) studied earlier in Ref. Chalker-1
. The other two circuits depicted in panels (b) and (c) of Fig. 1 can be
obtained from the brickwork circuit by reordering the local gates. Some
examples of the 2D generalization of such Floquet circuits are depicted in
Fig. 2. We prove that Floquet random quantum circuits of this type, subject to
either periodic or open boundary condition, are all described by the same
effective field theory as that of the CUE ensemble in the limit of large on-
site Hilbert space dimension, irrespective of the ordering of the local gates
and the dimensionality of the lattice of qudits.

We first consider the Floquet quantum circuits consisting of a $1$D lattice
chain of $L$ qudits, each of which contains $q\rightarrow\infty$ internal
states. The dimension of many-body Hilbert space is $N=q^{L}$. The time
evolutions of these Floquet circuits are discrete and time periodic, and the
evolutions over one period contain $2\leq M\leq L$ substeps (where $M$ is
model specific). For all of these Floquet quantum circuits, each qudit
(labeled by an integer $n=1,2,...,L$) is coupled to its neighbors on the left
hand side (at site $n-1$) and the right hand side (at site $n+1$) at two
different substeps (labeled by integers $s^{(n-1,n)}$ and $s^{(n,n+1)}$,
respectively). We consider both the periodic and open boundary conditions, and
the qudit label $n$ is defined modulo $L$ for periodic boundary condition. The
local gate that couples a pair of neighboring qudits at sites $n$ and $n+1$ is
given by a $q^{2}\times q^{2}$ random CUE matrix $w^{(n,n+1)}$, and is
represented diagrammatically by a blue box in Fig. 1. Unitary gates acting on
different pairs of neighboring qudits are independent and uncorrelated. We
consider all possible orderings of these quantum gates represented by
different configurations of $\left\\{1\leq s^{(n,n+1)}\leq
M|n=1,2,...,L^{\prime}\right\\}$ with constraint $s^{(n,n+1)}\neq
s^{(n-1,n)}.$ Here $L^{\prime}$ is defined as $L^{\prime}=L$ for periodic
boundary condition and $L^{\prime}=L-1$ for open boundary condition. The total
number of substeps $M$ is given by the total number of different integers in
the set $\left\\{s^{(n,n+1)}\right\\}$, and we consider all possible values of
$2\leq M\leq L^{\prime}$. For the staircase circuit in Fig. 1(c),
$M=L^{\prime}$ and $s^{(n,n+1)}=n$, while for the brickwork circuit in Fig.
1(a), $M=2$ and $s^{(n,n+1)}=2-(n\,\mathrm{mod}\,2)$ (for open boundary
condition or periodic boundary condition with even $L$).

For any of the Floquet random quantum circuits described above, the Floquet
operator can be expressed as

$\displaystyle\begin{aligned}
&U=W^{(\sigma(1),\sigma(1)+1)}W^{(\sigma(2),\sigma(2)+1)}...W^{(\sigma(L^{\prime}-1),\sigma(L^{\prime}-1)+1)}W^{(\sigma(L^{\prime}),\sigma(L^{\prime})+1)},\end{aligned}$
(20)

where $\sigma\in S_{L^{\prime}}$ represents a permutation of numbers
$1,2,...,L^{\prime}$. $W$ is defined as

$\displaystyle W^{(n,n+1)}=w^{(n,n+1)}\otimes 1^{(n,n+1)}.$ (21)

$1^{(n,n+1)}$ represents an identical matrix operating in the Hilbert space of
all sites except for $n$ and $n+1$. As mentioned earlier, $w^{(n,n+1)}$ acts
on the qudits at sites $n$ and $n+1$, and is drawn randomly and independently
from the CUE ensemble of dimension $q^{2}$. We consider all possible quantum
circuits whose Floquet operator can be expressed in the form of Eq. (20) for
arbitrary permutation $\sigma\in S_{L^{\prime}}$. For example, $\sigma$ is the
identical permutation for the staircase circuit in Fig. 1(c), while for the
brickwork circuit in Fig. 1(a) with periodic boundary condition and even $L$,
$\sigma$ is given by

$\displaystyle\sigma(2k-1)=k,\qquad\sigma(2k)=L/2+k,\qquad 1\leq k\leq L/2.$
(22)

In appendix A, we prove that the Floquet operator $U$ given by Eq. (20) with
arbitrary permutation $\sigma\in S_{L^{\prime}}$ obeys the conditions:

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=\frac{1}{N}\delta_{\bm{\mathrm{i}}\bm{\mathrm{i}}^{\prime}}\delta_{\bm{\mathrm{j}}\bm{\mathrm{j}}^{\prime}},\end{aligned}$
(23a) $\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle=\frac{1}{N^{2}}\left(\delta_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{1}^{\prime}}\delta_{\bm{\mathrm{j}}_{1}\bm{\mathrm{j}}_{1}^{\prime}}\delta_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{2}^{\prime}}\delta_{\bm{\mathrm{j}}_{2}\bm{\mathrm{j}}_{2}^{\prime}}+\delta_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}^{\prime}}\delta_{\bm{\mathrm{j}}_{1}\bm{\mathrm{j}}_{2}^{\prime}}\delta_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}^{\prime}}\delta_{\bm{\mathrm{j}}_{2}\bm{\mathrm{j}}_{1}^{\prime}}\right).\end{aligned}$
(23b)

Here the $L$-dimensional vector
$\bm{\mathrm{i}}=(i^{(1)},i^{(2)},...,i^{(L)})$ labels the many-body state of
the circuits, and its $n$-th component $i^{(n)}=1,2,...,q$ indexes the single-
particle state of the $n$-th qudit. We note that Eq. (23a) holds for arbitrary
$q$ while Eq. (23b) is derived in the limit of $q\rightarrow\infty$.

The above discussion can be straightforwardly generalized to higher-
dimensional Floquet quantum circuits with similar configurations. In
particular, let us consider now a $D$-dimensional cubic lattice of qudits,
with $L$ sites in each direction. The single-particle (many-body) Hilbert
space dimension is $q\rightarrow\infty$ ($N=q^{L^{D}}$). The time evolution is
again discrete and periodic, and is composed of local two-qudit gates that
couple separately all pairs of neighboring qudits. During one period, each
qudit is coupled to all of its neighboring qudits at different substeps by
different local unitary gates drawn randomly and independently from the CUE
ensemble of dimension $q^{2}$. For this type of $D$-dimensional Floquet
circuits with any possible ordering of these local gates, and with either
periodic or open boundary condition, we prove that the Floquet operator $U$
still satisfies Eq. (23). In this case, $\bm{\mathrm{i}}$ can be considered as
a $L^{D}$ dimensional vector whose component $i^{(\bm{\mathrm{n}})}=1,2,...,q$
labels the single-particle state of the qudit $\bm{\mathrm{n}}$ in the
$D$-dimensional lattice. The derivation is relegated to appendix B.

Eqs. (23a) and (23b) are also obeyed if $U$ is drawn randomly from a CUE
ensemble of dimension $N\rightarrow\infty$. As a result, the effective field
theory (Eq. (19)) for the Floquet quantum circuits under consideration is
equivalent to that of the CUE ensemble of dimension $N$, and the two-level
correlation functions for the current models are given by the CUE level
correlation function in Eq. (6).

###  4.2 Quadratic fluctuations

In the following, we present the detailed derivation of the two-level
correlation function $\bar{R}_{2}(\Delta\phi)$ for the quasi-energies of the
Floquet random quantum circuits described above (or equivalently the CUE
ensemble of the same dimension $N$) using the effective field theory Eq. (19)
and the moments of the Floquet operator Eq. (23).

Substituting Eq. (23a) into the quadratic order effective action $S_{\rm
eff}^{(2)}$ in Eq. (19b), we obtain

$\displaystyle\begin{aligned} S_{\rm
eff}^{(2)}=&(2RN+1)\sum_{\bm{\mathrm{i}}_{1},\bm{\mathrm{i}}_{2}}\operatorname{tr}\left(Z^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}Z_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}\right)-\frac{\alpha\beta}{N}\sum_{\bm{\mathrm{i}}_{1},\bm{\mathrm{i}}_{2}}\operatorname{tr}\left(Z^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{1}}Z_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{2}}\right).\end{aligned}$
(24)

We now divide $Z$ into the diagonal component
$X_{\bm{\mathrm{i}}\bm{\mathrm{i}}}\equiv Z_{\bm{\mathrm{i}}\bm{\mathrm{i}}}$
and the off-diagonal component $Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}\equiv
Z_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$ (for $\bm{\mathrm{i}}\neq\bm{\mathrm{j}}$)
in the Hilbert space, and Fourier transform the diagonal component
$X_{\bm{\mathrm{j}}\bm{\mathrm{j}}}$ with respect to $\bm{\mathrm{j}}$:

$\displaystyle\begin{aligned}
&X(\bm{\mathrm{k}})=\sum_{\bm{\mathrm{j}}}X_{\bm{\mathrm{j}}\bm{\mathrm{j}}}e^{-i2\pi\bm{\mathrm{k}}\cdot\bm{\mathrm{j}}/q}.\end{aligned}$
(25)

Here $\bm{\mathrm{k}}$ is an $L^{D}$ dimensional vector and the summation over
each component runs over $k^{(\bm{\mathrm{n}})}=0,1,...,q-1$. In terms of
$X(\bm{\mathrm{k}})$ and $Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$, the effective
action $S_{\rm eff}^{(2)}$ can be rewritten as

$\displaystyle\begin{aligned} S_{\rm
eff}^{(2)}=\,&(2RN+1)\frac{1}{N}\sum_{\bm{\mathrm{k}}}\operatorname{tr}\left(X^{\dagger}(\bm{\mathrm{k}})X(\bm{\mathrm{k}})\right)-\frac{\alpha\beta}{N}\operatorname{tr}\left(X^{\dagger}(0)X(0)\right)\\\
&+(2RN+1)\sum_{\bm{\mathrm{i}}_{1}\neq\bm{\mathrm{i}}_{2}}\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}Y_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}\right).\end{aligned}$
(26)

From the equation above, one finds that the bare propagator for
$X(\bm{\mathrm{k}}=0)$ acquires the form

$\displaystyle\begin{aligned} \left\langle
X^{uv}(0)(X^{\dagger})^{v^{\prime}u^{\prime}}(0)\right\rangle_{0}=\,\delta_{uu^{\prime}}\delta_{vv^{\prime}}\dfrac{N}{(2RN+1)-\alpha\beta}.\end{aligned}$
(27)

Here the angular bracket with subscript $0$ represents the averaging over the
Gaussian fluctuation of matrix $Z$ governed by the action $S_{\rm eff}^{(2)}$
(Eq. (26)). Taking the replica limit $R\rightarrow 0$ and setting
$\alpha=\beta$ to $e^{i\Delta\phi/2}$, the bare propagator for $X(0)$ becomes

$\displaystyle\begin{aligned} \left\langle
X^{uv}(0)(X^{\dagger})^{v^{\prime}u^{\prime}}(0)\right\rangle_{0}=\,\delta_{uu^{\prime}}\delta_{vv^{\prime}}\dfrac{N}{1-e^{i\Delta\phi}},\end{aligned}$
(28)

which diverges when $\Delta\phi\rightarrow 0$. The corresponding mode is
massless.

By contrast, the propagators for $X(\bm{\mathrm{k}}\neq 0)$ and
$Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$ are given by, respectively

$\displaystyle\begin{aligned} \left\langle
X^{uv}(\bm{\mathrm{k}})(X^{\dagger})^{v^{\prime}u^{\prime}}(\bm{\mathrm{k}}^{\prime})\right\rangle_{0}=\,&\delta_{uu^{\prime}}\delta_{vv^{\prime}}\delta_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\dfrac{N}{2RN+1},\\\
\left\langle
Y^{uv}_{\bm{\mathrm{i}}\bm{\mathrm{j}}}(Y^{\dagger})^{v^{\prime}u^{\prime}}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle_{0}=\,&\delta_{uu^{\prime}}\delta_{vv^{\prime}}\delta_{\bm{\mathrm{i}}\bm{\mathrm{i}}^{\prime}}\delta_{\bm{\mathrm{j}}\bm{\mathrm{j}}^{\prime}}\frac{1}{2RN+1}.\end{aligned}$
(29)

$X(\bm{\mathrm{k}}\neq 0)$ and $Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$ are
therefore massive, and give rise to a $\alpha$, $\beta$ independent
contribution to the generating function $\mathcal{Z}^{(R)}$ at the quadratic
order.

In summary, the $Z$ fluctuations can be divided into two categories: the
massless fluctuation $X(0)$ and the massive fluctuations
$X(\bm{\mathrm{k}}\neq 0)$ and $Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$. The
massive modes $X(\bm{\mathrm{k}}\neq 0)$ and
$Y_{\bm{\mathrm{i}}\bm{\mathrm{j}}}$ contribute a nonessential constant to the
generating function, while the soft mode $X(0)$ governs the spectral
statistics PRL ; Winer .

###  4.3 Quartic fluctuations

We now investigate the contribution to the self-energy from the quartic order
fluctuations (Eq. (19c)). Note that the last term in $S_{\rm eff}^{(4)}$ (Eq.
(19c)) does not contribute in the replica limit $R\rightarrow 0$ due to its
special structure in the replica space. We can therefore focus on the first
term in $S_{\rm eff}^{(4)}$ which can be expressed in terms of
$X(\bm{\mathrm{k}})$ and $Y$ as

$\displaystyle\begin{aligned} &S_{\rm
eff}^{(4-1)}=-\frac{(\alpha\beta)^{2}}{2N^{4}}\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{2}}\operatorname{tr}\left(X^{\dagger}(\bm{\mathrm{k}}_{1})X(\bm{\mathrm{k}}_{2})\right)\operatorname{tr}\left(X^{\dagger}(-\bm{\mathrm{k}}_{1})X(-\bm{\mathrm{k}}_{2})\right)\\\
&-\frac{(\alpha\beta)^{2}}{2N^{2}}\sum_{\bm{\mathrm{i}}_{1}\neq\bm{\mathrm{i}}_{2},\bm{\mathrm{i}}_{3}\neq\bm{\mathrm{i}}_{4}}\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}Y_{\bm{\mathrm{i}}_{3}\bm{\mathrm{i}}_{4}}\right)\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}Y_{\bm{\mathrm{i}}_{4}\bm{\mathrm{i}}_{3}}\right)\\\
&-\frac{(\alpha\beta)^{2}}{2N^{3}}\sum_{\bm{\mathrm{i}}_{1}\neq\bm{\mathrm{i}}_{2}}\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\left[\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}X(\bm{\mathrm{k}})\right)\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}X(-\bm{\mathrm{k}})\right)+\operatorname{tr}\left(X^{\dagger}(-\bm{\mathrm{k}})Y_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}\right)\operatorname{tr}\left(X^{\dagger}(\bm{\mathrm{k}})Y_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}\right)\right].\end{aligned}$
(30)

Here we have used Eq. (23b). Comparing with the quadratic action $S_{\rm
eff}^{(2)}$ in Eq. (26), one can see that the quartic action $S_{\rm
eff}^{(4-1)}$ is of higher order in $1/N$. Note that this counting does not
hold for nonergodic circuits and the higher order fluctuations become
important (see appendix C).

The self-energies for the $X$ and $Y$ components from the quartic interactions
are given by, respectively,

$\displaystyle\begin{aligned}
(\Sigma_{X})^{ab,ba}(\bm{\mathrm{k}},\bm{\mathrm{k}})=&\frac{(\alpha\beta)^{2}}{N^{4}}\left\langle
X^{ab}(-\bm{\mathrm{k}})(X^{\dagger})^{ba}(-\bm{\mathrm{k}})\right\rangle_{0}=\frac{(\alpha\beta)^{2}}{N^{3}}\frac{1}{1-\alpha\beta\delta_{\bm{\mathrm{k}},0}},\\\
(\Sigma_{Y})^{ab,ba}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2},\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}=&\frac{(\alpha\beta)^{2}}{N^{2}}\left\langle
Y_{\bm{\mathrm{i}}_{2}\bm{\mathrm{i}}_{1}}^{ab}(Y^{\dagger})^{ba}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}\right\rangle_{0}=\frac{(\alpha\beta)^{2}}{N^{2}}.\end{aligned}$
(31)

Here we have taken the replica limit and ignored the self-energy corrections
that vanish in the limit $R\rightarrow 0$. From the equation above, we can see
that the self-energy from the quartic order fluctuations is negligible in the
large $N$ limit for the massive modes $X(\bm{\mathrm{k}}\neq 0)$ and $Y$. For
the massless mode $X(0)$, the self-energy from the quartic interactions can be
ignored if we consider a not too small energy separation $\Delta\phi\gg 1/N$.
To recover the fine structure of the nearby quasi-energy levels, higher order
fluctuations of $X(0)$ are needed.

We note that, for integrable systems, fluctuations beyond the quadratic order
are no longer negligible, even when the quasi-energy separation $\Delta\phi$
being probed is much larger than the mean level spacing. In appendix C, we
consider a noninteracting Floquet model whose single-particle dynamics within
one period is generated by random CUE matrices. In particular, the single-
particle Floquet operator for each particle is independently drawn from the
CUE ensemble. We find that the action for the quadratic fluctuations of this
noninteracting model is identical to that of the Floquet random quantum
circuits considered in this section (or equivalently the CUE ensemble), and is
given by Eq. (26). However, the quartic fluctuations are governed by a
different action which, compared with its chaotic counterpart Eq. (30), leads
to a much larger contribution to the self-energy of the $Z$ matrix field.
Unlike the chaotic model, the higher order fluctuations become important for
this noninteracting model.

###  4.4 Two-level correlation function

For the chaotic Floquet random quantum circuits under consideration here, we
can focus on the quadratic fluctuation’s contribution and neglect the higher
order corrections. $\mathcal{C}(\alpha,\beta)$ defined in Eq. (11) is then
approximately given by

$\displaystyle\begin{aligned} \mathcal{C}(\alpha,\beta)=&\lim_{R\rightarrow
0}\frac{\alpha\beta}{R^{2}}c_{1}\left\langle\frac{1}{N}\operatorname{tr}\left(X^{\dagger}(0)X(0)\right)+\frac{\alpha\beta}{N^{2}}\operatorname{tr}^{2}\left(X^{\dagger}(0)X(0)\right)\right\rangle_{0}=\frac{\alpha\beta}{(1-\alpha\beta)^{2}}.\end{aligned}$
(32)

Here the overall coefficient $c_{1}$ contains the unessential contribution
from the massive modes and its value in the replica limit is determined from
the fact that $\lim_{R\rightarrow 0}Z^{(R)}(\alpha,\beta)=1$. Using Eq. (12),
we obtain the result for the two-level correlation function

$\displaystyle\begin{aligned}
\bar{R}_{2}(\Delta\phi)=-\frac{1}{8\pi^{2}}\dfrac{1}{\sin^{2}\left(\Delta\phi/2\right)}+\frac{N^{2}}{4\pi^{2}}.\end{aligned}$
(33)

It is easy to see that $\left\langle\operatorname{Tr}U\right\rangle=0$ for the
Floquet circuits under consideration. Therefore the average quasi-energy
density for any of these circuits is homogeneous and given by
$\bar{\rho}=N/2\pi$. This also means that the first term in Eq. (33)
corresponds to the connected part of the two-level correlation function
$R_{2}^{\rm con}(\Delta\phi\neq 0)$ defined in Eq. (7). We emphasize that this
equation applies to all of the Floquet quantum circuits under consideration.

Comparing Eq. (33) with the exact CUE result Eq. (6), one can see that the
smooth part of $\bar{R}_{2}(\Delta\phi)$ is recovered while an oscillatory
term proportional to $\cos(N\Delta\phi)$ is missing. To recover the
oscillatory term or to extract the behavior of $\bar{R}_{2}(\Delta\phi)$ at
small energy separation $\Delta\phi\lesssim 1/N$, nonperturbative information
about higher order fluctuations is needed. This may be obtained by
consideration of nonstandard saddle points (see Ref. Altland-rev ), similar to
the calculation of the level correlation function for the Hamiltonian systems
described by the Gaussian ensembles Andreev-Altshuler ; Kamenev-GUE ; Kamenev
; Kamenev-Keldysh .

##  5 Weingarten calculus

To show that the field theoretical approach described above provides access
not only to the spectral statistics but also to other generic properties of
quantum circuits, in this section, we rederive the known results for the
Weingarten calculus, which has been employed extensively in studies of quantum
circuits. Let us now consider a Haar integral of a product of the matrix
elements of a unitary matrix $U$ over the unitary group in $q$ dimensions
$\mathcal{U}(q)$:

$\displaystyle\begin{aligned}
I=&\int_{\mathcal{U}(q)}dUU_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{p}j_{p}}U^{\dagger}_{j_{1}^{\prime}i_{1}^{\prime}}U^{\dagger}_{j_{2}^{\prime}i_{2}^{\prime}}...U^{\dagger}_{j^{\prime}_{p^{\prime}}i^{\prime}_{p^{\prime}}}\\\
=&\left\langle
U_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{p}j_{p}}U^{\dagger}_{j_{1}^{\prime}i_{1}^{\prime}}U^{\dagger}_{j_{2}^{\prime}i_{2}^{\prime}}...U^{\dagger}_{j^{\prime}_{p^{\prime}}i^{\prime}_{p^{\prime}}}\right\rangle_{\rm
CUE},\end{aligned}$ (34)

where $dU$ denotes the normalized ($\int_{\mathcal{U}(q)}dU=1$) Haar measure.
This integral can be interpreted as the product of the unitary matrix elements
averaged over the CUE ensemble, denoted by the angular bracket with subscript
CUE. In this section and appendices D and E, we use $U$ to denote a random CUE
matrix (i.e., a Haar distributed random unitary matrix).

It has been found that this integral $I$ can be expressed in terms of the
Weingarten function Samuel ; Weingarten ; collins2003 ; collins2006 ;
collins2021 ; Weingarten-rev ; Beenakker :

$\displaystyle\begin{aligned} I=\delta_{p,p^{\prime}}\sum_{\tau,\sigma\in
S_{p}}\mathrm{Wg}(\tau^{-1}\sigma)\prod_{k=1}^{p}\delta_{i_{k},i_{\sigma(k)}^{\prime}}\delta_{j_{k},j_{\tau(k)}^{\prime}}.\end{aligned}$
(35)

Here the double summation runs over all permutations $\tau,\sigma\in S_{p}$ of
the integers $1,2,...,p$, and the Weingarten function
$\mathrm{Wg}(\tau^{-1}\sigma)$ depends only on the cycle structure of the
product $\tau^{-1}\sigma$. More specifically, the permutation
$\tau^{-1}\sigma$ can be divided into $m$ disjoint cycles:

$\displaystyle\begin{aligned} \left(P_{1}^{(1)}\rightarrow
P_{2}^{(1)}\rightarrow...\rightarrow
P_{c_{1}}^{(1)}\right)\left(P_{1}^{(2)}\rightarrow
P_{2}^{(2)}\rightarrow...\rightarrow
P_{c_{2}}^{(2)}\right)..\left(P_{1}^{(m)}\rightarrow
P_{2}^{(m)}\rightarrow...\rightarrow P_{c_{m}}^{(m)}\right).\end{aligned}$
(36)

Here $\left\\{P_{k}^{(l)}\right\\}$ are different integers from the set of
integers under the permutation $\left\\{1,2,...,p\right\\}$, and satisfy
$P_{k}^{(l)}=(\tau^{-1}\sigma)^{k-1}(P_{1}^{(l)})$ (for $k\leq c_{l}$) and
$P_{1}^{(l)}=(\tau^{-1}\sigma)P_{c_{l}}^{(l)}$. The Weingarten function
$\mathrm{Wg}(\tau^{-1}\sigma)$ depends only on the lengths of these disjoint
cycles $\left\\{c_{k}|k=1,2,...,m\right\\}$, which obey the constraint
$\sum_{k=1}^{m}c_{k}=p$, irrespective of their order. In the following, we
will denote the Weingarten function of the permutation $\tau^{-1}\sigma$ given
by Eq. (36) as $\mathrm{Wg}(\tau^{-1}\sigma)=V^{(p)}_{c_{1},c_{2},...,c_{m}}$.

The Weingarten function can be uniquely determined by the recursions relation
Samuel ; Beenakker :

$\displaystyle\begin{aligned}
qV_{c_{1},...,c_{m},1}^{(p+1)}+\sum_{s=1}^{m}c_{s}V_{c_{1},...,c_{s-1},c_{s}+1,c_{s+1},...,c_{m}}^{(p+1)}=V_{c_{1},...,c_{m}}^{(p)},\end{aligned}$
(37a) $\displaystyle\begin{aligned}
qV_{c_{1},...,c_{m}}^{(p)}+\sum_{c=1}^{c_{1}-1}V_{c_{1}-c,c,c_{2},...,c_{m}}^{(p)}+\sum_{s=2}^{m}c_{s}V_{c_{1}+c_{s},c_{2},...,c_{s-1},c_{s+1},...,c_{m}}^{(p)}=0.\end{aligned}$
(37b)

Here $\left\\{c_{i}\right\\}$ represent an arbitrary set of $m$ positive
integers that satisfy $\sum_{k=1}^{m}c_{k}=p$ in Eq. (37a), and in Eq. (37b)
an additional constraint $c_{1}\geq 2$ is imposed.

In the large $q\rightarrow\infty$ limit Beenakker , the Weingarten function of
the identical permutation $id$ is given by
$\mathrm{Wg}(id)=V^{(p)}_{1,1,...,1}=q^{-p}+O(q^{-p-2})$. For all the
remaining permutations, the Weingarten function
$V^{(p)}_{c_{1},c_{2},...,c_{m}}=O(q^{m-2p})$ is of higher order in $1/q$
(since the number of disjoint cycles $m<p$). Using these results, one can
easily see that Eq. (23) is obeyed by the CUE random matrix $U$ of dimension
$N$ in the large $N$ limit.

These results for Weingarten calculus have been obtained earlier using various
methods Weingarten-rev ; Samuel ; collins2003 ; collins2006 ; Beenakker ; Zee
). In the following, we rederive these results using a field theoretical
method similar to the one employed earlier in the analysis of the statistics
of quasi-energy spectra. We believe that this field theoretical approach can
be straightforwardly generalized and serves as a general framework to study
the quantum circuits.

###  5.1 Sigma model derivation for the Weingarten calculus

To begin with, we rewrite the integral $I$ in Eq. (34) as a fermionic path
integral

$\displaystyle\begin{aligned}
I=&\int_{\mathcal{U}(q)}dU\int\mathcal{D}(\bar{\psi},\psi)e^{-S_{U}[\bar{\psi},\psi]}F[\bar{\psi},\psi],\\\
S_{U}[\bar{\psi},\psi]=&\sum_{u=1}^{R^{\prime}}\sum_{l,k=1}^{q}\left(\bar{\psi}^{+,u}_{k}U_{kl}\psi^{+,u}_{l}+\bar{\psi}^{-,u}_{k}U^{\dagger}_{kl}\psi^{-,u}_{l}\right),\\\
F[\bar{\psi},\psi]=&\prod_{k=1}^{p}\left(\psi^{-,k}_{i_{k}}\bar{\psi}^{-,k}_{j_{k}}\right)\prod_{l=1}^{p^{\prime}}\left(\psi^{+,l}_{j^{\prime}_{l}}\bar{\psi}^{+,l}_{i^{\prime}_{l}}\right).\end{aligned}$
(38)

Here the fermionic field $\psi^{s,u}_{k}$ carries three indices that label,
respectively, the component associated with $U/U^{\dagger}$ ($s=+/-$), the
replica space ($u=1,2,...,R^{\prime}$) and the Hilbert space in which the
unitary matrix $U$ acts ($k=1,2,...,q$). Note here the replica number
$R^{\prime}$ is an integer given by $R^{\prime}=\max(p,p^{\prime})$, and,
unlike in the earlier calculation of the level correlation function, it dose
not need to be set to zero at the end.

In Eq. (38), we have made use of the Wick’s theorem and the following
identities

$\displaystyle\begin{aligned}
&\left\langle\psi^{-,u}_{i_{u}}\bar{\psi}^{-,v}_{j_{v}}\right\rangle_{S_{U}}=\delta_{uv}U_{i_{u}j_{u}},\qquad\left\langle\psi^{+,u}_{j_{u}^{\prime}}\bar{\psi}^{+,v}_{i_{v}^{\prime}}\right\rangle_{S_{U}}=\delta_{uv}U^{\dagger}_{j_{u}^{\prime}i_{u}^{\prime}},\\\
&\left\langle\psi^{-}\bar{\psi}^{+}\right\rangle_{S_{U}}=\left\langle\psi^{+}\bar{\psi}^{-}\right\rangle_{S_{U}}=0.\end{aligned}$
(39)

The angular bracket with the subscript $S_{U}$ represents the functional
averaging over the fermionic field $\psi$ with the weight
$e^{-S_{U}[\bar{\psi},\psi]}$. Note that in Eq. (38) fermions with different
replica indices ($u$) or $U/U^{\dagger}$ indices ($s=\pm$) are uncoupled. In
fact, the replica space is introduced here such that, when computing the
expectation value of $\left\langle F[\bar{\psi},\psi]\right\rangle_{S_{U}}$
using the Wick’s contraction, fermionic field $\psi^{\mp}$ with index $i_{k}$
($j_{k}^{\prime}$) has to pair with $\bar{\psi}^{\mp}$ with index $j_{k}$
($i_{k}^{\prime}$) as they share the same replica index $k$, leading to the
factor of $U_{i_{k}j_{k}}$ ($U_{j_{k}^{\prime}i_{k}^{\prime}}^{\dagger}$) in
the integrand of $I$.

Applying the color-flavor transformation Zirnbauer_1996 ; Zirnbauer_1998 ;
Zirnbauer_1999 ; Zirnbauer_2021 , the Haar integral over the unitary matrix
$U$ in Eq. (38) can be converted into an integral over a complex
$R^{\prime}\times R^{\prime}$ matrix field $Z$:

$\displaystyle\begin{aligned}
I=&z_{1}\int\mathcal{D}(Z,Z^{\dagger})\det(1+Z^{\dagger}Z)^{-(2R^{\prime}+q)}\int\mathcal{D}(\bar{\psi},\psi)e^{-S_{Z}[\bar{\psi},\psi]}F[\bar{\psi},\psi],\\\
=&z_{1}\int\mathcal{D}(Z,Z^{\dagger})\det(1+Z^{\dagger}Z)^{-(2R^{\prime}+q)}\det(-ZZ^{\dagger})^{q}\left\langle
F[\bar{\psi},\psi]\right\rangle_{S_{Z}},\end{aligned}$ (40)

where

$\displaystyle\begin{aligned}
S_{Z}[\bar{\psi},\psi]=&\sum_{u,v=1}^{R^{\prime}}\sum_{l=1}^{q}\left(-\bar{\psi}^{+,u}_{l}Z_{uv}\psi^{-,v}_{l}+\bar{\psi}^{-,u}_{l}Z^{\dagger}_{uv}\psi^{+,v}_{l}\right),\\\
z_{1}^{-1}=&\int\mathcal{D}(Z,Z^{\dagger})\det(1+Z^{\dagger}Z)^{-(2R^{\prime}+q)}\det(-ZZ^{\dagger})^{q}.\end{aligned}$
(41)

Note that the matrix field $Z$ here acts in the replica space only. In the
second equality of Eq. (40), we have integrated out the fermions governed by
the new action $S_{Z}[\bar{\psi},\psi]$, and the angular bracket with the
subscript $S_{Z}$ represents the functional averaging with the weight
$\exp(-S_{Z}[\bar{\psi},\psi])$. The normalization constant $z_{1}$ is
determined from the fact that $I=1$ if we set $F=1$.

After the color-flavor transformation, fermions with different replica indices
interact through the matrix field $Z$, while those with different Hilbert
space indices become uncoupled. The fermionic propagator now acquires the form

$\displaystyle\begin{aligned}
&\left\langle\psi^{-,u}_{i_{u}}\bar{\psi}^{+,v}_{i_{v}^{\prime}}\right\rangle_{S_{Z}}=-\delta_{i_{u}i_{v}^{\prime}}Z^{-1}_{uv},\qquad\left\langle\psi^{+,u}_{j_{u}^{\prime}}\bar{\psi}^{-,v}_{j_{v}}\right\rangle_{S_{Z}}=\delta_{j_{u}^{\prime}j_{v}}(Z^{\dagger})^{-1}_{uv},\\\
&\left\langle\psi^{-}\bar{\psi}^{-}\right\rangle_{S_{Z}}=\left\langle\psi^{+}\bar{\psi}^{+}\right\rangle_{S_{Z}}=0.\end{aligned}$
(42)

Making use of these results, we obtain

$\displaystyle\begin{aligned} \left\langle
F[\bar{\psi},\psi]\right\rangle_{S_{Z}}=&\delta_{p,p^{\prime}}\sum_{\tau,\sigma\in
S_{p}}\operatorname{sgn}(\tau^{-1}\sigma)\prod_{k=1}^{p}\delta_{i_{k},i_{\sigma(k)}^{\prime}}\delta_{j_{k},j_{\tau(k)}^{\prime}}Z^{-1}_{k\sigma(k)}(Z^{\dagger})^{-1}_{\tau(k)k}.\end{aligned}$
(43)

Inserting Eq. (43) into Eq. (40), and applying the transformation

$\displaystyle\begin{aligned} W=Z^{-1},\end{aligned}$ (44)

whose Jacobian leads to a contribution of
$2R^{\prime}\operatorname{tr}\ln(WW^{\dagger})$ to the action, we find

$\displaystyle\begin{aligned} &I=\delta_{p,p^{\prime}}\sum_{\tau,\sigma\in
S_{p}}g(\sigma,\tau)\prod_{k=1}^{p}\delta_{i_{k},i_{\sigma(k)}^{\prime}}\delta_{j_{k},j_{\tau(k)}^{\prime}}.\end{aligned}$
(45)

Here $g(\sigma,\tau)$ is given by

$\displaystyle\begin{aligned} &g(\sigma,\tau)=\left\langle
G[W,W^{\dagger};\sigma,\tau]\right\rangle_{W}\equiv\dfrac{\int\mathcal{D}(W,W^{\dagger})e^{-S_{w}[W,W^{\dagger}]}G[W,W^{\dagger};\sigma,\tau]}{\int\mathcal{D}(W,W^{\dagger})e^{-S_{w}[W,W^{\dagger}]}},\\\
&S_{w}[W,W^{\dagger}]=(2R^{\prime}+q)\operatorname{tr}\ln\left(1+WW^{\dagger}\right),\\\
&G[W,W^{\dagger};\sigma,\tau]=\operatorname{sgn}(\tau^{-1}\sigma)\prod_{k=1}^{p}W_{k\sigma(k)}W^{\dagger}_{\tau(k)k}.\end{aligned}$
(46)

The angular bracket with subscript $W$ represents the averaging over $W$ with
the action $S_{w}[W,W^{\dagger}]$.

Using this field theoretical approach, the Haar integral over the unitary
matrix $U$, whose matrix elements are highly correlated due to the constraint
$UU^{\dagger}=\bf{1}$, has been transformed to an integration over all complex
matrix $W$ governed by the action $S_{w}$. Note that Eq. (45) is equivalent to
the known result Eq. (35) if $g(\sigma,\tau)=\mathrm{Wg}(\tau^{-1}\sigma)$. In
the following, we’re going to show that $g(\sigma,\tau)$ defined by Eq. (46)
depends only on the cycle structure of $\tau^{-1}\sigma$ and is indeed the
Weingarten function $\mathrm{Wg}(\tau^{-1}\sigma)$.

###  5.2 General properties of the Weingarten function

If we apply the following transformation in Eq. (46)

$\displaystyle\begin{aligned} &W_{ij}\rightarrow W_{i\chi(j)},\qquad
W^{\dagger}_{ij}\rightarrow W^{\dagger}_{\chi(i)j},\end{aligned}$ (47)

where $\chi\in S_{p}$ represents an arbitrary permutation of $p$ numbers, the
action $S_{w}[W,W^{\dagger}]$ remains invariant, while
$G[W,W^{\dagger};\sigma,\tau]$ transforms to
$G[W,W^{\dagger};\chi\sigma,\chi\tau]$. This proves that $g(\sigma,\tau)$
obeys the condition

$\displaystyle g(\sigma,\tau)=g(\chi\sigma,\chi\tau),$ (48)

for an arbitrary permutation $\chi\in S_{p}$. Setting $\chi$ to $\tau^{-1}$,
one can immediately see that $g$ is a function of $\tau^{-1}\sigma$ only.

In an analogous way, one can also prove that $g(\sigma,\tau)$ depends only on
the cycle structure of $\tau^{-1}\sigma$:

$\displaystyle g(\tau,\sigma)=V^{(p)}_{c_{1},c_{2},...,c_{m}},$ (49)

with $m$ being the total number of disjoint cycles in $\tau^{-1}\sigma$ and
$c_{i}$ the length of the $i$-th cycle. Rearranging the $W$ and $W^{\dagger}$
terms in $G[W,W^{\dagger};\sigma,\tau]$, we can rewrite $g(\sigma,\tau)$ as

$\displaystyle\begin{aligned}
g(\sigma,\tau)=&\operatorname{sgn}(\tau^{-1}\sigma)\left\langle\prod_{k=1}^{m}\left(\prod_{l=1}^{c_{k}}W_{P_{l}^{(k)}\sigma(P_{l}^{(k)})}W^{\dagger}_{\sigma(P_{l}^{(k)})P_{l+1}^{(k)}}\right)\right\rangle_{W}\\\
=g(\tau^{-1}\sigma,id)=&(-1)^{p-m}\left\langle\prod_{k=1}^{m}\left(\prod_{l=1}^{c_{k}}W_{P_{l}^{(k)}P_{l+1}^{(k)}}W^{\dagger}_{P_{l+1}^{(k)}P_{l+1}^{(k)}}\right)\right\rangle_{W}.\end{aligned}$
(50)

Here the integers $\left\\{P_{l}^{(k)}\right\\}$ denote the cycle structure of
$\tau^{-1}\sigma$ (Eq. (36)) and they satisfy
$P_{l}^{(k)}=(\tau^{-1}\sigma)^{l-1}(P_{1}^{(k)})$. We introduced here the
notation that $P_{l}^{(k)}=P_{(l\,\mathrm{mod}\,c_{k})}^{(k)}$. In the second
equality, we have used $\operatorname{sgn}(\tau^{-1}\sigma)=(-1)^{p-m}$ as
well as Eq. (48).

Note that the action $S_{w}[W,W^{\dagger}]$ is invariant under the
transformation

$\displaystyle\begin{aligned} &W_{ij}\rightarrow W_{\chi(i)\chi(j)},\qquad
W^{\dagger}_{ij}\rightarrow W^{\dagger}_{\chi(i)\chi(j)},\end{aligned}$ (51)

for arbitrary $\chi\in S_{p}$. We now apply this transformation and choose the
permutation $\chi$ defined by

$\displaystyle\chi(P_{l}^{(k)})=b_{k}+l,$ (52)

where $b_{k}=\sum_{h=1}^{k-1}c_{h}$ for $k\geq 2$ and $b_{1}=0$. After this
transformation, $g(\sigma,\tau)$ becomes

$\displaystyle\begin{aligned}
g(\sigma,\tau)=&(-1)^{p-m}\left\langle\prod_{k=1}^{m}\left(\prod_{l=1}^{c_{k}}W_{b_{k}+l,b_{k}+(l+1\,\mathrm{mod}\,c_{k})}W^{\dagger}_{b_{k}+(l+1\,\mathrm{mod}\,c_{k}),b_{k}+(l+1\,\mathrm{mod}\,c_{k})}\right)\right\rangle_{W}.\end{aligned}$
(53)

This equation shows that $g(\sigma,\tau)$ depends only on the lengths
$\left\\{c_{i}|i=1,...,m\right\\}$ of the disjoint cycles of
$\tau^{-1}\sigma$.

Using these results, we prove in appendix E that the function $g(\sigma,\tau)$
defined in Eq. (46) satisfies the recursion relation Eq. (37) and therefore is
given by the Weingarten function $g(\sigma,\tau)=\mathrm{Wg}(\tau^{-1}\sigma)$
for $p\leq q$. See Ref. Samuel  for a discussion about extrapolating the
results for the case of $p\leq q$ to $p>q$. In appendix D, we also provide the
derivation for the asymptotic behavior of the Weingarten function
$g(\sigma,\tau)$ in the large $q\rightarrow\infty$ limit.

##  6 Conclusion

In summary, we derive an effective field theory which can be employed to
investigate the spectral statistics of the Floquet operator for a large class
of Floquet quantum systems. We apply it to a family of Floquet random quantum
circuits whose Floquet operators are composed of two-qudit random unitaries
acting on pairs of neighboring qudits in a $D$-dimensional lattice. Universal
RMT statistics has been found in the limit of infinite local Hilbert space
dimensions, irrespective of the ordering of local gates, dimensionality of the
qudit lattice, and the choice of the boundary condition (open or periodic).
This field theoretical approach has also been used to rederive the known
results for the Weingarten calculus, which is a method to compute Haar
integrals of polynomial functions of the matrix elements.

The wide applicability of the Weingarten calculus in quantum circuit studies
suggests that this field theoretical approach may be generalized to
investigate other fundamental features of quantum many-body systems, which is
a direction for future works. In particular, the field theory may be useful
for the investigation of phase transition in quantum circuits, such as the
measurement-induced phase transition Skinner ; Fisher-1 ; Fisher-2 ; potter
which arises from the competition between the unitary dynamics and the
projective measurements. A renormalization group analysis of the field theory
may be available to carefully examine this phase transition. Moreover, a
mapping has been found between the entanglement growth (or the operator
spreading in the high dimensional case) in random quantum circuits and the
classical statistical mechanics problem of the surface growth in the Kardar-
Parisi-Zhang universality class Nahum-2017 ; Nahum-2018 . The field
theoretical approach may provide some insight into this mapping in a more
generic setting.

In the present paper, we focus on the Weingarten calculus for the unitary
group having in mind random quantum circuits whose local gates are given by
the CUE random matrices. The derivation can be immediately generalized to
other compact groups collins2006 ; collins2009 ; Wg-COE ; Wg-SymSpace , making
use of the color-flavor transformation for the corresponding group
Zirnbauer_2021 . This generalization may have direct applications to random
quantum circuits with different symmetry classes HJ-2018 . Moreover, as a
simple example, the effective field theory given by Eq. (19) is employed above
to study the spectral statistics of Floquet quantum circuits whose local gates
are drawn from the CUE ensemble. In fact, it is applicable to quantum circuits
with random gates drawn from various ensembles of unitary matrices. Therefore
it is a useful analytical tool to study the Thouless energy - the energy
separation below which the RMT statistics appears, and also to look for the
putative transitions between the ergodic phase and many-body localization
phase in various Floquet random quantum circuit models Chalker-2 ; Chalker-5 ;
Chalker-6 ; Huse ; Huse-2 ; Prosen-Local .

Replica trick is used in the present paper to derive the sigma model for
generic Floquet quantum systems, which was originally constructed in the
supersymmetric formalism Zirnbauer_1996 ; Zirnbauer_1998 ; Zirnbauer_1999 ;
QKR . One advantage of using the replica method is that many terms in the
action for higher order fluctuations are irrelevant in the replica limit and
therefore can be omitted. However, unlike the supersymmetric calculation which
allows for a nonperturbative analysis Zirnbauer_1996 ; Efetov , the current
replica calculation is perturbative and can not recover the oscillating term
in the level correlation function. Consideration of nonstandard saddle points
Andreev-Altshuler ; Kamenev-GUE ; Kamenev ; Altland-rev  may recover the
missing oscillating term and is left for a future study.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences under Award No. DE-SC0001911. Y.L. acknowledges a post-
doctoral fellowship from the Simons Foundation “Ultra-Quantum Matter” Research
Collaboration.

##  A Derivation of the moments of Floquet operator

In this appendix, we prove that the Floquet operator $U$ given by Eq. (20)
with arbitrary permutation $\sigma\in S_{L^{\prime}}$ follows the conditions
Eq. (23a) and (23b), in the limit $q\rightarrow\infty$. In other words, the
moments of the Floquet operator $U$ for any Floquet random quantum circuit
that is related to the brickwork circuit studied in Ref. Chalker-1  by
reordering the local quantum gates within one period are given by Eq. (23). We
make use of the fact that the moments of the Floquet operator $U$ are given by
the product of the moments of the independent random CUE matrices
$w^{(n,n+1)}$, which also take the form of Eq. (23) after the replacement
$N\rightarrow q^{2}$ and $\bm{\mathrm{i}}\rightarrow(i^{(n)},i^{(n+1)})$. Here
$q^{2}$ and $(i^{(n)},i^{(n+1)})$ are, respectively, the dimension and label
of the two-qudit Hilbert space in which the CUE matrix $w^{(n,n+1)}$ operates.

For any Floquet operator $U$ which takes the form of Eq. (20), we can express
its second order moment $\left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle$
as

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{n=1}^{L^{\prime}}\left\langle
w^{(n,n+1)}_{a^{(n)},b^{(n+1)};c^{(n)},d^{(n+1)}}(w^{(n,n+1)})^{\dagger}_{c^{\prime(n)},d^{\prime(n+1)};a^{\prime(n)},b^{\prime(n+1)}}\right\rangle\\\
=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{n=1}^{L^{\prime}}\frac{1}{q^{2}}\delta_{a^{(n)},a^{\prime(n)}}\delta_{b^{(n+1)},b^{\prime(n+1)}}\delta_{c^{(n)},c^{\prime(n)}}\delta_{d^{(n+1)},d^{\prime(n+1)}}.\end{aligned}$
(54)

As earlier, vector $\bm{\mathrm{i}}$ labels the many-body state of all the
qudits in the circuit and its $n$-th element $i^{(n)}$ indexes the single-
particle state of the $n$-th qudit. $\bm{\mathrm{k}}=(k^{(1)},...,k^{(L)})$
($\bm{\mathrm{k}}=(k^{(2)},...,k^{(L-1)})$) is a $L$-dimensional
($(L-2)$-dimensional) vector for periodic (open) boundary condition. $k^{(n)}$
labels the single-particle state at site $n$ at an intermediate substep and
one has to sum over all possible states $k^{(n)}=1,...,q$.
$a^{(n)}\,(b^{(n+1)})$ and $c^{(n)}\,(d^{(n+1)})$ are respectively the row and
column indices the label the Hilbert space of the site $n$ ($n+1$) for the CUE
random matrix $w^{(n,n+1)}_{a^{(n)},b^{(n+1)};c^{(n)},d^{(n+1)}}$, and they
are defined as follows. Note first that only the unitaries $w^{(n,n+1)}$ and
$w^{(n-1,n)}$ act on the qudit at site $n$, and
$\left\\{a^{(n)},b^{(n)},c^{(n)},d^{(n)}\right\\}$ depend on the ordering of
these two unitaries. In particular, if $\sigma^{-1}(n)<\sigma^{-1}(n-1)$
(i.e., unitary $w^{(n,n+1)}$ is applied before $w^{(n-1,n)}$), we have

$\displaystyle\begin{aligned} a^{(n)}=i^{(n)},\quad c^{(n)}=k^{(n)},\quad
b^{(n)}=k^{(n)},\quad d^{(n)}=j^{(n)},\end{aligned}$ (55)

and otherwise,

$\displaystyle\begin{aligned} a^{(n)}=k^{(n)},\quad c^{(n)}=j^{(n)},\quad
b^{(n)}=i^{(n)},\quad d^{(n)}=k^{(n)}.\end{aligned}$ (56)

For periodic boundary condition, Eqs. (55) and (56) apply to arbitrary site
$n$ which is defined modulo $L$. By contrast, for open boundary condition,
these two equations hold for $1<n<L$. At the boundaries $n=1,L$, we instead
have

$\displaystyle a^{(1)}=i^{(1)},\qquad c^{(1)}=j^{(1)},\qquad
b^{(L)}=i^{(L)},\qquad d^{(L)}=j^{(L)}.$ (57)

The row and column indices for $(w^{(n,n+1)})^{\dagger}$ are defined in the
same way.
$\left\\{a^{\prime(n)},b^{\prime(n)},c^{\prime(n)},d^{\prime(n)}\right\\}$ are
related to $\left\\{i^{\prime(n)},j^{\prime(n)},k^{\prime(n)}\right\\}$ in
exactly the same way as $\left\\{a^{(n)},b^{(n)},c^{(n)},d^{(n)}\right\\}$ are
related to $\left\\{i^{(n)},j^{(n)},k^{(n)}\right\\}$.

For both the cases of $\sigma^{-1}(n)<\sigma^{-1}(n-1)$ and
$\sigma^{-1}(n)>\sigma^{-1}(n-1)$,
$\left\\{a^{(n)},b^{(n)},c^{(n)},d^{(n)}\right\\}$ are related with
$\left\\{i^{(n)},j^{(n)},k^{(n)},k^{(n)}\right\\}$ by a permutation, and the
following identity

$\displaystyle\delta_{a^{(n)},a^{\prime(n)}}\delta_{b^{(n)},b^{\prime(n)}}\delta_{c^{(n)},c^{\prime(n)}}\delta_{d^{(n)},d^{\prime(n)}}=\delta_{i^{(n)},i^{\prime(n)}}\delta_{j^{(n)},j^{\prime(n)}}\left(\delta_{k^{(n)},k^{\prime(n)}}\right)^{2},$
(58)

is always true for periodic boundary condition (and for open boundary
condition when $1<n<L$). Making use of this property in Eq. (54), we find that
the second order moment of the Floquet operator for the periodic boundary
condition is given by

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\frac{1}{q^{2L}}\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{n=1}^{L}\delta_{i^{(n)},i^{\prime(n)}}\delta_{j^{(n)},j^{\prime(n)}}(\delta_{k^{(n)},k^{\prime(n)}})^{2}=\frac{1}{q^{L}}\prod_{n=1}^{L}\delta_{i^{(n)},i^{\prime(n)}}\delta_{j^{(n)},j^{\prime(n)}}.\end{aligned}$
(59)

The second order moment of $U$ for the open boundary condition can be
evaluated in a similar way:

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\frac{1}{q^{2(L-1)}}\delta_{i^{(1)},i^{\prime(1)}}\delta_{j^{(1)},j^{\prime(1)}}\delta_{i^{(L)},i^{\prime(L)}}\delta_{j^{(L)},j^{\prime(L)}}\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{n=2}^{L-1}\delta_{i^{(n)},i^{\prime(n)}}\delta_{j^{(n)},j^{\prime(n)}}(\delta_{k^{(n)},k^{\prime(n)}})^{2}\\\
=&\frac{1}{q^{L}}\prod_{n=1}^{L}\delta_{i^{(n)},i^{\prime(n)}}\delta_{j^{(n)},j^{\prime(n)}}.\end{aligned}$
(60)

As a concrete example, consider the staircase circuit depicted in Fig. 1(c)
with periodic boundary condition. In this case,
$\left\\{a^{(n)},b^{(n)},c^{(n)},d^{(n)}\right\\}$ are given by Eq. (56), and
the second order moment for the Floquet operator $U$ is given by

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\left\langle
w^{(1,2)}_{i^{(1)},i^{(2)};k^{(1)},k^{(2)}}(w^{(1,2)})^{\dagger}_{k^{\prime(1)},k^{\prime(2)};i^{\prime(1)},i^{\prime(2)}}\right\rangle\\\
&\times\prod_{n=2}^{L-1}\left\langle
w^{(n,n+1)}_{k^{(n)},i^{(n+1)};j^{(n)},k^{(n+1)}}(w^{(n,n+1)})^{\dagger}_{j^{\prime(n)},k^{\prime(n+1)};k^{\prime(n)},i^{\prime(n+1)}}\right\rangle\\\
&\times\left\langle
w^{(L,1)}_{k^{(L)},k^{(1)};j^{(L)},j^{(1)}}(w^{(L,1)})^{\dagger}_{j^{\prime(L)},j^{\prime(1)};k^{\prime(L)},k^{\prime(1)}}\right\rangle\\\
&=\frac{1}{q^{L}}\delta_{\bm{\mathrm{i}},\bm{\mathrm{i^{\prime}}}}\delta_{\bm{\mathrm{j}},\bm{\mathrm{j^{\prime}}}}.\end{aligned}$
(61)

Reordering the two-qudit quantum gates only changes the row and column labels
of $w^{(n,n+1)}$ and $(w^{(n,n+1)})^{\dagger}$ in the equation above but does
not affect the final result. Note that this result for the second order moment
$\left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle$
is valid for arbitrary $q$.

We now evaluate the fourth order moment of $U$ and consider the
$q\rightarrow\infty$ limit for simplicity. It can be expressed as

$\displaystyle\begin{aligned} &\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle=\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\prod_{n=1}^{L^{\prime}}\left\langle
w^{(n,n+1)}_{a^{(n)}_{1},b^{(n+1)}_{1};c^{(n)}_{1},d^{(n+1)}_{1}}w^{(n,n+1)}_{a^{(n)}_{2},b^{(n+1)}_{2};c^{(n)}_{2},d^{(n+1)}_{2}}\right.\\\
&\times\left.(w^{(n,n+1)})^{\dagger}_{c^{\prime(n)}_{1},d^{\prime(n+1)}_{1};a^{\prime(n)}_{1},b^{\prime(n+1)}_{1}}(w^{(n,n+1)})^{\dagger}_{c^{\prime(n)}_{2},d^{\prime(n+1)}_{2};a^{\prime(n)}_{2},b^{\prime(n+1)}_{2}}\right\rangle\\\
=&\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\prod_{n=1}^{L^{\prime}}\frac{1}{q^{4}}\sum_{s^{(n)}=1}^{2}\left(\begin{aligned}
&\delta_{s^{(n)},1}\delta_{a^{(n)}_{1},a^{\prime(n)}_{1}}\delta_{b^{(n+1)}_{1},b^{\prime(n+1)}_{1}}\delta_{c^{(n)}_{1},c^{\prime(n)}_{1}}\delta_{d^{(n+1)}_{1},d^{\prime(n+1)}_{1}}\\\
&\times\delta_{a^{(n)}_{2},a^{\prime(n)}_{2}}\delta_{b^{(n+1)}_{2},b^{\prime(n+1)}_{2}}\delta_{c^{(n)}_{2},c^{\prime(n)}_{2}}\delta_{d^{(n+1)}_{2},d^{\prime(n+1)}_{2}}\\\
+&\delta_{s^{(n)},2}\delta_{a^{(n)}_{1},a^{\prime(n)}_{2}}\delta_{b^{(n+1)}_{1},b^{\prime(n+1)}_{2}}\delta_{c^{(n)}_{1},c^{\prime(n)}_{2}}\delta_{d^{(n+1)}_{1},d^{\prime(n+1)}_{2}}\\\
&\times\delta_{a^{(n)}_{2},a^{\prime(n)}_{1}}\delta_{b^{(n+1)}_{2},b^{\prime(n+1)}_{1}}\delta_{c^{(n)}_{2},c^{\prime(n)}_{1}}\delta_{d^{(n+1)}_{2},d^{\prime(n+1)}_{1}}\end{aligned}\right).\end{aligned}$
(62)

Here $\bm{\mathrm{k}}_{1,2}$ and $\bm{\mathrm{k}}^{\prime}_{1,2}$, similar to
$\bm{\mathrm{k}}$ and $\bm{\mathrm{k}}^{\prime}$, are $L$-dimensional
($(L-2)$-dimensional) vectors whose components label the intermediate single-
particle states for all sites (sites in the bulk of the lattice $1<n<L$) for
periodic (open) boundary condition.
$\left\\{a^{(n)}_{1/2},b^{(n)}_{1/2},c^{(n)}_{1/2},d^{(n)}_{1/2}\right\\}$ are
related to $\left\\{i^{(n)}_{1/2},j^{(n)}_{1/2},k^{(n)}_{1/2}\right\\}$ in the
same way as $\left\\{a^{(n)},b^{(n)},c^{(n)},d^{(n)}\right\\}$ are related to
$\left\\{i^{(n)},j^{(n)},k^{(n)}\right\\}$ (see Eqs. (55), (56) and (57)), and
similarly for
$\left\\{a^{\prime(n)}_{1/2},b^{\prime(n)}_{1/2},c^{\prime(n)}_{1/2},d^{\prime(n)}_{1/2}\right\\}$.
The summation over $s^{(n)}=1,2$ is simply introduced here to distinguish the
two terms within the bracket.

In the second equality in Eq. (62), terms involving
$\delta_{k^{(n)}_{1},k^{\prime(n)}_{1}}\delta_{k^{(n)}_{2},k^{\prime(n)}_{2}}\delta_{k^{(n)}_{1},k^{\prime(n)}_{2}}\delta_{k^{(n)}_{2},k^{\prime(n)}_{1}}$,
compared with those involving
$\left(\delta_{k^{(n)}_{1},k^{\prime(n)}_{1}}\delta_{k^{(n)}_{2},k^{\prime(n)}_{2}}\right)^{2}$
or
$\left(\delta_{k^{(n)}_{1},k^{\prime(n)}_{2}}\delta_{k^{(n)}_{2},k^{\prime(n)}_{1}}\right)^{2}$,
are of higher order in $1/q$, since the number of free summations is reduced.
For this reason, terms with $s^{(n)}\neq s^{(n+1)}$ for at least one $n$,
compared with those with all $s^{(n)}$ identical, are of higher order in the
large $q$ expansion and can be ignored. As a result, we obtain

$\displaystyle\begin{aligned} &\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle=\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\frac{1}{q^{4L^{\prime}}}\left(\begin{aligned}
\prod_{n=1}^{L^{\prime}}&\delta_{a^{(n)}_{1},a^{\prime(n)}_{1}}\delta_{b^{(n+1)}_{1},b^{\prime(n+1)}_{1}}\delta_{c^{(n)}_{1},c^{\prime(n)}_{1}}\delta_{d^{(n+1)}_{1},d^{\prime(n+1)}_{1}}\\\
&\times\delta_{a^{(n)}_{2},a^{\prime(n)}_{2}}\delta_{b^{(n+1)}_{2},b^{\prime(n+1)}_{2}}\delta_{c^{(n)}_{2},c^{\prime(n)}_{2}}\delta_{d^{(n+1)}_{2},d^{\prime(n+1)}_{2}}\\\
+\prod_{n=1}^{L^{\prime}}&\delta_{a^{(n)}_{1},a^{\prime(n)}_{2}}\delta_{b^{(n+1)}_{1},b^{\prime(n+1)}_{2}}\delta_{c^{(n)}_{1},c^{\prime(n)}_{2}}\delta_{d^{(n+1)}_{1},d^{\prime(n+1)}_{2}}\\\
&\times\delta_{a^{(n)}_{2},a^{\prime(n)}_{1}}\delta_{b^{(n+1)}_{2},b^{\prime(n+1)}_{1}}\delta_{c^{(n)}_{2},c^{\prime(n)}_{1}}\delta_{d^{(n+1)}_{2},d^{\prime(n+1)}_{1}}\end{aligned}\right)\\\
&=\frac{1}{q^{2L}}\left(\prod_{n=1}^{L}\delta_{i_{1}^{(n)},i_{1}^{\prime(n)}}\delta_{i_{2}^{(n)},i_{2}^{\prime(n)}}\delta_{j_{1}^{(n)},j_{1}^{\prime(n)}}\delta_{j_{2}^{(n)},j_{2}^{\prime(n)}}+\prod_{n=1}^{L}\delta_{i_{1}^{(n)},i_{2}^{\prime(n)}}\delta_{i_{2}^{(n)},i_{1}^{\prime(n)}}\delta_{j_{1}^{(n)},j_{2}^{\prime(n)}}\delta_{j_{2}^{(n)},j_{1}^{\prime(n)}}\right).\end{aligned}$
(63)

In the second equality, we used the fact that Eq. (54) reduces to Eq. (59)
(Eq. (60)) upon substituting Eqs. (55) and (56) (Eqs. (55), (56) and (57)).
This proves Eq. (23) for any Floquet operator described by Eq. (20) with
arbitrary permutation $\sigma\in S_{L^{\prime}}$.

##  B Generalization to higher dimensional Floquet quantum circuits

In this appendix, we generalize the calculation in appendix A to an arbitrary
high dimension $D$. More specifically, we consider now a $D$-dimensional
Floquet random quantum circuit which consists of a $D$-dimensional cubic
lattice of qudits. During one period, each qudit is coupled to all its
neighboring qudits at different time substeps by different two-qudit unitary
gates. All these unitary gates are drawn randomly and independently from the
CUE ensemble. We will prove that, for any Floquet quantum circuit of this type
with an arbitrary ordering of the two-qudit unitary gates, and for both
periodic and open boundary conditions, the second and fourth order moments of
the Floquet operator $U$ take the form Eq. (23), same as the CUE ensemble of
dimension $N=q^{L^{D}}$, with $q$ being the Hilbert space dimension of each
qudit and $L$ the number of lattice sites in each direction. For simplicity,
we consider the large $q\rightarrow\infty$ limit for the calculation of the
fourth order moment.

In the following, we employ the notation that each qudit is labeled by a
$D$-dimensional vector $\bm{\mathrm{n}}$ whose $\alpha$-th component
$n_{\alpha}$ represents the coordinate in the lattice in dimension
$\alpha=1,...,D$. The neighbors of qudit $\bm{\mathrm{n}}$ are denoted by
$\bm{\mathrm{n}}\pm\bm{\mathrm{e}}_{\alpha}$, with $\bm{\mathrm{e}}_{\alpha}$
being the unit vector in dimension $\alpha$. $n_{\alpha}$ is defined modulo
$L$ for periodic boundary condition. The many-body state is indexed by
$\bm{\mathrm{i}}$ and the single-particle state at site $\bm{\mathrm{n}}$ is
labeled by $i^{(\bm{\mathrm{n}})}$.

As in the 1D case, the second order moment $\left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle$
can be expressed as

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\left\langle
w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{a^{(\bm{\mathrm{n}},\alpha)},b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)};c^{(\bm{\mathrm{n}},\alpha)},d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}}(w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})})^{\dagger}_{c^{\prime(\bm{\mathrm{n}},\alpha)},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)};a^{\prime(\bm{\mathrm{n}},\alpha)},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}}\right\rangle\\\
=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\frac{1}{q^{2}}\delta_{a^{(\bm{\mathrm{n}},\alpha)},a^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}}\delta_{c^{(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}}.\end{aligned}$
(64)

Here $w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}$
represents the random CUE matrix that couples the qudits at sites
$\bm{\mathrm{n}}$ and $\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha}$.
$\prod_{\bm{\mathrm{n}}}^{\prime}$ stands for
$\prod_{n_{1}=1}^{L^{\prime}}...\prod_{n_{D}=1}^{L^{\prime}}$, where
$L^{\prime}$ is defined as $L$ ($L-1$) for periodic (open) boundary condition.
We use $k^{(\bm{\mathrm{n}},m)}$ to denote the single-particle state for qudit
$\bm{\mathrm{n}}$ at a intermediate substep, and it needs to be summed over
($\sum_{k^{(\bm{\mathrm{n}},m)}=1}^{q}$). $\sum_{\bm{\mathrm{k}}}$ is short
for the product
$\prod_{\bm{\mathrm{n}}}\prod_{m=1}^{M_{\bm{\mathrm{n}}}}\sum_{k^{(\bm{\mathrm{n}},m)}=1}^{q}$
where $m$ runs over $m=1,2,...,M_{\bm{\mathrm{n}}}$ and $\bm{\mathrm{n}}$ runs
over all sites in the lattice. $M_{\bm{\mathrm{n}}}+1$ is the number of
neighbors for site $\bm{\mathrm{n}}$. For periodic boundary condition
$M_{\bm{\mathrm{n}}}=2D-1$, while for open boundary condition

$\displaystyle M_{\bm{\mathrm{n}}}=2D-D_{\bm{\mathrm{n}}}-1,$ (65)

where $D_{\bm{\mathrm{n}}}$ denotes the number of components of vector
$\bm{\mathrm{n}}$ satisfying $n_{\alpha}=1$ or $L$.
$\sum_{\bm{\mathrm{k}}^{\prime}}$ is defined in the same manner. The row and
column indices
$\left\\{a^{(\bm{\mathrm{n}},\alpha)},b^{(\bm{\mathrm{n}},\alpha)},c^{(\bm{\mathrm{n}},\alpha)},d^{(\bm{\mathrm{n}},\alpha)}\right\\}$
and
$\left\\{a^{\prime(\bm{\mathrm{n}},\alpha)},b^{\prime(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)},d^{\prime(\bm{\mathrm{n}},\alpha)}\right\\}$
label the single-particle state at site $\bm{\mathrm{n}}$ at a certain time
substep and depend on the ordering of the quantum gates that couples sites
$\bm{\mathrm{n}}$ with all its neighbors.

Consider first the periodic boundary condition. For arbitrary ordering of the
two-qudit unitary gates, the $4D$ numbers in the set
$S_{\bm{\mathrm{n}}}^{(1)}=\left\\{a^{(\bm{\mathrm{n}},\alpha)},b^{(\bm{\mathrm{n}},\alpha)},c^{(\bm{\mathrm{n}},\alpha)},d^{(\bm{\mathrm{n}},\alpha)}|\alpha=1,2,...,D\right\\}$
are always related to those in the set
$S_{\bm{\mathrm{n}}}^{(2)}=\left\\{i^{(\bm{\mathrm{n}})},j^{(\bm{\mathrm{n}})},k^{(\bm{\mathrm{n}},m)},k^{(\bm{\mathrm{n}},m)}|m=1,2,...,2D-1\right\\}$
by a permutation. Row and column indices for $w^{\dagger}$ are defined in an
analogous way. In particular, the numbers in the set
$S_{\bm{\mathrm{n}}}^{\prime(1)}=\left\\{a^{\prime(\bm{\mathrm{n}},\alpha)},b^{\prime(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)},d^{\prime(\bm{\mathrm{n}},\alpha)}|\alpha=1,2,...,D\right\\}$
are related to those in the set
$S_{\bm{\mathrm{n}}}^{\prime(2)}=\left\\{i^{\prime(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})},k^{\prime(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}|m=1,2,...,2D-1\right\\}$
by the same permutation that connects $S_{\bm{\mathrm{n}}}^{(1)}$ to
$S_{\bm{\mathrm{n}}}^{(2)}$. That means, for arbitrary ordering of the two-
qudit gates, the following identity is always satisfied

$\displaystyle\prod_{\alpha=1}^{D}\delta_{a^{(\bm{\mathrm{n}},\alpha)},a^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{b^{(\bm{\mathrm{n}},\alpha)},b^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{c^{(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{d^{(\bm{\mathrm{n}},\alpha)},d^{\prime(\bm{\mathrm{n}},\alpha)}}=\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}\left(\prod_{m=1}^{2D-1}\delta_{k^{(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}}\right)^{2}.$
(66)

Inserting this equation into the second equality of Eq. (64), we obtain the
second order moment for the Floquet operator

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{\bm{\mathrm{n}}}\frac{1}{q^{2D}}\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}\left(\prod_{m=1}^{2D-1}\delta_{k^{(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}}\right)^{2}\\\
=&\frac{1}{q^{L^{D}}}\prod_{\bm{\mathrm{n}}}\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}.\end{aligned}$
(67)

For the open boundary condition, extra care must be taken for the sites
$\bm{\mathrm{n}}$ at the boundaries. For any site $\bm{\mathrm{n}}$ sits in
the bulk of the lattice, as in the case of periodic boundary condition, the
sets $S_{\bm{\mathrm{n}}}^{(1)}$ and $S_{\bm{\mathrm{n}}}^{(2)}$
($S_{\bm{\mathrm{n}}}^{\prime(1)}$ and $S_{\bm{\mathrm{n}}}^{\prime(2)}$)
defined above are always related by a $4D$-number permutation for any ordering
of the quantum gates. Moreover, the permutation that connects
$S_{\bm{\mathrm{n}}}^{(1)}$ and $S_{\bm{\mathrm{n}}}^{(2)}$ is identical to
the one that connects $S_{\bm{\mathrm{n}}}^{\prime(1)}$ and
$S_{\bm{\mathrm{n}}}^{\prime(2)}$. By contrast, at the boundary, for any
component of the vector $\bm{\mathrm{n}}$ that satisfies $n_{\alpha}=1$
($n_{\alpha}=L$), $b^{(\bm{\mathrm{n}},\alpha)}$ and
$d^{(\bm{\mathrm{n}},\alpha)}$ ($a^{(\bm{\mathrm{n}},\alpha)}$ and
$c^{(\bm{\mathrm{n}},\alpha)}$) have to be removed from the set
$S_{\bm{\mathrm{n}}}^{(1)}=\left\\{a^{(\bm{\mathrm{n}},\beta)},b^{(\bm{\mathrm{n}},\beta)},c^{(\bm{\mathrm{n}},\beta)},d^{(\bm{\mathrm{n}},\beta)}|\beta=1,2,...,D\right\\}$.
After applying this procedure, all the $2M_{\bm{\mathrm{n}}}+2$ numbers left
in the new set $S_{\bm{\mathrm{n}}}^{(1)}$ are related to those in the set
$S_{\bm{\mathrm{n}}}^{(2)}=\left\\{i^{(\bm{\mathrm{n}})},j^{(\bm{\mathrm{n}})},k^{(\bm{\mathrm{n}},m)},k^{(\bm{\mathrm{n}},m)}|m=1,2,...,M_{\bm{\mathrm{n}}}\right\\}$
by a permutation.
$\left\\{a^{\prime(\bm{\mathrm{n}},\alpha)},b^{\prime(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)},d^{\prime(\bm{\mathrm{n}},\alpha)}\right\\}$
are defined in the same way. They are related to
$\left\\{i^{\prime(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})},k^{\prime(\bm{\mathrm{n}},m)}\right\\}$
in exactly the same way as
$\left\\{a^{(\bm{\mathrm{n}},\alpha)},b^{(\bm{\mathrm{n}},\alpha)},c^{(\bm{\mathrm{n}},\alpha)},d^{(\bm{\mathrm{n}},\alpha)}\right\\}$
are related to
$\left\\{i^{(\bm{\mathrm{n}})},j^{(\bm{\mathrm{n}})},k^{(\bm{\mathrm{n}},m)}\right\\}$.
This leads to

$\displaystyle\begin{aligned} &\prod_{\alpha|n_{\alpha}\neq
L}\delta_{a^{(\bm{\mathrm{n}},\alpha)},a^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{c^{(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)}}\prod_{\beta|n_{\beta}\neq
1}\delta_{b^{(\bm{\mathrm{n}},\beta)},b^{\prime(\bm{\mathrm{n}},\beta)}}\delta_{d^{(\bm{\mathrm{n}},\beta)},d^{\prime(\bm{\mathrm{n}},\beta)}}\\\
&=\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}\left(\prod_{m=1}^{M_{\bm{\mathrm{n}}}}\delta_{k^{(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}}\right)^{2}.\end{aligned}$
(68)

Here the first (second) product $\prod_{\alpha|n_{\alpha}\neq L}$
($\prod_{\beta|n_{\beta}\neq 1}$) runs over all directions $\alpha$ ($\beta$)
that obey $n_{\alpha}\neq L$ ($n_{\beta}\neq 1$).

Using this result, we find that Eq. (64) reduces to

$\displaystyle\begin{aligned} \left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle=&\sum_{\bm{\mathrm{k}},\bm{\mathrm{k}}^{\prime}}\prod_{\bm{\mathrm{n}}}\left(\frac{1}{q^{M_{\bm{\mathrm{n}}}+1}}\\!\\!\\!\\!\prod_{\alpha|n_{\alpha}\neq
L}\\!\\!\\!\\!\delta_{a^{(\bm{\mathrm{n}},\alpha)},a^{\prime(\bm{\mathrm{n}},\alpha)}}\delta_{c^{(\bm{\mathrm{n}},\alpha)},c^{\prime(\bm{\mathrm{n}},\alpha)}}\\!\\!\\!\\!\prod_{\beta|n_{\beta}\neq
1}\\!\\!\\!\\!\delta_{b^{(\bm{\mathrm{n}},\beta)},b^{\prime(\bm{\mathrm{n}},\beta)}}\delta_{d^{(\bm{\mathrm{n}},\beta)},d^{\prime(\bm{\mathrm{n}},\beta)}}\right)\\\
=&\prod_{\bm{\mathrm{n}}}\left[\frac{1}{q^{M_{\bm{\mathrm{n}}}+1}}\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}\prod_{m=1}^{M_{\bm{\mathrm{n}}}}\left(\sum_{k^{(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}}\delta_{k^{(\bm{\mathrm{n}},m)},k^{\prime(\bm{\mathrm{n}},m)}}^{2}\right)\right]\\\
=&\frac{1}{q^{L^{D}}}\prod_{\bm{\mathrm{n}}}\delta_{i^{(\bm{\mathrm{n}})},i^{\prime(\bm{\mathrm{n}})}}\delta_{j^{(\bm{\mathrm{n}})},j^{\prime(\bm{\mathrm{n}})}}.\end{aligned}$
(69)

We have therefore proved Eq. (23) for the $D$-dimensional Floquet random
quantum circuit which has an arbitrary ordering of the two-qudit random
unitary gates and is subject to either periodic or open boundary condition.

The fourth order moment $\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle$,
in the limit of $q\rightarrow\infty$, can be expressed in a form similar to
the 1D expression Eq. (62):

$\displaystyle\begin{aligned} &\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle=\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\left\langle
w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{a^{(\bm{\mathrm{n}},\alpha)}_{1},b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1};c^{(\bm{\mathrm{n}},\alpha)}_{1},d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{a^{(\bm{\mathrm{n}},\alpha)}_{2},b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2};c^{(\bm{\mathrm{n}},\alpha)}_{2},d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\right.\\\
&\times\left.(w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})})^{\dagger}_{c^{\prime(\bm{\mathrm{n}})}_{1},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{1};a^{\prime(\bm{\mathrm{n}})}_{1},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{1}}(w^{(\bm{\mathrm{n}},\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})})^{\dagger}_{c^{\prime(\bm{\mathrm{n}})}_{2},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{2};a^{\prime(\bm{\mathrm{n}})}_{2},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha})}_{2}}\right\rangle\\\
=&\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\frac{1}{q^{4}}\sum_{s^{(\bm{\mathrm{n}},\alpha)}=1}^{2}\left(\begin{aligned}
&\delta_{s^{(\bm{\mathrm{n}},\alpha)},1}\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{1},a^{\prime(\bm{\mathrm{n}},,\alpha)}_{1}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{1},c^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\\\
&\times\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{2},a^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{2},c^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\\\
+&\delta_{s^{(\bm{\mathrm{n}},\alpha)},2}\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{1},a^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{1},c^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\\\
&\times\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{2},a^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{2},c^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\end{aligned}\right).\end{aligned}$
(70)

Here indices
$\left\\{a^{(\bm{\mathrm{n}},\alpha)}_{1/2},b^{(\bm{\mathrm{n}},\alpha)}_{1/2},c^{(\bm{\mathrm{n}},\alpha)}_{1/2},d^{(\bm{\mathrm{n}})}_{1/2}\right\\}$
are related to
$\left\\{i^{(\bm{\mathrm{n}})}_{1/2},j^{(\bm{\mathrm{n}})}_{1/2},k^{(\bm{\mathrm{n}},m)}_{1/2}\right\\}$
in the same way as
$\left\\{a^{(\bm{\mathrm{n}},\alpha)},b^{(\bm{\mathrm{n}},\alpha)},c^{(\bm{\mathrm{n}},\alpha)},d^{(\bm{\mathrm{n}})}\right\\}$
are related to
$\left\\{i^{(\bm{\mathrm{n}})},j^{(\bm{\mathrm{n}})},k^{(\bm{\mathrm{n}},m)}\right\\}$
(similarly for the indices of $w^{\dagger}$). As before,
$s^{(\bm{\mathrm{n}},\alpha)}$ is introduced to distinguish the two terms in
the bracket.

Similar to the 1D case, in the second equality in the equation above, terms
that contain the factor
$\delta_{k^{(\bm{\mathrm{n}},m)}_{1},k^{\prime(\bm{\mathrm{n}},m)}_{1}}\delta_{k^{(\bm{\mathrm{n}},m)}_{2},k^{\prime(\bm{\mathrm{n}},m)}_{2}}\delta_{k^{(\bm{\mathrm{n}},m)}_{1},k^{\prime(\bm{\mathrm{n}},m)}_{2}}\delta_{k^{(\bm{\mathrm{n}},m)}_{2},k^{\prime(\bm{\mathrm{n}},m)}_{1}}$
are of higher order in the large $q$ expansion, compared with those with the
factor
$\left(\delta_{k^{(\bm{\mathrm{n}},m)}_{1},k^{\prime(\bm{\mathrm{n}},m)}_{1}}\delta_{k^{(\bm{\mathrm{n}},m)}_{2},k^{\prime(\bm{\mathrm{n}},m)}_{2}}\right)^{2}$
or
$\left(\delta_{k^{(\bm{\mathrm{n}},m)}_{1},k^{\prime(\bm{\mathrm{n}},m)}_{2}}\delta_{k^{(\bm{\mathrm{n}},m)}_{2},k^{\prime(\bm{\mathrm{n}},m)}_{1}}\right)^{2}$
due to the reduced number of free summations. For this reason, we only need to
keep terms for which all $s^{(\bm{\mathrm{n}},\alpha)}$ are identical in the
leading order of the large $q$ expansion. This leads to

$\displaystyle\begin{aligned} &\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{1}\bm{\mathrm{i}}^{\prime}_{1}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{2}}\right\rangle\\\
=&\sum_{\bm{\mathrm{k}}_{1},\bm{\mathrm{k}}_{1}^{\prime},\bm{\mathrm{k}}_{2},\bm{\mathrm{k}}_{2}^{\prime}}\left(\begin{aligned}
&\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\frac{1}{q^{4}}\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{1},a^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{1},c^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\\\
&\times\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{2},a^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{2},c^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\\\
&+\prod_{\bm{\mathrm{n}}}^{\prime}\prod_{\alpha=1}^{D}\frac{1}{q^{4}}\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{1},a^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{1},c^{\prime(\bm{\mathrm{n}},\alpha)}_{2}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2}}\\\
&\times\delta_{a^{(\bm{\mathrm{n}},\alpha)}_{2},a^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{b^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},b^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\delta_{c^{(\bm{\mathrm{n}},\alpha)}_{2},c^{\prime(\bm{\mathrm{n}},\alpha)}_{1}}\delta_{d^{(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{2},d^{\prime(\bm{\mathrm{n}}+\bm{\mathrm{e}}_{\alpha},\alpha)}_{1}}\end{aligned}\right)\\\
&=\frac{1}{q^{2L^{D}}}\left(\prod_{\bm{\mathrm{n}}}\delta_{i_{1}^{(\bm{\mathrm{n}})},i_{1}^{\prime(\bm{\mathrm{n}})}}\delta_{i_{2}^{(\bm{\mathrm{n}})},i_{2}^{\prime(\bm{\mathrm{n}})}}\delta_{j_{1}^{(\bm{\mathrm{n}})},j_{1}^{\prime(\bm{\mathrm{n}})}}\delta_{j_{2}^{(\bm{\mathrm{n}})},j_{2}^{\prime(\bm{\mathrm{n}})}}+\prod_{\bm{\mathrm{n}}}\delta_{i_{1}^{(\bm{\mathrm{n}})},i_{2}^{\prime(\bm{\mathrm{n}})}}\delta_{i_{2}^{(\bm{\mathrm{n}})},i_{1}^{\prime(\bm{\mathrm{n}})}}\delta_{j_{1}^{(\bm{\mathrm{n}})},j_{2}^{\prime(\bm{\mathrm{n}})}}\delta_{j_{2}^{(\bm{\mathrm{n}})},j_{1}^{\prime(\bm{\mathrm{n}})}}\right).\end{aligned}$
(71)

In the second equality here, we made use of the intermediate results in the
derivation of the second order moment (i.e., the fact that Eq. (64) is
identical to Eq. (67) and Eq. (69) for periodic and open boundary condition,
respectively).

##  C Effective field theory of a noninteracting Floquet model

In this appendix, we compare the effective field theory of the Floquet random
quantum circuits studied earlier (or equivalently the CUE ensemble) with that
of a noninteracting Floquet model, and show that the quartic fluctuations
become important in the noninteracting case. We consider now a time-periodic
noninteracting system whose many-body Floquet operator is simply given by the
tensor product of all single-particle Floquet operator $w^{(n)}$:

$\displaystyle\begin{aligned} U=w^{(1)}\otimes w^{(2)}\otimes...\otimes
w^{(L)}.\end{aligned}$ (72)

All $\left\\{w^{(n)}|n=1,...,L\right\\}$ are $q\times q$ unitary matrices
drawn randomly and independent from the CUE ensemble. As before, we consider
the limit of large single-particle Hilbert space dimension
$q\rightarrow\infty$ for simplicity.

Since all single-particle Floquet operators $w^{(n)}$ are independent, the
many-body SFF $K(t)$ is given by the product of all single-particle SFFs
$k^{(n)}(t)$:

$\displaystyle\begin{aligned} K(t)=\prod_{n=1}^{L}k^{(n)}(t),\qquad
k^{(n)}(t)=\left\langle|\operatorname{Tr}(w^{(n)})^{t}|^{2}\right\rangle.\end{aligned}$
(73)

Here $k^{(n)}(t)$ is simply the SFF for the CUE ensemble of dimension $q$ (Eq.
(9) after the replacement $N\rightarrow q$). Therefore, one can easily see
that the many-body SFF now exhibits a fast growing ramp $K(t)=t^{L}$ which
plateaus at time $t=q$. The many-body quasi-energies are uncorrelated at the
energy separation of the order of many-body level spacing $\Delta\phi\sim
1/N$, with $N=q^{L}$ being the dimension of many-body Hilbert space PRL .

The many-body Floquet operator $U$ given by Eq. (72) satisfies Eq. (23a) but
not Eq. (23b). Its fourth order moment in the large $q$ limit takes the form,

$\displaystyle\begin{aligned} &\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}_{1}^{\prime}\bm{\mathrm{i}}_{1}^{\prime}}U^{\dagger}_{\bm{\mathrm{j}}_{2}^{\prime}\bm{\mathrm{i}}_{2}^{\prime}}\right\rangle=\prod_{n=1}^{L}\left\langle
w^{(n)}_{i^{(n)}_{1},j^{(n)}_{1}}w^{(n)}_{i^{(n)}_{2},j^{(n)}_{2}}(w^{(n)})^{\dagger}_{j^{\prime(n)}_{1},i^{\prime(n)}_{1}}(w^{(n)})^{\dagger}_{j^{\prime(n)}_{2},i^{\prime(n)}_{2}}\right\rangle\\\
=&\frac{1}{q^{2L}}\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{1},i^{\prime(n)}_{1}}\delta_{j^{(n)}_{1},j^{\prime(n)}_{1}}\delta_{i^{(n)}_{2},i^{\prime(n)}_{2}}\delta_{j^{(n)}_{2},j^{\prime(n)}_{2}}+\delta_{i^{(n)}_{1},i^{\prime(n)}_{2}}\delta_{j^{(n)}_{1},j^{\prime(n)}_{2}}\delta_{i^{(n)}_{2},i^{\prime(n)}_{1}}\delta_{j^{(n)}_{2},j^{\prime(n)}_{1}}\right).\end{aligned}$
(74)

Note that for the chaotic Floquet circuits described by Eq. (20), to have
nonzero moment $\left\langle
U_{\bm{\mathrm{i}}_{1}\bm{\mathrm{j}}_{1}}U_{\bm{\mathrm{i}}_{2}\bm{\mathrm{j}}_{2}}U^{\dagger}_{\bm{\mathrm{j}}_{1}^{\prime}\bm{\mathrm{i}}_{1}^{\prime}}U^{\dagger}_{\bm{\mathrm{j}}_{2}^{\prime}\bm{\mathrm{i}}_{2}^{\prime}}\right\rangle$,
if $i^{(n)}_{1}=i^{\prime(n)}_{1/2}$ is obeyed by one of the site, it must
holds for all other sites as well. Otherwise, the fourth order moment can be
ignored in the leading order of the large $q$ expansion. However, this is no
longer the case for the noninteracting model. It is this difference that leads
to the contrasting behaviors of the many-body level statistics for the two
models.

Since the second order moment of the Floquet operator $\left\langle
U_{\bm{\mathrm{i}}\bm{\mathrm{j}}}U^{\dagger}_{\bm{\mathrm{j}}^{\prime}\bm{\mathrm{i}}^{\prime}}\right\rangle$
for the current noninteracting model is identical to that of the chaotic
Floquet circuits discussed in section 4, the quadratic fluctuations in the
effective theory Eq. (19) are therefore also governed by the action in Eq.
(24). On the other hand, the quartic order effective action of the
noninteracting model is different from that of the chaotic model (Eq. (30))
and assumes the form

$\displaystyle\begin{aligned} S_{\rm
eff}^{(4-1)}=&-\frac{(\alpha\beta)^{2}}{2}\sum_{\bm{\mathrm{i}}_{1},...,\bm{\mathrm{i}}_{8}}\operatorname{tr}\left(Z^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}Z_{\bm{\mathrm{i}}_{3}\bm{\mathrm{i}}_{4}}\right)\operatorname{tr}\left(Z^{\dagger}_{\bm{\mathrm{i}}_{5}\bm{\mathrm{i}}_{6}}Z_{\bm{\mathrm{i}}_{7}\bm{\mathrm{i}}_{8}}\right)\\\
&\times\frac{1}{q^{2L}}\left[\begin{aligned}
&\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i^{(n)}_{6},i^{(n)}_{5}}\delta_{i^{(n)}_{3},i^{(n)}_{4}}\delta_{i^{(n)}_{7},i^{(n)}_{8}}+\delta_{i^{(n)}_{2},i^{(n)}_{5}}\delta_{i^{(n)}_{6},i^{(n)}_{1}}\delta_{i^{(n)}_{3},i^{(n)}_{8}}\delta_{i^{(n)}_{7},i^{(n)}_{4}}\right)\\\
&-\prod_{n=1}^{L}\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i^{(n)}_{6},i^{(n)}_{5}}\delta_{i^{(n)}_{3},i^{(n)}_{4}}\delta_{i^{(n)}_{7},i^{(n)}_{8}}\end{aligned}\right].\end{aligned}$
(75)

Here we have inserted Eq. (74) into $S_{\rm eff}^{(4)}$ in Eq. (19c), and, as
before, ignored the last term in $S_{\rm eff}^{(4)}$ which is not important in
the replica limit $R\rightarrow 0$.

Rewriting Eq. (75) in terms of $X$ and $Y$, which denote, respectively, the
diagonal and off-diagonal components of the matrix $Z$ in the Hilbert space,
we obtain

$\displaystyle\begin{aligned} &S_{\rm
eff}^{(4-1)}=-\frac{(\alpha\beta)^{2}}{2N^{2}}\sum_{\bm{\mathrm{k}}_{1},...,\bm{\mathrm{k}}_{4}}\operatorname{tr}\left(X^{\dagger}(\bm{\mathrm{k}}_{1})X(\bm{\mathrm{k}}_{2})\right)\operatorname{tr}\left(X^{\dagger}(\bm{\mathrm{k}}_{3})X(\bm{\mathrm{k}}_{4})\right)\\\
&\times\left[\prod_{n=1}^{L}\left(\delta_{k^{(n)}_{1},0}\delta_{k^{(n)}_{2},0}\delta_{k^{(n)}_{3},0}\delta_{k^{(n)}_{4},0}+\frac{1}{q^{2}}\delta_{k^{(n)}_{1},-k^{(n)}_{3}}\delta_{k^{(n)}_{2},-k^{(n)}_{4}}\right)-\prod_{n=1}^{L}\left(\delta_{k^{(n)}_{1},0}\delta_{k^{(n)}_{2},0}\delta_{k^{(n)}_{3},0}\delta_{k^{(n)}_{4},0}\right)\right]\\\
&-\frac{(\alpha\beta)^{2}}{2N^{2}}\sum_{\bm{\mathrm{i}}_{1}\neq\bm{\mathrm{i}}_{2},\bm{\mathrm{i}}_{3}\neq\bm{\mathrm{i}}_{4},\bm{\mathrm{i}}_{5}\neq\bm{\mathrm{i}}_{6},\bm{\mathrm{i}}_{7}\neq\bm{\mathrm{i}}_{8}}\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}Y_{\bm{\mathrm{i}}_{3}\bm{\mathrm{i}}_{4}}\right)\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{5}\bm{\mathrm{i}}_{6}}Y_{\bm{\mathrm{i}}_{7}\bm{\mathrm{i}}_{8}}\right)\\\
&\times\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i^{(n)}_{6},i^{(n)}_{5}}\delta_{i^{(n)}_{3},i^{(n)}_{4}}\delta_{i^{(n)}_{7},i^{(n)}_{8}}+\delta_{i^{(n)}_{2},i^{(n)}_{5}}\delta_{i^{(n)}_{6},i^{(n)}_{1}}\delta_{i^{(n)}_{3},i^{(n)}_{8}}\delta_{i^{(n)}_{7},i^{(n)}_{4}}\right)\\\
&-\frac{(\alpha\beta)^{2}}{2N^{2}}\sum_{\bm{\mathrm{i}}_{1}\neq\bm{\mathrm{i}}_{2},\bm{\mathrm{i}}_{5}\neq\bm{\mathrm{i}}_{6}}\sum_{\bm{\mathrm{k}}_{3},\bm{\mathrm{k}}_{4}}\left[\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}X(\bm{\mathrm{k}}_{3})\right)\operatorname{tr}\left(Y^{\dagger}_{\bm{\mathrm{i}}_{5}\bm{\mathrm{i}}_{6}}X(\bm{\mathrm{k}}_{4})\right)+\operatorname{tr}\left(X^{\dagger}(-\bm{\mathrm{k}}_{3})Y_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}\right)\operatorname{tr}\left(X^{\dagger}(-\bm{\mathrm{k}}_{4})Y_{\bm{\mathrm{i}}_{5}\bm{\mathrm{i}}_{6}}\right)\right]\\\
&\times\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i^{(n)}_{6},i^{(n)}_{5}}\delta_{k^{(n)}_{3},0}\delta_{k^{(n)}_{4},0}+\frac{1}{q}\delta_{i^{(n)}_{2},i^{(n)}_{5}}\delta_{i^{(n)}_{6},i^{(n)}_{1}}\delta_{k^{(n)}_{3},-k^{(n)}_{4}}\right).\end{aligned}$
(76)

Here $X(\bm{\mathrm{k}})$ represents the Fourier transform of the diagonal
component $X_{\bm{\mathrm{i}}\bm{\mathrm{i}}}$ with respect to
$\bm{\mathrm{i}}$ (Eq. (25)), and the $n$-th component of momentum
$\bm{\mathrm{k}}$ is denoted by $k^{(n)}$.

From Eq. (76), we obtain the self-energy for $X(\bm{\mathrm{k}})$ and $Y$ from
the quartic interactions, in the replica limit $R\rightarrow 0$,

$\displaystyle\begin{aligned}
(\Sigma_{X})^{ab;ba}(\bm{\mathrm{k}},\bm{\mathrm{k}})=&\frac{(\alpha\beta)^{2}}{N^{2}}\left\langle
X^{ab}(-\bm{\mathrm{k}})(X^{\dagger})^{ba}(-\bm{\mathrm{k}})\right\rangle_{0}\left[\prod_{n=1}^{L}\left(\delta_{k^{(n)},0}+\frac{1}{q^{2}}\right)-\delta_{\bm{\mathrm{k}},0}\right]\\\
=&\frac{(\alpha\beta)^{2}}{N}\frac{1}{1-\alpha\beta\delta_{\bm{\mathrm{k}},0}}\left[\left(1+\frac{1}{q^{2}}\right)^{r_{0}}\frac{1}{q^{2(L-r_{0})}}-\delta_{\bm{\mathrm{k}},0}\right],\\\
(\Sigma_{Y})^{ab;ba}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2};\bm{\mathrm{i}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{1}}=&\frac{(\alpha\beta)^{2}}{N^{2}}\sum_{\bm{\mathrm{i}}_{3}\neq\bm{\mathrm{i}}_{4}}\left\langle
Y^{ab}_{\bm{\mathrm{i}}_{3}\bm{\mathrm{i}}_{4}}(Y^{\dagger})^{ba}_{\bm{\mathrm{i}}_{4}\bm{\mathrm{i}}_{3}}\right\rangle_{0}\\\
&\times\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i^{(n)}_{3},i^{(n)}_{4}}\delta_{i_{1}^{\prime(n)},i_{2}^{\prime(n)}}+\delta_{i^{(n)}_{2},i^{(n)}_{3}}\delta_{i^{(n)}_{1},i^{(n)}_{4}}\delta_{i_{2}^{\prime(n)},i^{(n)}_{3}}\delta_{i_{1}^{\prime(n)},i^{(n)}_{4}}\right)\\\
=&\frac{(\alpha\beta)^{2}}{N}\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i_{1}^{\prime(n)},i_{2}^{\prime(n)}}+\frac{1}{q}\delta_{i^{(n)}_{2},i_{2}^{\prime(n)}}\delta_{i^{(n)}_{1},i_{1}^{\prime(n)}}\right)\\\
&-\frac{(\alpha\beta)^{2}}{N}\prod_{n=1}^{L}\left(\delta_{i^{(n)}_{2},i^{(n)}_{1}}\delta_{i_{1}^{\prime(n)},i_{2}^{\prime(n)}}+\frac{1}{q}\delta_{i^{(n)}_{2},i_{2}^{\prime(n)}}\delta_{i^{(n)}_{1},i_{1}^{\prime(n)}}\delta_{i^{(n)}_{1},i^{(n)}_{2}}\right).\end{aligned}$
(77)

Here we have used the expressions for the bare propagators Eqs.(27) and (29).
$r_{0}$ represents the number of zero components in the L-dimensional vector
$\bm{\mathrm{k}}$. We can therefore see that the self-energy for $X(0)$ is of
the order of $L/q^{2}N(1-\exp(i\Delta\phi))$, and the self-energy for
$X(\bm{\mathrm{k}})$ with $r_{0}<L$ zero components in $\bm{\mathrm{k}}$ is of
the order of $1/q^{2(L-r_{0})}N$. Furthermore, the self-energy
$(\Sigma_{Y})^{ab;ba}_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2};\bm{\mathrm{i}}^{\prime}_{2}\bm{\mathrm{i}}^{\prime}_{1}}$
is of the order of $1/(Nq^{L-r^{\prime}})$ where $r^{\prime}$ represents the
number of components satisfying $i_{2}^{(n)}=i_{1}^{(n)}$ and
$i_{2}^{\prime(n)}=i_{1}^{\prime(n)}$. Comparing Eqs. (77) with Eq. (31), we
find the self-energies for $X(\bm{\mathrm{k}})$ (with at least one zero
component in $\bm{\mathrm{k}}$) and
$Y_{\bm{\mathrm{i}}_{1}\bm{\mathrm{i}}_{2}}$ (when $i_{1}^{(n)}=i_{2}^{(n)}$
is satisfied for as least one site $n$) of the current noninteracting model
are much larger compared with their couterparts of the chaotic Floquet
circuits. Therefore, in contrast to the chaotic case, the quartic order
fluctuations for the current noninteracting model are essential.

##  D Large $q$ expansion of the Weingarten function

In this appendix, we study the asymptotic behavior of the Weingarten function
$g(\sigma,\tau)$ (Eq. (46)) for the unitary group $U(q)$ in the large
$q\rightarrow\infty$ limit. Expanding the action $S_{w}[W,W^{\dagger}]$ in
powers of $W$, we obtain

$\displaystyle\begin{aligned}
&S_{w}[W,W^{\dagger}]=q\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\operatorname{tr}\left(WW^{\dagger}\right)^{n},\end{aligned}$
(78)

where we have approximated the overall coefficient as $q+2R^{\prime}\approx
q$. One can see from rescaling $W$ by $W\rightarrow W/\sqrt{q}$ that the
higher order terms in the expansion above are also of higher order in $1/q$.

To the leading order in the large $q$ expansion, we can therefore apply the
Gaussian approximation and keep only the quadratic term in
$S_{w}[W,W^{\dagger}]$:

$\displaystyle
S_{w}^{(2)}[W,W^{\dagger}]=q\operatorname{tr}\left(WW^{\dagger}\right).$ (79)

The bare $W$ propagator is given by

$\displaystyle\begin{aligned} \left\langle
W_{uv}W^{\dagger}_{v^{\prime}u^{\prime}}\right\rangle_{0}=\delta_{uu^{\prime}}\delta_{vv^{\prime}}\frac{1}{q},\end{aligned}$
(80)

where the angular bracket with subscript $0$ denotes the functional averaging
with the Gaussian action $S_{w}^{(2)}$. Applying the Wick theorem to evaluate
the expectation value of $\left\langle
G[W,W^{\dagger};\sigma,\tau]\right\rangle_{0}$ (Eq. (46)), we find that, in
the Gaussian approximation, $g(\sigma,\tau)$ is nonvanishing only when
$\sigma=\tau$:

$\displaystyle\begin{aligned} &g(\sigma,\sigma)=\prod_{k=1}^{p}\left\langle
W_{k\sigma(k)}W^{\dagger}_{\sigma(k)k}\right\rangle_{0}=\frac{1}{q^{p}},\end{aligned}$
(81)

consistent with the known result for the Gaussian approximation to the
Weingarten function $\mathrm{Wg}(id)=g(\sigma,\sigma)$ Beenakker . The
Weingarten functions for all the remaining permutations $\tau^{-1}\sigma\neq
id$ are of higher order in $1/q$.

To deduce the leading order term in the large $q$ expansion of
$g(\sigma\neq\tau)$, we should go beyond the Gaussian approximation and keep
the higher order terms in $S_{w}[W,W^{\dagger}]$. Using Eq. (53), we find

$\displaystyle\begin{aligned}
&g(\sigma,\tau)\propto\prod_{k=1}^{m}\left\langle\left(\prod_{l=1}^{c_{k}}W_{b_{k}+l,b_{k}+(l+1\,\mathrm{mod}\,c_{k})}W^{\dagger}_{b_{k}+(l+1\,\mathrm{mod}\,c_{k}),b_{k}+(l+1\,\mathrm{mod}\,c_{k})}\right)\left[q\operatorname{tr}\left(WW^{\dagger}\right)^{c_{k}}\right]\right\rangle_{0}\\\
&\propto\prod_{k=1}^{m}\left(q\prod_{l=1}^{c_{k}}\left\langle
W_{b_{k}+l,b_{k}+(l+1\,\mathrm{mod}\,c_{k})}W^{\dagger}_{b_{k}+(l+1\,\mathrm{mod}\,c_{k}),b_{k}+l}\right\rangle_{0}\right.\\\
&\qquad\quad\times\left.\left\langle
W^{\dagger}_{b_{k}+(l+1\,\mathrm{mod}\,c_{k}),b_{k}+(l+1\,\mathrm{mod}\,c_{k})}W_{b_{k}+(l+1\,\mathrm{mod}\,c_{k}),b_{k}+(l+1\,\mathrm{mod}\,c_{k})}\right\rangle_{0}\right)\\\
&\propto\prod_{k=1}^{m}q^{1-2c_{k}},\end{aligned}$ (82)

which leads to the known result $g(\sigma,\tau)\propto q^{m-2p}$ Beenakker .
The extra factor of $q\operatorname{Tr}\left(WW^{\dagger}\right)^{c_{k}}$ in
the first line of the equation above represents the interaction vertex arising
from the $c_{k}$-th order term in the expansion of $S_{w}[W,W^{\dagger}]$ (Eq.
(78)).

##  E Recursion relation for the Weingarten function

In this appendix, we will show that the function $g(\sigma,\tau)$ given by Eq.
(46) satisfies the recursion relations for the Weingarten function, i.e., Eq.
(37). To derive the first recursion relation Eq. (37a), we make use of the
identity $UU^{\dagger}=\bf{1}$, which leads to

$\displaystyle\begin{aligned}
&\int_{\mathcal{U}(q)}dUU_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{p}j_{p}}U^{\dagger}_{j_{1}^{\prime}i_{1}^{\prime}}U^{\dagger}_{j_{2}^{\prime}i_{2}^{\prime}}...U^{\dagger}_{j^{\prime}_{p}i^{\prime}_{p}}\\\
=&\\!\\!\\!\\!\\!\\!\\!\\!\sum_{i_{p+1},i^{\prime}_{p+1},j_{p+1}=1}^{q}\\!\\!\\!\\!\\!\\!\\!\\!\delta_{i_{p+1},i^{\prime}_{p+1}}\delta_{j_{p+1},j^{\prime}_{p+1}}\int_{\mathcal{U}(q)}dUU_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{p}j_{p}}U_{i_{p+1}j_{p+1}}U^{\dagger}_{j_{1}^{\prime}i_{1}^{\prime}}U^{\dagger}_{j_{2}^{\prime}i_{2}^{\prime}}...U^{\dagger}_{j^{\prime}_{p}i^{\prime}_{p}}U^{\dagger}_{j^{\prime}_{p+1}i^{\prime}_{p+1}}.\end{aligned}$
(83)

Using Eq. (45), the left hand side (L.H.S.) of the equation above can be
expressed as

$\displaystyle\begin{aligned}
&\text{L.H.S}.=\sum_{\tau^{\prime},\sigma^{\prime}\in
S_{p}}\left(\prod_{k=1}^{p}\delta_{i_{k},i_{\sigma^{\prime}(k)}^{\prime}}\delta_{j_{k},j_{\tau^{\prime}(k)}^{\prime}}\right)g(\sigma^{\prime},\tau^{\prime}),\end{aligned}$
(84)

while the right hand side (R.H.S.) acquires the form

$\displaystyle\begin{aligned}
&\text{R.H.S.}=\\!\\!\\!\\!\\!\\!\\!\\!\sum_{i_{p+1},i^{\prime}_{p+1},j_{p+1}=1}^{q}\sum_{\tau^{\prime\prime},\sigma^{\prime\prime}\in
S_{p+1}}\delta_{i_{p+1},i^{\prime}_{p+1}}\delta_{j_{p+1},j^{\prime}_{p+1}}\left(\prod_{k=1}^{p+1}\delta_{i_{k},i_{\sigma^{\prime\prime}(k)}^{\prime}}\delta_{j_{k},j_{\tau^{\prime\prime}(k)}^{\prime}}\right)g(\sigma^{\prime\prime},\tau^{\prime\prime}).\end{aligned}$
(85)

We then choose an arbitrary pair of permutations $\sigma,\tau\in S_{p}$, and
set $i_{k}=i_{\sigma(k)}^{\prime},j_{k}=j_{\tau(k)}^{\prime}$, for
$k=1,2,...,p$. For simplicity, we consider the case where $p<q$ and all
integers in the set $\left\\{i^{\prime}_{k}|k=1,...,p\right\\}$
($\left\\{j^{\prime}_{k}|k=1,...,p+1\right\\}$) are different from each other.
In this case, the L.H.S. then becomes

$\displaystyle\begin{aligned}
&\text{L.H.S.}=g(\sigma,\tau)=V^{(p)}_{c_{1},c_{2},...,c_{m}},\end{aligned}$
(86)

where $\left\\{c_{k}|k=1,2,...,m\right\\}$ represent the lengths of disjoint
cycles for the permutation $\tau^{-1}\sigma$ and obey the constraint
$\sum_{k=1}^{m}c_{k}=p$. Here we made use of the result that $g(\sigma,\tau)$
depends only on $\left\\{c_{k}\right\\}$ and expressed it as
$V^{(p)}_{c_{1},c_{2},...,c_{m}}$ (see Sec. 5.2). On the other hand, the
R.H.S. becomes

$\displaystyle\begin{aligned}
\text{R.H.S.}=&\sum_{\tau^{\prime\prime},\sigma^{\prime\prime}\in
S_{p+1}}\sum_{i^{\prime}_{p+1}=1}^{q}\left(\prod_{k=1}^{p}\delta_{i_{\sigma(k)}^{\prime},i_{\sigma^{\prime\prime}(k)}^{\prime}}\delta_{j_{\tau(k)}^{\prime},j_{\tau^{\prime\prime}(k)}^{\prime}}\right)\delta_{i^{\prime}_{p+1},i_{\sigma^{\prime\prime}(p+1)}^{\prime}}\delta_{j^{\prime}_{p+1},j_{\tau^{\prime\prime}(p+1)}^{\prime}}g(\sigma^{\prime\prime},\tau^{\prime\prime}).\end{aligned}$
(87)

Let us now examine all terms in the summation in Eq. (87). The contribution is
nonzero for arbitrary $i_{p+1}^{\prime}$ if $\sigma^{\prime\prime}=\sigma_{p}$
and $\tau^{\prime\prime}=\tau_{p}$, where $\sigma_{p}$ and $\tau_{p}$ are
defined as

$\displaystyle\begin{aligned}
\sigma_{p}(k)=\begin{cases}\sigma(k),&k=1,2,...,p,\\\
p+1,&k=p+1,\end{cases}\qquad\tau_{p}(k)=\begin{cases}\tau(k),&k=1,2,...,p,\\\
p+1,&k=p+1.\end{cases}\end{aligned}$ (88)

The permutation $\tau^{-1}_{p}\sigma_{p}$, compared with $\tau^{-1}\sigma$,
contains an extra disjoint cycle of length $1$. More specifically, from
earlier result (Eq. (50)), we find that $g(\sigma_{p},\tau_{p})$ is given by

$\displaystyle\begin{aligned}
&g(\sigma_{p},\tau_{p})=(-1)^{p-m}\left\langle\left(W_{p+1,p+1}W^{\dagger}_{p+1,p+1}\right)\prod_{k=1}^{m}\left(\prod_{l=1}^{c_{k}}W_{P_{l}^{(k)},\sigma(P_{l}^{(k)})}W^{\dagger}_{\sigma(P_{l}^{(k)}),P_{l+1}^{(k)}}\right)\right\rangle_{W}\\\
=&V^{(p+1)}_{c_{1},c_{2},...c_{m},1}.\end{aligned}$ (89)

Compared with $g(\sigma,\tau)=V^{(p)}_{c_{1},c_{2},...c_{m}}$, it contains an
extra factor of $W_{p+1,p+1}W^{\dagger}_{p+1,p+1}$. Here the integers
$\left\\{P_{l}^{(k)}\right\\}$ denote the cycle structure of $\tau^{-1}\sigma$
as in Eq. (36). The total contribution from this pair of permutations
($\sigma_{p},\tau_{p}$) to Eq. (87) is

$\displaystyle\begin{aligned}
\sum_{i_{p+1}^{\prime}=1}^{q}g(\sigma_{p},\tau_{p})=qV^{(p+1)}_{c_{1},c_{2},...,c_{m},1},\end{aligned}$
(90)

where the factor $q$ comes from the summation over $i_{p+1}^{\prime}$.

Another nonvanishing contribution comes from the case when
$i_{p+1}^{\prime}=i_{\sigma(l)}^{\prime}$ for an arbitrary positive integer
$l\leq p$. In this case, the contribution is nonvanishing if the permutations
$\sigma^{\prime\prime}=\sigma_{l}$ and $\tau^{\prime\prime}=\tau_{p}$, where
$\tau_{p}$ is defined in Eq. (88) and $\sigma_{l}$ is given by

$\displaystyle\begin{aligned} &\sigma_{l}(k)=\begin{cases}\sigma(k),&k\neq
l,p+1,\\\ p+1,&k=l,\\\ \sigma(l),&k=p+1.\end{cases}\end{aligned}$ (91)

In this case, it is easy to see that
$\delta_{i^{\prime}_{p+1},i_{\sigma_{l}(p+1)}^{\prime}}=\delta_{i^{\prime}_{\sigma(l)},i_{\sigma_{l}(l)}^{\prime}}=\delta_{i^{\prime}_{p+1},i_{\sigma(l)}^{\prime}}=1$.

$g(\sigma_{l},\tau_{p})$ can be obtained from $g(\sigma,\tau)$ by replacing
the factor $W_{l,\sigma(l)}W^{\dagger}_{\sigma(l),\tau^{-1}\sigma(l)}$ in Eq.
(50) with
$W_{l,p+1}W^{\dagger}_{p+1,p+1}W_{p+1,\sigma(l)}W^{\dagger}_{\sigma(l),\tau^{-1}\sigma(l)}$.
This results in adding a number to the disjoint cycle which number $l$ belongs
(labeled by $j$):

$\displaystyle\begin{aligned} g(\sigma,\tau)=&\left\langle
W_{l,\sigma(l)}W^{\dagger}_{\sigma(l),\tau^{-1}\sigma(l)}...\right\rangle_{W}=V_{c_{1},...,c_{j},...,c_{m}}^{(p)},\\\
g(\sigma_{l},\tau_{p})=&\left\langle
W_{l,\sigma_{l}(l)}W^{\dagger}_{\sigma_{l}(l),\tau_{p}^{-1}\sigma_{l}(l)}W_{\tau_{p}^{-1}\sigma_{l}(l),\sigma_{l}\tau_{p}^{-1}\sigma_{l}(l)}W^{\dagger}_{\sigma_{l}\tau_{p}^{-1}\sigma_{l}(l),\tau_{p}^{-1}\sigma_{l}\tau_{p}^{-1}\sigma_{l}(l)}...\right\rangle_{W},\\\
=&\left\langle
W_{l,p+1}W^{\dagger}_{p+1,p+1}W_{p+1,\sigma(l)}W^{\dagger}_{\sigma(l),\tau^{-1}\sigma(l)}...\right\rangle_{W}=V_{c_{1},...,c_{j}+1,...,c_{m}}^{(p+1)}.\end{aligned}$
(92)

Here “$...$” represents all remaining terms that stay the same in
$g(\sigma_{l},\tau_{p})$ and $g(\sigma,\tau)$. Summing over all possibilities
of $l$, we find the total contribution of this type is

$\displaystyle\begin{aligned}
\sum_{l}g(\sigma_{l},\tau_{p})=\sum_{s=1}^{m}c_{s}V_{c_{1},...,c_{s-1},c_{s}+1,c_{s+1},...,c_{m}}^{(p+1)}.\end{aligned}$
(93)

Combining everything (Eqs. (90) and (93)), we find the R.H.S. of Eq. (83)

$\displaystyle\begin{aligned}
\text{R.H.S.}=qV_{c_{1},...,c_{m},1}^{(p+1)}+\sum_{s=1}^{m}c_{s}V_{c_{1},...,c_{s-1},c_{s}+1,c_{s+1},...,c_{m}}^{(p+1)}.\end{aligned}$
(94)

Equating this to the L.H.S. (Eq, (86)), we arrive at the recursion relation
Eq. (37a) which holds for arbitrary pair of $\sigma,\tau\in S_{p}$ (or
equivalently arbitrary configuration of $\left\\{c_{j}\right\\}$ that
satisfies $\sum_{j=1}^{m}c_{j}=p$) when $p<q$.

The recursion relation Eq. (37b) can be proved in a similar way using the
identity:

$\displaystyle\begin{aligned}
\sum_{j_{p},j_{p}^{\prime}=1}^{q}\delta_{j_{p},j_{p}^{\prime}}\int_{\mathcal{U}(q)}dUU_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{p-1}j_{p-1}}U_{i_{p}j_{p}}U^{\dagger}_{j_{1}^{\prime}i_{1}^{\prime}}U^{\dagger}_{j_{2}^{\prime}i_{2}^{\prime}}...U^{\dagger}_{j_{p-1}^{\prime}i_{p-1}^{\prime}}U^{\dagger}_{j_{p}^{\prime}i_{p}^{\prime}}=0,\end{aligned}$
(95)

for $i_{p}\neq i_{p}^{\prime}$. Using Eq. (45), one can rewrite this equation
as

$\displaystyle\begin{aligned}
&\sum_{j_{p},j_{p}^{\prime}=1}^{q}\sum_{\tau^{\prime},\sigma^{\prime}\in
S_{p}}\delta_{j_{p},j_{p}^{\prime}}\left(\prod_{k=1}^{p}\delta_{i_{k},i_{\sigma^{\prime}(k)}^{\prime}}\delta_{j_{k},j_{\tau^{\prime}(k)}^{\prime}}\right)g(\sigma^{\prime},\tau^{\prime})=0.\end{aligned}$
(96)

We now set $j_{k}=j_{k}^{\prime}$ and $i_{k}=i_{\chi(k)}^{\prime}$ for $1\leq
k\leq p$. $\chi\in S_{p}$ represents an arbitrary nonidentical permutation of
numbers $1,2,...,p$ which satisfies $\chi(p)\neq p$. It contains $m$ disjoint
cycles with lengths $\left\\{c_{k}|k=1,2,...,m\right\\}$. Since $\chi\neq id$
and $\chi(p)\neq p$, $p$ stays in a cycle (labeled by $1$) of length
$c_{1}\geq 2$. We consider the case $p\leq q$ and choose the set
$\left\\{i^{\prime}_{k}|k=1,...,p\right\\}$
($\left\\{j^{\prime}_{k}|k=1,...,p-1\right\\}$) whose elements are all
different from each other for simplicity. Eq. (96) can be rewritten as

$\displaystyle\begin{aligned}
&\sum_{j_{p}^{\prime}=1}^{q}\sum_{\tau^{\prime},\sigma^{\prime}\in
S_{p}}\left(\prod_{k=1}^{p}\delta_{i_{\chi(k)}^{\prime},i_{\sigma^{\prime}(k)}^{\prime}}\delta_{j_{k}^{\prime},j_{\tau^{\prime}(k)}^{\prime}}\right)g(\sigma^{\prime},\tau^{\prime})=0.\end{aligned}$
(97)

Let us now examine all terms in the summation on the L.H.S. of Eq. (97). For
arbitrary $j_{p}^{\prime}$, the contribution is nonvanishing when
$\sigma^{\prime}=\chi$ and $\tau^{\prime}=id$. The total contribution of this
type is

$\displaystyle\begin{aligned}
\sum_{j_{p}^{\prime}=1}^{q}g(\chi,id)=qg(\chi,id)=qV_{c_{1},c_{2},...,c_{m}}^{(p)}.\end{aligned}$
(98)

Another type of nonvanishing contributions comes from the case where
$j_{p}^{\prime}=j_{l}^{\prime}$ and $l$ also stays within the same cycle that
$p$ belongs (cycle $1$). In this case, the contribution is nonvanishing if
$\sigma^{\prime}=\chi$ and $\tau=\tau_{l}$, where $\tau_{l}$ is defined as

$\displaystyle\begin{aligned} \tau_{l}(k)=\begin{cases}k,&k\neq p,l,\\\
l,&k=p,\\\ p,&k=l.\end{cases}\end{aligned}$ (99)

$g(\chi,\tau_{l})$ can be obtained from $g(\chi,id)$ by replacing
$W^{\dagger}_{ll}$ and $W^{\dagger}_{pp}$ in $G[W,W^{\dagger};\chi,id]$ (Eq.
(46)) with $W^{\dagger}_{lp}$ and $W^{\dagger}_{pl}$. This breaks the cycle
$1$ of length $c_{1}$ into two disjoint cycles of lengths $c_{1}-c$ and $c$
($c$ obeys $0<c<c_{1}$ and depends on $l$):

$\displaystyle\begin{aligned} &g(\chi,id)=V^{(p)}_{c_{1},c_{2},...,c_{m}}\\\
=&\left\langle\left(W_{p,\chi(p)}W^{\dagger}_{\chi(p),\chi(p)}...W_{\chi^{-1}(l),l}W^{\dagger}_{l,l}W_{l,\chi(l)}W^{\dagger}_{\chi(l),\chi(l)}...W_{\chi^{-1}(p),p}W^{\dagger}_{p,p}\right)\prod_{k\neq
1}^{m}(...)\right\rangle_{W},\\\
&g(\chi,\tau_{l})=V^{(p)}_{c_{1}-c,c,c_{2},...,c_{m}}\\\
=&\left\langle\left(W_{p,\chi(p)}W^{\dagger}_{\chi(p),\chi(p)}...W_{\chi^{-1}(l),l}W^{\dagger}_{l,p}\right)\left(W_{l,\chi(l)}W^{\dagger}_{\chi(l),\chi(l)}...W_{\chi^{-1}(p),p}W^{\dagger}_{p,l}\right)\prod_{k\neq
1}^{m}(...)\right\rangle_{W}.\end{aligned}$ (100)

Here $\prod_{k\neq 1}^{m}(...)$ represents the contribution from all remaining
cycles other than $1$, and stays the same for both $g(\chi,id)$ and
$g(\chi,\tau_{l})$. Taking into account all possible $l$ belonging to the
cycle $1$, we find that the total contribution from such terms is

$\displaystyle\begin{aligned} \sum_{l\in\text{cycle
1}}g(\chi,\tau_{l})=\sum_{c=1}^{c_{1}-1}V_{c_{1}-c,c,c_{2},...,c_{m}}^{(p)}.\end{aligned}$
(101)

The last type of nonvanishing contribution arises from the case when
$j_{p}^{\prime}=j_{l}^{\prime}$ and $l$ does not belong to the cycle $1$ (we
label the cycle which $l$ belongs by $2$). In this case, the contribution is
nonvanishing if $\sigma^{\prime}=\chi$ and $\tau^{\prime}=\tau_{l}$ (Eq.
(99)). As mentioned earlier, $g(\chi,\tau_{l})$ can be obtained from
$g(\chi,id)$ by replacing $W^{\dagger}_{ll}$ and $W^{\dagger}_{pp}$ in
$G[W,W^{\dagger};\chi,id]$ with $W^{\dagger}_{lp}$ and $W^{\dagger}_{pl}$,
respectively. However, since now $l,p$ belong to different cycles, this
replacement combines the cycle $1$ and cycle $2$ into one cycle of length
$c_{1}+c_{2}$:

$\displaystyle\begin{aligned}
&g(\chi,id)=V_{c_{1},c_{2},c_{3},...,c_{m}}^{(p)}\\\
=&\left\langle\left(W_{p,\chi(p)}W^{\dagger}_{\chi(p),\chi(p)}...W_{\chi^{-1}(p),p}W^{\dagger}_{p,p}\right)\left(W_{l,\chi(l)}W^{\dagger}_{\chi(l),\chi(l)}...W_{\chi^{-1}(l),l}W^{\dagger}_{l,l}\right)\prod_{k>2}^{m}(...)\right\rangle,\\\
&g(\chi,\tau_{0})=V_{c_{1}+c_{2},c_{3},...,c_{m}}^{(p)}\\\