††institutetext: Laboratoire de Physique Théorique et Modélisation, CNRS UMR
8089, CY Cergy Paris Université, 95302 Cergy-Pontoise Cedex, France

# Frame potential of Brownian SYK model of Majorana and Dirac fermions

Anastasiia Tiutiakina 1 Andrea De Luca 1 Jacopo De Nardis
<EMAIL_ADDRESS>
###### Abstract

We consider the Brownian SYK, i.e. a system of $N$ Majorana (Dirac) fermions
with a white-noise $q$-body interaction term. We focus on the dynamics of the
Frame potentials, a measure of the scrambling and chaos, given by the moments
of the overlap between two independent realisations of the model. By means of
a Keldysh path-integral formalism, we compute its early and late-time value.
We show that, for $q>2$, the late time path integral saddle point correctly
reproduces the saturation to the value of the Haar frame potential. On the
contrary, for $q=2$, the model is quadratic and consistently we observe
saturation to the Haar value in the restricted space of Gaussian states
(gaussian Haar). The latter is characterised by larger system size corrections
that we correctly capture by counting the Goldstone modes of the Keldysh
saddle point.

##  1 Introduction and main results

The phenomenon of information scrambling in strongly interacting many-body
systems during unitary time evolution has gained significant attention due to
its deep connections with holography, black hole physics, and quantum
information Hayden2007 ; Shenker2014 ; Shenker2014-2 ; Maldacena2016 ;
PRXQuantum.2.010306 ; Almheiri2020 ; PhysRevX.9.031048 ; jian2022linear .
Recently, a related question has emerged regarding $k$-designs, i.e. ensembles
of quantum systems able to approximate, up to the $k$-th moment, uniformly
(Haar) distributed quantum states through unitary time evolution from a pure
and simple initial state Gross2007 ; Roberts2017 ; Hearth2023 ;
PhysRevX.13.011043 . Such a question is strictly related to the emergence of
eigenstate thermalization (ETH) in closed quantum systems PhysRevE.50.888 ;
DAlessio2016 ; PhysRevLett.129.170603 ; PhysRevE.57.7313 ; PhysRevA.43.2046 .
More recently, a related problem has been posed in the context of
deterministic Hamiltonian systems undergoing final read-out measurements: the
randomness induced by quantum measurements can approximate with different
accuracy Haar random states, a feature recently dubbed deep thermalization
Lucas2022 ; Hearth2023 ; PRXQuantum.4.010311 ; Choi2023 . Generic chaotic
quantum systems are indeed expected to thermalise at late time, i.e. the
reduced density matrices obtained by tracing out most of the degrees of
freedom are expected to converge to the standard Gibbs statistical density
matrix, which reduces to the identity matrix in the absence of conserved
quantities. However, the presence of different sources of stochasticity (e.g.
noise or measurements), raises the question whether, beyond the convergence of
the stochastic average, individual realisation (or trajectories) of the system
are exploring uniformly the accessible Hilbert space of pure states. However,
the study of strongly interacting system, in particular with the additional
ingredient of stochasticity, is generally analytically (and also numerically)
very challenging, except for a few specific systems made of random matrices or
dual unitary quantum circuits Ho2022 ; Ho2023 ; Ippoliti2022 . For this
reason, it is important to develop an intuition on some sufficiently-generic
but tractable cases. Among these, a prominent example is the Brownian SYK
model PhysRevLett.70.3339 ; 1806.06840 , describing $L$ Majorana (or Dirac)
fermions with all-to-all interactions (as for the standard SYK), with a time-
dependent white-noise interaction coupling. It serves as a paradigmatic model
for quantum scrambling and chaos Jian2021 ; 1806.06840 ; Snderhauf2019 ;
PhysRevE.78.021106 ; PhysRevX.9.031048 ; PhysRevB.100.064305 .

In this paper, we focus on investigating the degree of mixing in the time
evolution of the Brownian SYK model. To address this question, we employ a
measure of mixing called the Frame potential. Specifically, we consider two
identical copies of an initial density matrix $\rho_{0}$ and compute the
overlap between the time evolutions of these copies under two independent
realisations of the noise. The Frame potential, denoted as $F^{(k)}(T)$, is
given by the expression:

$F^{(k)}(T)=\left<{\rm
Tr}\left(U_{1}(T)\rho_{0}U_{1}^{\dagger}(T)U_{2}(T)\rho_{0}U_{2}^{\dagger}(T)\right)^{k}\right>,$
(1)

Here, $U_{1}$ and $U_{2}$ represent the (noisy) evolution operators of two
independent Brownian SYK systems with randomly chosen couplings and $F^{(k)}$
equals the $k$-th moment (averaged over both noise realizations) of their
overlap. In the case of sufficiently mixing time evolutions and starting from
a reference pure state $\rho_{0}=\ket{\Psi}\bra{\Psi}$, the Frame potential
exhibits an exponential decay from its initial value one. This decay manifests
the decreasing probability that the two density matrices remain close as time
progresses. Eventually, at infinite time $T\to\infty$, it is expected that the
density matrix $\rho(T)$ fully explores the entire Hilbert space in a uniform
manner, implying that its statistics can be effectively obtained replacing the
time-evolution operator with a Haar distributed unitary operator in the
Hilbert space $\mathcal{H}$ as

$\rho(T)\to U^{\dagger}\rho_{0}U,\qquad U\sim\mbox{Haar}(\mathcal{H}).$ (2)

Consequently, the Frame potentials for all $k$’s are expected to converge at
large times to their Haar-averaged value, i.e. where both $U_{1}(T)$ and
$U_{2}(T)$ are replaced by two independently drawn Haar-distributed unitary
matrices. We denote this limit case as $F_{\rm Haar}^{(k)}$ and it is easy to
prove (see for instance Ippoliti2022 ) that it corresponds to the minimum over
all possible distributions of unitary operators on $\mathcal{H}$. As detailed
in this paper, this convergence indeed occurs for any $q>2$ (with interactions
coupling more than 2 fermions), thus generalising the results of
jian2022linear . An exception arises in the Gaussian integrable case $q=2$,
with exactly two-fermions interaction 10.21468/SciPostPhys.12.1.042 ;
Bernard2021 ; swann2023spacetime . In this case, the evolution only exhibits
mixing behaviour within the space of Gaussian states. There, we can introduce
a Gaussian Haar ( gHaar) PhysRevE.104.014146  measure and demonstrate that the
Frame potentials converge for all $k$ to the latter at late times.
Intriguingly, the value of the gHaar-integrated Frame potential is equivalent
to the generic case in the limit of a large fermion number $L$, with
logarithmic corrections in $L$. As demonstrated in this work, these
corrections arise from fluctuations of Goldstone modes around the late-time
saddle point, which are only present in the Gaussian case where both the
Keldysh action and the saddle point’s symmetries are generated by continuous
groups.

Figure 1: The rescaled Frame potential in the logarithmic scale. Number of
fermions in both cases $L=20$. Brown colour corresponds to the number of
replicas $k=3$ , green colour corresponds to $k=2$. Dotted line is the
theoretical prediction of the Keldysh calculation at early time, and dashed
line is the prediction for the late time behaviour. From the numerical data,
we can see that the time of the saturation can be estimated directly from the
theoretical prediction. The saturation time for $q=2$: $T^{q=2}_{\rm
Dirac}=4\left(2\log(2)-\frac{k}{L}\log(L)\right)+O\left(\frac{1}{L}\right),~{}~{}T^{q=2}_{\rm
Majorana}=4\left(2\log(2)-\frac{k-1}{L}\log(L)\right)+O\left(\frac{1}{L}\right);$
and for $q>2$; $T^{q>2}_{\rm Dirac}=T^{q>2}_{\rm Majorana}=q^{2}\cdot
2\log(2)+O\left(\frac{1}{L}\right)$.

The structure of this paper is as follows. In Section 2, we derive the Frame
potential for the Brownian Majorana SYK system using the framework of the
thermofield double state. We outline the Keldysh path integral approach, which
facilitates averaging over the Brownian noise, and obtain an effective action
for the problem. Moving to Section 3, we explore the Swinger-Dyson equations
by performing a variation of the action. Initially, we consider a naive saddle
point solution that does not mix different Keldysh contours and replicas. The
resulting Frame potential exhibits exponential decay, as depicted by the
dotted lines in Fig. (1). Subsequently, we investigate a non-trivial solution
in the large-time limit, wherein different replicas and contours are mixed.
Notably, this solution exhibits distinct group symmetries for cases with $q>2$
and $q=2$. The continuous symmetry group characterizes the $q=2$ case, while
discrete symmetries are present for $q>2$. We then proceed to analyse the
fluctuations around the saddle points. For the $q>2$ case, these fluctuations
trivially cancel with the normalization. However, in the $q=2$ case,
fluctuations governed by the Goldstone theorem lead to massless modes in the
action. Consequently, the integration over these fluctuations deviates from
the normalization, resulting in logarithmic corrections ($\log(L)$) to the
Frame potential. During the evolution, we observe that the Frame potential
reaches a saturation value exponentially small with respect to the system
size. This behaviour is illustrated by the dashed lines in Fig.(1). In Section
4, we extend our analysis to the complex Brownian SYK model. Here, we repeat
the calculations performed in the Majorana case and find an analogous
exponential decay of the Frame potential (depicted by dotted lines in Fig.
(1)) and subsequent saturation (depicted by dashed lines in the same figure).
Finally, in Section 5, we conclude the paper by summarizing our findings. We
also discuss potential future directions and prospects for further exploration
within this research framework.

##  2 Brownian SYK model and Keldysh path integral representation of the Frame
potential

We consider the time evolution generated by the time-dependent Hamiltonian

$H(h,t)=i^{\frac{q}{2}}\sum_{1\leq i_{1}<i_{2}<..<i_{q}\leq
L}h_{i_{1}...i_{q}}(t)\hat{\chi}_{i_{1}}...\hat{\chi}_{i_{q}},$ (3)

where $\hat{\chi}_{i}$ are $L$ Majorana fermions and anticommutation relations
are applied $\\{\hat{\chi}_{i},\hat{\chi}_{j}\\}=\delta_{ij}$, here $q$ is an
even integer constant. This model has close analogies with the celebrated SYK
Sachdev:1992fk ; Gu:2019jub ; 1806.06840 , but with the important difference
that $h_{i_{1}\cdots i_{q}}$ are not constant in time here. On the contrary,
we take them to follow a white noise distribution,

$\left<h_{\textbf{i}}(t_{1})h_{\textbf{j}}(t_{2})\right>=\frac{2^{q-1}q!}{q^{2}L^{q-1}}\delta_{\textbf{i},\textbf{j}}\delta(t_{1}-t_{2})=\sigma^{2}\delta_{\textbf{i},\textbf{j}}\delta(t_{1}-t_{2}),$
(4)

where we denote the collective set of indices
$\textbf{i}=(i_{1},\ldots,i_{q})$ (similarly for j) and set
$\delta_{\textbf{i},\textbf{j}}=\delta_{i_{1},j_{1}}\ldots\delta_{i_{q},j_{q}}$.
For this reason, the current model is named Brownian SYK. We shall here assume
$L$ to be an even integer, so that the operators $\hat{\chi}_{i}$ can be
embedded into a $\mathcal{N}=2^{L/2}$ dimensional Hilbert space. We are
interested in studying the scrambling dynamics induced by the time evolution
(3). In order to do so, we focus on a specific initial condition, known as
thermofield double $\left(\left|{\rm TFD}\right>\right)$ Gu2017SpreadOE ;
1911.11977 ; Almheiri2020 , which has already been employed in the context of
Brownian SYK in Jian2021 . In practice, we consider two copies of the system,
prepared in a maximally entangled state. In the following, we shall address
the two copies as left (L) and right (R) halves. The selection of the
$\left|{\rm TFD}\right>$ state is not uniquely determined and relies on the
choice of a basis. However, different definitions yield states that are
connected through local unitary transformations. To provide a concrete
illustration, we offer an explicit characterization of $\left|{\rm
TFD}\right>$ in the context of the SYK model. Within this framework, the
Majorana fermion operators are denoted as $\hat{\chi}_{j,L}$ and
$\hat{\chi}_{j,R}$, where $j$ takes values from $1$ to $L$. (It should be
noted that $L$ must be an even number to ensure a well-defined Hilbert space).
A convenient specification for the state can be formulated using Dirac
fermions

$\left(\hat{c}_{j,L}-\hat{c}^{\dagger}_{j,R}\right)\left|{\rm
TFD}\right>=0~{}~{}\text{and}~{}~{}\left(\hat{c}^{\dagger}_{j,L}-\hat{c}_{j,R}\right)\left|{\rm
TFD}\right>=0,$ (5)

which can be expressed in the terms of Majorana fermions as

