aainstitutetext: The Laboratory for Quantum Gravity & Strings, Department of
Mathematics & Applied Mathematics, University of Cape Town, Cape Town, South
Africabbinstitutetext: National Institute for Theoretical and Computational
Sciences (NITheCS),
South Africa

# A Study of the SYK2 Model with Twisted Boundary Conditions

Jeff Murugan a Ruach Pillay Slayen a,b and Hendrik J. R. Van Zyl

###### Abstract

We study a version of the 2-body Sachdev-Ye-Kitaev (SYK2) model whose complex
fermions exhibit twisted boundary conditions on the thermal circle. As we
show, this is physically equivalent to coupling the fermions to a
1-dimensional external gauge field $A(t)$. In the latter formulation, the
gauge field itself can be thought of as arising from a radial symmetry
reduction of a $(2+1)$-dimensional Chern-Simons gauge field
$A_{\mu}(t,\mathbf{x})$. Using the diagnostic tools of the out-of-time-order
correlator (OTOC) and spectral form factor (SFF), which probe the sensitivity
to initial conditions and the spectral statistics respectively, we give a
detailed and pedagogical study of the integrable/chaotic properties of the
model. We find that the twisting has no effect on the OTOCs and, by extension,
the early-time chaos properties of the model. It does, however, have two
notable effects on the spectral form factor; an enhancement of the early-time
slope and the emergence of an explicit disorder scale needed for the
manifestation of zero modes. These zero modes are responsible for the late-
time exponential ramp in the quadratic SYK model.

††dedication: Dedicated to the memory of Roman Jackiw: A pioneer of many
dimensions, in few.

##  1 Introduction

Sachdev-Ye-Kitaev (SYK) models Sachdev:1992fk ; Maldacena:2016hyu ;
Polchinski:2016xgd ; Garcia-Garcia:2016mno ; Gross:2016kjj ; Gross:2017aos ;
Saad:2018bqo ; Rosenhaus:2018dtp ; Gu:2019jub  detail the properties of a
family of (0+1)-dimensional quantum mechanics of $N\gg 1$ Majorana fermions,
$\chi_{a}$, with all-to-all random $q$-fermi interactions. The $q=0$ free
theory with action,

$\displaystyle
I_{\mathrm{free}}=\sum_{a=1}^{N}\int\\!\\!dt\,\,\chi_{a}\dot{\chi}_{a}\,,$ (1)

is an exactly conformal quantum mechanics of $2^{N/2}$ zero-energy states and
a global $SO(N)$ flavour symmetry. The free theory admits $N$ relevant
deformations implemented by perturbing (1) with fermion monomial operators. In
this sense, the standard ($q=4$) SYK model Sachdev:1992fk ,

$\displaystyle
I_{\mathrm{SYK}}=\int\\!\\!dt\left(\sum_{a=1}^{N}\,\,\chi_{a}\dot{\chi}_{a}+\sum_{a,b,c,d=1}^{N}\frac{1}{4!}J_{abcd}\,\,\chi_{a}\chi_{b}\chi_{c}\chi_{d}\right)\,,$
(2)

can be considered a perturbation of (1) by a quartic monomial that preserves
the full $SO(N)$ flavour group as long as the fermion coupling $J_{abcd}$ is
disordered. Here the couplings are time-independent (and so the disorder is
quenched) and drawn randomly from a Gaussian distribution with
$\overline{J_{abcd}}=0$ and $\overline{J_{abcd}J^{abcd}}=3!J^{2}/N^{3}$. This
model has some well-documented and remarkable properties Maldacena:2016hyu :

  * •

It is solvable. At zero temperature and large $N$, the solution to the
Schwinger-Dyson equation for the fermion 2-point function is

$\displaystyle\langle\chi_{a}(t)\chi_{b}(0)\rangle=\left(\frac{1}{4\pi
J^{2}}\right)^{1/4}\frac{\mathrm{sgn(t)}}{|t|^{1/2}}\delta_{ab}\,,$ (3)

in the infrared and for $t\gg 1/J$. Dimensional analysis of this solution
shows that $[\chi]=1/4$ in the IR.

  * •

Indeed, the SYK model manifests an emergent conformal symmetry at low energies
and large $N$.

  * •

The model is famously maximally chaotic in the sense that the quantum Lyapunov
exponent, extracted from the out-of-time-ordered thermal fermion 4-point
function, saturates the Maldacena-Shenker-Stanford (MSS) bound
Maldacena:2015waa  $\lambda_{L}\leq 2\pi T$ at temperature $T$.

With all this going for it, one might very reasonably presume the SYK model to
furnish the CFT dual to an AdS2/CFT1 correspondence Sarosi:2017ykf ;
Almheiri:2014cka ; Jevicki:2016bwu ; Engelsoy:2016xyb ; Jevicki:2016ito ;
Cvetic:2016eiv . However, this hope was convincingly dashed in Jensen:2016pah
where it was pointed out that the SYK model cannot have a conventional weakly-
coupled, weakly-curved gravity dual for two crucial reasons:

  1. 1.

The large $N$ Majoranas $\chi_{a}(t)$ of the SYK model are dual to $N$
degenerate bulk fermions $\Psi_{a}(t,x)$ whose contribution to the bulk
partition function is comparable to the classical saddle, consequently
invalidating the saddle-point approximation.

  2. 2.

The model does not exhibit the large $N$ factorization expected of theories
with a conventional gravitational dual.

Interpreting these observations is tricky. In Jensen:2016pah  for example, it
was suggested that perhaps the correct statement is that the singlet sector of
the SYK model is dual to a 2-dimensional higher spin theory. To support this
assertion, it was shown there that the gravitational dynamics near an AdS2
throat can be re-written in terms of an effective quantum hydrodynamical
action that precisely matches with the Schwarzian action found in
Maldacena:2016hyu  from the low energy, large $N$ solution of the SYK model.
Another view point is that the ills of the SYK model could be cured by gauging
a large subgroup of the $SO(N)$ global flavour symmetry Jensen:2016pah . This
has, however, proven a difficult goal to realize.

With the goal of constructing such an AdS2/gauged-SYK correspondence in mind
then, we began a study of a generalization of the original SYK model, in which
the $N$ Majoranas are replaced by $N\gg 1$ spinless complex fermions - the so-
called charged SYK model Sachdev:2015efa ; Bulycheva:2017uqj . This system
possesses a global $U(1)$ symmetry and for the case $q=4$, for example, is
described by the action

$\displaystyle I_{\mathrm{CSYK}}=\int
dt\left[\sum_{a=1}^{N}\psi_{a}^{\dagger}\left(\frac{d}{dt}-\mu\right)\psi_{a}-\sum_{a,b,c,d=1}^{N}J_{ab;cd}\,\psi^{\dagger}_{a}\psi^{\dagger}_{b}\psi_{c}\psi_{d}\right]\,.$
(4)

Here, $\mu$ is a chemical potential that tunes the $U(1)$ charge
$\displaystyle Q=\sum_{a=1}^{N}\langle\psi^{\dagger}_{a}\psi_{a}\rangle$ and
the $J_{ab;cd}$ are now complex, independent, random couplings drawn from a
Gaussian distribution with zero mean and that satisfy

$\displaystyle J_{ab;cd}$ $\displaystyle=$ $\displaystyle-
J_{ba;cd}=-J_{ab;dc}=J^{*}_{cd;ab}\,,$
$\displaystyle\overline{|J_{ab;cd}|^{2}}$ $\displaystyle=$ $\displaystyle
J^{2}\,.$

In addition, the fermions close the algebra $\\{\psi_{a},\psi_{b}\\}=0$,
$\\{\psi_{a},\psi_{b}^{\dagger}\\}=\delta_{ab}$. This complex SYK model shares
many properties with its more well-understood real counterpart and that we
will discuss in more detail in section 3 below. Our original goal in this
article was to gauge the global $U(1)$ by coupling the complex fermions to a
$U(1)$ gauge field $A(t)$ and investigating its properties. In particular, we
set out with the aim of understanding if the resulting model,

  * •

possesses a sensible classical gravitational dual,

  * •

retains the chaos structure of the original (complex) SYK model and,

  * •

is sensitive to the breaking of the reparameterization invariance as the
theory is moved away from its conformal limit.

Importantly, in (0+1)-dimensions the gauge field is non-dynamical since a
Maxwell kinetic term cannot be added to the action but, being an odd-
dimensional spacetime means that a Chern-Simons term can. In this case it is
simply linear in $A(t)$. The result is a form of topological quantum mechanics
Dunne:1989hv ; Floreanini:1991iu ; Jackiw:1997kp  originally studied by Dunne
et.al. in the context of finite temperature Chern-Simons theory Dunne:1996yb ,
with an added layer of random interactions. However111We would like to thank
the anonymous referee for highlighting this point to us., as we will show
below, this coupling to an external gauge field - i.e. one that is not
integrated over in the path integral - is equivalent to a twisting of the
fermion boundary conditions on the thermal circle.

This twisting of the boundary conditions of random mass, and/or randomly
interacting fermions, is expected to introduce novel effects to the disorder
and interactions that can drastically change the behavior of the system.
Depending on the parameters of the model, such as the mass distribution, the
“Chern-Simons” coupling, and the chemical potential, the system can exhibit
different phases with different properties, such as localization,
delocalization, Luttinger liquid behavior, and gapping Chowdhury:2021qpy .
These phases can be characterized by various physical observables, such as the
energy levels, the charge density, the gap equation, and the density of
states. As a warmup to the problem we really want to study, viz the charged
SYK4 model with twisted boundary conditions in one spacetime dimension, in
this article we turn our attention to the simpler, integrable, quadratic SYK
model. Even in this much simplified setting, we show that there are some
interesting lessons to be learnt, lessons that we will port to the canonical
$q=4$ model in forthcoming work.

The SYK2 model, with a disordered 2-fermi interaction term, is equivalent to a
theory of free fermions with a random mass matrix. As we shall see, using the
standard SYK approach to formulate the theory in terms of the collective field
variables, the Schwinger-Dyson equations may be combined into a single
quadratic equation. This may then be solved for the large $N$ two-point
correlation function for arbitrary coupling strength. This is in contrast to
the $q>2$ cases of the model, in which one can only find an analytic
expression for the large $N$ two-point function in the strong-coupling/low-
energy conformal limit. This great simplification makes the SYK2 model
particularly amenable to analytic study, where it has been shown to have the
following basic features: firstly, it is an integrable model, as the OTOC does
not exhibit the exponential growth characteristic of quantum chaos
Gross:2016kjj . Secondly, it has an average spectral density given by the
Wigner semicircle law in the large $N$ limit Gross:2016kjj ; Cotler:2017jue ,
just like the Gaussian RMT ensembles. However, despite these seemingly trivial
features, the model has proven to be a remarkably rich object of study.

The rest of the article is structured as follows: In the interests of
pedagogy, in Section 2 we introduce the coupling of the complex SYK model to
the external “Chern-Simons” gauge field a la Dunne:1996yb ; Dunne:1989hv , and
show how the resulting model is equivalent to twisting the fermion boundary
conditions on the thermal circle. Section 3 is taken up with describing some
of the technicalities of the disorder averaging procedures we use in the rest
of the paper. The crux of the paper is contained in Sections 4 and 5 where we
compute the out-of-time-order correlator (or, more precisely, the double
commutator) and the spectral form factor of the model respectively. We
summarise our results and comment on future directions in Section 6. Finally,
we collect some technical details of our computations such as the derivation
of the 2-replica spectral form factor, the partition functions and the mass
matrices in the appendices.

##  2 Chern-Simons quantum mechanics

Our starting point is the following Lagrangian Dunne:1989hv ;
Floreanini:1991iu ; Dunne:1996yb ; Jackiw:1997kp , written in the imaginary
time formalism, is given by

$L=\sum_{j=1}^{N_{f}}\psi_{j}^{\dagger}(\partial_{\tau}-iA+m)\psi_{j}\,,$ (5)

where the fermions close the algebra

$\left\\{\psi_{i}^{\dagger},\psi_{j}\right\\}=\delta_{ij}.$ (6)

$U(1)$ gauge transformations act on the fields as

$\displaystyle\psi$ $\displaystyle\rightarrow$ $\displaystyle
e^{i\omega(\tau)}\psi$ $\displaystyle A$ $\displaystyle\rightarrow$
$\displaystyle A+\partial_{\tau}\omega.$ (7)

Note that the path integral corresponding to (5) is invariant under the gauge
transformation. In order to maintain the antiperiodicity of the gauge
transformed fermion fields, the phase must satisfy

$\displaystyle\omega(\beta)-\omega(0)=2\pi k,\quad k\in\mathbb{Z}.$ (8)

Transformations for which $k\neq 0$ are called large gauge transformations.
This gauge freedom means that we need only consider constant gauge field
configurations since, for any $A(\tau)$ we may always make the gauge
transformation 222Note that this is a small gauge transformation since
$\omega(\beta)=\omega(0).$

$\displaystyle\omega(\tau)=\frac{a}{\beta}\tau-\int_{0}^{\tau}d\tau^{\prime}A(\tau^{\prime}),$
(9)

where $a=\int_{0}^{\beta}d\tau A(\tau)$, to get to the constant configuration

$\displaystyle A(\tau)\rightarrow\frac{a}{\beta}.$ (10)

It should be emphasised that, although the gauge field coupling naively
resembles a constant chemical potential, its transformation properties are
different. In particular, all expressions in the imaginary time or real-time
formalism should be invariant under the large gauge transformations of the
form $a\rightarrow a+2\pi k$ for integer $k$.

Furthermore we emphasize that, since the gauge field is not being integrated
in the path integral, one could redefine the fermion fields as

$\psi^{\prime}(\tau)=\psi(\tau)e^{-i\int_{0}^{\tau}d\tau^{\prime}A(\tau^{\prime})}$

which eliminates the gauge field in (5). Where the fields $\psi$ satisfy the
usual antiperiodicity condition, the fields $\psi^{\prime}$ satisfy

$\psi^{\prime}(\beta)=e^{-ia}\psi^{\prime}(0).$

As such, the introduction of the gauge field is equivalent to introducing a
twisted boundary condition.

One of the many features of this model is that its Green’s function can be
computed analytically. In particular in the $m\rightarrow 0$ limit,