$\hat{c}_{j,L}=\frac{1}{\sqrt{2}}\left(\hat{\chi}_{2j-1,L}+i\hat{\chi}_{2j,L}\right)~{}~{}~{}\text{and}~{}~{}~{}\hat{c}_{j,R}=\frac{1}{\sqrt{2}}\left(\hat{\chi}_{2j,R}-i\hat{\chi}_{2j-1,R}\right).$
(6)

Then we have the standard definition in the terms of Majorana fermions:

$\left(\hat{\chi}_{j,L(R)}+i\hat{\chi}_{j,L(R)}\right)\left|{\rm
TFD}\right>=0,~{}~{}~{}\forall j=1,...,L,$ (7)

From the Eq. (5) we see that the ${\rm TFD}$ state in the eigenbasis of
$n_{j,L(R)}=\hat{c}^{\dagger}_{j,L(R)}\hat{c}_{j,L(R)}$ basis can be expressed
as a product of
$\left|0\right>_{L}\left|0\right>_{R}+\left|1\right>_{L}\left|1\right>_{R}$ on
each site. More explicitly, we consider the initial state:

$\left|{\rm
TFD}\right>=\frac{1}{\sqrt{\mathcal{N}}}\sum_{j=1}^{\mathcal{N}}\ket{j}\ket{j},~{}~{}~{}~{}~{}\rho_{0}=\left|{\rm
TFD}\right>\left<{\rm TFD}\right|,$ (8)

where each $\left|j\right>$ represents a possible string of zeros and ones.
Here $\mathcal{N}$ is the dimension of the Hilbert space. This state has the
following property,

$\operatorname{Tr}^{(2)}[A\otimes B\left|{\rm TFD}\right>\left<{\rm
TFD}\right|C\otimes D]=\frac{1}{\mathcal{N}}\operatorname{Tr}[AB^{t}D^{t}C],$
(9)

where the trace on the left is on the doubled Hilbert space, while the one on
the right is on the single copy. This identity holds for arbitrary operators
$A,B/C,D$ acting respectively on the first/second copy of the Hilbert space.
Except for the initial entanglement, the two halves evolve according to two
uncoupled unitary operators,

$\left|{\rm TFD}\right>(t)=U_{L}(t)\otimes U_{R}(t)\left|{\rm
TFD}\right>,\qquad\rho(t)=\left|{\rm TFD}(t)\right>\left<{\rm TFD}(t)\right|,$
(10)

where $U_{L}$ and $U_{R}$ are the unitary time evolution operators acting on
the two halves. We assume that the time evolution in each half is generated by
an independent realisation of the Brownian SYK

$U_{h^{\sigma}}(t+dt)=e^{-iH(h^{\sigma},t)dt}U_{h^{\sigma}}(t),$ (11)

where the subscript $\sigma=L,R$ labels the corresponding half and
$H_{L}(t),H_{R}(t)$ have the form of Eq. (3) in terms of two sets of Majoranas
$\hat{\chi}_{i,\sigma}$ and independently generated white noises
$h^{\sigma}_{\textbf{i}}(t)$. Starting with this initial state, we calculate
the $k$-th moments of the Frame potential averaged over the noise:

$F^{(k)}(T)=\int d\mu\left(h_{1}\right)d\mu\left(h_{2}\right){\rm
Tr}\left(U_{h_{1}}(T)\rho_{0}U^{\dagger}_{h_{1}}(T)U_{h_{2}}(T)\rho_{0}U^{\dagger}_{h_{2}}(T)\right)^{k},$
(12)

here, $U_{h_{1}}(T)=U_{h^{L}_{1}}(T)\otimes U^{t}_{h^{R}_{1}}(T)$, (and the
same for $h^{L(R)}_{2}$), where the couplings $h^{L(R)}_{1(2)}$ are
independent random variables with the distribution
$d\mu\left(h\right)=e^{-\frac{h^{2}}{2\sigma^{2}}}dh$ and variance defined in
Eq. (4). The initial condition $\rho_{0}$ is the double thermofield state that
we discussed earlier. Then, the Frame potential takes the form:

$F^{(k)}(T)=\frac{1}{\mathcal{N}^{2k}}\int d\mu(h)\left|{\rm
Tr}[U_{h}(T)]\right|^{2k},$ (13)

where $\mathcal{N}$ is the dimension of the Hilbert space, that comes from the
definition of the initial state Eq. (8). In order to get to this expression
for the Frame potebtial, we used the property Eq. (9) together with
$U_{h}=U_{h^{L}_{1}}(t)U_{h^{R}_{1}}(t)U_{h^{R}_{2}}^{{\dagger}}(t)U_{h^{L}_{2}}^{{\dagger}}(t)$.
Also, we introduced the measure over the noisy evolution
$d\mu(h)=d\mu(h_{1})d\mu(h_{2})$. The way of calculating this integral is the
following. First, we move to the path integral representation using Keldysh
technique. Then, we do average with respect to Gaussian random variables, and
finally, calculate the path integral using the saddle-point approximation.

In order to use Keldysh path integral formalism, we need to construct coherent
states Shankar:2017zag . To do so, we introduce an additional set of Majorana
operators $\eta_{i}$ and create fermionic operators in the form:

$\hat{c}_{i}=\frac{\hat{\chi}_{i}+i\hat{\eta}_{i}}{\sqrt{2}},~{}~{}~{}~{}\hat{c}^{\dagger}_{i}=\frac{\hat{\chi}_{i}-i\hat{\eta}_{i}}{\sqrt{2}}.$
(14)

With this notation, the coherent state can be defined as

$\left|\psi\right>=\left|0\right>-\psi\left|1\right>=e^{-\psi\hat{c}^{\dagger}\left|0\right>},$
(15)

where $\psi$ is a Grassmann variable that parameterizes our state. The set of
these states is not orthonormal, and their overlap is given by:

$\left<\psi|\psi^{\prime}\right>=\left(\left<0\right|-\left<1\right|\overline{\psi}\right)\left(\left|0\right>-\psi^{\prime}\left|1\right>\right)=1+\overline{\psi}\psi^{\prime}=e^{\overline{\psi}\psi^{\prime}}.$
(16)

Concerning the path integral formalism, we need to get two more useful
identities. First, using the Eq. (16) we can express an identity operator as
an integral over the Grasmann fields:

$\int d\overline{\psi}\int d\psi
e^{-\overline{\psi}\psi}\left|\psi\right>\left<\psi\right|=1,$ (17)

and secondly we can express the trace of an operator $\hat{O}$ in terms of
Grassmann variables as:

${\rm
Tr}(\hat{O})=\sum_{n_{1},n_{2},..=0,1}\left<n_{1}...n_{N}\right|\hat{O}\left|n_{1}...n_{N}\right>=\int
d\overline{\psi}d\psi
e^{-\overline{\psi}\psi}\left<\psi\right|\hat{O}\left|-\psi\right>,$ (18)

here $-1$ comes from commuting left and right coherent states. Using Eq. (17)
and Eq. (18) we can rewrite Eq. (13) as following

$F^{(k)}=\frac{1}{\mathcal{N}^{2}}\int
d\mu(h)d\overline{\psi}d\psi\prod_{l=1}^{k}\left<\psi^{-,l}_{1}\right|U_{1}(\epsilon)\left|\psi^{-,l}_{2}\right>\left<\psi^{-,l}_{2}\right|U_{2}(\epsilon)\left|\psi^{-,l}_{3}\right>...\left<\psi^{-,l}_{2N}|-\psi^{-,l}_{1}\right>\cdot\\\
\cdot\left<\psi^{+,l}_{2N}\right|U^{{\dagger}}_{2N}(\epsilon)\left|\psi^{+,l}_{2N-1}\right>\left<\psi^{+,l}_{2N-1}\right|U^{{\dagger}}_{2N-1}(\epsilon)\left|\psi^{+,l}_{2N-2}\right>...\left<\psi^{+,l}_{1}|-\psi^{+,l}_{2N}\right>e^{-\overline{\psi}^{s,l}_{i}\psi^{s,l}_{i}},$
(19)

where $s=\pm$ indicates one of the two Keldysh contours, $l$ counts replicas,
and $i$ labels time. Each matrix element here has the form,

$\left<\psi^{-,l}_{1}\right|U_{1}(\epsilon)\left|\psi^{-,l}_{2}\right>=\left<\psi^{-,l}_{1}\right|e^{-i^{\frac{q}{2}+1}\epsilon
h_{i_{1}i_{2}..i_{q}}^{1}(\hat{c}_{i_{1}}^{\dagger}+\hat{c}_{i_{1}})...(\hat{c}_{i_{q}}^{\dagger}+\hat{c}_{i_{q}})}\left|\psi^{-,l}_{2}\right>.$
(20)

Let us notice that for any matrix element of this form we can just replace
$\hat{c}^{\dagger}$ with $\overline{\psi}$ and $\hat{c}$ with $\psi$. We do
not need normal ordering here because $H$ Eq. (3) has one term in which the
order is opposite. We can reverse the order at the cost of a minus sign
(remember $i_{k}\neq i_{l}$ everywhere and all fermion operators anticommute),
obtain the action, and once again reverse the order of the corresponding
Grassmann variables $\psi,\overline{\psi}$ at the cost of an extra minus sign.
Therefore, we shall write the usual Keldysh action:

$S(\overline{\psi},\psi)=\sum_{l,s=\pm}s\int_{0}^{T}dt\left(\overline{\psi}_{j}^{s,l}\partial_{t}\psi_{j}^{s,l}+H(\hat{c}^{\dagger}\rightarrow\overline{\psi},\hat{c}\rightarrow\psi)\right).$
(21)

Here the continuous limit was taken $\epsilon\rightarrow 0$,
$N\rightarrow\infty$, and $N\epsilon\rightarrow T$. We can now move back to
Majorana fermions:

$\chi^{s,l}_{i}=\frac{\overline{\psi}^{s,l}_{i}+\psi^{s,l}_{i}}{\sqrt{2}},~{}~{}~{}~{}\eta^{s,l}_{i}=i\frac{\overline{\psi}^{s,l}_{i}-\psi^{s,l}_{i}}{\sqrt{2}},$
(22)

note that the variables we consider now are not operators but Grassmann
variables and $i$ is the index of the Fermionic degrees of freedom. In this
limit, actions for $\chi$ and $\eta$ fields separate Shankar:2017zag , and the
partition function factorizes. So, the part that depends on $\eta$ is a
spectator and can be integrated out Eq. (13) giving a contribution to the
overall normalisation. This can be fixed using the fact that the frame
potential equals one at $T=0$. Therefore, for the Frame potential, we obtain:

$F^{(k)}=\frac{F^{(k)}(T)}{F^{(k)}(T\rightarrow
0)}=\frac{Z(\eta)}{F^{(k)}(T\rightarrow 0)\mathcal{N}^{2}}\int d\mu(h)d\chi
e^{-S(\chi)},$ (23)

where we have

$S(\chi)=\sum_{l,s=\pm}s\int_{0}^{T}dt\left(\frac{1}{2}\chi_{j}^{s,l}\partial_{t}\chi_{j}^{s,l}+i^{\frac{q}{2}+1}\sum_{i_{1}<i_{2}<..<i_{q}}h_{i_{1}...i_{q}}(t)\chi^{s,l}_{i_{1}}...\chi^{s,l}_{i_{q}}\right).$
(24)

This action is a part of the standard path integral formalism for the SYK
model 1806.06840 ; Gu:2019jub ; Lunkin:2020tbq . Here, we are interested in
averaging the Frame potential over the random noisy couplings $h$. To do so,
we calculate the Gaussian integral over each of the variables and use the
property given by equation (4), which defines the variance $\sigma^{2}$:

$\int
dh_{i_{1}...i_{q}}e^{-\frac{1}{2\sigma^{2}}h^{2}_{i_{1}..i_{q}}+i^{\frac{q}{2}+1}h_{i_{1}..i_{q}}M_{i_{1}..i_{q}}}=e^{-\frac{(-1)^{\frac{q}{2}}\sigma^{2}}{2}M^{2}_{i_{1}...i_{q}}},$
(25)

where
$M_{i_{1}...i_{q}}=\sum_{l}(\chi^{-,l}_{i_{1}}...\chi^{-,l}_{i_{q}}-\chi^{+,l}_{i_{1}}...\chi^{+,l}_{i_{q}})$.
Introducing the bilocal fields $G(t_{1},t_{2})$, $\Sigma(t_{1},t_{2})$ and
performing Hubbard-Stratonovich transformation,

$\int
d\hat{G}\prod_{s,s^{\prime}}\delta\left(G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})-\frac{1}{L}\chi^{s,l}_{i}(t_{1})\chi^{s^{\prime},l^{\prime}}_{i}(t_{2})\right)=\int
d\hat{\Sigma}d\hat{G}e^{-\frac{L}{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}\left(G^{ss^{\prime}}_{ll^{\prime}}-\frac{1}{L}\chi^{s,l}_{i}\chi^{s^{\prime},l^{\prime}}_{i}\right)},$
(26)