$\langle\psi_{i}^{\dagger}(\tau_{1})\psi_{j}(\tau_{2})\rangle\equiv
G^{0}_{ij}(\tau_{1}-\tau_{2})=\left\\{\begin{array}[]{cc}\frac{\delta_{ij}}{2}\left(-1-i\tan\left(\frac{a}{2}\right)\right)e^{i\frac{a}{\beta}(\tau_{1}-\tau_{2})}&-\beta<\tau_{1}-\tau_{2}<0\\\
-\frac{i\delta_{ij}}{2}\tan\left(\frac{a}{2}\right)&\tau_{1}-\tau_{2}=0\\\
\frac{\delta_{ij}}{2}\left(1-i\tan\left(\frac{a}{2}\right)\right)e^{i\frac{a}{\beta}(\tau_{1}-\tau_{2})}&0<\tau_{1}-\tau_{2}<\beta\\\
\frac{i\delta_{ij}}{2}\tan\left(\frac{a}{2}\right)&(\tau_{1}-\tau_{2})=\beta.\end{array}\right.$
(11)

As explained in some detail in Dunne:1996yb , the Lagrangian (5) can be
derived as a reduction of (2+1)-dimensional Chern-Simons matter Jackiw:1997kp
, and provides a simplified toy model in which some of the subtleties of its
higher-dimensional avatar may be clarified. Inspired by the Sachdev-Ye-Kitaev
model Sachdev:1992fk , we will instead be interested in perturbing the model
with random $q$-body fermion operators. In this sense, this is an SYKq model
coupled to an external gauge field. As noted in the introduction, the SYKq
model has a number of interesting features such as its exact solvability at
large $N$, and the fact that it is maximally (quantum) chaotic
Maldacena:2016hyu ; Polchinski:2016xgd ; Garcia-Garcia:2016mno ; Gross:2016kjj
; Gross:2017aos ; Saad:2018bqo ; Rosenhaus:2018dtp . The latter suggests that
it should admit a holographic description as a black hole (in two dimensions)
Penington:2019kki ; Trunin:2020vwy , while the former allows (in principle)
for details of the holographic dictionary to be systematically unpacked.

A special case of this system is when $q=2$ where the model is non-chaotic and
may be solved exactly Gross:2016kjj  and so serves as an important test case
for the gauged SYK model. In particular, we would like to understand whether
the inclusion of the gauge field spoils its solvability properties and/or
introduces chaotic regions in its parameter space333In particular, it would be
interesting to explore whether integrable-chaotic transitions are present in
parameter space Sorokhaibam:2019qho .. To be explicit, we will focus on the
model,

$L=\sum_{j=1}^{N_{f}}\psi_{j}^{\dagger}(\partial_{\tau}-iA)\psi_{j}+\sum_{ij}\psi^{\dagger}_{i}m_{ij}\psi_{j}$
(12)

where the random mass matrix $m_{ij}=m_{ji}^{*}$ is drawn from a Gaussian
distribution with zero mean,

$P(m_{ij})=\left\\{\begin{array}[]{cc}(\frac{N}{\pi
m^{2}})^{\frac{1}{4}}e^{-\frac{N}{2m^{2}}|m_{ij}|^{2}}&i\neq j\\\
\sqrt{\frac{N}{2\pi
m^{2}}}e^{-\frac{N}{2m^{2}}m_{ii}^{2}}&i=j\end{array}\right.$ (13)

and variance,

$\overline{m_{ij}m_{ji}}=\overline{m_{ii}^{2}}=\frac{m^{2}}{N}\,.$ (14)

Here we have defined $\displaystyle\overline{F}=\int F\prod_{i\neq
j}P(m_{ij})dm_{ij}\prod_{i}P(m_{ii})dm_{ii}$, and use an overline to denote
the disorder averaging throughout. In addition, we need to specify whether we
consider an annealed or quenched averaging. In the former we perform the
disorder averaging at the beginning, obtaining a new post-averaging action
while in the latter we compute correlation functions for mixed choices of
$m_{ij}$ and perform the averaging afterwards. In this article we will utilise
both approaches. Our computation of the spectral form factor will use annealed
averaging while the out-of-time-ordered correlator will use quenched
averaging. At large $N$, we expect these results to coincide Baldwin2020 , but
hope that the development of both approaches will prove useful when different
perturbations of (5) are considered.

##  3 Disorder averaging

Before presenting our findings for the SFF and OTOC for the model described by
(12), it will be useful to take note of some results that are relevant for the
computation of general correlation functions of disordered system.

###  3.1 Quenched averaging

We begin by collecting results that will be useful for carrying out a quenched
averaging prescription. In this case, all Feynman diagrams consist of
integrals of the two-point function (11) with some collection of $m_{ij}$
terms inserted as vertices. As such, the convolution of $n$ two-point
functions is given by

$\displaystyle\int_{0}^{\beta}d\tau_{1}\int_{0}^{\beta}d\tau_{2}\cdots\int_{0}^{\beta}d\tau_{n}G^{0}_{k_{1}k_{1}}(t_{1}-\tau_{1})G^{0}_{k_{2}k_{2}}(\tau_{1}-\tau_{2})\cdots
G^{0}_{k_{n}k_{n}}(\tau_{n}-t_{2})$ $\displaystyle=$ $\displaystyle
e^{i\frac{a}{\beta}(t_{1}-t_{2})}\sum_{p=0}^{n}\frac{(t_{1}-t_{2})^{n-p}\beta^{p}(-1)^{p}}{2^{p+1}p!(n-p)!}\left.\partial^{\,p}_{x}\left(\tanh
x+\mathrm{sgn}(t_{1}-t_{2})\right)\right|_{x\rightarrow\tanh^{-1}(-i\tan(\frac{a}{2}))}$

where $G^{0}_{ij}$ is defined in (11). The mass averaging can be carried out
for arbitrary collections of mass insertions. We can do the integrals
separately for the diagonal couplings, and pairwise separately for the off-
diagonal ones. Consequently,

$\displaystyle\overline{m_{ij}^{s}m_{ji}^{r}}$ $\displaystyle=$ $\displaystyle
r\frac{m^{2r}}{N^{r}}\Gamma(r)\delta_{s,r}\ \ \ i\neq j$
$\displaystyle\overline{m_{ii}^{2r}}$ $\displaystyle=$ $\displaystyle
2^{r}\frac{m^{2r}}{N^{r}}\frac{\Gamma(r+\frac{1}{2})}{\sqrt{\pi}}\ \ \ r\in
N_{0}$ (15) $\displaystyle\overline{m_{ii}^{2r+1}}$ $\displaystyle=$
$\displaystyle 0\ \ \ r\in N_{0}$

Since the two-point function $G^{0}_{ij}$ is proportional to $\delta_{ij}$,
this forces the indices of the mass insertions to order into a sequence.
Consequently, for all diagrams one has to perform averages of the form

$\sum_{k_{i}=1}^{N}\overline{m_{ik_{1}}m_{k_{1}k_{2}}\cdots m_{k_{s}j}}$ (16)

In order for this quantity to be non-zero, the number of $m_{ij}$ terms must
match the number of $m_{ji}$ terms for all $i,j$. The leading large $N$
contribution comes from the so-called “rainbow diagrams” Gross:2016kjj . These
are illustrated in Fig. (1)

Figure 1: The diagrams relevant for four mass insertions. The fermion
propagators are indicated by solid lines and are directed from
$\psi^{\dagger}$ to $\psi$. The dashed lines are related to the mass-
averaging, which require the masses to be paired up to give a non-zero result.
The top two diagrams are rainbow diagrams for which the dashed lines do not
intersect. Diagrams such as this represent the leading large $N$ contribution
since the smallest number of flavor indices are contracted.

and result in

$\sum_{k_{i}=1}^{N}\overline{m_{ik_{1}}m_{k_{1}k_{2}}\cdots
m_{k_{2l-1}j}}=\frac{1}{l+1}\left(\begin{array}[]{c}2l\\\
l\end{array}\right)m^{2l}\delta_{ij}+\mathcal{O}(\frac{1}{N})$ (17)

The combinatoric factor can be computed by solving a recursive equation, see
Fig. (2).

Figure 2: The number of rainbow diagrams for $2l$ mass insertions can be
determined by a recursive relation, displayed here in diagrammatic form. The
double dashed lines sum over all rainbow diagrams involving the number of
insertions in the square box.

###  3.2 Imaginary time two-point function

Remarkably, these results may be combined to give a closed form expression for
the fully dressed two-point correlation function. The diagrams we need to sum
correspond to

$\displaystyle\overline{\langle\psi_{i}^{\dagger}(\tau)\psi_{j}\rangle}\equiv\bar{G}_{ij}(\tau)$
(18) $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\int_{0}^{\beta}d\tau_{1}\cdots\int_{0}^{\beta}d\tau_{2l}\
G^{0}_{ii}(\tau-\tau_{1})G^{0}_{k_{1}k_{1}}(\tau_{1}-\tau_{2})\cdots
G^{0}_{k_{2l-1}k_{2l-1}}(\tau_{2l-1}-\tau_{2l})G^{0}_{jj}(\tau_{2l})\times$
$\displaystyle\left(\sum_{k_{i}=1}^{N}\overline{m_{ik_{1}}m_{k_{1}k_{2}}\cdots
m_{k_{2l-1}j}}\right)\,.$

By taking the large $N$ limit and making use of an integral representation of
the hypergeometric function we obtain

$\displaystyle\overline{G}_{i_{1}j_{2}}(\tau)$ $\displaystyle=$
$\displaystyle\delta_{j_{2}i_{1}}e^{\frac{ia}{\beta}\tau}\frac{1}{\pi}\int_{-1}^{1}dy\sqrt{1-y^{2}}e^{-2my\tau}\left(\tanh\left(-i\frac{a}{2}+my\beta\right)+\textnormal{sgn}(\tau)\right)$
(19) $\displaystyle=$
$\displaystyle\delta_{j_{2}i_{1}}\textnormal{sgn}(\tau)e^{i\frac{a}{\beta}\tau}\frac{2}{\pi}\int_{-1}^{1}dy\sqrt{1-y^{2}}\frac{e^{-2my\tau}}{1+e^{-2\textnormal{sgn}(\tau)my\beta}e^{i\textnormal{sgn}(\tau)a}}$
$\displaystyle=$
$\displaystyle\delta_{j_{2}i_{1}}\textnormal{sgn}(\tau)\frac{2}{\pi}\int_{-1}^{1}dy\sqrt{1-y^{2}}\frac{e^{-\left(2my+i\frac{a}{\beta}\right)\tau}}{1+e^{-\textnormal{sgn}(\tau)\left(2my+i\frac{a}{\beta}\right)\beta}}\,.$

Taking $0<\tau<\beta$ and transforming the $y$-coordinate, we obtain

$\displaystyle\overline{G}_{i_{1}j_{2}}(\tau)$ $\displaystyle=$
$\displaystyle\delta_{j_{2}i_{1}}e^{i\frac{a}{\beta}\tau}\int_{-2m}^{2m}dy\frac{\sqrt{4m^{2}-y^{2}}}{2\pi
m^{2}}\frac{e^{-y\tau}}{1+e^{-y\beta}e^{ia}}\,,$ (20)

which clearly reproduces the correct limit for $a\rightarrow 0$, (see Appendix
C of Gross:2016kjj ) and is invariant under large gauge transformation444Note
that the overall phase accounts for the transformation of the fermion fields
under the gauge transformation. On the other hand, when $m\beta>>0$, the
integral may be computed analytically. The leading contribution to this
integral is given by

$\displaystyle\overline{G}_{i_{1}j_{2}}(\tau)$ $\displaystyle=$
$\displaystyle\delta_{j_{2}i_{1}}\textnormal{sgn}(\tau)e^{i\frac{a}{\beta}\tau}\frac{2}{\pi}\int_{0}^{1}dy\sqrt{1-y^{2}}e^{-2my|\tau|}$
(21) $\displaystyle=$
$\displaystyle\frac{\delta_{j_{2}i_{1}}}{2m\tau}e^{i\frac{a}{\beta}\tau}\left(I_{1}(2m|\tau|)-L_{1}(2m|\tau|)\right)\,,$

where $I_{1}$ is the modified Bessel function and $L_{1}$ is the modified
Struve function and differs from the expression C.7 in Gross:2016kjj  only by
an overall phase. For general $m\beta$ we can write (20) as

$\bar{G}_{i_{1}j_{2}}(\tau)=\delta_{j_{2}i_{1}}\frac{2}{\pi\beta}\int_{-1}^{1}dy\sqrt{1-y^{2}}\sum_{n=-\infty}^{\infty}\frac{e^{i(2n-1)\frac{\pi}{\beta}\tau}}{\frac{(2n-1)\pi
i}{\beta}+(2my+i\frac{a}{\beta})}$ (22)

after making use of an identity in Dunne:1996yb . This expression can also be
obtained by performing the inverse discrete Fourier transform and setting
$\omega_{n}=\frac{\pi}{\beta}(2n-1)$ to give,

$\displaystyle\bar{G}_{j_{2}i_{1}}(\omega_{n})=\frac{2}{\pi}\delta_{j_{2},i_{1}}\int_{-1}^{1}dy\sqrt{1-y^{2}}\frac{-1}{i\omega_{n}-(2my+i\frac{a}{\beta})}$

Note that while the above expression for an individual mode is not invariant
under the large gauge transformation $a\rightarrow a+2\pi k$, the set of all
modes is.

###  3.3 Real-time two-point function

Given $\overline{G}_{j_{2}j_{1}}(\omega_{n})$, the real-time retarded Green’s
function is usually obtained by analytically continuing
$i\omega_{n}\rightarrow\omega+i\epsilon$. In addition, however, we need to
Wick rotate the gauge field $a\rightarrow-ia$ as we are transforming from real
to imaginary time. To disambiguate between the real- and imaginary time
expressions we also replace $\beta\rightarrow T$. This step is essential for
preserving the gauge invariance of our expressions. Then, taking the Fourier
transform and carrying out the integral over $y$, we obtain

$\displaystyle\bar{G}_{j_{2}i_{1}}(t)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\delta_{j_{2},i_{1}}\int\frac{d\omega}{2\pi}\int_{-1}^{1}dy\sqrt{1-y^{2}}\frac{-1}{\omega+i\epsilon-(2my+\frac{a}{T})}e^{-i\omega
t}$ (24) $\displaystyle=$ $\displaystyle
ie^{-i\frac{a}{T}t}\frac{J_{1}(2mt)}{mt}\theta(t)\,,$

where $J_{1}$ is a Bessel function of the first kind. This agrees with the
corresponding expression in Michel:2016kwn  up to an additional phase factor
in $t$, which does not affect the late-time decay of the expression
$\frac{J_{1}(2mt)}{mt}\sim t^{-\frac{3}{2}}$. Without simultaneously Wick
rotating the gauge field, we would instead obtain an overall (decaying)
exponential factor which breaks the gauge invariance. The answer (24) may also
be obtained by deforming the integration contour Michel:2016kwn . By first
performing the integral over $y$ followed by the analytic continuation and
Wick rotation we obtain

$\bar{G}_{j_{2}i_{1}}(\omega)=-\delta_{j_{2}i_{1}}\left(\omega-\frac{a}{T}\right)\frac{1-\sqrt{1-\frac{4m^{2}}{(\omega-\frac{a}{T}+i\epsilon)^{2}}}}{2m^{2}}$
(25)

The square root factor introduces a branch cut that is shifted just off the
real axis. For positive time we close the contour in the lower complex plane
and get a non-vanishing result since the branch cut is contained in the
contour. The contour may be deformed to one that hugs the branch cut, yielding
(24). The combination $\frac{J_{1}(2mt)}{mt}$ decays as $t^{-\frac{3}{2}}$ at
late times.

###  3.4 Annealed averaging

The choice of a Gaussian distribution for the random coupling makes it
particularly convenient to carry out an annealed averaging. Depending on the
quantity of interest, a number of replicas may be constructed, the random
couplings for the replicas identified and the Gaussian integral may be
performed exactly. The general integrals we will be interested in are simple
Gaussians,

$\displaystyle\int\prod_{i\neq
j}P(m_{ij})\,dm_{ij}\prod_{k}P(m_{kk})\,dm_{ii}\exp\Biggl{[}m_{ij}\left(\sum_{\alpha}w_{\alpha}\int
d\tau\psi_{i_{\alpha}}^{\dagger}\psi_{j_{\alpha}}\right)$ (26)
$\displaystyle+$ $\displaystyle m_{kk}\left(\sum_{\alpha}w_{\alpha}\int
d\tau\psi_{k_{\alpha}}^{\dagger}\psi_{k_{\alpha}}\right)\Biggr{]}\,,$

where the label $\alpha$ runs over replicas. The weight factor $w_{\alpha}$ is
relevant when computing real-time correlators and complex conjugates.

In Appendix A we compute the disorder average for the $2$-replica spectral
form factor, while in Appendix B we compute the partition function, involving
a single replica. At large $N$ we find the expected agreement between the two-
point function computed using either the annealed or quenched prescriptions.

##  4 Double commutator

The out-of-time-ordered correlator (or the closely-related double commutator
correlation function) is a widely used contemporary diagnostic of quantum
chaos Larkin1969QuasiclassicalMI ; Almheiri:2013hfa ; Shenker:2013pqa ;
Shenker:2013yza ; Roberts:2014isa ; Roberts:2014ifa ; Stanford:2015owe ;
Maldacena:2015waa . It will serve a similar purpose in our study of the model
(12). Specifically, the correlator

$C(t)=\theta(t)\langle\sqrt{\rho}\left\\{\psi^{\dagger}(t)_{i},\psi(0)_{j}\right\\}^{\dagger}\sqrt{\rho}\left\\{\psi^{\dagger}(t)_{i},\psi(0)_{j}\right\\}\rangle$
(27)

is expected to exhibit exponential growth at early times for chaotic systems
Maldacena:2015waa . The rate of exponential growth is characterised by the
spectrum of Lyaponov exponents. The splitting of the factor $\rho$ into the
two factors of $\sqrt{\rho}$ regulates some divergences and places the
insertions on opposite sides of the thermal circle Stanford:2015owe .

To compute $C(t)$ (in the quenched averaging prescription) we need to perform
the standard summing over all Feynman diagrams. As emphasised in section 3,
this involves the sum of convolutions of the undressed two-point function,
which follows from the Wick contractions. The structure of the Wick
contractions for $C(t)$ has a special feature due to the anti-commutators.
Specifically, as a result of the fermionic statistics, contractions involving
an anti-commutator of fields vanish unless the fields in the anti-commutator
are contracted with each other. At tree level this implies factorisation of
(27) into two retarded propagators, a feature which persists even beyond tree
level where perturbative two-body insertions may always be combined into an
anti-commutators of fields. This implies that, for our model, Wick
contractions can only occur between fields that are inserted on the same side
of the thermal circle, dramatically simplifying (27) to

$\displaystyle C(t)$ $\displaystyle=$
$\displaystyle\theta(t)G_{ij}(t)(G_{ij}(t))^{*}$ $\displaystyle
G_{ij}(t_{1}-t_{2})$ $\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\int
dt_{1}^{\prime}dt_{2}^{\prime}\cdots
dt_{n}^{\prime}G^{0}_{ii}(t_{1}-t_{1}^{\prime})G^{0}_{k_{1}k_{1}}(t_{1}^{\prime}-t_{2}^{\prime})\cdots
G^{0}_{jj}(t_{n}^{\prime}-t_{2})\times$ (28)
$\displaystyle\sum_{k_{i}=1}^{N}m_{ik_{1}}m_{k_{1}k_{2}}\cdots
m_{k_{n-1}j}\,,$

where we have used the fact that $G^{0}_{ij}\propto\delta_{ij}$. We did not,
however, impose any further conditions on $G^{0}_{ij}$ nor did we need to use
the disorder averaging of the mass terms. The result (28) only makes use of
the fact that we have perturbative two-body interactions and the special form
of the Wick contractions for anti-commutators of fields.

The results of section 3 can be used directly to simplify the above
expression. However, a comment of caution is in order. At first glance the
above expression looks like the absolute-value-squared of the dressed two-
point function however, unlike the Wick contractions, the mass averaging need
not involve mass insertions from the same side of the thermal circle.

To leading order in $1/N$ then, we can write

$\bar{C}(t)=\theta(t)\bar{G}_{ij}(t)(\bar{G}_{ij}(t))^{*}+\frac{\theta(t)}{N}\bar{F}_{ijji}(t)+O(\frac{1}{N^{2}})$
(29)

so that at leading order in $N$, the double commutator factorises into the
product of dressed two-point functions. The $1/N$ corrections come in two
forms. The first are of the kind in (17) which will further dress the two-
point function. The second come from corrections involving mass averaging on
opposite sides of the thermal circle. We will focus on the latter since these
give rise to ladder diagrams, the proliferation of which is associated with
the spectrum of growth exponents Stanford:2015owe .

We may regroup the terms in (28) to isolate the masses that are involved in
averages on opposite sides of the thermal circle. The sum over mass insertions
for masses averaged on the same side of the thermal give rise to dressed two-
point functions. This means that

$\displaystyle\frac{1}{N}F_{ijji}(t)$ $\displaystyle=$
$\displaystyle\sum_{r=1}^{\infty}\left(\int dt_{1}dt_{2}\cdots
dt_{n}\bar{G}_{ii}(t-t_{1})\bar{G}_{k_{1}k_{1}}(t_{1}-t_{2})\cdots\bar{G}_{jj}(t_{r-1})\right)\times$
(30) $\displaystyle\left(\int dt^{\prime}_{1}dt^{\prime}_{2}\cdots
dt^{\prime}_{n}\bar{G}_{ii}(t-t^{\prime}_{1})\bar{G}_{l_{1}l_{1}}(t^{\prime}_{1}-t^{\prime}_{2})\cdots\bar{G}_{jj}(t^{\prime}_{r-1})\right)^{*}\times$
$\displaystyle\sum_{k_{i},l_{i}=1}^{N}(m_{ik_{1}}m_{k_{1}k_{2}}\cdots
m_{k_{r-1}j})(m_{il_{1}}m_{l_{1}l_{2}}\cdots m_{l_{r-1}j})^{*}\,,$

where we have isolated the $r$ masses involved in averaging on opposite sides
of the thermal circle. Taking the Fourier transform and isolating the leading
large $N$ contribution, we find

$\displaystyle\bar{F}_{ijji}(t)$ $\displaystyle=$ $\displaystyle
2\delta_{ij}\sum_{r=1}^{\infty}\int
d\omega_{1}d\omega_{2}e^{-i(\omega_{1}-\omega_{2})t}\left(\bar{G}(\omega_{1})(\bar{G}(\omega_{2}))^{*}\right)^{r+1}m^{2r}$
$\displaystyle=$ $\displaystyle 2\delta_{ij}\int
d\omega_{1}d\omega_{2}e^{-i(\omega_{1}-\omega_{2})t}\frac{m^{2}\left(\bar{G}(\omega_{1})(\bar{G}(\omega_{2}))^{*}\right)^{2}}{1-m^{2}\bar{G}(\omega_{1})(\bar{G}(\omega_{2}))^{*}}$
$\displaystyle=$ $\displaystyle 2\delta_{ij}\int
d\omega_{1}d\omega_{2}e^{-i(\omega_{1}-\omega_{2})t}\left(-\bar{G}(\omega_{1})(\bar{G}(\omega_{2}))^{*}+\frac{1}{1-m^{2}\bar{G}(\omega_{1})(\bar{G}(\omega_{2}))^{*}}\right)\,,$

in agreement with the results of Michel:2016kwn ; Gross:2016kjj . The Green’s
function in the above expression is the Fourier transform of the real-time
Green’s function (25).

By shifting the integration variables
$\omega_{j}\rightarrow\omega_{j}+\frac{a}{T}$, we obtain an integral of the
type considered in Gross:2016kjj . The phase factors from the two-point
functions cancel out and thus the gauge field has no effect on the double
commutator. As pointed out in Polchinski:2016xgd ; Gross:2016kjj , the
resulting expression does not give rise to exponential growth of the OTOC so
that we conclude that our model is not chaotic. This is not surprising since
neither is the underlying SYK2 model.

##  5 Spectral form factor

A hallmark of chaotic systems is that they exhibit spectral rigidity, in the
sense that the eigenvalues of the Hamiltonian repel each other i.e. they
forget as much information as possible, modulo what they are compelled to
remember by their symmetries. This heuristic picture is succinctly encoded in
the so-called Bohigas-Giannoni-Schmit (BGS) conjecture Bohigas:1983er  which
holds that any quantum chaotic system should exhibit the spectral statistics
of a random matrix ensemble consistent with the time-reversal symmetries of
the model. One very effective diagnostic tool for this random matrix theory
(RMT) behaviour is the spectral form factor (SFF), the Fourier transform of
the 2-point correlator of eigenvalues of the Hamiltonian. At the level of the
SFF, a signature of the level-repulsion indicating RMT behaviour, is the
emergence of a linear scaling with time - the so-called linear ramp - at
intermediate times Cotler:2016fpe ; Saad:2018bqo ; Liu:2018hlr ;
Cotler:2017jue .

The SFF also encodes specific non-chaotic features of disorder-averaged
integrable many-body systems, such as the (usual Majorana) SYK2 model
considered in Winer:2020mdc . There, the authors show that the 2-replica
action possesses an unusual, large $SU(2)$ symmetry, which gives rise to an
exponential, as opposed to a linear, ramp in the SFF. They further conjecture
that this exponential ramp should be a feature of all non-interacting
disordered systems, and attribute this behaviour of the SFF to the presence of
a high-dimensional manifold of saddle points resulting from the symmetry
group. The gauged SYK2 model which we consider in this article provides
another setting in which to test these claims.

We note in passing that the complex SYK2 model with chemical potential has
been considered previously in Liao:2020lac , largely using different methods
to the path integral approach utilised here as well as in Winer:2020mdc . In
real time - in which the SFF is calculated - our model is superficially
equivalent to theirs with the important distinction that the gauge field
transforms differently to a chemical potential and acquires a temperature
dependence in constant gauge (10).

Practically speaking, the spectral form factor can be computed by analytically
continuing the partition function $Z=Z(\beta)$ to $Z(\beta+iT)$. The SFF is
subsequently defined as

$g(\beta;T)\equiv\frac{\overline{Z(\beta+iT)Z^{*}(\beta+iT)}}{\overline{Z(\beta)}^{2}}\,.$
(31)

Its disconnected and connected components are defined respectively by

$\displaystyle g_{d}(\beta;T)$ $\displaystyle=$
$\displaystyle\frac{\overline{Z(\beta+iT)}\
\overline{Z^{*}(\beta+iT)}}{\overline{Z(\beta)}^{2}},$ (32) $\displaystyle
g_{c}(\beta;T)$ $\displaystyle=$ $\displaystyle g(\beta;T)-g_{d}(\beta;T)\,.$
(33)

Note that according to these definitions we are treating the disorder as
annealed, and not quenched. This allows us to perform the disorder average
over the denominator and numerator separately.

We will work primarily in the infinite temperature limit $g(\beta=0;T)\equiv
g(T)$. To calculate $Z(iT)$, we simply Wick rotate from Euclidean to real
time, under which $\tau\rightarrow it$, $\ A(\tau)\rightarrow-iA(t)$ and
integration runs from 0 to $T$. Working in the constant gauge
$A(\tau)=a/\beta$ which in real time becomes $A(t)=a/T$, our action,

$\displaystyle\tilde{I}=iS=-\int_{0}^{T}dt\left[\sum_{i=1}^{N}\psi_{i}^{\dagger}\left(\frac{d}{dt}-i\frac{a}{T}\right)\psi_{i}+i\sum_{i,j=1}^{N}m_{ij}\psi_{i}^{\dagger}\psi_{j}\right].$
(34)

As is usual in the path integral computation of the SFF, we consider the
physics of fermions in periodic real time, with period $T$ and with the
fermions satisfying anti-periodic boundary conditions.

Averaging over the disorder (see Appendix A for more details), yields the
following exact expression in terms of bilocal, collective fields $G$ and
$\Sigma$,

$\displaystyle g(T)$ $\displaystyle\equiv$
$\displaystyle\overline{Z(iT)Z(iT)^{*}}=\prod_{\alpha,\beta}\int\mathscr{D}\Sigma_{\alpha\beta}\mathscr{D}G_{\alpha\beta}e^{NI(G,\Sigma)}\,,$
(35)

where,

$\displaystyle I(G,\Sigma)$ $\displaystyle=$
$\displaystyle\log\det\left[-\delta(t-t^{\prime})\left(\delta_{\alpha\beta}\frac{d}{dt}-i\frac{a}{T}\sigma^{z}_{\alpha\beta}\right)+\Sigma_{\beta\alpha}(t,t^{\prime})\right]$
$\displaystyle-\int
dtdt^{\prime}\sum_{\alpha,\beta}\left[\Sigma_{\beta\alpha}(t^{\prime},t)G_{\alpha\beta}(t,t^{\prime})-(-1)^{\alpha+\beta}\frac{m^{2}}{2}G_{\alpha\beta}(t,t^{\prime})G_{\beta\alpha}(t^{\prime},t)\right],$

and $\sigma^{z}$ is the usual Pauli $z$-matrix. Greek indices run over the two
replicas. In order to tackle this path integral we will assume that time-
translation invariant configurations of $G$ and $\Sigma$ dominate
Winer:2020mdc . This allows us to take a discrete Fourier transform

$\displaystyle G_{\alpha\beta}(t)$ $\displaystyle=\frac{1}{T}\sum_{n\
\text{odd}}G_{\alpha\beta}(\omega_{n})e^{-i\omega_{n}t},\quad\Sigma_{\alpha\beta}(t)=\frac{1}{T}\sum_{n\
\text{odd}}\Sigma_{\alpha\beta}(\omega_{n})e^{-i\omega_{n}t},\qquad\omega_{n}=\frac{n\pi}{T},$
(37)

and ensures that variables with different Matsubara frequencies decouple. In
terms of these variables the path integral becomes

$g(T)=\prod_{n\ \textnormal{odd}}\prod_{\alpha,\beta}\left(\int
d\Sigma_{\alpha\beta}(\omega_{n})dG_{\alpha\beta}(\omega_{n})e^{NI_{n}(G(\omega_{n}),\Sigma(\omega_{n}))}\right)$
(38)

where now

$\displaystyle I_{n}(G,\Sigma)$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\text{Tr}\log\left[i\omega_{n}\delta_{\alpha\beta}+i\frac{a}{T}\sigma^{z}_{\alpha\beta}+\Sigma_{\alpha\beta}(\omega_{n})\right]$
(39)
$\displaystyle\qquad\quad-\Sigma_{\alpha\beta}(\omega_{n})G_{\beta\alpha}(\omega_{n})+(-1)^{\alpha+\beta}\frac{m^{2}}{2}G_{\alpha\beta}(\omega_{n})G_{\beta\alpha}(\omega_{n})\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\Bigg{\\{}\text{Tr}\log\left[\left(i\omega_{n}\sigma^{z}_{\alpha\gamma}+i\frac{a}{T}\delta_{\alpha\gamma}+\tilde{\Sigma}_{\alpha\gamma}(\omega_{n})\right)\sigma^{z}_{\gamma\beta}\right]$
$\displaystyle\qquad\quad-(-1)^{\alpha+\beta}\tilde{\Sigma}_{\alpha\beta}(\omega_{n})\tilde{G}_{\beta\alpha}(\omega_{n})+\frac{m^{2}}{2}\tilde{G}_{\alpha\beta}(\omega_{n})\tilde{G}_{\beta\alpha}(\omega_{n})\Bigg{\\}}\,,$

In this expression repeated indices are summed over, the trace is taken over
the replica indices, and in the second equality we have defined
$\tilde{G}=G\sigma^{z}$ and $\tilde{\Sigma}=\Sigma\sigma^{z}$. This may be
rewritten as

$\displaystyle I(G,\Sigma)$ $\displaystyle=$
$\displaystyle\text{Tr}\log\left[\left(i\sigma^{z}_{\alpha\gamma}\delta^{nn^{\prime}}\omega_{n}+i\frac{a}{T}\delta^{nn^{\prime}}_{\alpha\gamma}+\tilde{\Sigma}_{\alpha\gamma}^{nn^{\prime}}\right)\sigma^{z}_{\gamma\beta}\right]-\text{Tr}\left[\left((-1)^{\alpha+\gamma}\tilde{\Sigma}_{\alpha\gamma}^{np}-\frac{m^{2}}{2}\tilde{G}_{\alpha\gamma}^{np}\right)\tilde{G}_{\gamma\beta}^{pn^{\prime}}\right],$

with the traces now taken over both replica and (odd $n$) frequency spaces. In
addition, we have introduced the frequency space indices $n,n^{\prime},p$ and
defined the matrix
$\delta^{nn^{\prime}}_{\alpha\beta}\equiv\delta^{nn^{\prime}}\delta_{\alpha\beta}$.

###  5.1 Saddle point approximation

The overall factor of $N$ appearing in (35) makes the SFF well suited to a
saddle point analysis in the large $N$ limit. From (5) we find the
$(\tilde{G},\tilde{\Sigma})$ equations of motion,

$\displaystyle\frac{\delta
I}{\delta\tilde{\Sigma}_{\alpha\beta}^{nn^{\prime}}}=0\quad$
$\displaystyle\Rightarrow\quad\sigma^{z}_{\alpha\rho}\left[i\sigma^{z}_{\rho\gamma}\delta^{np}\omega_{n}+i\frac{a}{T}\delta_{\rho\gamma}^{np}+\tilde{\Sigma}_{\rho\gamma}^{np}\right]\sigma^{z}_{\gamma\nu}(-1)^{\nu+\beta}\tilde{G}_{\nu\beta}^{pn^{\prime}}=\delta_{\alpha\beta}^{nn^{\prime}},$
(41) $\displaystyle\frac{\delta
I}{\delta\tilde{G}_{\alpha\beta}^{nn^{\prime}}}=0\quad$
$\displaystyle\Rightarrow\quad\tilde{\Sigma}_{\alpha\beta}^{nn^{\prime}}=(-1)^{\alpha+\beta}m^{2}\tilde{G}_{\alpha\beta}^{nn^{\prime}},$
(42)

which can be written as a single equation of motion for $\tilde{\Sigma}$,

$\displaystyle\left(i\sigma^{z}\omega+i\frac{a}{T}+\tilde{\Sigma}\right)\tilde{\Sigma}=m^{2},$
(43)

suppressing the frequency and replica indices. Similarly, we can rewrite the
on-shell action (5),

$\displaystyle I(\tilde{\Sigma})$ $\displaystyle=$
$\displaystyle\text{Tr}\log\left[(i\sigma^{z}\omega+i\frac{a}{T}+\tilde{\Sigma})\sigma^{z}\right]-\frac{1}{2}\text{Tr}\tilde{\Sigma}^{2}\,,$
(44)

Note that without the $\sigma^{z}\omega$ term, the action (44) is invariant
under general unitary transformations $\tilde{\Sigma}\rightarrow
U^{\dagger}\tilde{\Sigma}U$. In the presence of the $\sigma^{z}\omega$ term,
this symmetry is explicitly broken to transformations,
$\tilde{\Sigma}\rightarrow\tilde{U}^{\dagger}\tilde{\Sigma}\tilde{U}$ for
unitaries $\tilde{U}$ satisfying
$\tilde{U}\sigma^{z}\omega\tilde{U}^{\dagger}=\sigma^{z}\omega$, where we
recall that $\sigma^{z}\omega$ is actually a matrix in replica and frequency
space given with matrix elements
$(\sigma^{z}\omega)_{\alpha\beta}^{nn^{\prime}}=\sigma^{z}_{\alpha\beta}\delta^{nn^{\prime}}\omega_{n}$.
The breaking of this symmetry and the associated Goldstone zero modes play a
crucial role in giving rise to the exponential ramp exhibited by the SFF in
this system Liao:2021ofk .

Solutions to (43) are diagonal in both replica and frequency indices, and are
given by

$\displaystyle\tilde{\Sigma}_{\alpha\beta}^{nn^{\prime}}=\frac{1}{2}\left[-i\left(\xi_{\alpha}\omega_{n}+\frac{a}{T}\right)+c_{\alpha}^{n}\sqrt{4m^{2}-\left(\xi_{\alpha}\omega_{n}+\frac{a}{T}\right)^{2}}\right]\delta^{nn^{\prime}}_{\alpha\beta},$
(45)

where $\xi_{1}=1,\ \xi_{2}=-1$ and the constants $c_{\alpha}^{n}=\pm 1$ can be
chosen independently for each replica and frequency index. Imposing the
condition that for large frequencies our solution must vanish Altland:2000fj
we obtain the following solutions

$\displaystyle\mathcal{S}_{\alpha\beta}^{nn^{\prime}}=\begin{cases}\frac{1}{2}\left[-i(\xi_{\alpha}x_{n}+b)+c_{\alpha}^{n}\sqrt{4-(\xi_{\alpha}x_{n}+b)^{2}}\right]\delta^{nn^{\prime}}_{\alpha\beta},&|\xi_{\alpha}x_{n}+b|\leq
2,\\\
\frac{1}{2}\left[-i(\xi_{\alpha}x_{n}+b)+i\text{sgn}(\xi_{\alpha}x_{n}+b)\sqrt{(\xi_{\alpha}x_{n}+b)^{2}-4}\right]\delta^{nn^{\prime}}_{\alpha\beta},&|\xi_{\alpha}x_{n}+b|>2.\end{cases}$
(46)

in terms of the dimensionless quantities
$\mathcal{S}\equiv\frac{1}{m}\tilde{\Sigma}$,
$x_{n}\equiv\frac{\omega_{n}}{m}=\frac{n\pi}{mT}$ and $b\equiv\frac{a}{mT}$.
Evidently, there are many saddles, one for each value of the constants
$c_{\alpha}^{n}$ for frequencies in the range $|\xi_{\alpha}x_{n}+b|\leq 2$.
It remains to check which of these contribute to the SFF at large $N$.

In what follows, we work with a form of on-shell action (44) that has been
regularized by subtracting off the action for the free theory $\log(i\omega)$,
to obtain

$\displaystyle I(\mathcal{S})$ $\displaystyle=$
$\displaystyle\text{Tr}\log\left[\left(\sigma^{z}+\frac{b}{x}-\frac{i}{x}\mathcal{S}\right)\sigma^{z}\right]-\frac{1}{2}\text{Tr}\mathcal{S}^{2}\,,$
(47)

where, from here on, we suppress the frequency dependence of $x_{n}$. To see
which saddle points solutions contribute to the SFF, we must analyse
fluctuations $\delta\mathcal{S}^{nn^{\prime}}_{\alpha\beta}$ about the various
saddles $\mathcal{S}^{nn^{\prime}}_{\alpha\beta}$. To this end, consider the
variation of the on-shell action (47) under
$\mathcal{S}\rightarrow\mathcal{S}+\delta\mathcal{S}$. After simplifying the
expression by using the equation of motion (43) we find that

$\displaystyle\delta I(\mathcal{S})$
$\displaystyle=-\frac{1}{2}\sum_{\alpha,\beta,n,n^{\prime}}M^{nn^{\prime}}_{\alpha\beta}\delta\mathcal{S}^{nn^{\prime}}_{\alpha\beta}\delta\mathcal{S}^{n^{\prime}n}_{\beta\alpha},$
(48)

with the mass matrix

$\displaystyle M^{nn^{\prime}}_{\alpha\beta}\equiv
1+\mathcal{S}^{nn}_{\alpha\alpha}\mathcal{S}^{n^{\prime}n^{\prime}}_{\beta\beta}.$
(49)

Note that here $\mathcal{S}^{nn}_{\alpha\alpha}$ denotes a single diagonal
component of the matrix $\mathcal{S}$ and not a trace (i.e. no summation over
repeated indices is implied in the above expression). Expressed in the above
terms, and normalising by dividing by $\overline{Z(0)}^{2}$, in the saddle
point approximation,

$g(T)=\mathcal{N}\sum_{\textnormal{saddles}}e^{NI(\mathcal{S}_{cl})}\prod_{\alpha,\beta}\int\left(\prod_{n,n^{\prime}\
\textnormal{odd}}\mathscr{D}\delta\mathcal{S}_{\alpha\beta}^{nn^{\prime}}\right)e^{N\delta
I(\mathcal{S}_{cl})}\,.$ (50)

###  5.2 The contributing saddles

For each frequency there are four possibilities, dictated by the values of
$|\xi_{1}x_{n}+b|$ and $|\xi_{2}x_{n}+b|$. If both are less (greater) than $2$
we refer to this region as $A$ ($D$). Regions $B$ and $C$ cover the other two
possibilities where either $|\xi_{1}x_{n}+b|\leq 2$ or $|\xi_{2}x_{n}+b|\leq
2$ respectively.

In Appendix C we verify that the mass matrices always have a non-negative real
part in all four regions, ensuring that all saddles do contribute. Moreover,
the real part of the mass matrices are all strictly positive except when
$\alpha\neq\beta$, $n^{\prime}=-n$ inside region $A$, where the vanishing mass
matrix signals the presence of zero modes which dominate the contribution to
the fluctuations. The region $A$ is, however, only present if
$|b|=|\frac{a}{mT}|<2$. For a given gauge field, this sets a natural scale for
the disorder strength necessary for zero modes to appear.

Furthermore, note that the mass matrices satisfy the following complex
conjugation relation

$\left(M^{nn^{\prime}}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\pm
1}\right)^{*}=M^{nn^{\prime}}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\mp
1}\,,$ (51)