we obtain the following expression for an effective action:

$S(\chi,G,\Sigma)=-\sum_{s,l}s\int_{0}^{T}\frac{1}{2}\chi_{j}^{s,l}\partial_{t}\chi_{j}^{s,l}+\frac{L}{4q^{2}}\sum_{l,l^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(2G_{ll^{\prime}}^{ss^{\prime}}\right)^{q}-\\\
-\frac{L}{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})\left(G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})-\frac{1}{L}\chi^{s,l}_{i}(t_{1})\chi^{s^{\prime},l^{\prime}}_{i}(t_{2})\right),$
(27)

where $c_{++}=c_{--}=-1$ and $c_{+-}=c_{-+}=1$. Notice that in order to form
bilocal fields in the expression Eq. (25), we performed $\frac{q(q-1)}{2}$
transpositions of the Grassman variables, which gave us the prefactor in the
interaction term $(-1)^{\frac{q(q-1)}{2}}$ together with $(-1)^{\frac{q}{2}}$
from the Gaussian integration, resulting in $(-1)^{\frac{q^{2}}{2}}$, which is
unity due to considering $q$ as an even integer. In order to proceed further,
we integrate out the fields $\chi_{i}$, and the action takes the form:

$\displaystyle S(G,\Sigma)=L\log({\rm
Pf}(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}))+\frac{L}{4q^{2}}\sum_{l,l^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(2G_{ll^{\prime}}^{ss^{\prime}}\right)^{q}-\frac{L}{2}\text{Tr}(\hat{\Sigma}^{T}\hat{G}),$
(28)

Here, the indexes $s$ and $s^{\prime}$ take two different values, $+$ and $-$,
which indicate one of the Keldysh contours. These contours are two loops, as
can be seen from equation (19). Notice that after averaging over the noise,
the fields from the different contours become mixed, and in the action, we
have not only $G^{++}$ and $G^{--}$, but also $G^{+-}$ and $G^{-+}$.
Therefore, we can think of these operators as matrices in the $+-$ space.
Also, let’s take into account that the derivative term
$\frac{1}{2}\chi\partial_{t}\chi$ has a positive sign for the $+$ contour and
a negative sign for the $-$ contour. Therefore, this term can also be
expressed as a $2\times 2$ matrix with Keldysh contour indexes in the
following way:

$-\sum_{s,l}s\int_{0}^{T}\frac{1}{2}\chi_{j}^{s,l}\partial_{t}\chi_{j}^{s,l}=-\frac{1}{2}\chi^{l}_{j}\left(\hat{\sigma}_{z}\partial_{t}\right)\chi^{l}_{j}.$
(29)

Here, we introduced the Pauli matrix $\hat{\sigma}_{z}$ in the $+-$ space, and
we assume a Kronecker delta in the replica space for this term. In the next
section, we are going to calculate this functional integral over the fields
$G$ and $\Sigma$ using the saddle-point approximation at large $L$.

##  3 Majorana SYK model and saddle points

We are in a position to obtain the equations of motion at the saddle points.
In order to do so, we consider the variation of the action with respect to the
bilocal fields:

$\displaystyle\frac{\delta
S}{\delta\hat{\Sigma}}:~{}~{}~{}\hat{G}^{-1}=\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(30) $\displaystyle\frac{\delta
S}{\delta\hat{G}}:~{}~{}~{}\Sigma^{ss^{\prime}}_{ll^{\prime}}=\frac{c_{ss^{\prime}}}{q}(2G^{ss^{\prime}}_{ll^{\prime}})^{q-1}\delta(t_{12}).$
(31)

Setting the variation to zero and summing over all possible solutions, we can
find the time dependence of the Frame potential.

###  3.1 Short-time behaviour

Let us now proceed naively by looking for a saddle point that is diagonal in
the $+-$ space. Then the obvious solution is:

$\hat{\Sigma}=0,~{}~{}~{}~{}~{}\hat{G}^{-1}=\hat{\sigma}_{z}\partial_{t}.$
(32)

With this solution, the Frame potential can be expressed as:

$F^{(k)}=\frac{1}{\mathcal{N}^{2k}}\int
d\mu\left(h_{1}\right)d\mu\left(h_{2}\right)\left({\rm Tr}(U){\rm
Tr}(U^{{\dagger}})\right)^{k}=\\\
=\frac{1}{\mathcal{N}^{2k}}\det(\hat{\sigma}^{z}\partial_{t})^{\frac{L}{2}}e^{-\frac{L}{4q^{2}}\sum_{ll^{\prime}}\int_{0}^{T}dt(2G_{ll^{\prime}}^{++})^{q}+(2G_{ll^{\prime}}^{--})^{q}}.$
(33)

This functional determinant can be calculated using diagonalisation of the
operator $\hat{\sigma}^{z}\partial_{t}$. By establishing antiperiodic boundary
conditions, we obtain eigenvalues $\lambda_{n}=\pm\frac{\pi i(2n+1)}{T}$,
where $n=0,\pm 1,\pm 2,..$. In the $+-$ space, we have two eigenvalues
$\lambda_{n}$ and $-\lambda_{n}$. However, coherent states are over complete,
so we consider $\lambda_{n}=\frac{\pi i}{T},\frac{3\pi
i}{T},\dots,\frac{(4N+1)\pi i}{4T}$, which prevents us from treating matrices

$\begin{pmatrix}-\lambda_{n}&0\\\
0&\lambda_{n}\end{pmatrix}~{}~{}~{}\text{and}~{}~{}~{}\begin{pmatrix}\lambda_{n}&0\\\
0&-\lambda_{n}\end{pmatrix},$ (34)

as different matrices. Notice that this determinant should also be taken in
the replica space, where we have $k$ copies of the system. Taking the limit
$N\rightarrow\infty$ and following the regularization technique from the book
Nakahara:2003nw  (notice that we have $\omega=0$), we get

$\det(\hat{\sigma}_{z}\partial_{t})^{\frac{L}{2}}=2^{\frac{2kL}{2}},$ (35)

where $2^{\frac{L}{2}}$ is the dimension of the Hilbert space for the initial
operator $U$. Notice that the matrix $G$ is also diagonal in the replica
space. Therefore,
$G^{++}_{ll^{\prime}}=\frac{1}{2}\delta_{ll^{\prime}}=-G_{ll^{\prime}}^{--}$,
and we have

$e^{-\frac{L}{4q^{2}}\sum_{ll^{\prime}}\int_{0}^{T}dt(2G_{ll^{\prime}}^{++})^{q}+(2G_{ll^{\prime}}^{--})^{q}}=e^{-\frac{LTk}{2q^{2}}}.$
(36)

The Frame potential on this solution takes the form

$F^{(k)}=\frac{\mathcal{N}^{2k}}{\mathcal{N}^{2k}}e^{-\frac{LTk}{2q^{2}}}=e^{-\frac{LTk}{2q^{2}}}.$
(37)

We can observe that on this solution, the Frame potential exhibits an
exponential decay with respect to time and all the replicas are factorised.
This behaviour implies that the contribution from this part of the Frame
potential is relevant during the early times of the system’s evolution but
becomes increasingly suppressed as time progresses.

###  3.2 Long-time behaviour

In this section, we consider two fundamentally different cases: $q>2$ and
$q=2$. The difference arises from the fact that in the case of $q=2$, the
action is invariant under the continuous symmetry group $SO(k)\times SO(k)$,
so the fields $\Sigma$ and $G$ can be rotated in the replica space by the
generators $O_{1}=e^{i\hat{A}_{1}\hat{\sigma}_{z}}$,
$O_{2}=e^{i\hat{A}_{2}\hat{1}}$ (which form $SO(k)\times SO(k)$ group), here
$\hat{A}_{i}$ is a skew-symmetric matrix of the size $k$ and by
$\hat{A}\hat{\sigma}_{z}$ we suppose the tensor product. Whereas for $q>2$,
the action is invariant under a discrete permutation group $P_{k}$. Therefore,
in the first case, the saddle point physics will be characterised by the
Goldstone theorem, and the integral over the fluctuations around the saddle
point solution will provide us with a non-trivial correction, which we also
aim to study. Another difference between these two examples is that the
solutions for $q=2$ will form a continuous manifold, while for $q>2$, we will
have just a discrete sum.

####  3.2.1 The non-Gaussian case $q>2$

Here, we begin by considering the non-diagonal solution in the $+-$ space. To
do so, let’s examine the following ansatz

$\hat{G}(t_{1},t_{2})=\frac{f(|t_{12}|)}{2}\begin{pmatrix}\text{{\rm
sgn}}(t_{12})&-\tau^{T}\\\ \tau&-\text{{\rm sgn}}(t_{12})\end{pmatrix}.$ (38)

Here, $f(|t_{12}|)$ is an arbitrary function that depends on the module of
difference between time $t_{1}$ and $t_{2}$, which will be found later. The
matrix $\tau$ is a permutation matrix in the replica space. Let’s discuss the
function $\rm sgn(t_{12})$ in the ansatz. In the discretised time one can see
that the action consists of the terms $G_{i,i-1}$, where $i$ labels time.
Because the solution depends on the difference between times $\rm
sgn(t_{i}-t_{i-1})=\rm sgn(\epsilon\times(i-i+1))=\rm sgn(\epsilon)$. So the
$\rm sgn(t_{12})$ is never equal to zero and defined as

$\rm sgn(t_{12})=\begin{cases}1,~{}~{}~{}t_{1}\geq t_{2}\\\
-1,~{}~{}~{}t_{1}<t_{2}\end{cases}.$ (39)

The choice of this ansatz is motivated by the fact that we can permute
replicas in the action without altering its form. Notice that in the case
$q=2$ we have continuous $SO(k)\times SO(k)$ symmetry of the action instead of
the permutation symmetry for general $q$. This case will be considered later
in this work. By using equation (30), we obtain:

$\hat{\Sigma}(t_{1},t_{2})=\frac{1}{q}\delta(t_{12})f(0)^{q-1}\begin{pmatrix}0&-\tau^{T}\\\
\tau&0\end{pmatrix}.$ (40)

Now, our objective is to determine the function $f(t)$. To achieve this, we
perform the Fourier transformation $\hat{\Sigma}(t)=\int
d\omega\hat{\Sigma}(\omega)e^{-i\omega t}$ and substitute this ansatz into the
first equation of motion (30). This procedure gives:

$\hat{G}(\omega)=\frac{1}{\omega^{2}+\left(\frac{f(0)^{q-1}}{q}\right)^{2}}\begin{pmatrix}i\omega&-\frac{f(0)^{q-1}}{q}\tau^{T}\\\
\frac{f(0)^{q-1}}{q}\tau&-i\omega\end{pmatrix}.$ (41)

By performing the inverse Fourier transform and comparing with the ansatz, we
arrive at the following result:

$\displaystyle\hat{G}(t_{1},t_{2})=\frac{e^{-\frac{|t_{12}|}{q}}}{2}\begin{pmatrix}\text{{\rm
sgn}}(t_{12})&-\tau^{T}\\\ \tau&-\text{{\rm
sgn}}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}(t_{1},t_{2})=\frac{1}{q}\delta(t_{12})\begin{pmatrix}0&-\tau^{T}\\\
\tau&0\end{pmatrix}.$ (42)

Here, we again should make a remark on the diagonal elements of this solution.
Usually, for Grassmann variables, ${\rm{\rm sgn}}(0)=\chi(t)\chi(t)=0$.
However, in the discretised time the product of Majorana fermions involves
fermions at slightly different time arguments as we explained in Eq. (39), and
thus conventionally $\chi(t)\chi(t+\epsilon)=\frac{1}{2}$.

Our objective is to understand the behaviour of the action on this solution.
Firstly, let’s note that the second term in the action, Eq. (28), cancels due
to the $c_{ss^{\prime}}$ prefactor and the fact that at zero time ${\rm{\rm
sgn}}(0+\epsilon)=1$

$\frac{L}{4q^{2}}\sum_{l,l^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(2G_{ll^{\prime}}^{ss^{\prime}}\right)^{q}=0.$
(43)

On the other hand the third term will give us a non-trivial contribution

$-\frac{L}{2}\int_{0}^{T}dt\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})=-\frac{L}{2}{\rm
Tr}(\hat{\Sigma}^{T}\hat{G})=-\frac{LTk}{2q}.$ (44)

Finally let’s proceed with the calculation of the first term in the action. In
fact, the Pfaffian can be understood as a normal thermal partition function
with the real inverse temperature equal to $T$:

$e^{S_{1}}=({\rm Pf}[\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}])^{L}={\rm
Tr}(e^{-T\hat{H}})^{L},$ (45)

where

$\hat{H}=-\frac{1}{2}\begin{pmatrix}\hat{\chi}^{+}&i\hat{\chi}^{-}\par\end{pmatrix}\hat{\Sigma}_{0}\begin{pmatrix}\hat{\chi}^{+}\\\
i\hat{\chi}^{-}\end{pmatrix},$ (46)