in each of the four regions, a fact that will be important in demonstrating
that we obtain a real expression for the SFF. Using (47) we define
$I_{J}(x)\equiv I(\mathcal{S}^{(J)})$ where $J=A,B,C,D$, so the classical
action takes the form,

$\displaystyle I(S_{cl})=$ $\displaystyle
2\Bigg{[}\sum_{\begin{subarray}{c}n\pi>0\\\ n\
\text{odd}\end{subarray}}^{2mT-|a|}I_{A}^{(c,\tilde{c})}(x_{n})+\frac{1}{2}\sum_{\begin{subarray}{c}n\pi>2mT-|a|\\\
n\
\text{odd}\end{subarray}}^{2mT+|a|}\left(I_{C}^{(c)}(\text{sgn}(a)x)+I_{C}^{(\tilde{c})}(\text{sgn}(a)x)\right)+\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\pi>2mT+|a|\\\
n\ \text{odd}\end{subarray}}^{\infty}\\!\\!\\!\\!\\!\\!I_{D}(x_{n})\Bigg{]},$
(52)

for strong disorder $|a|/T<2m$, while for weak disorder $|a|/T>2m$,

$\displaystyle I(S_{cl})=$ $\displaystyle
2\Bigg{[}\sum_{\begin{subarray}{c}n\pi>0\\\ n\
\text{odd}\end{subarray}}^{|a|-2mT}I_{D}(x_{n})+\frac{1}{2}\sum_{\begin{subarray}{c}n\pi>|a|-2mT\\\
n\
\text{odd}\end{subarray}}^{|a|+2mT}\left(I_{C}^{(c)}(\text{sgn}(a)x)+I_{C}^{(\tilde{c})}(\text{sgn}(a)x)\right)+\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\pi>2mT+|a|\\\
n\ \text{odd}\end{subarray}}^{\infty}\\!\\!\\!\\!\\!\\!I_{D}(x_{n})\Bigg{]}.$
(53)