and

$\hat{\Sigma}_{0}=\frac{1}{q}\begin{pmatrix}0&-\tau^{T}\\\
\tau&0\end{pmatrix}.$ (47)

Let us explain the form of the Hamiltonian. The derivative in the $--$ sector
Eq. (45) comes with an extra minus, so the Majorana fields in the partition
function should come with an extra $i$. Finally, we obtain

$({\rm Pf}[\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}])^{L}={\rm
Tr}(e^{i\frac{T}{q}\hat{\chi}^{-}\tau\hat{\chi}^{+}})^{L}.$ (48)

Similar expression was calculated in ReplicesSymmetry , and here we use the
analogous method. Let us observe that any arbitrary permutation can be
expressed as a direct sum of disjoint cycles, $\tau=\bigoplus_{i}\tau^{c}_{i}$
which implies,

$\tau=\left(\begin{array}[]{c|c|c|c}\tau^{c}_{1}&0&0&0\\\ \hline\cr
0&\tau^{c}_{2}&0&0\\\ \hline\cr 0&0&...&0\\\ \hline\cr
0&0&0&\tau^{c}_{m}\end{array}\right),$ (49)

then ${\rm Tr}(e^{i\frac{T}{q}\hat{\chi}^{-}\tau\hat{\chi}^{+}})={\rm
Tr}(e^{i\frac{T}{q}\sum_{i}\hat{\chi}_{i}^{-}\tau^{c}_{i}\hat{\chi}_{i}^{+}})$
, where the sum goes over the cycles. Also notice that each cycle can be
represented in the canonical form, therefore we have

$\tau^{c}_{i}=\begin{cases}\delta^{\alpha+1,\beta},~{}~{}~{}n_{c}~{}~{}~{}\text{odd}\\\
{\rm{\rm
sgn}}(\beta-\alpha)\delta^{\alpha+1,\beta},~{}~{}~{}n_{c}~{}~{}~{}\text{even},\end{cases}$
(50)

where $n_{c}$ represents the length of a cycle, and
$\delta^{n_{c}+1,\beta}=\delta^{1,\beta}$. Here, ${\rm{\rm
sgn}}(\beta-\alpha)$ merely affects the permutation on the boundaries. The
boundary conditions have been chosen in a way to keep the parity $P$ equal to
$1$. It can be observed that for any $\tau$, this solution does not respect
the parity of the initial state (the initial state is a thermofield double
state between, where for each site, the Fermi parity is $P=1$). Therefore, the
boundary conditions Eq. (50) arise from this fact. Then, for each cycle, we
have:

$-T\hat{H}_{c}=i\frac{T}{q}\left(\hat{\chi}^{-}_{\alpha}\delta^{\alpha+1,\beta}\hat{\chi}^{+}_{\beta}+(-1)^{n_{c}+1}\hat{\chi}^{-}_{n_{c}}\hat{\chi}^{+}_{1}\right).$
(51)

Introducing complex fermions

$c_{\alpha}=\frac{\hat{\chi}^{+}_{\alpha}+i\hat{\chi}^{-}_{\alpha}}{\sqrt{2}},~{}~{}~{}~{}c^{\dagger}_{\alpha}=\frac{\hat{\chi}^{+}_{\alpha}-i\hat{\chi}^{-}_{\alpha}}{\sqrt{2}},$
(52)

and inserting this to the Hamiltonian we obtain

$-T\hat{H}(\tau_{c})=\frac{T}{2q}\sum_{\alpha}(c^{\dagger}_{\alpha}c_{\alpha+1}+c_{\alpha}c_{\alpha+1}+h.c.)+\frac{T}{2q}((-1)^{n_{c}+1}c^{\dagger}_{n}c_{1}+h.c.),$
(53)

this is the Kitaev chain Kitaev2001-da  with periodic/antiperiodic boundary
conditions, Eq. (50) and odd/even length of chain respectively. The
diagonalisation of the chain is readerly done, giving

${\rm Tr}(e^{-T\hat{H}(\tau_{c})})\sim
e^{\frac{n_{c}T}{2q}},~{}~{}~{}T\rightarrow\infty.$ (54)

Now we can finally sum over all the cycles, noticing that the size of the
permutation matrix $k$ should be equal to the sum of the length of cycles of
this permutation $k=\sum_{c}n_{c}$

${\rm Tr}(e^{-T\hat{H}})=\prod_{c}e^{\frac{n_{c}T}{2q}}=e^{\frac{kT}{2q}}.$
(55)

Consequently we find the following expression for the Pfaffian:

$e^{S_{1}}=({\rm
Pf}[\hat{\sigma}_{z}\partial_{t}-\Sigma])^{L}=e^{\frac{LkT}{2q}},$ (56)

and the Frame potential on this solution can be expressed as:

$F^{(k)}(\tau)=\frac{1}{\mathcal{N}^{2k}}e^{S_{1}+S_{2}+S_{3}}=\frac{1}{\mathcal{N}^{2k}}(e^{\frac{LkT}{2q}+0-\frac{LkT}{2q}})=\frac{1}{\mathcal{N}^{2k}}.$
(57)

We shall now perform a summation over all the possible saddle points after we
discuss the fluctuations in the system. Because the fluctuations can give a
correction that depends on the permutation matrix, so the summation may not be
trivial Eq. (96).

####  3.2.2 The Gaussian case $q=2$

In the case $q=2$, the action has a continuous symmetry of $SO(k)\times SO(k)$
as we discussed in the beginning of the section. Therefore, the ansatz that we
examine here consists of orthogonal matrices instead of permutations. Using
the same procedure as in the previous subsection, one can see that for $q=2$,
the solution is:

$\displaystyle\hat{G}(t_{1},t_{2})=\frac{e^{-\frac{|t_{12}|}{q}}}{2}\begin{pmatrix}\text{{\rm
sgn}}(t_{12})&-\theta^{T}\\\ \theta&-\text{{\rm
sgn}}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}(t_{1},t_{2})=\frac{1}{q}\delta(t_{12})\begin{pmatrix}0&-\theta^{T}\\\
\theta&0\end{pmatrix},$ (58)

where $\theta$ is an orthogonal matrix in the replica space. Again, we are
interested in the calculation of the action on this solution.Notice that the
second and third terms in the action are the same as in the case of
permutation matrices. Therefore, our objective now is to calculate the first
term, namely the Pfaffian.

Let us now apply the same procedure to the orthogonal matrices. Any arbitrary
orthogonal matrix can be expressed in block diagonal form:

$\displaystyle\theta_{2m}=\left(\begin{array}[]{c|c|c|c|c}R_{1}&0&0&0&0\\\
\hline\cr 0&R_{2}&0&0&0\\\ \hline\cr 0&0&...&0&0\\\ \hline\cr
0&0&0&R_{m-1}&0\\\ \hline\cr
0&0&0&0&R_{m}\end{array}\right),~{}~{}~{}\text{and}~{}~{}~{}\theta_{2m+1}=\left(\begin{array}[]{c|c|c|c|c}R_{1}&0&0&0&0\\\
\hline\cr 0&R_{2}&0&0&0\\\ \hline\cr 0&0&...&0&0\\\ \hline\cr 0&0&0&R_{m}&0\\\
\hline\cr 0&0&0&0&\pm 1\\\ \end{array}\right),$ (69)

where, $m$ is a positive integer, and $R_{i}$ is an orthogonal matrix in two-
dimensional space. Then, ${\rm
Tr}(e^{i\frac{T}{q}\hat{\chi}^{-}\theta\hat{\chi}^{+}})={\rm
Tr}(e^{i\frac{T}{q}\sum_{i}\hat{\chi}_{i}^{-}R_{i}\hat{\chi}_{i}^{+}})$. Here,
we should also take into account that not all orthogonal matrices will
preserve the parity of the initial state. Therefore, $R_{i}$ should be
rotation matrices. Also, in the case of an odd size, the last element on the
diagonal is positive unity because we want to preserve parity. Now, let’s
consider the Hamiltonian for the two-dimensional orthogonal submatrix:

$-T\hat{H}_{R}=i\frac{T}{q}\left(\cos(\phi)\hat{\chi}^{-}_{1}\hat{\chi}^{+}_{1}+\sin(\phi)\hat{\chi}^{-}_{1}\hat{\chi}^{+}_{2}-\sin(\phi)\hat{\chi}^{-}_{2}\hat{\chi}^{+}_{1}+\cos(\phi)\hat{\chi}^{-}_{2}\hat{\chi}^{+}_{2}\right).$
(70)

Here $\phi$ is an angle which parametrizes the rotational matrix, and we can
again move to the complex fermions Eq. (52), which gives a simple two-fermions
Hamiltonian:

$-T\hat{H}_{R}=\frac{T}{2q}(2\cos(\phi)(c_{1}^{\dagger}c_{1}+c_{2}^{\dagger}c_{2}-1)+2\sin(\phi)(c_{1}c_{2}-c_{1}^{\dagger}c_{2}^{\dagger}))=\\\
=\frac{T}{q}\left(c_{1}^{\dagger}~{}~{}c_{2}\right)\begin{pmatrix}\cos{\phi}&-\sin{\phi}\\\
-\sin{\phi}&-\cos{\phi}\end{pmatrix}\begin{pmatrix}c_{1}\\\
c^{\dagger}_{2}\end{pmatrix},$ (71)

after the diagonalisation the trace can be easily calculated. Using the well
known formula $\rm
Tr(e^{c^{\dagger}_{i}\Gamma_{ij}c_{j}})=\det(1+e^{\Gamma_{ij}})$ we obtain:

${\rm
Tr}(e^{-T\hat{H}_{R}})=\det\left(1+\exp\left(\frac{T}{q}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\right)\right)=2+e^{\frac{T}{q}}+e^{-\frac{T}{q}}\rightarrow
e^{\frac{T}{q}},~{}~{}T\rightarrow\infty.$ (72)

Therefore the Frame potential on the particular solution which depends on an
orthogonal matrix $\theta$ is

$F^{(k)}(\theta)=\frac{1}{\mathcal{N}^{2k}}e^{\frac{LkT}{2q}+0-\frac{LkT}{2q}}=\frac{1}{\mathcal{N}^{2k}}.$
(73)

] As in the previous subsection, the Frame potential for the particular
rotational matrix does not depend on this matrix. Therefore, we should first
consider the fluctuations around the saddle point, and after sum over all
possible solutions. As we shall see in the coming section, the fluctuations
are responsible for difference between $q>2$ and $q=2$ cases. In the $q=2$
example they bring the finite $L$ corrections $\frac{\log(L)}{L}$, while in
the case $q>2$ fluctuations are not relevant and leave the next correction of
the order $\frac{1}{L}$.

###  3.3 Symmetries and Goldstone modes

In the next subsection we shall calculate fluctuations around the saddle point
which we discussed earlier. For the case $q=2$, we need to understand how many
modes will be massive and carry non-trivial corrections to the Frame
potential. To do so, we can use the Goldstone theorem. Let $B$ be a continuous
group of a global symmetry of the action and $H$ is a subgroup of $B$ which
leaves the solution of the equations of motion unchanged. Then the number of
massless modes is equal to $R_{ml}={\rm dim}(B)-{\rm dim}(H)$. Therefore, the
number of massive modes is $R_{m}=d-R_{ml}$ where $d$ is the number of degrees
of freedom of our system. Let us start with identification of the degrees of
freedom. By the definition:

$G_{ll^{\prime}}=\begin{pmatrix}G^{++}_{ll^{\prime}}&G^{+-}_{ll^{\prime}}\\\
G^{-+}_{ll^{\prime}}&G^{--}_{ll^{\prime}}\end{pmatrix},$ (74)

and
$G^{++}_{ll^{\prime}}=\frac{1}{L}\sum^{L}_{i=1}(\chi^{+})^{i}_{l}(\chi^{+})^{i}_{l^{\prime}}$
,where $L$ is the number of fermionic modes, first we notice that the matrix
$G_{ll^{\prime}}$ is a skew-symmetric matrix because of the anticommutation
relations on fields $\chi$, which contains $2k(2k-1)/2$ degrees of freedom.
Let us check the commutation relations for this matrix

$[G_{ll^{\prime}},G_{mm^{\prime}}]=\frac{1}{L^{2}}\sum_{i}\sum_{j}[\chi^{i}_{l}\chi^{i}_{l^{\prime}},\chi^{j}_{m}\chi^{j}_{m^{\prime}}],$
(75)

expanding the commutator and using anticommutation relations
$\\{\chi^{i}_{l},\chi^{j}_{m}\\}=\frac{1}{2}\delta^{ij}\delta^{lm}$ we have,