where the sums run over positive values of $x$ only and we have used the
symmetries of $I_{D}(x)=I_{D}(-x)$, $I^{(c)}_{B}(x)=I^{(c)}_{C}(-x)$ and
$I^{(c,\tilde{c})}_{A}(x)=I^{(\tilde{c},c)}_{A}(-x)$. Note that the factors of
$I_{A}$, $I_{B}$ and $I_{C}$ all depend on the values $c\equiv
c_{1}^{n},\tilde{c}\equiv c_{2}^{n}$ labelling the various saddles while
$I_{D}$ is the same for all saddle points. The sum over saddle points then
involves a sum over all allowed choices of $c_{\alpha}^{n}$. To proceed we
will evaluate these expressions in the early- and late-time limits.

###  5.3 SFF for $mT\ll 1$

In the early-time limit, all the contributions to the Matsubara sums in
(52),(53) come from $I_{D}$, unless there is an odd integer $n$ within $2mT$
of $|a|/\pi$.555In this case there is a single-frequency contribution from
$I_{C}$ constituting four saddles points, see Appendix E. When $a$ is not too
close to being an odd multiple of $\pi$ we may simply neglect this
contribution. Since $I_{D}$ does not involve the constants $c_{\alpha}^{n}$
there is a single saddle in this limit and we find that

$\displaystyle e^{NI(\mathcal{S}_{cl})}$ $\displaystyle=$
$\displaystyle\exp\Bigg{\\{}2N\sum_{n>0,\text{odd}}I_{D}(x_{n})\Bigg{\\}}$
(54) $\displaystyle\approx$
$\displaystyle\prod_{n>0,\text{odd}}\left[\frac{(n\pi-a)(n\pi+a)}{n^{2}\pi^{2}}\right]^{2N}\exp\left\\{-Nm^{2}T^{2}\sum_{n>0,\text{odd}}\frac{2\left(n^{2}\pi^{2}+a^{2}\right)}{(n\pi-a)^{2}(n\pi+a)^{2}}\right\\}$
$\displaystyle=$
$\displaystyle\cos\left(\frac{a}{2}\right)^{2N}\exp\Bigg{\\{}-\frac{Nm^{2}T^{2}}{4}\sec^{2}\left(\frac{a}{2}\right)\Bigg{\\}},$

where, in the second equality, we have neglected $\mathcal{O}(m^{4}T^{4})$
terms in the exponent. The fluctuation integral does not involve zero modes
and we may approximate it as a constant at early times. After normalising by
dividing by $g(0)=|Z(0)|^{2}$, we obtain the following expression for the
spectral form factor

$\displaystyle g(T)=$
$\displaystyle\exp\left\\{-\frac{Nm^{2}T^{2}}{4}\sec^{2}\left(\frac{a}{2}\right)\right\\}.$
(55)

In the limit $a\rightarrow 0$ the argument of the exponent reduces to
$-Nm^{2}T^{2}/4$, in agreement with the Majorana result in Winer:2020mdc
666Note that our result is a factor of 2 larger, arising from the fact that
complex fermions have twice as many degrees of freedom as Majorana fermions..
We note also, in passing, that (55) is periodic in $a$ with period $2\pi$
ensuring that the early-time SFF is indeed gauge invariant under large gauge
transformations i.e. transformations for which $a\rightarrow a+2\pi k$,
$k\in\textbf{Z}$.

The exponent of (55) blows up to negative infinity as $a\to\pi$. However,
recall that our result is valid only up to values of $mT$ for which $a$ is a
distance $2mT$ away from $\pi$. We refer the reader to Appendix D where we
plot the analytic results up to $\mathcal{O}(m^{8}T^{8})$ for various values
of $a$ and Appendix E for a more careful treatment in the region $\pi-a<2mT\ll
1$.

###  5.4 SFF for $mT\gg 1$

At late times (and assuming strong disorder), we use the expression (52) for
the SFF and continue $x_{n}$ to a continuous variable so that the Matsubara
sums are approximated by integrals. This gives

$\displaystyle S_{cl}$ $\displaystyle=$
$\displaystyle\frac{mTN}{\pi}\Bigg{[}\int_{0}^{2-|b|}dx\,I_{A}^{(c,\tilde{c})}(x)+\frac{1}{2}\int_{2-|b|}^{2+|b|}dx\left(I_{C}^{(c)}(\text{sgn}(b)x)+I_{C}^{(\tilde{c})}(\text{sgn}(b)x)\right)$
(56) $\displaystyle+$
$\displaystyle\int_{2+|b|}^{\infty}dx\,I_{D}(x)\Bigg{]}\,.$

Summing over saddles labelled by all possible configurations of
$(c,\tilde{c})$, results in a vanishing SFF (see Appendix F). To see a
nonvanishing result we need to consider fluctuations about these saddles. Of
particular interest to us in the expression (56) is the inclusion of the
region $A$ where zero modes are possible. As we will see, the effect of these
zero modes is to give rise to an exponential ramp at late times. Moreover, the
presence of the gauge field $a$ sets a natural scale
$\lambda\equiv\frac{|a|}{2m}$ for the disorder in order for zero modes to be
supported. To elaborate, let’s rewrite the expression for the SFF as,

$\displaystyle g(T)$
$\displaystyle=\frac{\int\mathscr{D}\mathcal{S}e^{NI(\mathcal{S})}}{\int\mathscr{D}\mathcal{S}e^{-\frac{N}{2}\text{Tr}\mathcal{S}^{2}}},$
(57)

where we now normalise by dividing by the free theory Winer:2020mdc ;
Liao:2020lac

$\displaystyle\int\mathscr{D}\mathcal{S}e^{-\frac{N}{2}\text{Tr}\mathcal{S}^{2}}=\prod_{\alpha,\beta,n,n^{\prime}}\sqrt{\frac{2\pi}{N}}.$
(58)

We can then write

$\displaystyle g(T)=\sum_{c,\tilde{c}}e^{NI(\mathcal{S})}F(T),$ (59)

using the large $N$ saddle point approximation up to quadratic fluctuations
and summing over all saddles labelled by $c\equiv
c_{\alpha}^{n},\tilde{c}=c_{\beta}^{n^{\prime}}$. The quadratic fluctuation
path integral, containing the normalisation, is given by

$\displaystyle F(T)$
$\displaystyle=\prod_{\alpha,\beta,n,n^{\prime}}\sqrt{\frac{N}{2\pi}}\int\mathscr{D}\delta\mathcal{S}_{\alpha\beta}^{nn^{\prime}}\exp\left\\{-\frac{N}{2}M_{\alpha\beta}^{nn^{\prime}}\delta\mathcal{S}_{\alpha\beta}^{nn^{\prime}}\delta\mathcal{S}_{\beta\alpha}^{n^{\prime}n}\right\\}.$
(60)

using (48). Note that the mass matrix is dependent on the particular saddle
point solution considered. We have already verified that
$M_{\alpha\beta}^{nn^{\prime}}\geq 0$ for all saddles and frequency regions.
From the cases where $M_{\alpha\beta}^{nn^{\prime}}>0$, we expect a finite
contribution to the path integral, provided that the products over frequencies
converge. We also have zero modes in region $A$ where
$M_{\alpha\beta}^{nn^{\prime}}=0$ for some saddles.

Let us first consider massive modes and their fluctuations. For these modes
$M_{\alpha\beta}^{nn^{\prime}}>0$ and the quadratic fluctuation integral is
easily carried out in each case to obtain

$\displaystyle
F(T)=\prod_{\alpha,\beta,n,n^{\prime}}\sqrt{\frac{N}{2\pi}}\sqrt{\frac{2\pi}{NM_{\alpha\beta}^{nn^{\prime}}}}=\exp\left\\{-\frac{1}{2}\sum_{\alpha,\beta,n,n^{\prime}}\log
M_{\alpha\beta}^{nn^{\prime}}\right\\}\,.$ (61)

Note that the contribution from these massive fluctuations is independent of
$N$. While we will not evaluate it explicitly, we will verify that it is
convergent and results in a real SFF.

First, let’s prove convergence. The regions $A$, $B$ and $C$ are finite, so
that the contribution from the massive modes from this regions is clearly
finite. In the infinite region $D$, however, we must entertain the possibility
of UV divergences. For very large frequencies $x_{n},x_{n^{\prime}}\gg 1$, the
mass matrix given explicitly in (C) can be approximated by

$\displaystyle M_{\alpha\beta}^{nn^{\prime}}$ $\displaystyle\approx
1-\frac{4}{\xi_{\alpha}\xi_{\beta}x_{n}x_{n^{\prime}}}\quad\Rightarrow\quad\sum_{\alpha,\beta,n,n^{\prime}}\log
M_{\alpha\beta}^{nn^{\prime}}\sim\sum_{\alpha,\beta,n,n^{\prime}}\frac{1}{\xi_{\alpha}\xi_{\beta}x_{n}x_{n^{\prime}}}.$
(62)

The sum is not absolutely convergent, but by imposing any finite frequency
cutoff we can ensure that cancellations between either $n^{\prime}=-n$ or
$\alpha\neq\beta$ are sufficient to deal with any UV divergences Winer:2020mdc
.

Second, the realness of the SFF follows from the complex conjugation relations
of the mass matrices (51), and that the classical action (56) exhibits the
same complex conjugate relations. Consequently, after summing over the
saddles, one obtains a real answer.

We now turn our attention to the zero modes. As discussed, these occur in the
region $A$ where mass matrices with $\alpha\neq\beta$ and $n^{\prime}=-n$
vanish when we are considering a saddle with $c_{1}^{n}=-c_{2}^{-n}$. Due to
the vanishing of the mass matrix we need to compute the contribution of the
zero modes differently. The zero modes are Goldstone modes that have their
origin in the explicit breaking of the symmetry of the equations of motion
(43) under general unitary transformations $\tilde{\Sigma}\rightarrow
U^{\dagger}\tilde{\Sigma}U$ by the $\sigma^{z}\omega$ term. The full set of
zero modes are obtained by performing unitary transformations with the
residual symmetry and summing over solutions that are not related in this way.
Zero modes that are related by such unitary transformations are degenerate. To
be explicit, for each non-degenerate zero mode pair, from the normalisation
with respect to the free contribution, we obtain a factor of $\frac{N}{2\pi}$.
This should be multiplied by the volume of the factor group representing the
residual symmetry Kamenev:1999zz . This volume does not depend on $N$. Since
the massive modes also give a contribution that is $N$-independent, the
leading large $N$ SFF is obtained by considering the saddle points that
maximize the number of non-degenerate zero modes. The classical action for
these saddle points vanishes and the total contribution to the SFF for such
saddles is given by

$\displaystyle
2N^{2\left(2mT-|a|\right)/\pi}\left(\prod_{|x_{n}|<2-\frac{|a|}{mT}}\frac{\mathcal{V}^{n}_{12}\mathcal{V}_{21}^{n}}{2\pi}\right)\exp\left\\{-\frac{1}{2}\sum_{\alpha,\beta,n,n^{\prime}}^{{}^{\prime}}\log
M_{\alpha\beta}^{nn^{\prime}}\right\\},$ (63)

where the summation $\sum_{\alpha,\beta,n,n^{\prime}}^{{}^{\prime}}$ excludes
the zero modes. The first factor above is the leading large-$N$ contribution
from maximising the number of non-degenerate zero modes; the second is the
product of the $N$-independent factor group volumes and the third is the
contribution from the massive modes for the saddle. Note that, though this
expression for the ramp is not gauge-invariant under large gauge
transformations $a\rightarrow a+2\pi k$, we expect gauge-invariance to be
restored when higher order fluctuations are taken into account.

Up to quadratic order in fluctuations, all $N$ dependence is contained in the
prefactor $N^{2\left(2mT-|a|\right)/\pi}$, from which we observe the
exponential ramp typical of noninteracting disordered systems Winer:2020mdc .
Note that both $mT$ and $\frac{|a|}{T}T$ are important and we require strong
disorder in order for the ramp to be present. We expect the contribution from
the remaining modes contained in the exponential term to become important at
times $mT\gg N$, where the soft modes combine with the zero mode fluctuations
to produce the late time plateau chacteristic of the SFF Winer:2020mdc . In
addition, we expect that the coefficient in the $N$ dependent ramp may receive
corrections from higher order fluctuations than the quadratic ones we have
considered here Liao:2020lac . Note also that we have neglected saddle point
solutions which break time translation invariance, which may slightly alter
the results obtained.

For $mT\ll 1$ on the other hand, we had no contribution from region $A$.
Consequently, there are no zero mode fluctuations to consider. The remaining
massive mode fluctuations are subleading in $N$ and so may be neglected. Thus
our results for the $mT\ll 1$ SFF slope in (55) and (111) are in harmony with
the late time results.

##  6 Discussion

In this paper, we have studied the properties of random mass fermions in one
dimension with twisted boundary conditions. The introduction of random masses
(for the fermions) furnishes a promising toy model to study aspects of an
$AdS_{2}$/gauged-SYK correspondence. This is due to two main features that
this model inherits: i. the random $q$-body interaction characteristic of SYKq
models and ii. the gauge field term that has been shown Dunne:1996yb  to endow
the $0+1$ dimensional model with features of $2+1$-dimensional Chern-Simons
theory, as well as the solvability of both. As discussed in the introduction,
our original goal was to construct a gauged SYKq model with an eye on the
canonical $q=4$ case largely because of its known holographic connection to JT
gravity. However as we hope to have convinced any reader to have made it to
this point, the $q=2$ (quadratic) model is also rather rich and certainly
merits the detour.

As emphasized in the main text, one may naively expect that coupling the
fermions to a gauge field has the same effect as turning on a complex chemical
potential. However, the gauge field transforms differently to a chemical
potential. An important check is that all physical quantities be invariant
under large gauge transformations. Though the complex SYK model (for large
$N$) is solvable, it is not immediately obvious that the introduction of a
gauge field (or equivalently applying twisted boundary conditions) will retain
this crucial property. Happily, we found that it does, at least in this
quadratic case.

As we have noted, the complex SYK2 model with chemical potential has been
considered in Liao:2020lac , largely using different methods to the path
integral approach utilised here and in Winer:2020mdc . Their model is
equivalent to ours under the replacement of the chemical potential $\mu$ by
the gauge field $iA(\tau)$. What distinguishes our model is that in periodic
real time with period $T$, relevant for the calculation of the SFF, the
constant gauge configuration777Recall that in Euclidean time, the constant
gauge configuration is $A(\tau)=a/\beta$. becomes $A(t)=a/T$. Since $T$ is the
time argument of the SFF, the external gauge field in our model therefore
introduces additional “time” dependence into the SFF calculation that is
absent in the case of a constant chemical potential. Indeed, our SFF results
are not equivalent to those of Liao:2020lac  under the analytic continuation
$\mu\rightarrow ia/T$. Under this replacement, the result for the SFF
slope888We have correctly normalised their result and taken the limit $T\ll 0$
and $m=1$ in Liao:2020lac ,

$\displaystyle g(T)=\exp\left\\{-\frac{NT^{2}}{2}\cos(\mu
T)\right\\}\rightarrow\exp\left\\{-\frac{NT^{2}}{2}\cosh(a)\right\\},$ (64)

has qualitatively different behaviour and is not invariant under large gauge
transformations.

The primary foci of this article were the out-of-time-ordered correlator and
spectral form factor of this model, using a quenched and annealed averaging
protocol respectively. From a technical viewpoint many of the established
tools for solving the SYK2 model can be readily adapted to the present case.
In particular, the fully dressed two-point correlation function can be
computed exactly at large $N$. When transforming to real time and Wick
rotating the gauge field, the gauge field gives rise to an overall phase
factor for the two-point function. When considering expressions involving a
balance of the Green’s function and its conjugate, these phases may cancel. In
particular, the out-of-time-order correlator is unaffected with no exponential
growth for the out-of-time-order correlator, confirming that the SYK2 model’s
integrablity persists at early times even with the gauge field turned on.

The spectral form factor is more interesting. Here we observe two noteworthy
features introduced by the gauge field. The first is an enhancement of the
early time slope. The second is the introduction of an explicit scale for
strong disorder necessary for the presence of zero modes. These zero modes
give rise to an exponential ramp at late times, consistent with existing
results for the quadratic SYK model.

Our calculation results in an SFF slope controlled by the parameter $a$, which
can be thought of as taking values in the interval $[0,2\pi)$ due to the
invariance under large gauge transformations. For $a=0$ we obtain the result
of the ungauged theory studied in Winer:2020mdc . As we increase $a$ from 0 to
$\pi$, the slope decays faster and faster until the exponent diverges to
$-\infty$ for $a=\pi$.

The fact that we may use the parameter $a$ to tune the decay of the slope
presents an interesting possibility. In a quantum chaotic system there is a
timescale which signals the onset of RMT universality; the Thouless time
Nosaka:2018iat . From the perspective of the SFF, it is the time at which the
ramp begins, but identifying it for a given model is nontrivial, because the
beginning of the ramp is usually obscured by the slope. In order to find the
Thouless time, some method is then required to isolate the ramp from the
slope.

In a Gaussian Unitary Ensemble (GUE) this can be achieved by simply
calculating the connected and disconnected pieces of the SFF separately; the
disconnected part contains the slope while the connected part contains the
ramp. However, for the SYK4 model, this technique fails due to the fact that
the connected piece contains its own slope Cotler:2017jue . Various ways of
dealing with this have been proposed Nosaka:2018iat  involving variations of
the SFF. The “Gaussian-filtered SFF” Gharibyan:2018jrp  eliminates the slope
by softening the hard edge of the spectral density from which it arises, thus
allowing for the ramp to be visible at earlier times. The “connected unfolded
SFF” Garcia-Garcia:2018ruf  is the SFF computed using the unfolded spectrum
and removing the disconnected piece.

Our SFF results show that for SYK2 model, the twisting of the fermion boundary
conditions doesn’t significantly affect the large $T$ exponential ramp, while
at small $T$ we can kill off the slope by taking $a\rightarrow\pi$. If it were
the case that taking $a\rightarrow\pi$ had a similar effect on the slope of
the “gauged” complex SYK4 model, while leaving the ramp invariant, then this
limit could be used as another method of controlling the slope independently
of the ramp, for the purposes of identifying the Thouless time. Further
analysis of the gauged SYK4 model would be necessary in order to make precise
statements in this direction.

We briefly note the different mechanisms underlying the linear and exponential
ramps in the Majorana SYK4 model and the (Majorana and complex) SYK2 models,
respectively. For the Majorana SYK4 model, the contributions of the replica-
diagonal saddles to the full SFF just give the disconnected SFF Saad:2018bqo ,
which describes the slope. The ramp can only be seen by including the late
time contributions of the replica-non-diagonal saddles, which spontaneously
break time translation symmetry. In the SYK2 cases considered here and in
Liao:2020lac ; Winer:2020mdc , we see that the ramp has different origins:
considering only replica-diagonal solutions still leads to symmetry breaking
and an exponential ramp, which is not present in the disconnected SFF999The
replica-non-diagonal solutions then appear to play less of a role in these
models..

Specifically, in the Majorana SYK2 case studied in Winer:2020mdc , the $SU(2)$
conjugation symmetry of the two-replica action is spontaneously broken by the
form of the diagonal saddle point solutions above a critical frequency. The
associated Goldstone zero modes give rise to an exponential ramp. In the
complex SYK2 case studied here and in Liao:2020lac , the the two-replica
action contains a term which explicitly breaks the unitary conjugation
symmetry that it would otherwise possess. A subset of the associated Goldstone
soft modes are indeed zero modes (with vanishing action) and these give rise
to an exponential ramp.

In our model, we found that the presence of the gauge field $a$ sets a natural
timescale $\lambda\equiv\frac{|a|}{2m}$ in order for zero modes to be
supported. However, as this scale is not gauge invariant, we expect that
higher order corrections (beyond the quadratic fluctuations considered in this
thesis) will have an effect.

In summary, our results show that coupling the SYK2 model to a “topological”
gauge field gives rise to a solvable model that retains the non-chaos
properties of the former but which exhibits some interesting new features.
This is a promising start for further development of this direction,
specifically a gauged SYKq model. Evidently, the introduction of twisted
boundary conditions with SYKq models do not present systems that are
technically much more difficult than the SYKq models themselves. Of course,
the obvious test of this statement will be the maximally chaotic $q=4$ model.

In particular, a natural follow up to our work would be to study a twisted
SYK2 model perturbed by an SYK4 interaction. In this case, with the “gauge
field” turned off we obtain a complex version of the mass-deformed (Majorana)
SYK model studied in Garcia-Garcia:2018ruf ; Nosaka:2018iat , which exhibits a
chaotic-integrable transition. It would then be interesting to see whether
this transition is affected by the presence of a nontrivial twist of the
fermion boundary conditions. In Bhattacharya:2017vaz  it was found that the
Lyapunov exponent of the complex SYKq model in the $q\rightarrow\infty$ limit
is suppressed by a nonzero chemical potential, with the chemical potential
having an exponentially large dominance over the $q$-body interaction coupling
strength in determining the chaos behaviour of the system. In light of this it
would be interesting to see whether a nonzero twisting parameter has this same
effect of pushing the mass-deformed complex SYK model into an integrable phase
for all values of the SYK4 coupling, from the perspective of spectral
statistics (the SFF), as opposed to the OTOC. An analytic treatment of the SFF
for the mass deformed SYK model is still outstanding101010A qualitative
discussion from an analytic perspective is given in Winer:2020mdc ., so it is
likely that the prospective calculations discussed above would have to be
performed numerically. However, steps in this direction have already been
taken for the complex SYK2 model when perturbed by quartic, non-SYK4
interactions. Namely, the SFF of the SYK2 model perturbed by a non-random
4-fermi interaction term Liao:2021ofk ; Lau:2020qnl  and by a random 4-fermi
SYK${}_{2}^{2}$ interaction (the square of an SYK2 term) Liao:2021uwx  have
been studied. In both cases the solution for the SYK2 SFF serves as the
conceptual and analytical solution about which the full solutions are
constructed as perturbations, with the result that the quartic interactions
induce a mass for the SYK2 soft modes which would otherwise give rise to the
exponential ramp111111While the model in Liao:2021ofk  is chaotic, the model
in Liao:2021uwx  has chaotic and integrable regimes. In the integrable regime
no mass is induced for the soft modes and the exponential ramp survives.. The
exponential ramp is then suppressed, a necessary prerequisite for the
emergence of RMT statistics.

Another natural continuation of our work would then be to repeat our
calculation for the “gauged” SYK2 model perturbed by such non-SYK four-fermi
interactions, and to see if this has an effect on the induced mass which is
responsible for the chaotic-integrable transition. As we have emphasised
repeatedly, the effect of the external gauge field is equivalent to
introducing twisted boundary conditions for the theory which, at least based
on our results here, would impact quantities computed by utilising several
replicas.

Finally, we circle back to our initial motivation and point out that a gauged
SYK model may have a critical role to play in the formulation of a more
conventional weakly-coupled, small curvature gravity dual. We hope that our
work will stimulate further investigations of these topics in the future.

###### Acknowledgements.

We would like to thank Jonathan Shock for his collaboration during earlier
stages of this work, as well as his invaluable Mathematica expertise
throughout. We would also like to extend our thanks to the analymous referee
whose detailed constructive comments and observations regarding the
equivalence of the external gauge field to a twisting of the fermion boundary
conditions led to a much better understanding of the problem. JM acknowledges
support from the “Quantum Technologies for Sustainable Development” grant from
the National Institute for Theoretical and Computational Sciences of South
Africa (NITHECS) and a Simons Associateship at the Abdus Salam International
Center for Theoretical Physics, Trieste where some of this work was carried
out. RPS is supported by a graduate fellowship from the National Institute for
Theoretical and Computational Sciences and by the Shuttleworth Postgraduate
Scholarship Programme. H.J.R.vZ is supported by the “Quantum Technologies for
Sustainable Devlopment" grant from the National Institute for Theoretical and
Computational Sciences (NITHECS).