$[G_{ll^{\prime}},G_{mm^{\prime}}]=\frac{1}{2L}(\delta_{l^{\prime}m}G_{lm^{\prime}}-\delta_{lm^{\prime}}G_{ml^{\prime}}+\delta_{lm}G_{m^{\prime}l^{\prime}}-\delta_{l^{\prime}m^{\prime}}G_{lm}),$
(76)

which gives us the commutation relations of rotation generators in $2k$
dimensions. So we see that by the definition our matrix $G_{ll^{\prime}}$ is a
generator of $so(2k)$ algebra. In Fava:2023tgg , it was shown that
$S_{ll^{\prime}}=iG_{ll^{\prime}}$ satisfy orthogonality relation in the large
$L$ limit:

$SS^{T}=\mathbb{I},$ (77)

for our matrix $G$ it also means

$GG^{T}=\mathbb{I}.$ (78)

Therefore this matrix should be not only skew-symmetric but also orthogonal.
Space of these matrices isomorphic to the $SO(2k)/U(k)$ baker , which gives us
$d=\frac{2k(2k-1)}{2}-k^{2}=k(k-1)$ degrees of freedom.

Now let’s consider the group of the global symmetry of the action in more
details. Remind that at finite times, the action looks like Eq. (28). So after
the integration over the fields $G$ we have

$S=L\log({\rm
Pf}(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}))-\frac{L}{4}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(\Sigma^{ss^{\prime}}_{ll^{\prime}}\right)^{2}.$
(79)

One can see that the problematic term here is $L\log({\rm
Pf}(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}))$, which as we already noticed
has a rotation symmetry in replica space
$O_{1}=e^{i\hat{A}_{1}\hat{\sigma}_{z}}$, $O_{2}=e^{i\hat{A}_{2}\hat{1}}$ with
$\hat{A}_{i}$ skew-symmetric matrix of the size $k$ and by
$\hat{A}\hat{\sigma}_{z}$ we suppose the tensor product. Notice that matrices
$O_{0}=e^{i\hat{A}}$ form a $SO(k)$ symmetry group, therefore the space of the
symmetries is $SO(k)\times SO(k)$. So

$\det(\hat{\sigma}_{z}\partial-
O_{i}\hat{\Sigma}O_{i}^{T})=\det(O_{i}O_{i}^{T})\det(O_{i}^{-1}\hat{\sigma_{z}}\partial_{t}(O_{i}^{T})^{-1}-\hat{\Sigma}),$
(80)

since matrices $O_{i}$ commute with $\hat{\sigma}_{z}$, we have

$\det(\hat{\sigma}_{z}\partial-
O_{i}\hat{\Sigma}O_{i}^{T})=\det(\hat{\sigma}_{z}\partial-\hat{\Sigma}).$ (81)

Other possible transformations are not symmetries of the action, so the group
is $B=SO(k)\times SO(k)$ and ${\rm dim}(B)=k(k-1)$. Remind that the saddle
point solution has the form

$\hat{G}(t_{1},t_{2})=\frac{e^{-\frac{|t_{12}|}{2}}}{2}\begin{pmatrix}{\rm
sgn}(t_{12})&-\theta^{T}\\\ \theta&-{\rm
sgn}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}=\frac{1}{2}\delta(t_{12})f(0)\begin{pmatrix}0&-\theta^{T}\\\
\theta&0\end{pmatrix},$ (82)

where $\theta$ is an orthogonal matrix $k\times k$. From the form of the
solution, we see that the group that allows us to transform one solution into
another one is $B/H=SO(k)$, which means that $H=SO(k)$ too. So $R_{ml}={\rm
dim}(B)-{\rm dim}(H)=k(k-1)/2$. Therefore, we find

$R_{m}=d-R_{ml}=k(k-1)-k(k-1)/2=k(k-1)/2,$ (83)

giving the number of massive modes, which will allow us to calculate the
integral over the relevant fluctuations.

###  3.4 Fluctuations around the saddle points

Let us consider fluctuations in our system. In order to do so, we expand the
effective action around the saddle point solution
$S=S_{0}(G_{0},\Sigma_{0})+S^{(2)}(\delta G,\delta\Sigma)$ and integrate
fluctuations out, which suppose to give us the next order correction in the
fermionic degrees of freedom. We remind that the object of our interest is the
frame potential which we normalise to be equal to unity at zero time:

$F^{(k)}=\frac{F^{(k)}(T)}{F^{(k)}(0)}=\frac{1}{F^{(k)}(0)\mathcal{N}^{2k}}\int
DGD\Sigma e^{\left[L\log({\rm
Pf}(\hat{\sigma}_{z}\partial-\hat{\Sigma}))+\frac{L}{4}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(G^{ss^{\prime}}_{ll^{\prime}}\right)^{2}\right]}\cdot\\\
\cdot e^{-\frac{L}{2}\int
dt_{1}dt_{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})}.$
(84)

And the saddle point approximation gives us the following expression for the
Frame potential

$F^{(k)}(T)=\frac{1}{\mathcal{N}^{2k}}e^{S_{0}}\left<\int D\delta
GD\delta\Sigma e^{\frac{1}{2}S^{(2)}(\delta G,\delta\Sigma)}\right>,$ (85)

where the on-shell action $S_{0}$ was found in previous chapters. Now let’s
consider the quadratic fluctuation term $S^{(2)}(\delta G,\delta\Sigma)$:

$S^{(2)}(\delta
G,\delta\Sigma)=\frac{\delta^{2}S}{\delta\hat{G}\delta\hat{G}}\delta\hat{G}^{2}+\frac{\delta^{2}S}{\delta\hat{\Sigma}\delta\hat{\Sigma}}\delta\hat{\Sigma}^{2}+2\frac{\delta^{2}S}{\delta\hat{G}\delta\hat{\Sigma}}\delta\hat{\Sigma}\delta\hat{G},$
(86)

doing variation we can find

$\displaystyle\frac{\delta^{2}S}{\delta\hat{G}\delta\hat{G}}\delta\hat{G}^{2}=\frac{L}{2}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dtc^{ss^{\prime}}(\delta
G^{ss^{\prime}}_{ll^{\prime}})^{2},$ (87)
$\displaystyle\frac{\delta^{2}S}{\delta\hat{G}\delta\hat{\Sigma}}\delta\hat{\Sigma}\delta\hat{G}=-L\int
dt_{1}dt_{2}\delta\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})\delta
G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2}).$ (88)

Let’s consider the second term in Eq. (86), to do so we can express:

$L\log({\rm
Pf}(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}-\delta\hat{\Sigma}))=\frac{L}{2}{\rm
Tr}(\log(1-(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma})^{-1}\delta\hat{\Sigma}))+\frac{L}{2}Tr\log(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma})=\\\
=\frac{L}{2}Tr\log(\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma})-\frac{L}{2}Tr(\hat{G}\delta\hat{\Sigma})-\frac{L}{4}Tr(\hat{G}\delta\hat{\Sigma}\hat{G}\delta\hat{\Sigma}).$
(89)

The first term here belongs to $S_{0}$, the second term vanishes because of
the equations of motion, and the third term is the quadratic fluctuation which
we are looking for. Therefore,

$\frac{\delta^{2}S}{\delta\hat{\Sigma}\delta\hat{\Sigma}}\delta\hat{\Sigma}^{2}=-\frac{L}{4}\int
dt_{1}dt_{2}G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})\delta\Sigma^{s^{\prime}\alpha}_{l^{\prime}m}(t_{1},t_{2})G^{\alpha\alpha^{\prime}}_{mm^{\prime}}(t_{1},t_{2})\delta\Sigma^{\alpha^{\prime}s}_{m^{\prime}l}(t_{1},t_{2}),$
(90)

where the operator $G$ is a solution of the saddle point equation. Integrating
out the fluctuations, we obtain:

$\int D\delta GD\delta\Sigma e^{\frac{1}{2}S^{(2)}(\delta
G,\delta\Sigma)}\sim\left(\frac{4}{L}\right)^{k(k-1)/2}\left(\frac{1}{Det(Q(\theta))}\right)^{k(k-1)}.$
(91)

Here $Q(\theta)$ is a matrix that depends on the orthogonal matrices from the
saddle point solution and does not depend on the size of the system. We also
need to calculate the normalisation term, which is a standard Gaussian
integral

$F^{(k)}(0)=\int DGD\Sigma
e^{\left[\frac{L}{4}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dtc_{ss^{\prime}}\left(G^{ss^{\prime}}_{ll^{\prime}}\right)^{2}-\frac{L}{2}\int
dt_{1}dt_{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})\right]}=\\\
=\left(\frac{4}{L}\right)^{k(k-1)/2}\int D\delta\Sigma
e^{\frac{L}{4}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dtc^{ss^{\prime}}(\Sigma^{ss^{\prime}}_{ll^{\prime}})^{2}(t)}=\left(\frac{4}{L}\right)^{k(k-1)}.$
(92)

Putting it all together, we obtain for $q=2$:

$F^{(k)}_{q=2}=\frac{e^{S_{0}}}{\mathcal{N}^{2k}}\left(\frac{L}{4}\right)^{k(k-1)/2}\left<\frac{1}{Det(Q(\theta))^{k(k-1)}}\right>_{\theta},$
(93)

where $\left<...\right>_{\theta}$ means the summation over all saddle points.
Notice that here we are interested in the large $L$ expansion the expression
and term $\left<\frac{1}{Det(Q(\theta))^{k(k-1)}}\right>_{\theta}$ does not
depend on $L$. Therefore, taking the logarithm of the Frame potential we can
neglect this part:

$\lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left<F^{k}\right>_{q=2}\right)\sim-k\log(2)+\frac{k(k-1)}{2}\frac{1}{L}\log(L),$
(94)

which coincides with the result for fermionic Haar calculation, see Appendix
(133). Notice that for the generic case of $q$ interaction the symmetry of the
action and the solution is discrete, so there are no Goldstone modes in this
case. Due to this fact the integration over the fluctuation must be carried
over all the degrees of freedom which are ${\rm dim}(SO(2k)/U(k))$, also, we
have the same degrees of freedom at zero and at finite time. Therefore, the
fluctuation part from the finite time action will be cancelled by the
normalisation. Hence, for the generic $q>2$ case we have

$F^{(k)}_{q>2}=\frac{e^{S_{0}}}{\mathcal{N}^{2k}}\left<\frac{1}{Det(Q(\tau))^{k(k-1)}}\right>_{\tau},$
(95)

together with,

$\lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left<F^{k}\right>_{q>2}\right)\sim-k\log(2).$
(96)

One should notice that this observation also coincide with the calculation for
the Haar measure Eq. (125).

##  4 Dirac SYK model and Kelsdysh saddle points

In this part, we consider the Brownian evolution given by the Hamiltonian

$H(t)=\sum_{\begin{subarray}{c}1\leq i_{1}\leq i_{2}...\leq i_{q/2}\leq L\\\
1\leq i_{q/2+1}\leq...\leq i_{q}\leq
L\end{subarray}}h_{i_{1}...i_{q}}\hat{c}^{\dagger}_{i_{1}}..\hat{c}^{\dagger}_{i_{q/2}}\hat{c}_{i_{q/2+1}}...\hat{c}_{i_{q}}.$
(97)

Here, $\hat{c}^{\dagger}$ and $\hat{c}$ are Dirac Fermions, and standard anti-
commutation relations are applied
$\\{\hat{c}_{k},\hat{c}^{\dagger}_{l}\\}=\delta_{kl}$. As in the previous
chapter, $h_{i_{1}...i_{q}}$ are normally distributed random variables, but
this time, they are complex and have the variance:

$\left<h_{i_{1}...i_{q}}(t_{1})h^{*}_{i_{1}...i_{q}}(t_{2})\right>=\frac{2^{q}(q/2)!^{2}}{q^{2}L^{q-1}}\delta(t_{1}-t_{2}).$
(98)

In the next sections we are going to repeat the calculation for Dirac fermions
as we have done in the previous section for Majorana fermions.

###  4.1 Keldysh path integral and saddle points

We again choose the TFD state as initial state, and the object of our interest
has the form given in the Eq. (13), which can be expressed as two Keldysh
contours Eq. (19). The evolution operator is given by
$U_{i}(\epsilon)=e^{-i\epsilon H(t_{i})}$, where the Hamiltonian is from the
Eq. (97). Now we want to derive an effective action of the theory. We do so in
the similar manner as in the Majorana case. Firstly we average over the noise
and secondly perform the Hubbard-Stratonovich transformation:

$\int
d\hat{G}\prod_{s,s^{\prime}}\delta\left(G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})-\frac{1}{L}\overline{\psi}^{s,l}_{i}(t_{1})\psi^{s^{\prime},l^{\prime}}_{i}(t_{2})\right)=\int
d\hat{G}d\hat{\Sigma}e^{-\Sigma^{ss^{\prime}}_{ll^{\prime}}\left(G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})-\frac{1}{L}\overline{\psi}^{s,l}_{i}(t_{1})\psi^{s^{\prime},l^{\prime}}_{i}(t_{2})\right)}.$
(99)