## References

  * (1) S. Sachdev and J. Ye, Gapless spin fluid ground state in a random, quantum Heisenberg magnet, Phys. Rev. Lett. 70 (1993) 3339, [cond-mat/9212030].
  * (2) J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D 94 (2016), no. 10 106002, [arXiv:1604.07818].
  * (3) J. Polchinski and V. Rosenhaus, The Spectrum in the Sachdev-Ye-Kitaev Model, JHEP 04 (2016) 001, [arXiv:1601.06768].
  * (4) A. M. García-García and J. J. M. Verbaarschot, Spectral and thermodynamic properties of the Sachdev-Ye-Kitaev model, Phys. Rev. D 94 (2016), no. 12 126010, [arXiv:1610.03816].
  * (5) D. J. Gross and V. Rosenhaus, A Generalization of Sachdev-Ye-Kitaev, JHEP 02 (2017) 093, [arXiv:1610.01569].
  * (6) D. J. Gross and V. Rosenhaus, All point correlation functions in SYK, JHEP 12 (2017) 148, [arXiv:1710.08113].
  * (7) P. Saad, S. H. Shenker, and D. Stanford, A semiclassical ramp in SYK and in gravity, arXiv:1806.06840.
  * (8) V. Rosenhaus, An introduction to the SYK model, J. Phys. A 52 (2019) 323001, [arXiv:1807.03334].
  * (9) Y. Gu, A. Kitaev, S. Sachdev, and G. Tarnopolsky, Notes on the complex Sachdev-Ye-Kitaev model, JHEP 02 (2020) 157, [arXiv:1910.14099].
  * (10) J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, JHEP 08 (2016) 106, [arXiv:1503.01409].
  * (11) G. Sárosi, AdS2 holography and the SYK model, PoS Modave2017 (2018) 001, [arXiv:1711.08482].
  * (12) A. Almheiri and J. Polchinski, Models of AdS2 backreaction and holography, JHEP 11 (2015) 014, [arXiv:1402.6334].
  * (13) A. Jevicki, K. Suzuki, and J. Yoon, Bi-Local Holography in the SYK Model, JHEP 07 (2016) 007, [arXiv:1603.06246].
  * (14) J. Engelsöy, T. G. Mertens, and H. Verlinde, An investigation of AdS2 backreaction and holography, JHEP 07 (2016) 139, [arXiv:1606.03438].
  * (15) A. Jevicki and K. Suzuki, Bi-Local Holography in the SYK Model: Perturbations, JHEP 11 (2016) 046, [arXiv:1608.07567].
  * (16) M. Cvetič and I. Papadimitriou, AdS2 holographic dictionary, JHEP 12 (2016) 008, [arXiv:1608.07018]. [Erratum: JHEP 01, 120 (2017)].
  * (17) K. Jensen, Chaos in AdS2 Holography, Phys. Rev. Lett. 117 (2016), no. 11 111601, [arXiv:1605.06098].
  * (18) S. Sachdev, Bekenstein-Hawking Entropy and Strange Metals, Phys. Rev. X 5 (2015), no. 4 041025, [arXiv:1506.05111].
  * (19) K. Bulycheva, A note on the SYK model with complex fermions, JHEP 12 (2017) 069, [arXiv:1706.07411].
  * (20) G. V. Dunne, R. Jackiw, and C. A. Trugenberger, Topological (Chern-Simons) Quantum Mechanics, Phys. Rev. D 41 (1990) 661.
  * (21) R. Floreanini, R. Percacci, and E. Sezgin, Infinite dimensional algebras in Chern-Simons quantum mechanics, Phys. Lett. B 261 (1991) 51–56, [hep-th/9111053].
  * (22) R. Jackiw and S. Y. Pi, Reducing the Chern-Simons term by a symmetry, Phys. Lett. B 423 (1998) 364–368, [hep-th/9712087].
  * (23) G. V. Dunne, K.-M. Lee, and C.-h. Lu, On the finite temperature Chern-Simons coefficient, Phys. Rev. Lett. 78 (1997) 3434–3437, [hep-th/9612194].
  * (24) D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, Sachdev-Ye-Kitaev models and beyond: Window into non-Fermi liquids, Rev. Mod. Phys. 94 (2022), no. 3 035004, [arXiv:2109.05037].
  * (25) J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, Chaos, Complexity, and Random Matrices, JHEP 11 (2017) 048, [arXiv:1706.05400].
  * (26) G. Penington, S. H. Shenker, D. Stanford, and Z. Yang, Replica wormholes and the black hole interior, JHEP 03 (2022) 205, [arXiv:1911.11977].
  * (27) D. A. Trunin, Pedagogical introduction to the Sachdev–Ye–Kitaev model and two-dimensional dilaton gravity, Usp. Fiz. Nauk 191 (2021), no. 3 225–261, [arXiv:2002.12187].
  * (28) N. Sorokhaibam, Phase transition and chaos in charged SYK model, JHEP 07 (2020) 055, [arXiv:1912.04326].
  * (29) C. Baldwin and B. Swingle, Quenched vs annealed: Glassiness from sk to syk, Physical Review X 10 (08, 2020).
  * (30) B. Michel, J. Polchinski, V. Rosenhaus, and S. J. Suh, Four-point function in the IOP matrix model, JHEP 05 (2016) 048, [arXiv:1602.06422].
  * (31) A. I. Larkin and Y. N. Ovchinnikov, Quasiclassical method in the theory of superconductivity, Journal of Experimental and Theoretical Physics (1969).
  * (32) A. Almheiri, D. Marolf, J. Polchinski, D. Stanford, and J. Sully, An Apologia for Firewalls, JHEP 09 (2013) 018, [arXiv:1304.6483].
  * (33) S. H. Shenker and D. Stanford, Black holes and the butterfly effect, JHEP 03 (2014) 067, [arXiv:1306.0622].
  * (34) S. H. Shenker and D. Stanford, Multiple Shocks, JHEP 12 (2014) 046, [arXiv:1312.3296].
  * (35) D. A. Roberts, D. Stanford, and L. Susskind, Localized shocks, JHEP 03 (2015) 051, [arXiv:1409.8180].
  * (36) D. A. Roberts and D. Stanford, Two-dimensional conformal field theory and the butterfly effect, Phys. Rev. Lett. 115 (2015), no. 13 131603, [arXiv:1412.5123].
  * (37) D. Stanford, Many-body chaos at weak coupling, JHEP 10 (2016) 009, [arXiv:1512.07687].
  * (38) O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52 (1984) 1–4.
  * (39) J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, Black Holes and Random Matrices, JHEP 05 (2017) 118, [arXiv:1611.04650]. [Erratum: JHEP 09, 002 (2018)].
  * (40) J. Liu, Spectral form factors and late time quantum chaos, Phys. Rev. D 98 (2018), no. 8 086026, [arXiv:1806.05316].
  * (41) M. Winer, S.-K. Jian, and B. Swingle, An exponential ramp in the quadratic Sachdev-Ye-Kitaev model, Phys. Rev. Lett. 125 (2020) 250602, [arXiv:2006.15152].
  * (42) Y. Liao, A. Vikram, and V. Galitski, Many-body level statistics of single-particle quantum chaos, Phys. Rev. Lett. 125 (2020) 250601, [arXiv:2005.08991].
  * (43) Y. Liao and V. Galitski, Emergence of many-body quantum chaos via spontaneous breaking of unitarity, Phys. Rev. B 105 (2022), no. 14 L140202, [arXiv:2104.05721].
  * (44) A. Altland and A. Kamenev, Wigner-Dyson statistics from the Keldysh sigma-model, Phys. Rev. Lett. 85 (2000) 5615–5618, [cond-mat/0006363].
  * (45) A. Kamenev and M. Mezard, Level correlations in disordered metals: The replica sigma model, Phys. Rev. B 60 (1999) 3944–3954, [cond-mat/9903001].
  * (46) T. Nosaka, D. Rosa, and J. Yoon, The Thouless time for mass-deformed SYK, JHEP 09 (2018) 041, [arXiv:1804.09934].
  * (47) H. Gharibyan, M. Hanada, S. H. Shenker, and M. Tezuka, Onset of Random Matrix Behavior in Scrambling Systems, JHEP 07 (2018) 124, [arXiv:1803.08050]. [Erratum: JHEP 02, 197 (2019)].
  * (48) A. M. García-García, Y. Jia, and J. J. M. Verbaarschot, Universality and Thouless energy in the supersymmetric Sachdev-Ye-Kitaev Model, Phys. Rev. D 97 (2018), no. 10 106003, [arXiv:1801.01071].
  * (49) R. Bhattacharya, S. Chakrabarti, D. P. Jatkar, and A. Kundu, SYK Model, Chaos and Conserved Charge, JHEP 11 (2017) 180, [arXiv:1709.07613].
  * (50) P. H. C. Lau, C.-T. Ma, J. Murugan, and M. Tezuka, Correlated disorder in the SYK2 model, J. Phys. A 54 (2021), no. 9 095401, [arXiv:2003.05401].
  * (51) Y. Liao and V. Galitski, Universal dephasing mechanism of many-body quantum chaos, Phys. Rev. Res. 4 (2022), no. 1 L012037, [arXiv:2110.02976].

##  Appendix A Derivation of the 2-replica SFF

Defining

$\displaystyle
I_{\alpha}=-\int_{0}^{T}dt\sum_{i=1}^{N}\psi_{i_{\alpha}}^{\dagger}\frac{d}{dt}\psi_{i_{\alpha}},$
(65)

we will explicitly perform the disorder average using the probability
distributions (13), to obtain the 2-replica SFF. Note that we take each
replica to have the same gauge field configuration, so the Chern-Simons term
cancels out.

$\displaystyle\overline{Z(iT)Z(iT)^{*}}$ (67) $\displaystyle=$
$\displaystyle\int\left[\prod_{\alpha=1}^{2}\prod_{i_{\alpha}=1}^{N}\mathscr{D}\psi_{i_{\alpha}}^{\dagger}\mathscr{D}\psi_{i_{\alpha}}\right]e^{I_{1}+I_{2}}\left(\sqrt{\frac{N}{2\pi
m^{2}}}\right)^{N}\left(\sqrt{\frac{N}{\pi
m^{2}}}\right)^{N(N-1)}\exp\left\\{i\frac{a}{T}\int
dt\left(\psi^{\dagger}_{i_{1}}\psi_{i_{1}}-\psi^{\dagger}_{i_{2}}\psi_{i_{2}}\right)\right\\}$
$\displaystyle\prod_{i}\int dm_{R,ii}\exp\left\\{-i\int dt\
m_{R,ii}(\psi_{i_{1}}^{\dagger}\psi_{i_{1}}-\psi_{i_{2}}^{\dagger}\psi_{i_{2}})\right\\}e^{-\frac{N}{2m^{2}}m_{R,ii}^{2}}$
$\displaystyle\prod_{i<j}\int dm_{R,ij}\exp\left\\{-i\int dt\
m_{R,ij}(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}+\psi_{j_{1}}^{\dagger}\psi_{i_{1}}-\psi_{j_{2}}^{\dagger}\psi_{i_{2}}-\psi_{i_{2}}^{\dagger}\psi_{j_{2}})\right\\}e^{-\frac{N}{m^{2}}m_{R,ij}^{2}}$
$\displaystyle\prod_{i<j}\int dm_{I,ij}\exp\left\\{\int dt\
m_{I,ij}(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}-\psi_{j_{1}}^{\dagger}\psi_{i_{1}}+\psi_{j_{2}}^{\dagger}\psi_{i_{2}}-\psi_{i_{2}}^{\dagger}\psi_{j_{2}})\right\\}e^{-\frac{N}{m^{2}}m_{I,ij}^{2}}$
$\displaystyle=$
$\displaystyle\int\left[\prod_{\alpha=1}^{2}\prod_{i_{\alpha}=1}^{N}\mathscr{D}\psi_{i_{\alpha}}^{\dagger}\mathscr{D}\psi_{i_{\alpha}}\right]e^{I_{1}+I_{2}}\exp\left\\{i\frac{a}{T}\int
dt\left(\psi^{\dagger}_{i_{1}}\psi_{i_{1}}-\psi^{\dagger}_{i_{2}}\psi_{i_{2}}\right)\right\\}$
$\displaystyle\exp\left\\{-\frac{m^{2}}{2N}\int
dtdt^{\prime}\sum_{i}(\psi_{i_{1}}^{\dagger}\psi_{i_{1}}(t)-\psi_{i_{2}}^{\dagger}\psi_{i_{2}}(t))(\psi_{i_{1}}^{\dagger}\psi_{i_{1}}(t^{\prime})-\psi_{i_{2}}^{\dagger}\psi_{i_{2}}(t^{\prime}))\right.$
$\displaystyle\left.-\frac{m^{2}}{4N}\int
dtdt^{\prime}\sum_{i<j}(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}(t)+\psi_{j_{1}}^{\dagger}\psi_{i_{1}}(t)-\psi_{i_{2}}^{\dagger}\psi_{j_{2}}(t)-\psi_{j_{2}}^{\dagger}\psi_{i_{2}}(t))\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}(t^{\prime})+\psi_{j_{1}}^{\dagger}\psi_{i_{1}}(t^{\prime})-\psi_{i_{2}}^{\dagger}\psi_{j_{2}}(t^{\prime})-\psi_{j_{2}}^{\dagger}\psi_{i_{2}}(t^{\prime}))\right.$
$\displaystyle+\frac{m^{2}}{4N}\int
dtdt^{\prime}\sum_{i<j}(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}(t)-\psi_{j_{1}}^{\dagger}\psi_{i_{1}}(t)-\psi_{i_{2}}^{\dagger}\psi_{j_{2}}(t)+\psi_{j_{2}}^{\dagger}\psi_{i_{2}}(t))$
$\displaystyle\qquad\qquad\qquad\qquad(\psi_{i_{1}}^{\dagger}\psi_{j_{1}}(t^{\prime})-\psi_{j_{1}}^{\dagger}\psi_{i_{1}}(t^{\prime})-\psi_{i_{2}}^{\dagger}\psi_{j_{2}}(t^{\prime})+\psi_{j_{2}}^{\dagger}\psi_{i_{2}}(t^{\prime}))\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\int\left[\prod_{\alpha=1}^{2}\prod_{i_{\alpha}=1}^{N}\mathscr{D}\psi_{i_{\alpha}}^{\dagger}\mathscr{D}\psi_{i_{\alpha}}\right]e^{I_{1}+I_{2}}\exp\Bigg{\\{}-\frac{m^{2}}{2N}\int
dtdt^{\prime}\sum_{\alpha,\beta}\sum_{i,j}(-1)^{\alpha+\beta}\psi_{i_{\alpha}}^{\dagger}\psi_{j_{\alpha}}(t)\psi_{j_{\beta}}^{\dagger}\psi_{i_{\beta}}(t^{\prime})\Bigg{\\}}$
$\displaystyle\exp\left\\{i\frac{a}{T}\int
dt\left(\psi^{\dagger}_{i_{1}}\psi_{i_{1}}-\psi^{\dagger}_{i_{2}}\psi_{i_{2}}\right)\right\\}$

We now define the collective, bilocal field

$\displaystyle
G_{\alpha\beta}(t^{\prime},t)=\sum_{i_{\alpha}=1}^{N}\frac{1}{N}\psi^{\dagger}_{i_{\alpha}}(t)\psi_{i_{\beta}}(t^{\prime}),$
(68)

which we introduce to the partition function via the functional delta function

$\displaystyle\delta\left(G_{\alpha\beta}(t^{\prime},t)-\sum_{i=1}^{N}\frac{1}{N}\psi^{\dagger}_{i_{\alpha}}(t)\psi_{i_{\beta}}(t^{\prime})\right)$
(69)
$\displaystyle=\prod_{\alpha,\beta}\int\mathscr{D}\Sigma_{\alpha\beta}\exp\left\\{-N\Sigma_{\beta\alpha}(t,t^{\prime})\left(G_{\alpha\beta}(t^{\prime},t)-\sum_{i=1}^{N}\frac{1}{N}\psi^{\dagger}_{i_{\alpha}}(t)\psi_{i_{\beta}}(t^{\prime})\right)\right\\}.$
(70)

Note that we have introduced the auxilary field $\Sigma_{\beta\alpha}$ as a
Lagrange multiplier. This allows us to rewrite the quartic term in (67) as

$\displaystyle\exp\left\\{-\frac{m^{2}}{2N}\int
dtdt^{\prime}\sum_{\alpha,\beta}\sum_{i,j}(-1)^{\alpha+\beta}\psi_{i_{\alpha}}^{\dagger}\psi_{j_{\alpha}}(t)\psi_{j_{\beta}}^{\dagger}\psi_{i_{\beta}}(t^{\prime})\right\\}$
$\displaystyle=$
$\displaystyle\prod_{\alpha,\beta}\int\mathscr{D}\Sigma_{\alpha\beta}\mathscr{D}G_{\alpha\beta}\exp\Bigg{\\{}\int
dtdt^{\prime}\bigg{[}(-1)^{\alpha+\beta}\frac{Nm^{2}}{2}G_{\alpha\beta}(t^{\prime},t)G_{\beta\alpha}(t,t^{\prime})$
$\displaystyle\qquad\qquad\qquad\qquad-N\Sigma_{\beta\alpha}(t,t^{\prime})\left(G_{\alpha\beta}(t^{\prime},t)-\sum_{i=1}^{N}\frac{1}{N}\psi^{\dagger}_{i_{\alpha}}(t)\psi_{i_{\beta}}(t^{\prime})\right)\bigg{]}\Bigg{\\}},$
(71)

so that (67) becomes

$\displaystyle\int\left[\prod_{\alpha,\beta}\mathscr{D}\Sigma_{\alpha\beta}\mathscr{D}G_{\alpha\beta}\right]\left[\prod_{\alpha}\prod_{i_{\alpha}=1}^{N}\mathscr{D}\psi_{i_{\alpha}}^{\dagger}\mathscr{D}\psi_{i_{\alpha}}\right]$
$\displaystyle\exp\Bigg{\\{}\int
dtdt^{\prime}\sum_{\alpha,\beta}\sum_{i}\psi_{i_{\alpha}}^{\dagger}(t)\left[-\delta(t-t^{\prime})\left(\delta_{\alpha\beta}\frac{d}{dt}-i\sigma^{z}_{\alpha\beta}\frac{a}{T}\right)+\Sigma_{\beta\alpha}(t,t^{\prime})\right]\psi_{i_{\beta}}(t^{\prime})$
$\displaystyle-N\int
dtdt^{\prime}\sum_{\alpha,\beta}\left[\Sigma_{\beta\alpha}(t^{\prime},t)G_{\alpha\beta}(t,t^{\prime})-(-1)^{\alpha+\beta}\frac{m^{2}}{2}G_{\alpha\beta}(t,t^{\prime})G_{\beta\alpha}(t^{\prime},t)\right]\Bigg{\\}},$
(72)

where $\sigma^{z}$ is the usual Pauli $z$-matrix. Performing the Gaussian
integral over the fermions gives the following expression for the SFF

$\displaystyle\overline{Z(iT)Z(iT)^{*}}$ $\displaystyle=$
$\displaystyle\prod_{\alpha,\beta}\int\mathscr{D}\Sigma_{\alpha\beta}\mathscr{D}G_{\alpha\beta}e^{NI(G,\Sigma)}\,,$
(73)

where

$\displaystyle I(G,\Sigma)$ $\displaystyle=$
$\displaystyle\log\det\left[-\delta(t-t^{\prime})\left(\delta_{\alpha\beta}\frac{d}{dt}-i\sigma^{z}_{\alpha\beta}\frac{a}{T}\right)+\Sigma_{\beta\alpha}(t,t^{\prime})\right]$
$\displaystyle-\int
dtdt^{\prime}\sum_{\alpha,\beta}\left[\Sigma_{\beta\alpha}(t^{\prime},t)G_{\alpha\beta}(t,t^{\prime})-(-1)^{\alpha+\beta}\frac{m^{2}}{2}G_{\alpha\beta}(t,t^{\prime})G_{\beta\alpha}(t^{\prime},t)\right].$

##  Appendix B Partition function

To calculate the partition function (in Euclidean time) for our model, we
simply restrict the results (73) and (A) to a single replica, and analytically
continue $t\rightarrow-i\tau$, $a\rightarrow ia$ and $m\rightarrow im$, to
obtain

$\displaystyle\overline{Z(\beta)}$ $\displaystyle=$
$\displaystyle\int\mathscr{D}\Sigma\mathscr{D}Ge^{NI},$ $\displaystyle I$
$\displaystyle=$
$\displaystyle\log\det\left[-\delta(\tau-\tau^{\prime})\left(\frac{d}{d\tau}-i\frac{a}{\beta}\right)+\Sigma(\tau,\tau^{\prime})\right]$
(75) $\displaystyle-\int d\tau
d\tau^{\prime}\left[\Sigma(\tau^{\prime},\tau)G(\tau,\tau^{\prime})+\frac{m^{2}}{2}G(\tau,\tau^{\prime})G(\tau^{\prime},\tau)\right],$

with equations of motion given in the large-$N$ saddle point approximation by

$\displaystyle\delta(\tau-\tau^{\prime\prime})$ $\displaystyle=\int
d\tau^{\prime}\left[-\delta(\tau-\tau^{\prime})\left(\frac{d}{d\tau}-i\frac{a}{T}\right)+\Sigma(\tau-\tau^{\prime})\right]G(\tau^{\prime}-\tau^{\prime\prime}),$
(76) $\displaystyle G(\tau)$ $\displaystyle=-m^{2}\Sigma(\tau).$ (77)

We first consider the free theory where $m=0$. The operator
$\partial_{\tau}-iA(\tau)$ may be explicitly diagonalised as follows
Dunne:1996yb

$\displaystyle\left(\frac{d}{d\tau}-iA(\tau)\right)\psi_{n}(\tau)$
$\displaystyle=$
$\displaystyle\Lambda_{n}\psi_{n}(\tau)\quad\Rightarrow\quad\psi_{n}(\tau)=\exp\left\\{\int_{0}^{\tau}d\tau^{\prime}(iA(\tau^{\prime})+\Lambda_{n})\right\\}.$

Making use of the anti-periodicity of the fermionic field
$\psi(0)=-\psi(\beta)$, it follows that

$\displaystyle\Lambda_{n}\beta+i\int_{0}^{\beta}A(\tau)d\tau$ $\displaystyle=$
$\displaystyle(2n-1)\pi
i\quad\Rightarrow\quad\Lambda_{n}=-\frac{ia}{\beta}+\frac{(2n-1)\pi
i}{\beta},$ (79)

where $n\in\textbf{Z}$ and $a=\int_{0}^{\beta}A(\tau)d\tau$ as before. The
free partition is then given by

$\displaystyle Z_{0}(\beta)$
$\displaystyle=\left[\det\left(-\partial_{\tau}+iA(\tau)\right)\right]^{N}=\prod_{n=-\infty}^{\infty}\left(\frac{ia}{\beta}-\frac{(2n-1)\pi
i}{\beta}\right)^{N},$ (80)

which is formally infinite. We regulate it by dividing by the analogous
expression in the absence of the gauge field to obtain

$\displaystyle Z_{0}(\beta)$
$\displaystyle=\left[\frac{\det\left(-\partial_{\tau}+iA(\tau)\right)}{\det\left(-\partial_{\tau}\right)}\right]^{N}=\prod_{n=-\infty}^{\infty}\left[\frac{-\frac{ia}{\beta}+\frac{(2n-1)\pi
i}{\beta}}{\frac{(2n-1)\pi
i}{\beta}}\right]^{N}=\prod_{n=-\infty}^{\infty}\left[1-\frac{a}{(2n-1)\pi}\right]^{N}$
$\displaystyle=\left[\cos\left(\frac{a}{2}\right)\right]^{N}.$ (81)

This means that the free theory with $m=0$ has a free energy of

$\displaystyle F_{0}(\beta)=-\frac{1}{\beta}\log
Z_{0}(\beta)=-\frac{N}{\beta}\log\left[\cos\left(\frac{a}{2}\right)\right],$
(82)

which diverges whenever $a$ is an odd multiple of $\pi$.

Generalising to nonzero $m$, direct evaluation of the functional determinant
in (75) is best avoided in favour of a different method: formulating a
differential equation for the free energy Maldacena:2016hyu ; Winer:2020mdc
by exploiting the fact that the only dimensionful quantities in the problem
are $m$ and $\beta$, so all $\beta$ dependence must be through $m\beta$. It
follows that partial derivatives acting on the free energy are related as
$\displaystyle\beta\frac{d}{d\beta}=m\frac{d}{dm}$, so that

$\displaystyle m^{2}\beta\frac{dI}{d\beta}$
$\displaystyle=-2\frac{dI}{d(m^{-2})}=-\int d\tau
d\tau^{\prime}\Sigma(\tau^{\prime}-\tau)\Sigma(\tau-\tau^{\prime}).$ (83)

Here we have used that the on-shell action (B) only has explicit $m$
dependence since any dependence through the fields doesn’t contribute because
they satisfy the equations of motion $\frac{\delta I}{\delta\Sigma}=0$. The
equations of motion (76) may be manipulated to obtain

$\displaystyle\int d\tau
d\tau^{\prime}\Sigma(\tau-\tau^{\prime})\Sigma(\tau^{\prime}-\tau)=\left[\beta\frac{d}{d\tau}\Sigma(\tau)-ia\Sigma(\tau)\right]_{\tau\rightarrow
0^{+}}.$ (84)

Our result for the real-time propagator (20) and (77) gives

$\displaystyle\Sigma(\tau)$
$\displaystyle=m^{2}e^{i\frac{a}{\beta}\tau}\int_{-2m}^{2m}dy\frac{\sqrt{4m^{2}-y^{2}}}{2\pi
m^{2}}\frac{e^{-y\tau}}{1+e^{-y\beta}e^{ia}}\,,$ (85)

which combined with (83) with (84) yields the desired differential equation

$\displaystyle\frac{dI}{d\beta}$
$\displaystyle=\int_{-2m}^{2m}dy\frac{y\sqrt{4m^{2}-y^{2}}}{2\pi
m^{2}(1+e^{ia-\beta y})}.$ (86)

Integrating w.r.t. $\beta$ (swapping the order of integration) yields

$\displaystyle I$
$\displaystyle=\int_{-2m}^{2m}dy\frac{\sqrt{4m^{2}-y^{2}}\log\left(e^{\beta
y}+e^{ia}\right)}{2\pi m^{2}}+C.$ (87)

where $C$ is a constant of integration. The integrand may be expanded as a
polynomial in $\beta$ using

$\displaystyle\sqrt{4m^{2}-y^{2}}\log\left(e^{\beta y}+e^{ia}\right)$
$\displaystyle=\sum_{n=0}^{\infty}y^{n}\sqrt{4m^{2}-y^{2}}\frac{f_{n}(e^{ia})}{n!(1+e^{ia})^{n}}\beta^{n}$
(88)

where $f_{n}$ is a polynomial, given for the first few terms by

$\displaystyle f_{0}(x)$ $\displaystyle=\log\left(1+x\right)$ $\displaystyle
f_{1}(x)$ $\displaystyle=1,$ $\displaystyle f_{2}(x)$ $\displaystyle=x,$
$\displaystyle f_{3}(x)$ $\displaystyle=x^{2}-x,$ $\displaystyle f_{4}(x)$
$\displaystyle=x^{3}-4x^{2}+x,$ $\displaystyle f_{5}(x)$
$\displaystyle=x^{4}-11x^{3}+11x^{2}-x.$

We may then perform the integral term by term, to obtain

$\displaystyle I$
$\displaystyle=-\frac{1}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{2^{n}}{n!}\frac{\Gamma\left(\frac{n+1}{2}\right)}{\Gamma\left(\frac{n}{2}+2\right)}\frac{f_{n}(e^{ia})}{(1+e^{ia})^{n}}m^{n}\beta^{n}$
$\displaystyle=\log\cos\left(\frac{a}{2}\right)+\frac{1}{8}\sec^{2}\left(\frac{a}{2}\right)m^{2}\beta^{2}+\mathcal{O}(m^{4}\beta^{4}),$
(89)

where we have fixed the integration constant $C=\log(e^{-ia/2}/2)$ to ensure
that we recover (B) in the limit $m\rightarrow 0$. Some remarks are in order.
Note that this expression is invariant under large gauge transformations. It
diverges when $a$ is equal to an odd multiple of $\pi$. We verify that it is
equal, order by order, to the exponent of the $mT\ll 0$ SFF (D) up to
normalisation (unlike the partition function, the SFF must satisfy $g(0)=1$),
the replacement $m\rightarrow im$ (to account for the Wick rotation to real
time) and an overall factor of $2$ (since the SFF is the square of a partition
function). It therefore also serves as a useful check of the early time SFF
result in the main text, which is given by the disconnected component of the
SFF (32) as expected Cotler:2016fpe .

##  Appendix C SFF fluctuations: Mass matrices

For each frequency, only saddle point solutions whose mass matrices have
positive real parts (for both replica indices) will contribute in the large
$N$ saddle point approximation. In this appendix we carefully check this
condition for each saddle point in each of the four regions.

Region $A$: $|x_{n}\pm b|\leq 2$, from which it follows that $|b|<2$ and
$|x_{n}|<2-|b|$. Note that for $|b|>2$, this saddle point $A$ does not
contribute at all. The solution in this region is given by

$\displaystyle\mathcal{S}_{\alpha\beta}^{nn^{\prime}(A)}=\frac{1}{2}\left[-i(\xi_{\alpha}x_{n}+b)+c_{\alpha}^{n}\sqrt{4-(\xi_{\alpha}x_{n}+b)^{2}}\right]\delta^{nn^{\prime}}_{\alpha\beta},$
(90)

We can independently choose $c_{1}^{n}$ and $c_{2}^{n}$, giving four possible
solutions for each $n$. The corresponding mass matrices are given by

$\displaystyle
M^{nn^{\prime}(A)}_{\alpha\beta}=1+\frac{1}{4}\left[-i(\xi_{\alpha}x_{n}+b)+c_{\alpha}^{n}\sqrt{4-(\xi_{\alpha}x_{n}+b)^{2}}\right]\left[-i(\xi_{\beta}x_{n^{\prime}}+b)+c_{\beta}^{n^{\prime}}\sqrt{4-(\xi_{\beta}x_{n^{\prime}}+b)^{2}}\right]$
(91)

where $|\xi_{\alpha}x_{n}+b|\leq 2$ and $|\xi_{\beta}x_{n^{\prime}}+b|\leq 2$.
It can be verified that Re$\left[M^{nn^{\prime}}_{\alpha\beta}\right]\geq 0$
for all choices of $c_{\alpha}^{n}$ and $c_{\beta}^{n^{\prime}}$. We note that
when $c_{\alpha}^{n}=-c_{\beta}^{n^{\prime}}$ and
$\xi_{\alpha}x_{n}=\xi_{\beta}x_{n^{\prime}}$ ($\alpha\neq\beta$ and
$n=-n^{\prime}$)121212The other case $\alpha=\beta$, $n=n^{\prime}$ which
satisfies $\xi_{\alpha}x_{n}=\xi_{\beta}x_{n^{\prime}}$ is ruled out by the
condition $c_{\alpha}^{n}=-c_{\beta}^{n^{\prime}}$. the corresponding mass
matrices vanish: $M^{n,-n}_{12}=M^{n,-n}_{21}=0$. This signals the presence of
zero modes, whose significance is discussed in the main text. Since the mass
matrix is always non-negative, all four saddle point solutions contribute to
the SFF in this region. Note also that the mass matrices for the two saddles
with $c_{1}=c_{2}$ are complex conjugates of each other

$\displaystyle\left(M^{nn^{\prime}(A)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\pm
1}\right)^{*}=M^{nn^{\prime}(A)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\mp
1},$ (92)

a fact that will become relevant when evaluating the late time SFF later.

Region $D$: $|x_{n}\pm b|>2$. In this case the solution is completely
determined and there is only one choice for each frequency. This regime can be
written $|x_{n}|>2+|b|$ and $|x_{n}|<|b|-2$ (when $|b|<2$). The solution here
is given by

$\displaystyle\mathcal{S}_{\alpha\beta}^{nn^{\prime}(D)}=\frac{1}{2}\left[-i(\xi_{\alpha}x_{n}+b)+i\text{sgn}(\xi_{\alpha}x_{n}+b)\sqrt{(\xi_{\alpha}x_{n}+b)^{2}-4}\right]\delta^{nn^{\prime}}_{\alpha\beta},$
(93)

so we have

$\displaystyle M^{nn^{\prime}(D)}_{\alpha\beta}=$ $\displaystyle
1+\frac{1}{4}\left[-i(\xi_{\alpha}x_{n}+b)+i\text{sgn}(\xi_{\alpha}x_{n}+b)\sqrt{(\xi_{\alpha}x_{n}+b)^{2}-4}\right]$
$\displaystyle\qquad\times\left[-i(\xi_{\beta}x_{n^{\prime}}+b)+i\text{sgn}(\xi_{\beta}x_{n^{\prime}}+b)\sqrt{(\xi_{\beta}x_{n^{\prime}}+b)^{2}-4}\right]$
(94)