After we integrate over Grassmann fields we obtain the action:

$S=L\log(\det(\hat{\sigma}_{z}\partial-\hat{\Sigma}))+\frac{L}{2q^{2}}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int
dt_{1}dt_{2}\delta(t_{12})c_{ss^{\prime}}(2G^{ss^{\prime}}_{ll^{\prime}})^{q/2}(-2G^{s^{\prime}s}_{l^{\prime}l})^{q/2}\\\
-L\int
dt_{1}dt_{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2}).$
(100)

Here we are again interested in the equations of motions. Repeating the
procedure from the section with Majorana fermions, one can find:

$\displaystyle-(\hat{G}^{-1})^{T}=\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma},$
(101)
$\displaystyle\Sigma^{ss^{\prime}}_{ll^{\prime}}=\frac{c_{ss^{\prime}}}{q}(2G^{ss^{\prime}}_{ll^{\prime}})^{\frac{q}{2}-1}(-2G^{s^{\prime}s}_{l^{\prime}l})^{\frac{q}{2}}\delta(t_{12}).$
(102)

Notice that for Majorana fermions, by definition, matrices $\hat{G}$ and
$\hat{\Sigma}$ were real-valued. Now, they are allowed to be complex. The
early-time behaviour can also be calculated in the similar way so the
corresponding Frame potential is

$F^{(k)}(T\rightarrow 0)=e^{-\frac{LTk}{q^{2}}}.$ (103)

In the case $q>2$, we once again have permutation symmetry in the replica
space. Therefore, we choose an ansatz that consists of permutation matrices,
as it is given in the Eq. (38). Therefore, we arrive at the solution:

$\hat{G}(t_{1},t_{2})=\frac{e^{-\frac{|t_{12}|}{q}}}{2}\begin{pmatrix}\text{{\rm
sgn}}(t_{12})&-\tau^{T}\\\ \tau&-\text{{\rm
sgn}}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}(t_{1},t_{2})=\frac{1}{q}\delta(t_{12})\begin{pmatrix}0&-\tau^{T}\\\
\tau&0\end{pmatrix}.$ (104)

Let us calculate the on-shell action using this solution. We start with the
first term in the action:

$e^{S_{1}}=(\det[\hat{\sigma}_{z}\partial_{t}-\hat{\Sigma}])^{L}={\rm
Tr}(e^{-T\hat{H}})^{L}={\rm
Tr}(e^{iT\hat{c}^{\dagger}\hat{\Sigma}\hat{c}})^{L}=\det(1+e^{iT\hat{\Sigma}})^{L},$
(105)

notice that we can expand this exponent and use the fact that
$\hat{\Sigma}^{2}=-\frac{1}{q^{2}}\mathcal{I}$ as follows

$e^{iT\hat{\Sigma}}=1+iT\hat{\Sigma}-\frac{T^{2}(\hat{\Sigma})^{2}}{2}-i\frac{T^{3}(\hat{\Sigma})^{3}}{3!}+\frac{T^{4}(\hat{\Sigma})^{4}}{4!}+...=\\\
=\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k}\frac{1}{(2k)!}+i\hat{\Sigma}_{0}\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k+1}\frac{1}{(2k+1)!},$
(106)

where we have introduced the matrix

$\hat{\Sigma}_{0}=\begin{pmatrix}0&-\tau^{\dagger}\\\ \tau&0\end{pmatrix}.$
(107)

Inserting this result to the Eq. (105) we obtain

$e^{S_{1}}=\det\begin{pmatrix}1+\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k}\frac{1}{(2k)!}&-i\tau^{\dagger}\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k+1}\frac{1}{(2k+1)!}\\\
i\tau\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k+1}\frac{1}{(2k+1)!}&1+\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k}\frac{1}{(2k)!}\end{pmatrix}^{L}=\det\begin{pmatrix}A&B\\\
C&D\end{pmatrix}^{L}.$ (108)

The determinant of this block matrix can be calculated if we notice that
$[C,D]=0$ , because $D$ is a diagonal matrix in the replica space.Then
$\det\begin{pmatrix}A&B\\\ C&D\end{pmatrix}=\det(AD-BC)$ that gives:

$e^{S_{1}}=\det\left(\left(1+\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k}\frac{1}{(2k)!}\right)^{2}-\left(\sum\limits_{k=0}^{\infty}\left(\frac{T}{q}\right)^{2k+1}\frac{1}{(2k+1)!}\right)^{2}\right)^{L}=\\\
=\det\left(\left(1+e^{-\frac{T}{q}}\right)\left(1+e^{\frac{T}{q}}\right)\right)^{L}\rightarrow
e^{\frac{kTL}{q}},~{}~{}~{}T\rightarrow\infty.$ (109)

Now let’s move to the second term in the action:

$S_{2}=\frac{L}{2q^{2}}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int
dt_{1}dt_{2}\delta(t_{12})c_{ss^{\prime}}(2G^{ss^{\prime}}_{ll^{\prime}})^{q/2}(-2G^{s^{\prime}s}_{l^{\prime}l})^{q/2}=0,$
(110)

it is equal to zero on the solution because of the forward and backward
evolution, which comes with the different signs. Notice that at zero time the
expression under the integral is also equal to zero because the solution is
again defined with the normalised $\rm sgn(t_{12})$ Eq. (39). The third term
in the action gives,

$S_{3}=-L{\rm Tr}(\Sigma^{T}G)=-\frac{LkT}{q}.$ (111)

And the Frame potential on this solution takes the form:

$F^{(k)}(\tau)=\frac{1}{\mathcal{N}^{2k}}e^{S_{1}+S_{2}+S_{3}}=\frac{1}{\mathcal{N}^{2k}}(e^{\frac{LkT}{q}+0-\frac{LkT}{q}})=\frac{1}{\mathcal{N}^{2k}},$
(112)

where $\mathcal{N}=2^{L}$ is the dimension of the Hilbert space. Again, we
will do a summation over all the possible solutions after we consider the
fluctuations. Now let’s consider the case $q=2$. Notice that here the action
is symmetric under unitary transformations in the replica space. Therefore, we
consider an ansatz that consists of unitary matrices:

$\hat{G}(t_{1},t_{2})=\frac{f(|t_{12}|)}{2}\begin{pmatrix}{\rm
sgn}(t_{12})&-\tilde{u}^{\dagger}\\\ \tilde{u}&-{\rm
sgn}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}=\frac{1}{2}\delta(t_{12})f(0)\begin{pmatrix}0&-\tilde{u}^{\dagger}\\\
\tilde{u}&0\end{pmatrix},$ (113)

where $\tilde{u}$ is unitary matrix in $k\times k$ space. Inserting this
matrix to the Swinger-Dyson equations, we find:

$\hat{G}(t_{1},t_{2})=\frac{e^{-\frac{|t_{12}|}{2}}}{2}\begin{pmatrix}{\rm
sgn}(t_{12})&-\tilde{u}^{\dagger}\\\ \tilde{u}&-{\rm
sgn}(t_{12})\end{pmatrix},~{}~{}~{}\hat{\Sigma}=\frac{1}{2}\delta(t_{12})\begin{pmatrix}0&-\tilde{u}^{\dagger}\\\
\tilde{u}&0\end{pmatrix}.$ (114)

The action can be calculated exactly in the same way as for permutation
matrices. Hence, we have:

$F^{(k)}(\tilde{u})=\frac{1}{\mathcal{N}^{2k}}.$ (115)

As in the case of Majorana fermions, the difference between the two cases
$q>2$ and $q=2$ will come from the consideration of the fluctuations around
the saddle point solution and the presence of the Goldstone modes.

###  4.2 Degrees of freedom, symmetries and fluctuations

Let us start with the consideration of the Goldstone theorem for the case
$q=2$. The vector fields $\psi$ and $\overline{\psi}$ have $2k$ components, so
the group acting in this space is $U(2k)$, with $d=4k^{2}$ generators. This is
the number of degrees of freedom for our system. Now let’s consider the group
of the global symmetry of the action. Remember that at finite times, the
action looks like:

$S=L\log(\det(\hat{\sigma}_{z}\partial-\hat{\Sigma}))+\frac{L}{2}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dt(G^{++}_{ll^{\prime}}G^{++}_{l^{\prime}l}+G^{--}_{ll^{\prime}}G^{--}_{l^{\prime}l}-2G^{+-}_{ll^{\prime}}G^{-+}_{l^{\prime}l})\\\
-L\int
dt_{1}dt_{2}\Sigma^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2})G^{ss^{\prime}}_{ll^{\prime}}(t_{1},t_{2}).$
(116)

Again, the symmetries of the action are restricted by the term
$L\log(\det(\hat{\sigma}_{z}\partial-\hat{\Sigma}))$, which has rotation
symmetries in the replica space $U_{1}=e^{i\hat{A}_{1}\hat{\sigma}_{z}}$,
$U_{2}=e^{i\hat{A}_{2}\hat{1}}$, where $\hat{A}_{i}$ is a Hermitian matrix of
size $k$. By $\hat{A}\hat{\sigma}_{z}$, we again mean the tensor product.
Notice that matrices $U_{0}=e^{i\hat{A}}$ form the $U(k)$ group. Therefore,
the group of symmetries of the action is $B=U(k)\times U(k)$.

Remember that the saddle point solution has the form Eq. (114) where
$\tilde{u}$ is a unitary matrix. From the form of the solution, we see that
the group that allows us to transform one solution to another is $B/H=U(k)$,
which means that $H=U(k)$ as well. So the number of massless modes is
$R_{ml}=\text{dim}(B)-\text{dim}(H)=k^{2}$. Then we have for the number of
massive modes:

$R_{m}=d-R_{ml}=4k^{2}-k^{2}=3k^{2}.$ (117)

Let us consider fluctuations in our system. To do so, we again expand the
effective action around the saddle point solution
$S=S_{0}(G_{0},\Sigma_{0})+S^{(2)}(\delta G,\delta\Sigma)$ and integrate out
the fluctuations, which are expected to give us the next order correction in
the fermionic degrees of freedom. Remember that the object of our interest is
the frame potential, which we normalize to be equal to unity at zero time:

$\left<F^{k}\right>=\frac{F^{(k)}(T)}{F^{(k)}(0)}=\frac{1}{F^{(k)}(0)\mathcal{N}^{2k}}\int
DGD\Sigma
e^{\left[L\log(\det(\hat{\sigma}_{z}\partial-\hat{\Sigma}))\right]}\cdot\\\
\cdot
e^{\frac{L}{2}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dt(G^{++}_{ll^{\prime}}G^{++}_{l^{\prime}l}+G^{--}_{ll^{\prime}}G^{--}_{l^{\prime}l}-2G^{+-}_{ll^{\prime}}G^{-+}_{l^{\prime}l})-LTr(\hat{\Sigma}^{T}\hat{G})},$
(118)

the saddle point approximation again gives us the expression Eq. (85). And the
integration over fluctuation in the highest order in $L$ gives us

$F^{(k)}\sim\frac{1}{\mathcal{N}^{2k}}e^{S_{0}}\left(\frac{1}{\sqrt{L}}\right)^{2\cdot
3k^{2}}.$ (119)

We also need to calculate the normalisation term, which is just Gaussian
integral in fields $G$ and $\Sigma$ the number of degrees of freedom was found
before $d=4k^{2}$

$F^{(k)}(0)=\int DGD\Sigma
e^{\left[\frac{L}{2}\sum\limits_{l,l^{\prime}s,s^{\prime}}\int_{0}^{T}dt(G^{++}_{ll^{\prime}}G^{++}_{l^{\prime}l}+G^{--}_{ll^{\prime}}G^{--}_{l^{\prime}l}-2G^{+-}_{ll^{\prime}}G^{-+}_{l^{\prime}l})-LTr(\hat{\Sigma}^{T}\hat{G})\right]}\sim\left(\frac{1}{\sqrt{L}}\right)^{2\cdot
4k^{2}}.$ (120)

Therefore we get

$\left<F^{k}\right>=\frac{F^{(k)}(T)}{F^{(k)}(0)}\sim\frac{1}{\mathcal{N}^{2k}}e^{S_{0}}L^{k^{2}},$
(121)

and

$\lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left(\left<F^{k}\right>\right)\right)\sim-2k\log(2)+k^{2}\frac{1}{L}\log(L)+O\left(\frac{1}{L}\right),$
(122)

which coincides with the result for fermionic Haar calculation, Eq. 141. As in
the case of Majorana fermions, for the generic case of $q$ interactions, there
are no Goldstone modes. Due to this fact, the integration over the
fluctuations extends over all the degrees of freedom, which are
$\text{dim}(U(2k))$. This coincides with the number of degrees of freedom at
both zero and finite time. Therefore, the fluctuation part from the finite-
time action will be cancelled by the normalization. Hence, for the generic
case of $q>2$, we have:

$\lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left<F^{k}\right>_{q>2}\right)\sim-2k\log(2)+O\left(\frac{1}{L}\right).$
(123)

Notice that this observation also coincide with the calculation for the Haar
measure Eq. (137).

##  5 Conclusion and discussion

In conclusion, this paper explores the degree of mixing in the time evolution
of the Brownian SYK model, which serves as a paradigmatic model for quantum
scrambling and chaos. The investigation focuses on the Frame potential, a
quantifier of mixing behaviour. By considering two identical copies of an
initial density matrix and calculating the overlap between their independent
time evolutions, the Frame potential, denoted as $F^{(k)}(T)$, is evaluated.
From a saddle point analysis of the Keldysh action, we find that the Frame
potential exhibits exponential decay at early time and saturates to an
exponentially small, with the size of the system value at a timescale only
dependent on the number $q$ of interacting fermions. Our analysis reveals that
in the case of non-Gaussian SYK models with $q>2$, the Frame potential
converges to its Haar-averaged value, denoted as $F_{\rm Haar}$, at late
times. Such a late time behaviour manifests that the systems fully explores
the available Hilbert space in a uniform manner. However, a difference arises
in the Gaussian integrable case with two-fermion interactions. In this
scenario, the evolution exhibits mixing behaviour solely within the space of
Gaussian states. Here, a Gaussian Haar (gHaar) measure is introduced, and it
is demonstrated that the Frame potential converges to the gHaar measure at
late times. Remarkably, the value of the gHaar-integrated Frame potential is
found to be equivalent to the generic case in the limit of a large fermion
number $L$, albeit with logarithmic corrections in $L$. The presence of these
corrections arises from Goldstone fluctuations around the late-time saddle
point. Notably, these fluctuations are only present in the Gaussian case,
where both the Keldysh action and the saddle point have a continuous group
symmetry. Our work therefore extends recent findings on the Frame potential in
SYK-related models jian2022linear .

The study of scrambling dynamics in the context of continuous monitoring has
received significant attention in recent years, with a particular focus on
measurement-induced phase transitions in fermionic systems 2108.11973 ;
2302.12820 ; PhysRevX.9.031009 ; Antonini2022 . Exploring for example the
effects of continuous monitoring on the Frame potential in the SYK model
represents a promising avenue for extending the current work. Additionally,
investigating the behaviour of the Frame potential in non-Brownian dynamics
for the SYK model presents another compelling extension.

## Acknowledgement

We thank Lorenzo Piroli, Guido Giachetti and Adam Nahum for useful
discussions. This work has been partially funded by the ERC Starting Grant
101042293 (HEPIQ) (J.D.N. and A.T.) and by the ANR JCJC grant ANR-21-CE47-0003
(TamEnt) (A.D.L.).

##  Appendix A Haar-averaged Frame potential

As we discussed during our work, it is expected for the Frame potential to
converge to its Haar value during Brownian evolution. In this appendix, we
calculate the Haar values for the Frame potential. First, we consider the case
when the evolution matrix $U$ is selected from the corresponding group using
the Haar measure. The groups we consider here are the orthogonal group for
Majorana fermions and the unitary group for Dirac fermions. These Haar values
correspond to the evolution with an arbitrary number of interacting fermions
($q>2$). Second, we consider the case where the evolution matrix is quadratic
in fermions ($q=2$), i.e., $U=e^{c^{\dagger}hc}$. In this case, the matrix
values are again chosen with respect to the Haar measure, but the space of
integration is smaller compared to the full Haar average.

### Majorana Fermions

#### Haar averages

We are interested in the calculation of the Frame potential in the case when
$U$ is an orthogonal matrix taken from the Haar measure:

$F^{(k)}_{\rm Haar}=\frac{1}{\mathcal{N}^{2k}}\int_{\rm Haar}dU\left|{\rm
Tr}\left[U\right]\right|^{2k}.$ (124)

Average of the traces with respect to the Haar measure was intensively studied
in the Random Matrix theory. It was shown that for sufficiently large size of
matrices, averaging over orthogonal group can be expressed as an averaging
over standard normal variables Diaconis1994OnTE ; Diaconis2000LinearFO ;
Vasilchuk:2003xz ; Ram1995CharactersOB . This gives as the result

$\frac{1}{\mathcal{N}^{2k}}\int_{\rm Haar}dU\left|{\rm
Tr}\left[U\right]\right|^{2k}=\frac{1}{\mathcal{N}^{2k}}\mathbb{E}(Z^{2k})=\frac{1}{\mathcal{N}^{2k}}(2k)!!,$
(125)

where $Z$ is a normal variable. Fair to notice that these formulae are true
only for $k\leq\frac{L}{4}$, where $L$ is the matrix size.

#### Gaussian-Haar averages

If we assume that the unitary evolution $U$ is Gaussian, i.e.
$U=e^{\frac{1}{2}h_{ij}\chi_{i}\chi_{j}}$ , where $\chi_{j}$ are Majorana
fields, then we can express ${\rm Tr}(U)={\rm Pf}(1+e^{h_{ij}})$(it was shown
in Klich2014ANO ) where $e^{h_{ij}}=u$ is an orthogonal matrix with the size
$L\times L$ ($L$ is even). We are interested in the case when matrices $u$ are
random with respect to the Haar measure. Notice that we are averaging over
matrices with $\det(u)=1$ because $u=e^{\frac{1}{2}h}$ where $h$ is
antisymmetric matrix (all the eigenvalues of $h$ have conjugated pairs and lie
on the imaginary axis, it means that eigenvalues of $u$ are
$\lambda_{1}=e^{ib_{1}},\lambda^{*}_{1}=e^{-ib_{1}},...$ and their product is
1). This gives us the following averaging

$\left<F^{k}\right>_{\rm gHaar}=\left<{\rm
Tr}^{(2)}[\rho_{1}(t)\rho_{2}(t)]^{k}\right>_{\rm
gHaar}=\frac{1}{2^{Lk}}\int\left|\det(1+u)\right|^{k}d\mu(u)=\\\
=\frac{1}{2^{Lk}}\int\left|\prod^{L}_{j=1}(1+e^{i\theta_{j}})\right|^{k}\prod_{i<j}\left|2(\cos\theta_{i}-\cos\theta_{j})\right|^{2}\frac{d\theta_{1}..d\theta_{L/2}}{\left(\frac{L}{2}\right)!(2\pi)^{\frac{L}{2}}}=\\\
=\frac{1}{2^{Lk}}\int\left|\prod^{L/2}_{j=1}(2+2\cos(\theta_{j}))\right|^{k}\prod_{i<j}\left|2(\cos\theta_{i}-\cos\theta_{j})\right|^{2}\frac{d\theta_{1}..d\theta_{L/2}}{\left(\frac{L}{2}\right)!(2\pi)^{\frac{L}{2}}}.$
(126)

In the third line we used the fact that all eigenvalues have its conjugated
pair.The measure in this expression comes from the fact that we integrate
other orthogonal matrices with the determinant equal to 1 and even size(see
chapter 2.6 in Forrester ). Doing change of variables $\cos(\theta)=t$ and
calculating the Jacobian of this transformation we can get

$\left<F^{k}\right>_{\rm
gHaar}=\frac{1}{2^{Lk}}\frac{1}{F^{0}}\int_{-1}^{1}dt_{1}..\int_{-1}^{1}dt_{L/2}\prod^{L/2}_{j=1}|2(1+t_{j})|^{k}\prod_{i<j}|t_{i}-t_{j}|^{2}\prod_{k=1}^{L/2}(1-t_{k}^{2})^{-\frac{1}{2}}.$
(127)

Now we need to calculate these integrals. For simplicity, we can introduce the
function to name these integrals :

$F=2^{\frac{Lk}{2}}\int_{-1}^{1}dt_{1}..\int_{-1}^{1}dt_{L/2}\prod^{L/2}_{j=1}|1+t_{j}|^{k}\prod_{i<j}|t_{i}-t_{j}|^{2}\prod_{k=1}^{L/2}(1-t_{k}^{2})^{-\frac{1}{2}}.$
(128)

Here we can notice that it is just a special case of Selberg integral (see
chapters $3.6,3.7,4.1,4.7$ in Forrester )

$F=2^{\frac{Lk}{2}}2^{\frac{L}{2}+\frac{L(k-1)}{2}+\frac{L}{2}(\frac{L}{2}-1)}S_{\frac{L}{2}}\left(k-\frac{1}{2},-\frac{1}{2},1\right)=2^{Lk}2^{\frac{L}{2}(\frac{L}{2}-1)}S_{\frac{L}{2}}\left(k-\frac{1}{2},-\frac{1}{2},1\right).$
(129)

The Selberg integral has well known form in the terms of the Gamma functions,
which is:

$S_{L}(\lambda_{1},\lambda_{2},\lambda)=\int_{0}^{1}dt_{1}...\int_{0}^{1}dt_{L}\prod_{l=1}^{L}t_{l}^{\lambda_{1}}(1-t_{l})^{\lambda_{2}}\prod_{k<j}|t_{k}-t_{j}|^{2\lambda}=\\\
=\prod_{j=0}^{L-1}\frac{\Gamma(\lambda_{1}+1+j\lambda)\Gamma(\lambda_{2}+1+j\lambda)\Gamma(1+(j+1)\lambda)}{\Gamma(\lambda_{1}+\lambda_{2}+2+(L+j-1)\lambda)\Gamma(\lambda+1)}.$
(130)

Inserting this solution into the Eq. (127) we can have the final expression

$\left<F^{k}\right>_{\rm
gHaar}=\frac{1}{2^{Lk}}\frac{2^{Lk}2^{\frac{L}{2}(\frac{L}{2}-1)}S_{\frac{L}{2}}\left(k-\frac{1}{2},-\frac{1}{2},1\right)}{2^{\frac{L}{2}(\frac{L}{2}-1)}S_{\frac{L}{2}}\left(-\frac{1}{2},-\frac{1}{2},1\right)}=\frac{S_{\frac{L}{2}}\left(k-\frac{1}{2},-\frac{1}{2},1\right)}{S_{\frac{L}{2}}\left(-\frac{1}{2},-\frac{1}{2},1\right)}=\\\
=\prod_{j=0}^{\frac{L}{2}-1}\frac{\Gamma(k+\frac{1}{2}+j)}{\Gamma(\frac{1}{2}+j)}\frac{\Gamma(\frac{L}{2}+j)}{\Gamma(k+\frac{N}{2}+j)}.$
(131)

Therefore for $k=1$ we have

$\left<F^{1}\right>_{\rm gHaar}=\frac{2}{2^{L}}=\left<F^{1}\right>_{\rm
Haar},$ (132)

which is equal to the Frame potential in the case when the evolution matrix
was taken from the Haar ensemble of Orthogonal matrices. We are interested in
the large $L$ limit expansion of the Frame potential for arbitrary k,
therefore we can use ”Stirling like” formula for the Barnes functions, which
gives:

$\\\
lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left(\left<F^{k}\right>_{\rm
gHaar}\right)\right)\sim-k\log(2)+\frac{k(k-1)}{2}\frac{1}{L}\log(L).$ (133)

Let us compare this value to the full Haar average. Remember that for the full
Haar we have:

$\\\
lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left(\left<F^{k}\right>_{\rm
Haar}\right)\right)\sim-k\log(2)+\log((2k)!!)\frac{1}{L}.$ (134)

Therefore for large $L$ and fixed $k\leq\frac{L}{4}$:

$\left<F^{k}\right>_{\rm gHaar}\geq\left<F^{k}\right>_{\rm Haar},$ (135)

so we always have full Haar average being smaller than the Fermionic one.

### Dirac Fermions

#### Haar averages

Now let’s do a similar calculation in the Dirac fermions case. We are
interested in the calculation of the Frame potential in the case when $U$ is a
unitary matrix taken from the Haar measure:

$F^{(k)}_{\rm Haar}=\frac{1}{\mathcal{N}^{2k}}\int_{\rm Haar}dU\left|{\rm
Tr}\left[U\right]\right|^{2k}.$ (136)

In this case we also can use the result from the Random matrix theory. Which
gives

$\frac{1}{\mathcal{N}^{2k}}\int_{\rm Haar}dU\left|{\rm
Tr}\left[U\right]\right|^{2k}=\frac{1}{\mathcal{N}^{2k}}\mathbb{E}(Z^{2k})=\frac{1}{\mathcal{N}^{2k}}(k)!,$
(137)

where $Z$ is a standard complex normal variable. Again, we notice that these
formulae are true only for $k\leq\frac{L}{2}$, where $L$ is the matrix size.

#### Gaussian-Haar averages