where $|\xi_{\alpha}x_{n}+b|>2$ and $|\xi_{\beta}x_{n^{\prime}}+b|>2$. In this
case
Re$\left[M^{nn^{\prime}}_{\alpha\beta}\right]=M^{nn^{\prime}}_{\alpha\beta}>0$
for all allowed frequencies and replica indices. Of course since
$M^{nn^{\prime}}_{\alpha\beta}$ is real and independent of the constants
$c_{\alpha}^{n}$, we have that

$\displaystyle\left(M^{nn^{\prime}(D)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\pm
1}\right)^{*}=M^{nn^{\prime}(D)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\mp
1}.$ (95)

Region $B$: $|x_{n}+b|\leq 2$, $|x_{n}-b|>2$. We can choose the single
constant $c_{1}^{n}$, giving two possible solutions for each $n$. We can write
the region as

$\displaystyle|b|>2:$ $\displaystyle\qquad|x_{n}+b|<2,$ (96)
$\displaystyle|b|\leq 2:$ $\displaystyle\qquad|x_{n}+\text{sgn}(b)2|<|b|.$
(97)

The solutions here is given by

$\displaystyle\mathcal{S}_{\alpha\beta}^{nn^{\prime}(B)}=$
$\displaystyle=\frac{1}{2}\left[-i(\xi_{\alpha}x_{n}+b)+c_{\alpha}^{n}\sqrt{4-(\xi_{\alpha}x_{n}+b)^{2}}\right]\delta^{nn^{\prime}}_{\alpha
1}$
$\displaystyle+\frac{1}{2}\left[-i(\xi_{\beta}x_{n}+b)+i\text{sgn}(\xi_{\beta}x_{n}+b)\sqrt{(\xi_{\beta}x_{n}+b)^{2}-4}\right]\delta^{nn^{\prime}}_{2\beta}.$
(98)

The first diagonal component of the mass matrix is given by

$\displaystyle M^{nn^{\prime}(B)}_{11}=$ $\displaystyle
M^{nn^{\prime}(A)}_{11},$ (99)

where $|x_{n}+b|\leq 2$ and $|x_{n^{\prime}}+b|\leq 2$. It follows that
Re$\left[M^{nn^{\prime}}_{11}\right]>0$. Note that the zero modes from $A$ do
note feature here since they only exist for $\alpha\neq\beta$. The second
diagonal component is given by

$\displaystyle M^{nn^{\prime}(B)}_{22}=$ $\displaystyle
M^{nn^{\prime}(D)}_{22},$ (100)

where $|x_{n}-b|>2$ and $|x_{n^{\prime}}-b|>2$. In this case
Re$\left[M^{nn^{\prime}}_{22}\right]>0$ as above. Finally, we have the replica
non-diagonal components given by

$\displaystyle M^{nn^{\prime}(B)}_{\alpha\neq\beta}=$ $\displaystyle
1+\frac{1}{4}\left[-i(\xi_{\alpha}x_{n}+b)+i\text{sgn}(\xi_{\alpha}x_{n}+b)\sqrt{(\xi_{\alpha}x_{n}+b)^{2}-4}\right]$
$\displaystyle\qquad\times\left[-i(\xi_{\beta}x_{n^{\prime}}+b)+c_{\beta}^{n^{\prime}}\sqrt{(\xi_{\beta}x_{n^{\prime}}+b)^{2}-4}\right],$
(101)

where $|\xi_{\alpha}x_{n}+b|\leq 2$ and $|\xi_{\beta}x_{n^{\prime}}+b|>2$. In
all of the allowed frequency regimes, we have that
Re$\left[M^{nn^{\prime}}_{\alpha\neq\beta}\right]>0$. Thus the saddles for
both choices of $c_{1}^{n}$ contribute at each frequency in this regime. In
all three cases we have that

$\displaystyle\left(M^{nn^{\prime}(B)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\pm
1}\right)^{*}=M^{nn^{\prime}(B)}_{\alpha\beta}\big{|}_{c_{1}^{n}=c_{2}^{n^{\prime}}=\mp
1}.$ (102)

Region $C$: $|x_{n}-b|\leq 2$, $|x_{n}+b|>2$. We can choose the single
constant $c_{2}^{n}$, giving two possible solutions for each $n$. The region
can be written as

$\displaystyle|b|>2:$ $\displaystyle\qquad|x_{n}-b|<2,$ (103)
$\displaystyle|b|\leq 2:$ $\displaystyle\qquad|x_{n}-\text{sgn}(b)2|<|b|.$
(104)

Similarly to the above case, we find that the saddles for both choices of
$c_{2}^{n}$ contribute at each allowed frequency, since

$\displaystyle M^{nn^{\prime}(C)}_{11}$
$\displaystyle=M^{nn^{\prime}(D)}_{11}>0,$ (105) $\displaystyle
M^{nn^{\prime}(C)}_{22}$ $\displaystyle=M^{nn^{\prime}(A)}_{22}>0,$ (106)
$\displaystyle M^{nn^{\prime}(B)}_{\alpha\neq\beta}$ $\displaystyle>0.$ (107)

Also,

$\displaystyle\left(M^{nn^{\prime}(C)}_{\alpha\beta}\big{|}_{c_{1}^{n},c_{2}^{n^{\prime}}=\pm
1}\right)^{*}=M^{nn^{\prime}(C)}_{\alpha\beta}\big{|}_{c_{1}^{n},c_{2}^{n^{\prime}}=\mp
1}.$ (108)

The above constitutes an extensive list of all mass matrices one may encounter
when performing the saddle point approximation for the SFF.

##  Appendix D $mT\ll 1$ SFF

$\displaystyle g(T)=$
$\displaystyle\exp\Bigg{\\{}N\Bigg{(}-\frac{1}{4}m^{2}T^{2}\sec^{2}\left(\frac{a}{2}\right)+\frac{1}{48}m^{4}T^{4}(\cos
a-2)\sec^{4}\left(\frac{a}{2}\right)$
$\displaystyle\qquad-\frac{1}{2304}m^{6}T^{6}(-26\cos a+\cos
2a+33)\sec^{6}\left(\frac{a}{2}\right)$
$\displaystyle\qquad-\frac{m^{8}T^{8}(1191\cos
a-120\cos(2a)+\cos(3a)-1208)\sec^{8}\left(\frac{a}{2}\right)}{184320}+\mathcal{O}(m^{10}T^{10})\Bigg{)}\Bigg{\\}}.$
(109)

Note that this expression is invariant under large gauge transformations
$a\rightarrow a+2k\pi$, $k\in\textbf{Z}$. As $a$ approaches $\pi$, our result
becomes trustworthy only for smaller and smaller times $mT<\frac{\pi-a}{2}$,
since we have assumed this in order to neglect the contribution from the
saddle $C$. For $mT<\frac{\pi-a}{2}$ we always see the characteristic slope of
the SFF, see figure 3.

Figure 3: SFF (D) for $mT\ll 1$ for various $a$, to $\mathcal{O}(m^{2}T^{2})$
(blue) and order $\mathcal{O}(m^{8}T^{8})$ (yellow) in the exponent. Grey
vertical lines mark $mT=\frac{\pi-a}{2}$. Here we have set $N=1$, for large
$N$ the decay is much more rapid.

##  Appendix E SFF for $\pi-a<2mT\ll 1$

That gauge field configurations for which $a=\int_{0}^{\beta}dtA(t)$ is an odd
multiple of $\pi$ are unphysical can be motivated by the fact that in such
cases, we are effectively shifting all fermionic Matsubara frequencies to
bosonic Matsubara frequencies, for which $\omega_{n}=\frac{2n\pi}{T}$ for
integer $n$. This in turn violates the antiperiodicity requirement on our
fields across the thermal circle. Indeed, the partition function diverges in
such cases, as is shown in Appendix (B).

Figure 4: Log-log plots of the SFF (111) for $a=\pi-\epsilon$ and
$0<\epsilon<2mT\ll 1$, with $N=1$.

In evaluating the early time SFF (53), as $a$ approaches $\pi$ we are
eventually forced to take the contribution from frequencies $x_{n}$ in region
$C$, where $|n\pi-a|<2mT$, into account. To this end, we take
$\epsilon\equiv\pi-a>0$ where $0<\epsilon<2mT\ll 1$, so that the $n=1$
frequency contribution now comes from region $C$. We then evaluate the SFF
(53) the same as before, except that the $n=1$ contribution now comes from
region $C$ and not region $D$

$\displaystyle\overline{Z(iT)Z^{*}(iT)}=$ $\displaystyle\sum_{c,\tilde{c}=\pm
1}\exp\Bigg{\\{}N\left[I_{C}^{(c)}\left(\frac{\pi}{mT}\right)+I_{C}^{(\tilde{c})}\left(\frac{\pi}{mT}\right)+2\sum_{n>1,\text{odd}}I_{D}(x_{n})\right]\Bigg{\\}}.$
(110)

Expanding to $\mathcal{O}(m^{2}T^{2})$, performing the sums and implementing
the normalisation (including dividing by the factor of 4 since we sum over 4
saddles) we obtain

$\displaystyle g(T)=$
$\displaystyle\frac{1}{4}\left(1+\frac{m^{2}T^{2}}{e\epsilon^{2}}\exp\left\\{\frac{m^{2}T^{2}}{\epsilon^{2}}+\frac{\epsilon^{2}}{2m^{2}T^{2}}\right\\}\right)^{2N}\exp\left\\{N\frac{m^{2}T^{2}}{2\cos\epsilon-2}\right\\}.$
(111)

This result is valid to all orders in $\epsilon$ for times $\epsilon<2mT\ll 1$
(indeed it diverges as $mT\rightarrow 0$ for any finite $\epsilon$). We verify
that the same result is obtained had we set $a=3\pi-\epsilon$ at the start of
the calculation and accordingly evaluated the $C$ saddles at $x_{n}=3\pi/mT$.
We expect that the result is invariant under all large gauge transformations
$a\rightarrow a+2\pi k$, $k\in Z$ in this sense.

We plot our result (approximating $\cos\epsilon\approx 1-\epsilon^{2}/2$ to
obtain a result which is a function only of the combination $mT/\epsilon$) in
figure 4. On the left we zoom into the region $\epsilon<2mT<2\epsilon$ and we
see the end of the slope and the beginning of a region of increasing
behaviour: this is the characteristic “dip” of the SFF. We may numerically
estimate the dip time as

$\displaystyle mT_{\text{dip}}\approx 0.73(\pi-a),$ (112)

valid for values of the the gauge field strength $a$ such that $0<\pi-a\ll 1$.
If we simply ignore this range of validity and set $a=0$, we obtain a dip time
(for $m=1$) for the SYK2 model of

$\displaystyle T_{\text{dip}}\approx 0.73\pi.$ (113)

Using equations (17) and (33) of Winer:2020mdc  for the SFF131313Note that we
are here comparing the dip times of the complex and Majorana SYK2 models. This
is justified since we expect these results to differ only by a power of 2
(accounting for the difference between the number of degrees of freedom in the
models) which does not affect the dip time., we obtain a dip time of

$\displaystyle T_{\text{dip}}\approx 0.76\pi.$ (114)

Using equation (5b) of Liao:2020lac  for the SFF, with the chemical potential
switched off, we obtain a dip time of

$\displaystyle T_{\text{dip}}\approx 0.8\pi.$ (115)

It thus appears that our result for the $0<\pi-a\ll 1$ dip time remains valid
for values of $a$ outside of these bounds. It would be interesting to
investigate the reasons for this further, which we leave to future work.

Note that the increase we see after the dip time is exponential in $N$,
whereas we know that the late time ramp is exponential in $\log N$. This
exponential in $N$ behaviour therefore likely corresponds to
fluctuations141414See figure 1 of Winer:2020mdc  or Liao:2020lac . which are
known to exist at the start of the SYK2 exponential ramp, and cannot be
trusted up to the values of $mT$ plotted on the right of figure 4.

##  Appendix F Vanishing of the $mT\gg 1$ sum over saddles

We may use the following features of the real and imaginary parts of the on-
shell actions

$\displaystyle[I_{A}^{(1,1)}(x)]^{*}$
$\displaystyle=I_{A}^{(-1,-1)}(x),\qquad[I_{A}^{(1,-1)}(x)]^{*}=I_{A}^{(-1,1)}(x)$
$\displaystyle\text{Re}[I_{A}^{(1,1)}(x)]$
$\displaystyle=\text{Re}[I_{A}^{(1,-1)}(x)],\qquad[I_{C}^{(1)}(x)]^{*}=I_{C}^{(-1)}(x),$
(116)

to perform the sum over saddles $\\{c\\},\\{\tilde{c}\\}$ in (52) and rewrite
the SFF valid for large times $2mT>|a|$ as

$\displaystyle\overline{Z(iT)Z^{*}(iT)}$ $\displaystyle=$
$\displaystyle\prod_{\begin{subarray}{c}n\pi>0\\\ n\
\text{odd}\end{subarray}}^{2mT-|a|}\left[e^{2NI_{A}^{(1,1)}(x_{n})}+e^{2N[I_{A}^{(1,1)}(x_{n})]^{*}}+e^{2NI_{A}^{(1,-1)}(x_{n})}+e^{2N[I_{A}^{(1,-1)}(x_{n})]^{*}}\right]$
$\displaystyle\prod_{\begin{subarray}{c}n\pi>2mT-|a|\\\ n\
\text{odd}\end{subarray}}^{2mT+|a|}\left[e^{2NI_{C}^{(1)}(\text{sgn}(a)x_{n})}+e^{2N[I_{C}^{(1)}(\text{sgn}(a)x_{n})]^{*}}\right]\prod_{\begin{subarray}{c}n\pi>2mT+|a|\\\
n\ \text{odd}\end{subarray}}^{\infty}e^{2NI_{D}(x_{n})}$ $\displaystyle=$
$\displaystyle\prod_{\begin{subarray}{c}n\pi>0\\\ n\
\text{odd}\end{subarray}}^{2mT-|a|}2e^{2N\text{Re}[I_{A}^{(1,1)}(x_{n})]}\left[\cos\left(2N\text{Im}[I_{A}^{(1,1)}(x_{n})]\right)+\cos\left(2N\text{Im}[I_{A}^{(1,-1)}(x_{n})]\right)\right]$
$\displaystyle\prod_{\begin{subarray}{c}n\pi>2mT-|a|\\\ n\
\text{odd}\end{subarray}}^{2mT+|a|}4\cos^{2}\left(N\text{Im}[I_{C}^{(1)}(\text{sgn}(a)x_{n})]\right)e^{2N\text{Re}[I_{C}^{(1)}(\text{sgn}(a)x_{n})]}\prod_{\begin{subarray}{c}n\pi>2mT+|a|\\\
n\ \text{odd}\end{subarray}}^{\infty}e^{2NI_{D}(x_{n})}.$

In the limit $mT\gg 1$, the variable $x_{n}=\frac{n\pi}{mT}$ becomes
continuous and the late time expression for the SFF is proportional to

$\displaystyle\prod_{x=0}^{2-|b|}\left[\cos\left(2N\text{Im}[I_{A}^{(1,1)}(x)]\right)+\cos\left(2N\text{Im}[I_{A}^{(1,-1)}(x)]\right)\right]$
$\displaystyle=$
$\displaystyle\prod_{x=0}^{2-|b|}\cos\left[N\left(\text{Im}[I_{A}^{(1,1)}(x)+I_{A}^{(1,-1)}(x_{n})]\right)\right]\cos\left[N\left(\text{Im}[I_{A}^{(1,1)}(x_{n})-I_{A}^{(1,-1)}(x)]\right)\right]$

This factor, and so the SFF, vanishes when the argument of either cosine is
equal to $k\pi/2$, $k$ odd. For $k$ such that $k=\pm 2N$ this condition
becomes

$\displaystyle\text{Im}[I_{A}^{(1,1)}(x)\pm I_{A}^{(1,-1)}(x)]=\pm\pi.$ (119)

It is easily verified numerically that this condition is always satisfied for
some $x\in[0,2-|b|)$ when $|b|\leq 1$. The SFF (F) is also proportional to

$\displaystyle\prod_{x=2-|b|}^{2+|b|}\cos^{2}\left(N\text{Im}[I_{C}^{(1)}(\text{sgn}(b)x)]\right),$
(120)

and it can similarly be verified that this factor vanishes for some
$x\in(2-|b|,2+|b|)$ when $|b|>1$. Thus the SFF obtained purely from the saddle
point contributions vanishes at late times for any $b=\frac{a}{mT}$.