Again we assume that the unitary evolution $U$ is Gaussian, i.e.
$U=e^{ih_{ij}c^{\dagger}_{i}c_{j}}$ , where $c_{j}$ are Dirac fermions, then
we can express ${\rm Tr}(U)={\rm Pf}(1+u)$ where $e^{ih_{ij}}=u$ is a unitary
matrix with the size $L\times L$ ($L$ is even). The matrices $u$ are again
random with respect to Haar measure. This gives us the following averaging,

$\left<F^{k}\right>_{\rm gHaar}=\left<{\rm
Tr}^{(2)}[\rho_{1}(t)\rho_{2}(t)]^{k}\right>_{\rm
gHaar}=\frac{1}{2^{2Lk}}\int\left|\det(1+u)\right|^{k}d\mu(u)=\\\
=\frac{1}{2^{2Lk}}\int\prod^{L}_{j=1}\left|(1+e^{i\theta_{j}})\right|^{k}\prod_{i<j}\left|e^{i\theta_{i}}-e^{i\theta_{j}}\right|^{2}\frac{d\theta_{1}..d\theta_{N}}{N!(2\pi)^{N}}.$
(138)

The measure in this expression comes from the fact that we integrate other
unitary matrices(see chapter 2.6 in Forrester ). This integral is a Morris
integral which can be expressed in terms of Gamma functions:

$M_{L}(b,b,\lambda)=\int\prod^{L}_{j=1}\left|(1+e^{i\theta_{j}})\right|^{2b}\prod_{i<j}\left|e^{i\theta_{i}}-e^{i\theta_{j}}\right|^{2\lambda}\frac{d\theta_{1}..d\theta_{N}}{N!(2\pi)^{N}}=\\\
=\prod_{j=0}^{L-1}\frac{\Gamma(\lambda
j+2b+1)\Gamma(\lambda(j+1)+1)}{\Gamma(\lambda j+b+1)\Gamma(\lambda
j+2b+1)\Gamma(1+\lambda)}.$ (139)

Inserting this into the original expression for the Frame potential we can
get:

$\left<F^{k}\right>_{\rm
gHaar}=\frac{1}{2^{2Lk}}\frac{M_{L}(k,k,1)}{M_{L}(0,0,1)}=\frac{1}{2^{2Lk}}\prod^{L-1}_{j=0}\frac{\Gamma(j+2k+1)\Gamma(j+1)}{\Gamma(j+k+1)^{2}}.$
(140)

Again we are interested in the large $L$ limit expansion:

$\\\
lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left(\left<F^{k}\right>_{\rm
gHaar}\right)\right)\sim-2k\log(2)+k^{2}\frac{1}{L}\log(L).$ (141)

Comparing this to the full Haar average

$\\\
lim_{L\rightarrow\infty}\left(\frac{1}{L}\log\left(\left<F^{k}\right>_{\rm
Haar}\right)\right)\sim-2k\log(2)+\log(k!)\frac{1}{L},$ (142)

we again have

$\left<F^{k}\right>_{\rm gHaar}\geq\left<F^{k}\right>_{\rm Haar},$ (143)

so the full Haar average is again smaller than the fermionic one.

## References

  * [1] Patrick Hayden and John Preskill.  Black holes as mirrors: quantum information in random subsystems.  Journal of High Energy Physics, 2007(09):120–120, September 2007\.
  * [2] Stephen H. Shenker and Douglas Stanford.  Black holes and the butterfly effect.  Journal of High Energy Physics, 2014(3), March 2014.
  * [3] Stephen H. Shenker and Douglas Stanford.  Multiple shocks.  Journal of High Energy Physics, 2014(12), December 2014.
  * [4] Juan Maldacena, Stephen H. Shenker, and Douglas Stanford.  A bound on chaos.  Journal of High Energy Physics, 2016(8), August 2016.
  * [5] Akram Touil and Sebastian Deffner.  Information scrambling versus decoherence—two competing sinks for entropy.  PRX Quantum, 2:010306, Jan 2021.
  * [6] Ahmed Almheiri, Thomas Hartman, Juan Maldacena, Edgar Shaghoulian, and Amirhossein Tajdini.  Replica wormholes and the entropy of hawking radiation.  Journal of High Energy Physics, 2020(5), May 2020.
  * [7] Shenglong Xu and Brian Swingle.  Locality, quantum fluctuations, and scrambling.  Phys. Rev. X, 9:031048, Sep 2019.
  * [8] Shao-Kai Jian, Gregory Bentsen, and Brian Swingle.  Linear growth of circuit complexity from brownian dynamics, 2022.
  * [9] D. Gross, K. Audenaert, and J. Eisert.  Evenly distributed unitaries: On the structure of unitary designs.  Journal of Mathematical Physics, 48(5):052104, May 2007.
  * [10] Daniel A. Roberts and Beni Yoshida.  Chaos and complexity by design.  Journal of High Energy Physics, 2017(4), April 2017.
  * [11] Sumner N. Hearth, Michael O. Flynn, Anushya Chandran, and Chris R. Laumann.  Unitary k-designs from random number-conserving quantum circuits.  2023\.
  * [12] Alexios Christopoulos, Pierre Le Doussal, Denis Bernard, and Andrea De Luca.  Universal out-of-equilibrium dynamics of 1d critical quantum systems perturbed by noise coupled to energy.  Phys. Rev. X, 13:011043, Mar 2023.
  * [13] Mark Srednicki.  Chaos and quantum thermalization.  Phys. Rev. E, 50:888–901, Aug 1994.
  * [14] Luca D'Alessio, Yariv Kafri, Anatoli Polkovnikov, and Marcos Rigol.  From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics.  Advances in Physics, 65(3):239–362, May 2016.
  * [15] Silvia Pappalardi, Laura Foini, and Jorge Kurchan.  Eigenstate thermalization hypothesis and free probability.  Phys. Rev. Lett., 129:170603, Oct 2022.
  * [16] Sanjay Hortikar and Mark Srednicki.  Random matrix elements and eigenfunctions in chaotic systems.  Phys. Rev. E, 57:7313–7316, Jun 1998.
  * [17] J. M. Deutsch.  Quantum statistical mechanics in a closed system.  Phys. Rev. A, 43:2046–2049, Feb 1991.
  * [18] Maxime Lucas, Lorenzo Piroli, Jacopo De Nardis, and Andrea De Luca.  Generalized deep thermalization for free fermions.  Phys. Rev. A, 107:032215, Mar 2023.
  * [19] Jordan S. Cotler, Daniel K. Mark, Hsin-Yuan Huang, Felipe Hernández, Joonhee Choi, Adam L. Shaw, Manuel Endres, and Soonwon Choi.  Emergent quantum state designs from individual many-body wave functions.  PRX Quantum, 4:010311, Jan 2023.
  * [20] Joonhee Choi, Adam L. Shaw, Ivaylo S. Madjarov, Xin Xie, Ran Finkelstein, Jacob P. Covey, Jordan S. Cotler, Daniel K. Mark, Hsin-Yuan Huang, Anant Kale, Hannes Pichler, Fernando G. S. L. Brandão, Soonwon Choi, and Manuel Endres.  Preparing random states and benchmarking with many-body quantum chaos.  Nature, 613(7944):468–473, January 2023.
  * [21] Wen Wei Ho and Soonwon Choi.  Exact emergent quantum state designs from quantum chaotic dynamics.  Phys. Rev. Lett., 128:060601, Feb 2022.
  * [22] Harshank Shrotriya and Wen Wei Ho.  Nonlocality of deep thermalization.  2023\.
  * [23] Matteo Ippoliti and Wen Wei Ho.  Solvable model of deep thermalization with distinct design times.  Quantum, 6:886, December 2022.
  * [24] Subir Sachdev and Jinwu Ye.  Gapless spin-fluid ground state in a random quantum heisenberg magnet.  Phys. Rev. Lett., 70:3339–3342, May 1993.
  * [25] Phil Saad, Stephen H. Shenker, and Douglas Stanford.  A semiclassical ramp in syk and in gravity, 2018.
  * [26] Shao-Kai Jian and Brian Swingle.  Note on entropy dynamics in the brownian SYK model.  Journal of High Energy Physics, 2021(3), March 2021.
  * [27] Christoph Sünderhauf, Lorenzo Piroli, Xiao-Liang Qi, Norbert Schuch, and J. Ignacio Cirac.  Quantum chaos in the brownian SYK model with large finite n : OTOCs and tripartite information.  Journal of High Energy Physics, 2019(11), November 2019.
  * [28] Pedro Ribeiro, Julien Vidal, and Rémy Mosseri.  Exact spectrum of the lipkin-meshkov-glick model in the thermodynamic limit and finite-size corrections.  Phys. Rev. E, 78:021106, Aug 2008.
  * [29] Xiao Chen and Tianci Zhou.  Quantum chaos dynamics in long-range power law interaction systems.  Phys. Rev. B, 100:064305, Aug 2019.
  * [30] Denis Bernard, Fabian H. L. Essler, Ludwig Hruza, and Marko Medenjak.  Dynamics of fluctuations in quantum simple exclusion processes.  SciPost Phys., 12:042, 2022.
  * [31] Denis Bernard and Tony Jin.  Solution to the quantum symmetric simple exclusion process: The continuous case.  Communications in Mathematical Physics, 384(2):1141–1185, Jun 2021\.
  * [32] Tobias Swann, Denis Bernard, and Adam Nahum.  Spacetime picture for entanglement generation in noisy fermion chains, 2023.
  * [33] Denis Bernard and Lorenzo Piroli.  Entanglement distribution in the quantum symmetric simple exclusion process.  Phys. Rev. E, 104:014146, Jul 2021.
  * [34] Subir Sachdev and Jinwu Ye.  Gapless spin fluid ground state in a random, quantum Heisenberg magnet.  Phys. Rev. Lett., 70:3339, 1993.
  * [35] Yingfei Gu, Alexei Kitaev, Subir Sachdev, and Grigory Tarnopolsky.  Notes on the complex Sachdev-Ye-Kitaev model.  JHEP, 02:157, 2020.
  * [36] Yingfei Gu, Andrew Lucas, and Xiao-Liang Qi.  Spread of entanglement in a sachdev-ye-kitaev chain.  Journal of High Energy Physics, 2017:1–44, 2017.
  * [37] Geoff Penington, Stephen H. Shenker, Douglas Stanford, and Zhenbin Yang.  Replica wormholes and the black hole interior.  2019\.
  * [38] Ramamurti Shankar.  Quantum Field Theory and Condensed Matter.  Cambridge University Press, 8 2017.
  * [39] A. V. Lunkin, A. Yu. Kitaev, and M. V. Feigel’man.  Perturbed Sachdev-Ye-Kitaev Model: A Polaron in the Hyperbolic Plane.  Phys. Rev. Lett., 125(19):196602, 2020.
  * [40] M. Nakahara.  Geometry, topology and physics.  2003\.
  * [41] Shao-Kai Jian and Brian Swingle.  Phase transition in von neumann entanglement entropy from replica symmetry breaking.  2021\.
  * [42] A Yu Kitaev.  Unpaired majorana fermions in quantum wires.  Phys.–Usp., 44(10S):131–136, October 2001.
  * [43] Michele Fava, Lorenzo Piroli, Tobias Swann, Denis Bernard, and Adam Nahum.  Nonlinear sigma models for monitored dynamics of free fermions.  2 2023.
  * [44] Andrew Baker.  Matrix Groups An Introduction to Lie Group Theory.
  * [45] Shao-Kai Jian and Brian Swingle.  Phase transition in von neumann entanglement entropy from replica symmetry breaking, 2021.
  * [46] Michele Fava, Lorenzo Piroli, Tobias Swann, Denis Bernard, and Adam Nahum.  Nonlinear sigma models for monitored dynamics of free fermions, 2023.
  * [47] Brian Skinner, Jonathan Ruhman, and Adam Nahum.  Measurement-induced phase transitions in the dynamics of entanglement.  Phys. Rev. X, 9:031009, Jul 2019.
  * [48] Stefano Antonini, Gregory Bentsen, ChunJun Cao, Jonathan Harper, Shao-Kai Jian, and Brian Swingle.  Holographic measurement and bulk teleportation.  Journal of High Energy Physics, 2022(12), December 2022.
  * [49] Persi Diaconis and Mehrdad Shahshahani.  On the eigenvalues of random matrices.  Journal of Applied Probability, 31:49 – 62, 1994.
  * [50] Persi Diaconis and Steven N. Evans.  Linear functionals of eigenvalues of random matrices.  2000\.
  * [51] V. Vasilchuk.  On the Moments of the Traces of Unitary and Orthogonal Random Matrices.  eConf, C0306234:1207–1213, 2003.
  * [52] Arun Ram.  Characters of brauer’s centralizer algebras.  Pacific Journal of Mathematics, 169:173–200, 1995.
  * [53] Israel Klich.  A note on the full counting statistics of paired fermions.  Journal of Statistical Mechanics: Theory and Experiment, 2014, 2014\.
  * [54] Forrester.  Log-gases and random matrices.