# Quantum chaos and phase transition in the Yukawa-SYK model

Andrew Davis  Yuxuan Wang  Department of Physics, University of Florida,
Gainesville, Florida, 32601, USA

###### Abstract

We analyze the relation between quantum chaotic behavior and phase transition
of the Yukawa-SYK model as a function of filling and temperature, which
describes random Yukawa interactions between $N$ complex fermions and $M$
bosons in zero spatial dimensions, for both the non-Fermi liquid and
insulating states at finite temperature and chemical potential. We solve the
ladder equations for the out-of-time-order correlator (OTOC) for both the
bosons and fermions. Despite the appearance of the chemical potential in the
Hamiltonian, which explicitly introduces an additional energy scale, the OTOCs
for the fermions and bosons in the non-Fermi liquid state turn out to be
unaffected, and the Lyapunov exponents that diagnose chaos remain maximal. As
the chemical potential increases, the system is known to experience a first-
order transition from a critical phase to a gapped insulating phase. We
postulate that the boundary of the region in parameter space where each phase
is (meta)stable coincides with the curve on which the Lyapunov exponent is
maximal. By calculating the exponent in the insulating phase and comparing to
numerical results on the boundaries of stability, we show that this is
plausible.

##  I Introduction

The Sachdev-Ye-Kitaev (SYK) model [1] has attracted great interest in both
high-energy and condensed matter communities [2, 3, 4]. It consists of $N$
flavors of fermions in $(0+1)$-d coupled by a random, Gaussian, all-to-all
interaction, and it exhibits several remarkable properties [5]. First, it is
exactly solvable in the large-$N$ limit, despite the strong coupling. Second,
it possesses a (near) conformal symmetry in the infrared and has a holographic
dual, corresponding to a black hole in near-AdS2, making it a tractable
platform for the study of AdS/CFT correspondence [6]. Finally, it is a fast
scrambler – it is maximally chaotic as measured by the Lyapunov exponent
$\lambda_{L}=2\pi T$, a feature shared by black holes in Einstein gravity [5].
This makes it noteworthy as an example of a solvable system for quantum
gravity.

On the condensed matter side, the SYK model and its variants are useful
platforms for the understanding of non-Fermi liquid (nFL) physics [7, 8, 9],
which is relevant to the strange metal behavior of unconventional
superconductors [10, 11, 12, 13, 14, 15]. Non-Fermi liquids are characterized
by a power-law frequency dependence of the fermion self-energy and a resulting
spectral function without quasiparticle peaks, and occurs, for example, when
massless bosonic modes destroy the coherence of the fermions [16]. While the
fermions in the SYK model display nFL behavior in the near-conformal regime,
perhaps more relevant to condensed matter physics is a variant of the SYK
model, dubbed the Yukawa-SYK (Yukawa-SYK) model, in which the random four-
fermion interaction is replaced with a random Yukawa interaction between $N$
flavors of complex fermions and $M$ flavors of bosons with bare mass $m_{0}$.
The Yukawa-SYK model has been studied in the context of superconducting
instabilities of nFLs [17, 18, 19, 20] and has been generalized to finite
dimensions [21]. Interestingly, it has been shown that the pairing transition
is dual to that of a holographic superconductor in AdS2 spacetime [22].

At a generic filling factor, the Yukawa-SYK model has a nFL phase as well as
an incompressible, insulating phase. [23] The nFL state exhibits power-law
self-energies for both the fermions and bosons with $M/N$-dependent exponents,
and also has the interesting property that the boson mass flows to zero even
if the bare mass is large (this “self-tuning” also occurs in $(1+1)$-d but not
higher dimensions [21]). For small chemical potential, the nFL is stable, but
at larger $\mu$ there are nFL solutions with negative compressibility,
indicating an unstable phase. At a critical value of $\mu$, there is a first-
order transition to the insulating phase. The phase diagram of the model has
been studied analytically and numerically in [24, 25]. It is worth noting that
the original SYK model with complex fermions also displays a similar phase
diagram [26, 27].

At half filling with $\mu=0$, the Yukawa-SYK model is known to be maximally
chaotic [28]. In this work, we focus on the quantum chaotic properties of the
Yukawa-SYK model in $(0+1)$-d away from half filling as a function of chemical
potential $\mu$ and their relation to the phase transition. It is typical to
quantify the rate at which chaos develops by the Lyapunov exponent
$\lambda_{L}$ [29, 28, 21, 30]. In classical chaotic systems, even slight
differences in initial conditions result in trajectories $q(t)$ that diverge
exponentially: $\delta q(t)/\delta q(0)\propto e^{\lambda_{L}t}$. In quantum
many-body systems, the analogous object is an out-of-time-order correlator
(OTOC) [29] which also grows exponentially and defines a quantum Lyapunov
exponent (correlators of this form had been considered earlier by [31]). We
will compute the following correlator for the fermions:

$\displaystyle
F_{c}(t_{1},t_{2})\sim\Tr[\rho^{\frac{1}{2}}\\{c(t_{1}),c^{\dagger}(0)\\}\rho^{\frac{1}{2}}\\{c(t_{2}),c^{\dagger}(0)\\}^{\dagger}],$
(1)

When $t_{1}=t_{2}=t$, this defines the OTOC, which involves the trace of a
positive semidefinite operator and ensures a real Lyapunov spectrum. The
Lyapunov exponent can then be extracted from the leading exponential growth of
this OTOC, $F_{c}(t_{1},t_{2})\sim\exp[\lambda_{L}(t_{1}+t_{2})/2]$. It has
been argued [32] that in quantum systems at temperature $T$ there is an upper
bound [32] on the Lyapunov exponent $\lambda_{L}\leq 2\pi k_{B}T/\hbar$ (we
set $k_{B}=\hbar=1$ henceforth), which is known to be saturated by black holes
and SYK-like models [5].

The key result of this work is that $\lambda_{L}=2\pi T$ to leading order for
the entire nFL phase; on the other hand, $\lambda_{L}$ is exponentially
suppressed in the insulating phase. As was mentioned, the two phases are
separated by a first-order phase transition, which means there is a hysteresis
regime in $(\mu,T)$-space where both the nFL and insulating solutions are
stable or metastable [25], and the transition occurs somewhere in this region
when the global minimum of the free energy switches from one solution to the
other. From the upper bound for $\lambda_{L}$, we conjecture that the
boundaries of this hysteresis region are curves on which the $\lambda_{L}$ for
one of the phases attains its maximum value. Indeed, we numerically verify
that the curve on which the insulating solution disappears matches very well
with the curve on which $\lambda_{L}$ approaches $2\pi T$ for the insulating
solution at low temperatures. For the nFL solution, since $\lambda_{L}$
universally approaches the upper bound to leading order, one needs to go
beyond the conformal limit to compute the correction to $\lambda_{L}$ and to
verify the conjecture. We leave this to a future study.

The rest of this paper is organized as follows. In Section II we present the
Yukawa-SYK model and review the known results at zero temperature for both the
nFL and insulating states. We then map the zero temperature nFL Green’s
functions to finite temperature using the conformal symmetry of the model, and
we calculate the various Green’s functions that will be required for the
perturbative expansion of the OTOC. In Section III we show how the OTOC can be
computed by summing a series of ladder diagrams and solve the resulting
eigenfunction problem in detail. Finally, in Section IV we discuss the
implications of our results for the OTOC on the phase diagram of the Yukawa-
SYK model and compare with some known numerical results.

##  II Model and preliminaries

The Yukawa-SYK model [17, 18, 19, 20, 23] consists of $M$ flavors of complex
fermions and $N$ flavors of bosons interacting via a Yukawa term with a random
coupling. The bosons have bare mass $m_{0}$; the fermions are dispersionless,
and the fermion density is controlled by the chemical potential $\mu$; see
Refs. [33, 34] for related works on a similar model in string-inspired quiver
quantum mechanics. The Lagrangian of the system in imaginary time is

$\displaystyle L=$
$\displaystyle\sum_{i}^{M}\quantity[c_{i}^{\dagger}(\partial_{\tau}-\mu)c_{i}]+\sum_{\alpha}^{N}\quantity[\frac{1}{2}(\partial_{\tau}\phi_{\alpha})^{2}+\frac{1}{2}m_{0}^{2}\phi_{\alpha}^{2}]$
$\displaystyle+\frac{i}{\sqrt{MN}}\sum_{ij\alpha}t^{\alpha}_{ij}c_{i}^{\dagger}c_{j}\phi_{\alpha}.$
(2)

For each realization of the system, the coupling is drawn from a Gaussian
distribution with $\expectationvalue{t^{\alpha}_{ij}}=0$ and
$\expectationvalue{t^{\alpha}_{ij}t^{\beta}_{kl}}=\omega_{0}^{3}(\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk})\delta_{\alpha\beta}$,
with $\omega_{0}>0$. For the rest of this work, we assume the “weak-coupling
limit” $\omega_{0}\ll m_{0}$, in which, despite the name, the theory remains
non-perturbative. Without loss of generality we can restrict our attention to
$\mu>0$ because of particle-hole symmetry.

This model has already been analyzed at zero temperature in [17, 18, 19, 20,
23] and at finite temperature in [25]. We assume a replica-diagonal solution
and define the fermion and boson self-energies through
$G^{-1}(\omega)=i\omega+\mu+\Sigma(\omega)$ and
$D^{-1}(\Omega)=\Omega^{2}+m_{0}^{2}+\Pi(\Omega)$. The Schwinger-Dyson (SD)
equations in the imaginary time domain are then

$\displaystyle\Sigma(\tau,\tau^{\prime})=-\omega_{0}^{3}D(\tau,\tau^{\prime})G(\tau,\tau^{\prime})$
(3)
$\displaystyle\Pi(\tau,\tau^{\prime})=\frac{M}{N}\omega_{0}^{3}G(\tau,\tau^{\prime})G(\tau^{\prime},\tau)$
(4)

It was shown that in the weak-coupling limit, the only important energy scales
are $\mu$ and the ratio

$\displaystyle\omega_{F}\equiv\omega_{0}^{3}/m_{0}^{2}.$ (5)

For $\mu>\omega_{F}/2$, the SD equation has a (meta)stable solution
corresponding to an incompressible state: the filling $\nu=1$ regardless of
$\mu$, the bosonic self-energy vanishes, and the fermionic self-energy is a
constant. On the other hand, for small $\mu$, the SD equations admit a non-
Fermi liquid solution in which the Green’s functions have power-law forms. We
will compute the OTOC in both of these phases, beginning by deriving the
various finite-temperature Green’s functions that will appear in the expansion
of the OTOC.

###  II.1 Non-Fermi liquid state

At zero temperature, for sufficiently small chemical potential with
$\mu\lesssim\omega_{F}$, a solution of the Schwinger-Dyson equation, which is
stable for $\mu\ll\omega_{F}$ and metastable for $\mu\sim\omega_{F}$, is

$\displaystyle\Sigma(\omega)=-\mu+\omega_{f}^{1-x}|\omega|^{x}(\alpha+i\operatorname{sgn}\omega)\equiv-\mu+\tilde{\Sigma}(\omega)$
(6) $\displaystyle\Pi(\Omega)=-m_{0}^{2}+\beta
m_{0}^{2}\omega_{f}^{2x-1}|\Omega|^{1-2x}\equiv-m_{0}^{2}+\tilde{\Pi}(\Omega)$
(7)

where $\omega_{f}$ (related to but not to be confused with $\omega_{F}$) is a
dynamically generated energy scale below which these solutions are valid. The
parameters $\omega_{f}$, $\alpha$, $\beta$, and $x$ are related to the values
of $\omega_{0}$, $m_{0}$, $\omega_{F}$, $\mu$, and $M/N$. The expressions for
these parameters have been obtained in [23], but, as we shall see, they will
not enter the Lyapunov exponent. In particular we will make use of the
relations

$\displaystyle\frac{M}{N}=\frac{1-2x}{x}\frac{\cos(\pi x)}{\cos(\pi
x)+\frac{1-\alpha^{2}}{1+\alpha^{2}}},$
$\displaystyle\frac{\omega_{0}^{3}}{4\pi\beta
m_{0}^{2}\omega_{f}(1+\alpha^{2})}=\frac{-\Gamma(-2x)}{\Gamma^{2}(-x)}.$ (8)

In fact, the self-energies always have this form at sufficiently low
frequencies regardless of $\omega_{F}$ and $m_{0}$ (even in the strong
coupling limit $\omega_{0}\gg m_{0}$). In the conformal limit in which
$i\omega\ll\Sigma(\omega)$ and $\Omega^{2}\ll\Pi(\Omega)$, the Green’s
functions have the power-law forms $G(\omega)=\tilde{\Sigma}(\omega)^{-1}$ and
$D(\Omega)=\tilde{\Pi}(\Omega)^{-1}$, which can be recast in the time domain
as

$\displaystyle\int\differential{\tau^{\prime}}G(\tau,\tau^{\prime}){\Sigma}(\tau^{\prime},\tau^{\prime\prime})=$
$\displaystyle\;\delta(\tau-\tau^{\prime\prime})$ (9)
$\displaystyle\int\differential{\tau^{\prime}}D(\tau,\tau^{\prime}){\Pi}(\tau^{\prime},\tau^{\prime\prime})=$
$\displaystyle\;\delta(\tau-\tau^{\prime\prime}).$ (10)

We have further dropped the tildes on the self-energies. In doing so, we are
neglecting a contribution on the LHS of (9) which looks like
$\int\differential{\tau^{\prime}}G(\tau,\tau^{\prime})\mu\delta(\tau^{\prime}-\tau^{\prime\prime})$
(and an analogous contribution in (10)). This is permissible because it is
only nonzero at infinitely short time scales, whereas we are interested in the
low-energy (long-time) physics.

In the frequency domain, the Green’s functions are

$\displaystyle G(\omega)=$
$\displaystyle\frac{\omega_{f}^{x-1}}{\alpha^{2}+1}\frac{\alpha-i\operatorname{sgn}\omega}{|\omega|^{x}}$
(11) $\displaystyle D(\Omega)=$
$\displaystyle\beta^{-1}m_{0}^{-2}\omega_{f}^{1-2x}\frac{1}{|\Omega|^{1-2x}}.$
(12)

Using the Fourier transform identities

$\displaystyle\int\frac{\differential{\omega}}{2\pi}e^{i\omega\tau}\frac{1}{|\omega|^{2\Delta}}=$
$\displaystyle\frac{\sin(\pi\Delta)}{\pi}\Gamma(1-2\Delta)|\tau|^{2\Delta-1}$
(13)
$\displaystyle\int\frac{\differential{\omega}}{2\pi}e^{i\omega\tau}\frac{\operatorname{sgn}\omega}{|\omega|^{2\Delta}}=$
$\displaystyle
i\frac{\cos(\pi\Delta)}{\pi}\Gamma(1-2\Delta)|\tau|^{2\Delta-1}\operatorname{sgn}\tau,$
(14)

we find in the time domain

$\displaystyle
G(\tau-\tau^{\prime})=-A_{G}\operatorname{sgn}(\tau-\tau^{\prime})$
$\displaystyle\times\cos(\frac{\pi
x}{2}+\operatorname{sgn}(\tau-\tau^{\prime})\arctan\alpha)\frac{1}{|\tau-\tau^{\prime}|^{1-x}}$
(15) $\displaystyle
D(\tau-\tau^{\prime})=A_{D}\frac{1}{|\tau-\tau^{\prime}|^{2x}}$ (16)

where $A_{G}=\omega_{f}^{x-1}(\alpha^{2}+1)^{-1/2}\pi^{-1}\Gamma(1-x)$ and
$A_{D}\equiv\beta^{-1}m_{0}^{-2}\omega_{f}^{1-2x}\pi^{-1}\cos(\pi
x)\Gamma(2x)$ are constants.

Taken together, equations (3), (4), (9), and (10) are invariant under the
transformation $\tau\to f(\tau)$ and

$\displaystyle
G(\tau,\tau^{\prime})\rightarrow[f^{\prime}(\tau)f^{\prime}(\tau^{\prime})]^{\Delta}\frac{g(\tau)}{g(\tau^{\prime})}G(f(\tau),f(\tau^{\prime}))$
(17)
$\displaystyle\Sigma(\tau,\tau^{\prime})\rightarrow[f^{\prime}(\tau)f^{\prime}(\tau^{\prime})]^{1-\Delta}\frac{g(\tau)}{g(\tau^{\prime})}\Sigma(f(\tau),f(\tau^{\prime}))$
(18) $\displaystyle
D(\tau,\tau^{\prime})\rightarrow[f^{\prime}(\tau)f^{\prime}(\tau^{\prime})]^{1-2\Delta}D(f(\tau),f(\tau^{\prime}))$
(19)
$\displaystyle\Pi(\tau,\tau^{\prime})\rightarrow[f^{\prime}(\tau)f^{\prime}(\tau^{\prime})]^{2\Delta}\Pi(f(\tau),f(\tau^{\prime}))$
(20)

where $f$ and $g$ are arbitrary functions and

$\Delta=\frac{1-x}{2}$ (21)

is a scaling exponent. As discussed in [35], $f$ corresponds to the
reparametrization (conformal) symmetry, and $g$ corresponds to an emergent
symmetry of the complex fermions, which is a $U(1)$ gauge symmetry when
$|g|=1$.

Using the emergent symmetries, we now obtain the finite-temperature,
imaginary-time propagators as well as some other propagators that will be
necessary to compute the OTOC, using the method of [5] and [36]. Focusing
first on the bosons, we make the choice

$\displaystyle f(\tau)=\tan(\frac{\pi\tau}{\beta})$ (22)

in equation (19) to obtain

$\displaystyle
D(\tau-\tau^{\prime})=A_{D}\quantity(\frac{\pi}{\beta|\sin\frac{\pi(\tau-\tau^{\prime})}{\beta}|})^{2x},$
(23)

which solves the Schwinger-Dyson equations, is translationally invariant, and
has the required periodicity in $\beta$. This propagator is defined as
$D(\tau-\tau^{\prime})\equiv\expectationvalue{T\phi(\tau)\phi(\tau^{\prime})}$
where the $\tau$’s are in general complex and, importantly, $T$ orders the
fields according to the real part of $\tau$ (corresponding to imaginary time)
and ignores the imaginary part (corresponding to real time). This means that
the retarded propagator
$-iD^{R}(t-t^{\prime})\equiv\expectationvalue{[\phi(it),\phi(it^{\prime})]}\theta(t-t^{\prime})$
can be built from (23) by adding an infinitesimal real part to the arguments
to ensure the correct ordering:

$\displaystyle
D^{R}(t-t^{\prime})=-i[D(\varepsilon+it,it^{\prime})-D(it,\varepsilon+it^{\prime})]\theta(t-t^{\prime})$
$\displaystyle=-2A_{D}\sin(\pi
x)\quantity(\frac{\pi}{\beta\sinh\frac{\pi(t-t^{\prime})}{\beta}})^{2x}\theta(t-t^{\prime}).$
(24)

We will also need the Wightman propagator which connects two times separated
by half the thermal circle:

$\displaystyle
D^{lr}(t,t^{\prime})\equiv\expectationvalue{\phi(\beta/2+it)\phi(it^{\prime})}$
$\displaystyle=A_{D}\quantity(\frac{\pi}{\beta\cosh\frac{\pi(t-t^{\prime})}{\beta}})^{2x}.$
(25)

Because this propagator is automatically time-ordered, it is simply
$D(\beta/2+it,it^{\prime})$.

Now we do the same for the fermions, applying the transformation (17) with the
aforementioned choice of $f(\tau)$ and using time translation invariance to
set $\tau^{\prime}=0$. We have an additional freedom to choose $g(\tau)$ in
(17). The finite-temperature Green’s function is then, with $g(\tau)$ still
undetermined,

$\displaystyle G(\tau)=-A_{G}g(\tau)\operatorname{sgn}[\tan(\pi\tau/\beta)]$
$\displaystyle\times\cos(\frac{\pi
x}{2}+\operatorname{sgn}[\tan(\pi\tau/\beta)]\arctan\alpha)\quantity(\frac{\pi}{\beta|\sin\frac{\pi\tau}{\beta}|})^{1-x}.$
(26)

To determine $g(\tau)$, we use the short-time information of the Green’s
function which enforces the antiperiodicity condition
$G(\tau-\beta)=-G(\tau)$. The result is the constraint

$\displaystyle\frac{g(\tau-\beta)}{g(\tau)}=\frac{\cos(\frac{\pi
x}{2}+\arctan\alpha)}{\cos(\frac{\pi x}{2}-\arctan\alpha)}$ (27)

which is solved by

$\displaystyle g(\tau)=\left(\frac{\cos(\frac{\pi
x}{2}+\arctan\alpha)}{\cos(\frac{\pi
x}{2}-\arctan\alpha)}\right)^{-\tau/\beta},$ (28)

normalized such that $g(\tau)=1$ at $T=0$.

The propagator (II.1), with $g(\tau)$ given by (28), is
$G(\tau-\tau^{\prime})\equiv\expectationvalue{Tc(\tau)c^{\dagger}(\tau^{\prime})}$,
from which we can build the retarded propagator

$\displaystyle G^{R}(t-t^{\prime})\equiv$
$\displaystyle\expectationvalue{\\{c(it),c^{\dagger}(it^{\prime})\\}}\theta(t-t^{\prime})$
$\displaystyle=$
$\displaystyle[G(\varepsilon+it,it^{\prime})-G(it,\varepsilon+it^{\prime})]\theta(t-t^{\prime})$
$\displaystyle=$ $\displaystyle-A_{G}e^{i\arctan\alpha}\sin(\pi
x)g(i(t-t^{\prime}))$
$\displaystyle\times\quantity(\frac{\pi}{\beta\sinh\frac{\pi(t-t^{\prime})}{\beta}})^{1-x}\theta(t-t^{\prime})$
(29)

and (note the argument – these propagators are defined so that the later time
is always called $t$)

$\displaystyle
G^{A}(t-t^{\prime})\equiv\expectationvalue{\\{c(it^{\prime}),c^{\dagger}(it)\\}}\theta(t-t^{\prime})$
$\displaystyle=[G^{R}(t-t^{\prime})]^{*}.$ (30)

The Wightman propagators are given by

$\displaystyle
G^{r}(t-t^{\prime})\equiv\expectationvalue{c(\beta/2+it)c^{\dagger}(it^{\prime})}=G(\beta/2+it,it^{\prime})$
$\displaystyle=$ $\displaystyle-A_{G}\sqrt{\cos(\frac{\pi
x}{2}-\arctan\alpha)\cos(\frac{\pi x}{2}+\arctan\alpha)}$ $\displaystyle\times
g(i(t-t^{\prime}))\quantity(\frac{\pi}{\beta\cosh\frac{\pi(t-t^{\prime})}{\beta}})^{1-x}$
(31)

and

$\displaystyle
G^{l}(t-t^{\prime})\equiv\expectationvalue{c(it^{\prime})c^{\dagger}(\beta/2+it)}=[G^{r}(t-t^{\prime})]^{*}.$
(32)

###  II.2 Insulating state

At zero temperature and in the weak coupling limit $\omega_{0}\ll m_{0}$, the
system has a (meta)stable insulating phase for $\mu>\omega_{F}/2$ [23]. In
this phase, the Green’s functions are simply those of free bosons (with the
mass unchanged) and free fermions with a renormalized chemical potential

$\displaystyle\tilde{\mu}\equiv\mu-\frac{\omega_{F}}{2}>0.$ (33)

Fourier transforming to imaginary time yields

$\displaystyle G(\tau)=e^{-\tilde{\mu}\tau}\theta(\tau),\quad\quad
D(\tau)=\frac{e^{-m_{0}|\tau|}}{2m_{0}}.$ (34)

Because of the gapped nature of the phase, the Green’s functions are
approximately independent of temperature as long as $T\ll\tilde{\mu},m_{0}$.
The retarded and advanced Green’s functions are constructed as before (this
can also be done from the frequency domain with the usual analytic
continuation $i\omega\rightarrow\omega\pm i\delta$):

$\displaystyle G^{R,A}(t-t^{\prime})=e^{\mp
i\tilde{\mu}(t-t^{\prime})}\theta(t-t^{\prime}),$ $\displaystyle
D^{R}(t-t^{\prime})=\frac{\sin[m_{0}(t-t^{\prime})]}{m_{0}}\theta(t-t^{\prime}).$
(35)

Finally, we obtain the right- and left-pointing Wightman propagators

$\displaystyle G^{r,l}(t-t^{\prime})=e^{-\tilde{\mu}[\frac{\beta}{2}\pm
i(t-t^{\prime})]},$ $\displaystyle
D^{r,l}(t-t^{\prime})=\frac{1}{2m_{0}}e^{-m_{0}[\frac{\beta}{2}\pm
i(t-t^{\prime})]}.$ (36)

However, the Wightman propagators for the bosons contain a factor of
$\exp(-m_{0}\beta/2)$ and are exponentially suppressed, since we assume
$m_{0}\gg T$. (The fermion propagators are not similarly suppressed because we
make no assumption about the ratio $\tilde{\mu}/T$.) This means that we can
approximate $D^{r,l}\approx 0$, and we will entirely neglect the diagrams in
the ladder series with this type of rung.

##  III Solution of the ladder equation

For the Yukawa-SYK model, the OTOC in Eq. (1) is precisely defined as

$F_{c}(t_{1},t_{2})\\!=\\!\sum_{ij}\frac{\Tr}{M^{2}}[\rho^{\frac{1}{2}}\\{c_{i}(t_{1}),c_{j}^{\dagger}(0)\\}\rho^{\frac{1}{2}}\\{c_{i}(t_{2}),c_{j}^{\dagger}(0)\\}^{\dagger}],$
(37)

We evaluate the OTOC by the usual technique which is described in [36] (see
also [29, 37]) via a path integral defined on a complex time (Keldysh)
contour, which encircles the compactified imaginary time direction and also
has two real time folds (“rails”) at $\tau=0$ and $\tau=\beta/2$ (see Fig. 1).
Two operators are placed at the two ends ($\imaginary\tau=0$ and
$\imaginary\tau=t_{1,2}$) of each rail. Since imaginary time is strictly
increasing along this contour, contour-ordering is the same as $\tau$-ordering
and we can use the $\tau$-ordered Green’s functions of the previous section.
The four possible orderings of the operators in the double-anticommutator are
accounted for by displacing the operators at $t=0$ infinitesimally forward or
backward in the imaginary time direction so that they lie just after or just
before the real time folds; summing over the four different contour-orderings
yields the desired correlator.

Figure 1: The Keldysh contour used for the expansion of the OTOC, with two
real-time folds (note that real-time evolution corresponds to the imaginary
part of $\tau$). The real part of $\tau$ has period $\beta$; the ends of the
contour should be identified. The arrow indicates the direction of contour-
ordering. One ordering of the operators, from the expansion of the
anticommutators, is depicted.

At zeroth order in the interaction, each of the orderings can be computed by
contracting the operators using Wick’s theorem in the path integral along the
complex contour. The terms with contractions between operators across
different rails cancel after summing over different orderings, and the
survivors conspire to produce a product of the (free) retarded and advanced
Green’s functions.

In the interaction picture, we expand the interaction-induced time evolution
operator perturbatively, inserting copies of the interaction vertex along the
real time folds and using Wick’s theorem. Insertions along the imaginary time
axis, on the other hand, cancel upon summing over orderings. When operators on
the same rail are contracted, the resulting contributions are just part of the
self-energy, and these are accounted for by using the fully dressed
propagators rather than the non-interacting ones. On the contrary, the
contractions within interaction vertices between opposite rails are not
captured by the self-energy. Taking into account these diagrams (with
insertions on both sides of each rail), to leading order in $M$ and $N$, we
find that the OTOC can be expressed as a series of uncrossed ladder diagrams,
which can be written recursively as in Fig. 2. The fact that only ladder
diagrams need to be included at leading order in large-$N$ is a common theme
in related models and has been shown, for example, for the SYK model [5] and
for a weak coupling $\phi^{4}$ theory [29]. For the details of the large-$M,N$
analysis for the current model, see Appendix A. According to a modified
Feynman rule (cf. [29]), the vertical lines are retarded (advanced)
propagators and the horizontal lines are Wightman propagators. For
completeness, in Appendix B, we prove the Feynman rule using the lowest order
diagram of the ladder series.

The goal is to compute the correlator on the LHS of the upper equation in Fig.
2. As is indicated by the last diagram on the RHS, its full expression
involves a second correlator, defined as

$F_{\bar{c}}=\frac{1}{M^{2}}\sum_{ij}\Tr[\rho^{\frac{1}{2}}\\{c_{i}^{\dagger}(t_{1}),c_{j}^{\dagger}(0)\\}\rho^{\frac{1}{2}}\\{c_{i}^{\dagger}(t_{2}),c_{j}^{\dagger}(0)\\}^{\dagger}],$
(38)

given by the lower equation, to account for the two possible directions of
fermion loops (a complication that does not arise in models with Majorana
fermions). Even though $F_{\bar{c}}=0$ for free fermions, it is generated by
interaction effects at second order in perturbation theory. To see this, one
can plug the first term on the right hand side of the upper equation in Fig. 2
to the third term of the right hand side of the lower equation. More detail on
how this diagram follows from the Feynman rules can be found in Appendix C. In
turn, $F_{\bar{c}}$ contributes to $F_{c}$ via the last diagram on the RHS of
the upper equation in Fig. 2. Following the Feynman rules in a similar way, we
show the leading order contribution at fourth order in Appendix C.

As we shall see, including these diagrams is important for getting the correct
Lyapunov exponent. We note that there is some variation in the literature as
to whether this type of diagram is included, which we leave for future
investigations. In agreement with Refs. [21] and [38], we find that for
complex fermions such diagrams are indeed generated.

Figure 2: Feynman diagrams for the fermionic OTOC to leading order in $1/N$
and $1/M$. All propagators are fully renormalized by self-energy contributions
to leading order in $1/N$ and $1/M$. The disorder average involves averaging
products of independent Gaussian variables, which can be reduced by Isserlis’s
theorem to a sum of products of averages of all pairings of the variables.
Hence the disorder average, indicated by dashed lines, acts like another field
with a constant propagator which requires the fermion and boson flavors at one
end to match those at the other end. Diagrams with crossed rungs, and other
types of pairings in the disorder average, are suppressed in the large $M,N$
limit.

We express the ladder equation for the OTOC in terms of integral kernels

$\displaystyle K_{cc}$
$\displaystyle=\omega_{0}^{3}G^{R}(t_{13})G^{A}(t_{24})D^{lr}(t_{43})$
$\displaystyle K_{\phi c}$
$\displaystyle=\omega_{0}^{3}D^{R}(t_{13})D^{R}(t_{24})G^{r}(t_{43})$
$\displaystyle K_{c\phi}$
$\displaystyle=\frac{M}{N}\omega_{0}^{3}G^{R}(t_{13})G^{A}(t_{24})G^{l}(t_{43})$
$\displaystyle K_{\bar{c}\bar{c}}$
$\displaystyle=\omega_{0}^{3}G^{A}(t_{13})G^{R}(t_{24})D^{lr}(t_{43})$
$\displaystyle K_{\phi\bar{c}}$
$\displaystyle=\omega_{0}^{3}D^{R}(t_{13})D^{R}(t_{24})G^{l}(t_{43})$
$\displaystyle K_{\bar{c}\phi}$
$\displaystyle=\frac{M}{N}\omega_{0}^{3}G^{A}(t_{13})G^{R}(t_{24})G^{r}(t_{43})$
(39)

where $t_{ij}\equiv t_{i}-t_{j}$. The notation $K_{\phi c}$, for example,
indicates that it corresponds to portions of the diagrams which create bosons
at later times out of incoming fermions. The $\bar{c}$ refers to fermion lines
with arrows that are reversed compared to the OTOC (that is, an advanced
propagator on the left rail and a retarded one on the right). The kernels
involving $\bar{c}$ are simply the complex conjugates of the the corresponding
kernels with $c$. Denoting the fermion OTOC by $F_{c}(t_{1},t_{2})$, the
ladder equation is then

$\displaystyle F_{c}=$ $\displaystyle\frac{1}{M}G^{R}(t_{1})G^{A}(t_{2})$
$\displaystyle+K_{cc}\star F_{c}+K_{c\phi}\star K_{\phi c}\star
F_{c}+K_{c\phi}\star K_{\phi\bar{c}}\star F_{\bar{c}}$ $\displaystyle
F_{\bar{c}}=$ $\displaystyle K_{\bar{c}\bar{c}}\star
F_{\bar{c}}+K_{\bar{c}\phi}\star K_{\phi\bar{c}}\star
F_{\bar{c}}+K_{\bar{c}\phi}\star K_{\phi c}\star F_{c}$ (40)

The $\star$ denotes an integral convolution, that is,

$\displaystyle(K\star
F)(t_{1},t_{2})\equiv\int_{-\infty}^{\infty}\differential{t_{3}}\differential{t_{4}}K(t_{1},t_{2},t_{3},t_{4})F(t_{3},t_{4}).$
(41)

The inhomogeneous first term $G^{R}(t_{1})G^{A}(t_{2})/M$ in the series does
not increase exponentially like $F_{c}$ and quickly becomes negligible. The
inhomogeneous term acts as a source and ensures that the exponential growth is
suppressed by $1/M$, i.e., $F_{c}\propto e^{\lambda_{L}t}/M$. The value of
$\lambda_{L}$ can be determined by dropping the inhomogeneous term. After
doing this, the two equations in Fig. 2 can be solved by
$F_{c}=F_{\bar{c}}^{*}$, so that one only needs to solve one of them.

The strategy for solving this eigenfunction problem is to find functions
$f_{c}(t_{1},t_{2})$ (possibly complex) and $f_{\phi}(t_{1},t_{2})$ (real)
which satisfy

$\displaystyle K_{cc}\star f_{c}=k_{cc}f_{c},~{}~{}$ $\displaystyle
K_{\bar{c}\bar{c}}\star f_{c}^{*}=k_{cc}f_{c}^{*}$ $\displaystyle K_{\phi
c}\star f_{c}=k_{\phi c}f_{\phi},~{}~{}$ $\displaystyle K_{\phi\bar{c}}\star
f_{c}^{*}=k_{\phi c}f_{\phi}$ $\displaystyle K_{c\phi}\star
f_{\phi}=k_{c\phi}f_{c},~{}~{}$ $\displaystyle K_{\bar{c}\phi}\star
f_{\phi}=k_{c\phi}f_{c}^{*},$ (42)

where $k_{cc},k_{\phi c},k_{c\phi}$ are real constants. If we can find such
functions, then we can choose $F_{c}\sim f_{c}/M$ and $F_{\bar{c}}\sim
f_{c}^{*}/M$ in (III). Both equations reduce to the same algebraic requirement
that

$k_{cc}+2k_{c\phi}k_{\phi c}=1,$ (43)

where the factor of 2 comes from considering both $F_{c}$ and $F_{\bar{c}}$ on
the right hand side of the upper equation of Fig. 2.

###  III.1 Non-Fermi liquid state

We begin by writing out $K_{cc}$ in full:

$\displaystyle
K_{cc}=C_{cc}e^{-it_{13}L/\beta}\quantity(\frac{\pi}{\beta\sinh\frac{\pi
t_{13}}{\beta}})^{1-x}\theta(t_{13})e^{+it_{24}L/\beta}\quantity(\frac{\pi}{\beta\sinh\frac{\pi
t_{24}}{\beta}})^{1-x}\theta(t_{24})\quantity(\frac{\pi}{\beta\cosh\frac{\pi
t_{34}}{\beta}})^{2x}$ (44)

where

$C_{cc}\equiv\omega_{0}^{3}A_{G}^{2}A_{D}\sin^{2}(\pi x),\textrm{ and
}L=\log(\frac{\cos(\frac{\pi x}{2}+\arctan\alpha)}{\cos(\frac{\pi
x}{2}-\arctan\alpha)}).$ (45)

We first make the change of variable $t_{i}\equiv\frac{\beta}{2\pi}\phi_{i}$
and use the step functions to restrict the region of integration. The first
equation in (III) becomes

$\displaystyle\int_{-\infty}^{\phi_{1}}\differential{\phi_{3}}\int_{-\infty}^{\phi_{2}}\differential{\phi_{4}}C_{cc}e^{-i\phi_{13}L/2\pi}e^{+i\phi_{24}L/2\pi}\quantity(\frac{1}{2\sinh\frac{\phi_{13}}{2}})^{1-x}\quantity(\frac{1}{2\sinh\frac{\phi_{24}}{2}})^{1-x}\quantity(\frac{1}{2\cosh\frac{\phi_{34}}{2}})^{2x}f_{c}(\phi_{3},\phi_{4})=k_{cc}f_{c}(\phi_{1},\phi_{2})$
(46)

Now we perform another change of variable: $z_{1,3}\equiv
e^{-\phi_{1,3}},z_{2,4}\equiv-e^{-\phi_{2,4}}$ (note the minus sign). The
kernel equation becomes

$\displaystyle
C_{cc}\int_{z_{1}}^{\infty}\frac{\differential{z_{3}}}{|z_{3}|}\int_{-\infty}^{z_{2}}\frac{\differential{z_{4}}}{|z_{4}|}\quantity|\frac{z_{1}z_{4}}{z_{2}z_{3}}|^{\frac{iL}{2\pi}}\frac{|z_{1}z_{3}|^{\frac{1}{2}(1-x)}}{|z_{13}|^{1-x}}\frac{|z_{2}z_{4}|^{\frac{1}{2}(1-x)}}{|z_{24}|^{1-x}}\frac{|z_{3}z_{4}|^{x}}{|z_{34}|^{2x}}f_{c}(z_{3},z_{4})=k_{cc}f_{c}(z_{1},z_{2}),$
(47)

where $z_{ij}\equiv z_{i}-z_{j}$. We solve this by making the ansatz

$\displaystyle
f_{c}(z_{1},z_{2})=\quantity|\frac{z_{1}}{z_{2}}|^{\frac{iL}{2\pi}}\frac{|z_{1}z_{2}|^{\frac{1}{2}(1-x)}}{|z_{12}|^{1-x-h}}$
(48)

which will make the integrand only a function of the difference between the
$z$’s. With this choice the integral in (47) becomes

$\displaystyle
C_{cc}\quantity(\int_{z_{1}}^{\infty}\\!\\!\\!\\!\\!\differential{z_{3}}\int_{-\infty}^{z_{2}}\\!\\!\\!\\!\\!\differential{z_{4}}\frac{|z_{12}|^{1-x-h}}{|z_{13}|^{1-x}|z_{24}|^{1-x}|z_{34}|^{1+x-h}})f_{c}(z_{1},z_{2}).$
(49)

Here $h$ is a parameter that can be later tuned to satisfy (43). As we shall
see, it is directly related to $\lambda_{L}$. The expression in parentheses,
as promised, is simply a constant, despite appearing to be a function of
$z_{1}$ and $z_{2}$. It is invariant under the affine transformation
$z_{i}\rightarrow az_{i}+b$; in particular the integration variables can be
shifted and scaled to set $z_{1}=1,z_{2}=0$ and the integral can be evaluated:

$\displaystyle k_{cc}=$ $\displaystyle
C_{cc}\int_{1}^{\infty}\differential{z_{3}}\int_{-\infty}^{0}\differential{z_{4}}\frac{1}{|1-z_{3}|^{1-x}|z_{4}|^{1-x}|z_{34}|^{1+x-h}}$
$\displaystyle=$ $\displaystyle
C_{cc}\frac{\Gamma^{2}(x)\Gamma(1-x-h)}{\Gamma(1+x-h)}.$ (50)

Now we repeat this procedure for the second equation of (III), feeding in our
expression for $f_{c}$. Upon making the additional ansatz

$\displaystyle
f_{\phi}(z_{1},z_{2})=\frac{|z_{1}z_{2}|^{x}}{|z_{12}|^{2x-h}},$ (51)

again to make the integrand dependent only on the $z_{ij}$’s, we find

$\displaystyle k_{\phi c}$ $\displaystyle=C_{\phi
c}\int_{1}^{\infty}\differential{z_{3}}\int_{-\infty}^{0}\differential{z_{4}}\frac{1}{|1-z_{3}|^{2x}|z_{4}|^{2x}|z_{34}|^{2-2x-h}}$
$\displaystyle=C_{\phi c}\frac{\Gamma^{2}(1-2x)\Gamma(2x-h)}{\Gamma(2-2x-h)}$
(52)

where

$\displaystyle C_{\phi c}\equiv-4\omega_{0}^{3}A_{G}A_{D}^{2}\sin^{2}(\pi
x)(2\pi/\beta)^{3x-1}$ $\displaystyle\times\sqrt{\cos(\frac{\pi
x}{2}-\arctan\alpha)\cos(\frac{\pi x}{2}+\arctan\alpha)}.$ (53)

The final two equations of (III) are compatible with the same choice of
$f_{c}$ and $f_{\phi}$, and the result is

$\displaystyle k_{c\phi}$
$\displaystyle=C_{c\phi}\int_{1}^{\infty}\differential{z_{3}}\int_{-\infty}^{0}\differential{z_{4}}\frac{1}{|1-z_{3}|^{1-x}|z_{4}|^{1-x}|z_{34}|^{1+x-h}}$
$\displaystyle=C_{c\phi}\frac{\Gamma^{2}(x)\Gamma(1-x-h)}{\Gamma(1+x-h)}$ (54)

where

$\displaystyle C_{c\phi}\equiv-\frac{M}{N}\omega_{0}^{3}A_{G}^{3}\sin^{2}(\pi
x)(2\pi/\beta)^{1-3x}$ $\displaystyle\times\sqrt{\cos(\frac{\pi
x}{2}-\arctan\alpha)\cos(\frac{\pi x}{2}+\arctan\alpha)}.$ (55)

With these results, Eq. (43) reduces to an algebraic equation for $h$,

$\displaystyle\frac{\Gamma^{2}(x)\Gamma(1-x-h)}{\Gamma(1+x-h)}$
$\displaystyle\times\quantity[C_{cc}+2C_{c\phi}C_{\phi
c}\frac{\Gamma^{2}(1-2x)\Gamma(2x-h)}{\Gamma(2-2x-h)}]=1.$ (56)

The constants $C_{cc},C_{c\phi},C_{\phi c}$ involve a complicated mix of
$\alpha,\beta,\omega_{0},m_{0},\omega_{f},\frac{n_{s}M}{N},\mbox{ and }x$. We
can use the relations (II.1), along with the gamma function reflection
identity and the fact $\cos(\pi x/2-\arctan\alpha)\cos(\pi
x/2+\arctan\alpha)=\frac{1}{2}(\cos(\pi
x)+\frac{1-\alpha^{2}}{1+\alpha^{2}})$, to write (III.1) purely in terms of
$x$:

$\displaystyle
1=\frac{-1}{\Gamma(-x)\Gamma(x)}\frac{\Gamma^{2}(x)\Gamma(1-x-h)}{\Gamma(1+x-h)}$
$\displaystyle\times\quantity[1+2(1-2x)\frac{1}{\Gamma(2x)\Gamma(1-2x)}\frac{\Gamma^{2}(1-2x)\Gamma(2x-h)}{\Gamma(2-2x-h)}].$
(57)

Remarkably, this is exactly solved for all $x$ by

$h=-1.$ (58)

This can be verified analytically by making repeated use of the identity
$\Gamma(z+1)=z\Gamma(z)$. Transforming back to the original time variables, we
have

$\displaystyle
f_{c}(t_{1},t_{2})=e^{\frac{-iL}{\beta}(t_{1}-t_{2})}\frac{e^{-h\frac{2\pi}{\beta}\quantity(\frac{t_{1}+t_{2}}{2})}}{(2\cosh{\frac{2\pi}{\beta}\quantity(\frac{t_{1}-t_{2}}{2})})^{1-x-h}}.$
(59)

The OTOC grows in the same way as $f_{c}(t,t)$ at long times, and the Lyapunov
exponent can be read off from this expression:

$\lambda_{L}=\frac{-2\pi}{\beta}h={2\pi}T.$ (60)

Notice that the only effect of the chemical potential, which enters through
$L$, is to introduce an oscillatory factor in $f_{c}$ which does not affect
the OTOC or $\lambda_{L}$. This means that in any regime where the conformal
limit is valid, the spectral asymmetry does not disturb the OTOC from its form
at half-filling, which corresponds to $\mu=\alpha=L=0$. Furthermore,
$\lambda_{L}$ is independent of the ratio $M/N$.

Figure 3: Feynman diagrams for the bosonic OTOC. Notation is the same as in
Fig. 2.

Interestingly, by following the argument in [37], we can conclude that the
bosons are maximally chaotic as well with little additional computation, since
the boson ladder equation involves the same kernels as for the fermions. The
boson correlator $F_{\phi}(t_{1},t_{2})$ is defined as

$F_{\phi}(t_{1},t_{2})\\!=\\!\\!\sum_{\alpha\beta}\\!\frac{\Tr}{N^{2}}[\rho^{\frac{1}{2}}[\phi_{\alpha}(t_{1}),\phi_{\beta}(0)]\rho^{\frac{1}{2}}[\phi_{\alpha}(t_{2}),\phi_{\beta}(0)]^{\dagger}],$
(61)

and can be expressed diagrammatically via the coupled equations of Fig. 3 by
introducing the correlators $F_{c\phi}$ and $F_{\bar{c}\phi}$, which involve
the processes that convert two incoming bosons to two outgoing fermions.
Algebraically, the equations read

$\displaystyle F_{\phi}=F_{\phi}^{0}+K_{\phi c}\star
F_{c\phi}+K_{\phi\bar{c}}\star F_{\bar{c}\phi}$ $\displaystyle
F_{c\phi}=K_{cc}\star F_{c\phi}+K_{c\phi}\star F_{\phi}$ $\displaystyle
F_{\bar{c}\phi}=K_{\bar{c}\bar{c}}\star F_{\bar{c}\phi}+K_{\bar{c}\phi}\star
F_{\phi},$ (62)

where $F_{c\phi}$, $F_{\bar{c}\phi}$, and etc., represented diagrammatically
in Fig. 3, are “mixed” OTOC’s defined as, e.g.,

$F_{c\phi}(t_{1},t_{2})\\!=\\!\\!\sum_{i\alpha}\\!\frac{\Tr}{MN}\quantity[\rho^{\frac{1}{2}}[c_{i}(t_{1}),\phi_{\alpha}(0)]\rho^{\frac{1}{2}}[c_{i}(t_{2}),\phi_{\alpha}(0)]^{\dagger}].$
(63)

We neglect the inhomogeneous term in (III.1) as usual and substitute the first
equation into the second and third to eliminate $F_{\phi}$. We see that
$F_{c\phi}$ and $F_{\bar{c}\phi}$ obey precisely the same ladder equation as
$F_{c}$ and $F_{\bar{c}}$ respectively (cf. Eq. (III)), and thus they too are
proportional to $f_{c}$ and $f_{c}^{*}$ with the same exponent $h=-1$. Hence
the mixed OTOCs such as $F_{c\phi}$ also have the same Lyapunov exponent.

Substituting these solutions back into the first equation of (III.1) and using
the properties of the kernels (III), we find that the boson OTOC $F_{\phi}$ is
proportional to $f_{\phi}$, which in the time domain is

$\displaystyle
f_{\phi}(t_{1},t_{2})=\frac{e^{-h\frac{2\pi}{\beta}\quantity(\frac{t_{1}+t_{2}}{2})}}{(2\cosh{\frac{2\pi}{\beta}\quantity(\frac{t_{1}-t_{2}}{2})})^{2x-h}}.$
(64)

This has exactly the same long-time growth and Lyapunov exponent as the
fermion OTOC, only without the oscillation in $t_{1}-t_{2}$.

Intuitively, since the system is strongly interacting among both boson and
fermions, the exponential growth of chaos does not depend on which field
operators are used to perturb and probe the system.

###  III.2 Insulating state

We repeat the procedure for the (meta)stable insulating state. As mentioned
earlier (cf. Eq. (II.2)), the boson left-right propagators are exponentially
suppressed, so we drop the inhomogeneous terms, containing $K_{cc}$ and
$K_{\bar{c}\bar{c}}$, on the right-hand sides of (III):

$\displaystyle F_{c}=K_{c\phi}\star K_{\phi c}\star F_{c}+K_{c\phi}\star
K_{\phi\bar{c}}\star F_{\bar{c}}$ $\displaystyle
F_{\bar{c}}=K_{\bar{c}\phi}\star K_{\phi\bar{c}}\star
F_{\bar{c}}+K_{\bar{c}\phi}\star K_{\phi c}\star F_{c}.$ (65)

The kernels have the same structure as before, but the Green’s functions are
now the insulating ones. Again we solve this by finding $f_{c}$ and $f_{\phi}$
satisfying (III).

$\displaystyle K_{\phi c}\star f_{c}=k_{\phi c}f_{\phi}$ $\displaystyle
K_{c\phi}\star f_{\phi}=k_{c\phi}f_{c}$ (66)

Then the choice $F_{c}=f_{c},F_{\bar{c}}=f_{c}^{*}$ solves the insulating
ladder equation when

$1=2k_{c\phi}k_{\phi c}.$ (67)

By the same argument as in the previous subsection, the boson OTOC $F_{\phi}$
will again be proportional to the function $f_{\phi}$. Since we will see that
$f_{c}$ and $f_{\phi}$ both have the same exponential growth, we conclude that
the fermions and bosons have the same Lyapunov exponent in the insulating
phase.

We make the ansatz

$\displaystyle
f_{c}(t_{1},t_{2})=e^{\frac{\lambda_{L}}{2}(t_{1}+t_{2})-i\gamma(t_{1}-t_{2})}$
(68)

where $\lambda_{L}$ and $\gamma$ are real, and $\lambda_{L}$ is of course the
Lyapunov exponent. Since we obtain the OTOC, which must be real, by setting
$t_{1}=t_{2}$, $\gamma$ will drop out of the final result, but it is
nonetheless necessary to keep it as a parameter. As will be shown, in the
absence of the inhomogeneous term, $\lambda_{L}$ and $\gamma$ parametrize a
whole family of functions with the properties (III), and the requirement
$1=2k_{c\phi}k_{\phi c}$ constrains the solutions for $\lambda_{L}$ and
$\gamma$. The multiplicity of solutions is an artifact of neglecting the
inhomogeneous term in the ladder equation (the OTOC should be unique), but
physically the ladder equation represents the amplification of $F_{c}^{0}$ by
repeated applications of a kernel, and only the fastest-growing component will
survive this procedure. Therefore we will tune $\gamma$ to make $\lambda_{L}$
as large as possible.

For the first equation in (III.2), we first calculate the action of $K_{\phi
c}$ on our ansatz, using Eqs. (III), (II.2), (II.2), and (68) the result is

$\displaystyle K_{\phi c}\star f_{c}=$
$\displaystyle\omega_{0}^{3}\int_{-\infty}^{t_{1}}\differential{t_{3}}\int_{-\infty}^{t_{2}}\differential{t_{4}}\frac{\sin
m_{0}t_{13}}{m_{0}}\frac{\sin m_{0}t_{24}}{m_{0}}\times
e^{-\tilde{\mu}(\beta/2+it_{43})}e^{\frac{\lambda_{L}}{2}(t_{3}+t_{4})-i\gamma
t_{34}}$ $\displaystyle=$
$\displaystyle\frac{\omega_{0}^{3}}{m_{0}^{2}}e^{-\tilde{\mu}\beta/2}\frac{1}{(2i)^{2}}\sum_{\sigma,\sigma^{\prime}=\pm}\frac{\sigma\sigma^{\prime}}{(-i\sigma
m_{0}+i\tilde{\mu}+\lambda_{L}/2-i\gamma)(-i\sigma^{\prime}m_{0}-i\tilde{\mu}+\lambda_{L}/2+i\gamma)}e^{-i\tilde{\mu}t_{12}}e^{\frac{\lambda_{L}}{2}(t_{1}+t_{2})-i\gamma
t_{12}}$ (69)

from which we can read off $k_{\phi c}$, and we find that

$f_{\phi}(t_{1},t_{2})=e^{+i\tilde{\mu}t_{12}}e^{\frac{\lambda_{L}}{2}(t_{1}+t_{2})-i\gamma
t_{12}}.$ (70)

Now we act on this with $K_{c\phi}$:

$\displaystyle K_{c\phi}\star
f_{\phi}=\frac{M}{N}\omega_{0}^{3}\int_{-\infty}^{t_{1}}\differential{t_{3}}\int_{-\infty}^{t_{2}}\differential{t_{4}}e^{-i\tilde{\mu}t_{13}}e^{i\tilde{\mu}t_{24}}e^{-\tilde{\mu}(\beta/2-it_{43})}e^{+i\tilde{\mu}t_{34}}e^{\frac{\lambda_{L}}{2}(t_{3}+t_{4})-i\gamma
t_{34}}$
$\displaystyle=\frac{M}{N}\omega_{0}^{3}e^{-\tilde{\mu}\beta/2}\frac{1}{(i\tilde{\mu}+\lambda_{L}/2-i\gamma)(-i\tilde{\mu}+\lambda_{L}/2+i\gamma)}e^{\frac{\lambda_{L}}{2}(t_{1}+t_{2})-i\gamma
t_{12}},$ (71)

We recover $f_{c}$ as required, and the prefactor is $k_{c\phi}$. Combining
the results of (III.2) and (III.2) with the requirement $1=2k_{c\phi}k_{\phi
c}$, we find that the ansatz solves the ladder equation when

$\displaystyle
1=\frac{2M}{N}\frac{\omega_{0}^{6}}{m_{0}^{2}}e^{-\tilde{\mu}\beta}\quantity(\frac{1}{2i})^{2}\sum_{\sigma\sigma^{\prime}}\frac{\sigma\sigma^{\prime}}{(-i\sigma
m_{0}+i\tilde{\mu}+\lambda_{L}/2-i\gamma)(-i\sigma^{\prime}m_{0}-i\tilde{\mu}+\lambda_{L}/2+i\gamma)(i\tilde{\mu}+\lambda_{L}/2-i\gamma)(-i\tilde{\mu}+\lambda_{L}/2+i\gamma)}.$
(72)

Since we have assumed the boson mass is the largest energy scale in the
problem and since $\lambda_{L}$ cannot exceed $2\pi T$, we can neglect
$\tilde{\mu}$ and $\lambda_{L}$ relative to $m_{0}$ in the denominator of Eq.
(72). We will further assume and then verify that $\gamma$ is similarly
negligible. With these assumptions the eigenvalue equation reduces to

$\displaystyle\frac{2M}{N}\omega_{F}^{2}e^{-\tilde{\mu}\beta}=\quantity(\frac{\lambda_{L}}{2})^{2}+(\gamma-\tilde{\mu})^{2}.$
(73)

Clearly the optimal choice is $\gamma=\tilde{\mu}$, consistent with the
assumption $m_{0}\gg\gamma$; any other choice drags the Lyapunov exponent
down, so

$\displaystyle\lambda_{L}(\tilde{\mu},T)=\sqrt{\frac{8M}{N}}\omega_{F}e^{-\tilde{\mu}/2T}\textrm{,
for $T\ll\tilde{\mu}$}.$ (74)

Note that this choice also ensures that the bosonic OTOC is real, just like in
the nFL case. We see that for the insulating state $\lambda_{L}$ is
exponentially suppressed, consistent with the understanding that the system is
in a trivial gapped state.

Observe that the Lyapunov exponent is an increasing function of temperature
and a decreasing function of chemical potential. This is consistent with the
intuition that chaos is suppressed as the system goes deeper into the
insulating phase, where the filling is close to unity and there is little
phase space available for scrambling.

##  IV Discussion and conclusion

In this work, we generalized the Green’s functions of the Yukawa-SYK model at
weak coupling to finite temperatures and computed the OTOCs and Lyapunov
exponents of the nFL state in the conformal limit and the insulating state.
The fermions are maximally chaotic in the nFL phase, even away from half-
filling. The chemical potential merely gives a phase to the correlator
$F_{c}(t_{1},t_{2})$ which does not affect the OTOC or the Lyapunov exponent.
The bosons also saturate the chaos bound, as do the other “mixed” correlators
described above such as $F_{c\phi}$, etc.

For the insulating state, we have obtained an expression for $\lambda_{L}$ at
low temperature in Eq. (74). As the temperature increases, $\lambda_{L}$ for
the insulating solution increases. As $\lambda_{L}$ has an upper bound of
$2\pi T$, it is interesting to see where in the $(\mu,T)$ plane $\lambda_{L}$
for the insulating state approaches this value. Heuristically this should
correspond to the point when the metastable insulating state becomes unstable.
To determine the boundary in the $(\tilde{\mu},T)$-plane, we set
$\lambda_{L}(T)\sim 2\pi T$ and solve for the critical value $\mu_{c}$ as a
function of temperature in the low-temperature limit:

$\displaystyle\mu_{c}(T)=\frac{\omega_{F}}{2}+2T\log(\sqrt{\frac{M}{N}}\frac{\omega_{F}}{T}).$
(75)

We expect Eq. (74) to be correct for
$T\ll\tilde{\mu}\equiv\mu_{c}-{\omega_{F}}/{2}$, so $T\ll\omega_{F}$. This
curve is plotted in blue in Fig. 4.

Figure 4: The Lyapunov exponent, calculated from the insulating side in the
low-$T$ limit, saturates the chaos bound along the solid blue curves. The red
curves with data points come from numerical results [25] for the leftmost
boundary of the hysteresis region. For completeness, the numerical results for
the rightmost boundary of the hysteresis region are also shown in orange,
though we did not attempt to find this curve analytically. The $T=0$ data
points are analytic results – see [23]. The upper and lower panel have $M/N=1$
and $M/N=4$ respectively, and both have $\omega_{F}=1/4$. The agreement is
good at low $T$ but, unsurprisingly, fails when $T\sim\omega_{F}$, noticeable
in the lower panel.

To the left of this curve, the exponent violates this bound, so the insulating
state cannot be stable, demarcating the leftmost boundary of the hysteresis
region (to the right of this curve, we cannot say that the insulating state is
energetically favorable to the nFL state, only that it is not forbidden by the
chaos bound).

In [25], the nFL and insulating states were placed on a phase diagram as a
function of $(\mu,T)$, which we reproduce for the present model in Fig. 4. The
data were obtained by iteratively solving the Schwinger-Dyson equations as
described therein. (The index structure of the random coupling is slightly
different in the model of [25], but the Schwinger-Dyson equations have
essentially the same structure. The only change is that $M/N$ in this paper
corresponds to $4M/N$ in [25].) The two phases are separated by a first-order
transition (not shown), located in the hysteresis region bounded by the red
and orange curves, in which both states are local minima of the free energy.
On the left branch of the red curve, the insulating state becomes unstable,
and on the right, the nFL state becomes unstable.

The boundary of stability of the insulating state determined by the Lyapunov
exponent appears to agree with the numerical results where expected. Because
it is difficult to compute the boundary of the hysteresis region numerically
by the iterative method described in [25], especially at very low
temperatures, and because our equation (75) is a low-temperature
approximation, it is difficult to meaningfully quantify this agreement, but
nonetheless the connection between the chaos bound and the stability of the
phase is plausible. The rightmost boundary of the hysteresis region, beyond
which the nFL state is unstable, will require more work to calculate because
our analysis assumes the conformal limit which yields an exponent of $2\pi T$
regardless of $\mu$ and $T$. In reality there should be small corrections to
this result, but calculating these will require going beyond the conformal
limit. We surmise that the rightmost boundary of the hysteresis region
corresponds to the curve where these small corrections switch from negative to
positive and violate the chaos bound. In this sense the chaos bound may also
be interpreted as a stability bound [39]. We leave a detailed analysis of this
to a future study.

###### Acknowledgements.

We thank A.V. Chubukov and J. Schmalian for useful discussions. This work is
supported by startup funds at the University of Florida and by NSF under award
number DMR-2045781.

##  Appendix A Large-$M$,$N$ counting

In this appendix we justify the claim that the leading order diagrams in the
large-$M$,$N$ expansion of the OTOC (37) are those depicted in Fig. 2 by
explicitly computing the factors of $M$ and $N$ up to second nonvanishing
order in the interaction. At zeroth order, the fermion operators are free, so
the first term is

$\displaystyle\frac{1}{M^{2}}$
$\displaystyle\sum_{ij}\Tr[\rho^{\frac{1}{2}}\\{c_{i}^{0}(t_{1}),c_{j}^{0\dagger}(0)\\}\rho^{\frac{1}{2}}\\{c_{i}^{0}(t_{2}),c_{j}^{0\dagger}(0)\\}^{\dagger}]$
(76)
$\displaystyle=\frac{1}{M^{2}}\Tr[\rho^{\frac{1}{2}}G^{R,0}(t_{1})\delta_{ij}\rho^{\frac{1}{2}}G^{A,0}(t_{2})\delta_{ij}]$
(77) $\displaystyle=\frac{1}{M}G^{R,0}(t_{1})G^{A,0}(t_{2}).$ (78)

The superscript zeroes indicate free fermions and their free Green’s
functions. The summation over fermion flavors has produced one factor of $M$
for an overall factor of $1/M$ at this order. This term is represented by
diagram (a) in Fig. 5. The fermion flavors on each rail are constrained to be
the same, so the single factor of $M$ corresponds to the sum over the single
fermion index $i$ in the diagram.

Figure 5: Some diagrams at leading order in the large-$M$,$N$ expansion of the
fermion OTOC. All propagators in these diagrams are bare, unlike in the
figures in the main text.

At first order in the interaction, all diagrams vanish because every term is
proportional to the random coupling $t^{\alpha}_{ij}$ which has zero mean. All
diagrams with an odd number of interaction vertices vanish similarly. This
disorder averaging can be viewed as an additional field with a constant
propagator which forces the two fermion flavors and the boson flavor at one
end to match those at the other.

Figure 6: A diagram with two Type-I rungs, included in the ladder series.

At first nonvanishing order, the two vertices can either be placed on the same
rail or with one on each rail. The former option, however, merely forms part
of the self-energy for the fermion on that rail and is already included in
Fig. 2 – since those propagators are all dressed, the very first term in the
sum contains this contribution. All the self-energy diagrams, by construction,
are of leading order in $M$ and $N$, so these contributions to the OTOC are as
well. The latter option is depicted in (b) of Fig. 5 (the “Type-I rung”).
Compared to the zeroth order term, it contains one additional factor of $M$
from the new fermion sum, one additional factor of $N$ from the new boson sum,
and a factor of $1/MN$ from the interaction. Hence this diagram is also
leading in $M$ and $N$. In fact, the addition of each Type-I rung at any order
has the same effect and does not change the $M$ and $N$ factors. This is true
at all orders, so the diagram of Fig. 6 for example is also dominant, and it
is included in our ladder equation.

At second nonvanishing order, there are several ways of contracting the
vertices, and the disorder average is over a product of four $t$’s, which can
be decomposed by Isserlis’s theorem into a sum of averages of all pairings of
the $t$’s (diagrammatically, all the ways of pairing the vertices by dotted
lines). Again many of these diagrams contain pieces which form part of the
self-energy, such as (d) in Fig. 5, which is accounted for in Fig. 2 by the
second term in the sum. Another possibility is two Type-I rungs in sequence.
The important new diagram at this order is the “Type-II” rung (or “box”) of
(c) in Fig. 5. Compared to the zeroth order diagram, it has three more fermion
sums (a factor of $M^{3}$), one more boson sum (a factor of $N$), and four
more vertices (a factor of $1/M^{2}N^{2}$) for a combined factor of $M/N$.
Hence this diagram is also of leading order, as are any diagrams obtained by
appending more of these Type-II rungs.

Notice that the pairing of vertices in the disorder average is important. For
example, the Type-II rung must have the disorder average done as in (c). The
other two possible pairings of the vertices force too many constraints on the
flavors, and it can be easily verified that such terms are subleading in $M$
and/or $N$. Another example is shown in (a) of Fig. 7: if the disorder average
connected vertices on the same rung, this would be a relevant self-energy
contribution, but with the inter-rung pairing as shown, it is suppressed by
$1/N$.

Figure 7: Examples of subleading diagrams in the large-$M$,$N$ expansion of
the fermion OTOC. All propagators are bare, unlike in the figures in the main
text.

There are a few other diagrams to check, including crossed rungs, other
pairings in the disorder average, and diagram (b) of Fig. 7, but it is not
difficult to verify that they are all subleading. The dominant diagrams, at
all orders, are precisely those achieved by iteratively tacking on Type-I and
Type-II rungs to the zeroth order diagram as our ladder equation depicts.

##  Appendix B Feynman rule for the ladder diagrams in the OTOC

In the interaction picture, the leading correction to the fermion OTOC
$F_{c}(t_{1},t_{2})$ beyond self-energy effects is given by, using Wick’s
theorem on the path integral along the complex time contour,

$\displaystyle F_{1}=$
$\displaystyle\frac{1}{M}\int_{0}^{t_{1}}\differential{t}^{\prime}\int_{0}^{t_{2}}\differential{t}^{\prime\prime}\wick$
$\displaystyle=$
$\displaystyle\frac{\omega_{0}^{3}}{M}\int_{0}^{t_{1}}\differential{t}^{\prime}\int_{0}^{t_{2}}\differential{t}^{\prime\prime}\wick{\\{[\c{1}c^{\dagger}(it^{\prime})\c{2}\phi(it^{\prime})\c{3}c(it^{\prime}),\c{3}c^{\dagger}(it_{1})],\c{1}c(0)\\}\\{[\c{4}c^{\dagger}(\beta/2+it^{\prime\prime})\c{2}\phi(\beta/2+it^{\prime\prime})\c{5}c(\beta/2+it^{\prime\prime}),\c{5}c^{\dagger}(\beta/2+it_{2})],\c{4}c(\beta/2)\\}^{\dagger},}$
(79)

which corresponds to a ladder diagram with one rung, which we show in Fig. 8.
Here $I(t)\propto tc^{\dagger}(t)\phi(t)c(t)$ is the interaction vertex (with
the flavor summations suppressed for brevity), which we used in the second
step. The contraction is defined as the complex-time-contour-ordered
correlator

$\wick{\c{\mathcal{O}_1}(t_{1})\c{\mathcal{O}_2}(t_{2})}=\langle T_{\rm
cont}\mathcal{O}_{1}(t_{1})\mathcal{O}_{2}(t_{2})\rangle,$ (80)

where $T_{\rm cont}$ is the ordering operator on the complex time contour.
Therefore for contractions across two real-time folds, we get Wightman
propagators, e.g.,

$\wick{\c{2}\phi(0)\c{2}\phi(\beta/2+it)=D^{lr}(t),}$ (81)

Figure 8: The ladder diagram corresponding to (79).

Importantly, as we mentioned in the main text, the (anti)-commutators in Eq.
(79) are not defined in the usual way, but via altering the order the
operators on the complex time contour, i.e.,

$\displaystyle\wick{\left[{\mathcal{O}_{1}}(it_{1}),{\mathcal{O}_{2}}(it_{2})\right]}_{\pm}$
(82) $\displaystyle\equiv$ $\displaystyle\lim_{\epsilon\to
0}\mathcal{O}_{1}(\epsilon+it_{1})\mathcal{O}_{2}(it_{2})-\mathcal{O}_{1}(it_{1})\mathcal{O}_{2}(\epsilon+it_{2}),$

where $[\cdots]_{+}\equiv\\{\cdots\\}$. Therefore in this notation

$\wick{\\{\c{1}c(0),\c{1}c^{\dagger}(it)\\}\equiv
G^{R}(t),~{}~{}~{}\\{\c{1}c^{\dagger}(0),\c{1}c(it)\\}\equiv G^{A}(t).}$ (83)

The nested (anti)-commutators in Eq. (79) are defined using multiple
$\epsilon_{i}\to 0$, making sure the limits are taken such that $\epsilon_{i}$
corresponding to inner (anti)-commutators tend to zero first. It is easy to
verify that the (anti)-commutator defined this way satisfies the same
algebraic properties as the original version.

The expression in Eq. (79) can be simplified by successively extracting
contractions (c-numbers) from the nested (anti)-commutator. Importantly, when
contracting (anti)-commutators, the only nonvanishing contribution comes from
contractions that match the (anti)-commutation, i.e., between operators
separated by a comma. In particular,

$\wick{[{\bar{a}},\c{3}b]_{\pm}\c{3}{\bar{b}}=0,}$ (84)

because the order between $b$ and $\bar{b}$ is not altered by the
$\epsilon$-prescription above, and the two terms forming the commutator
defined in Eq. (82) cancel. As an example for nested commutators,

$\wick{\\{\c{1}a,[\c{1}{\bar{a}}\c{2}{\bar{b}},\c{2}b]\\}=-\\{\c{4}a,\\{\c{4}{\bar{a}},\c{3}b\\}\c{3}{\bar{b}}\\}+\\{\c{4}a,\c{4}{\bar{a}}\\{\c{3}b,\c{3}{\bar{b}}\\}\\}=\\{\c{4}a,\c{4}{\bar{a}}\\}\\{\c{3}b,\c{3}{\bar{b}}\\},}$
(85)

where in the first step we used the mathematical identity
$[\bar{a}\bar{b},b]\equiv-\\{\bar{a},b\\}\bar{b}+\bar{a}\\{\bar{b},\bar{b}\\}$,
and the first term in the middle expression vanishes because of (84).
Combining this recipe with Eq. (81), it can be shown in general that for an
OTOC with nested commutators, one can successively replace contractions
separated by a comma with the retarded or advanced Green’s functions, and
replace contractions across real-time folds via Wightman propagators. This is
the analog of Feynman rules for OTOCs.

Applying the Feynman rules, we get

$\displaystyle F_{1}=$
$\displaystyle\frac{\omega_{0}^{3}}{M}\int_{-\infty}^{\infty}\differential{t}^{\prime}\int_{-\infty}^{\infty}\differential{t}^{\prime\prime}G^{R}(t^{\prime})G^{R}(t_{1}-t^{\prime})D^{lr}(t^{\prime\prime}-t^{\prime})$
$\displaystyle\times G^{A}(t^{\prime\prime})G^{A}(t_{2}-t^{\prime\prime}).$
(86)

Higher-order ladder diagrams are obtained by replacing
$\int_{0}^{t_{1}}\differential{t}^{\prime}[I(it^{\prime}),c^{\dagger}(it_{1})]$
with

$\int_{0}^{t_{1}}\differential{t}_{n}^{\prime}\cdots\int_{0}^{\differential{t}_{2}^{\prime}}\differential{t}_{1}^{\prime}[I(it_{n}^{\prime}),[\cdots,[I(it_{1}^{\prime}),c^{\dagger}(it_{1})],$
(87)

which can then be contracted using the same procedure. It is straightforward
to see that the Feynman rule for the ladder diagrams in the OTOC is to assign
a retarded (advanced) Green’s function to all vertical lines, and a Wightman
correlator to the horizontal lines.

##  Appendix C Leading order diagrams involving $F_{\bar{c}}$ from Feynman
rules

At low orders in the expansion of the OTOC, one finds retarded Green’s
functions on one rail of the ladders and advanced ones on the other. However,
it is possible for both retarded and advanced Green’s functions to appear on
either rail. This first becomes relevant to the “anomalous” OTOC $F_{\bar{c}}$
at second order in perturbation theory, and relevant to the OTOC $F_{c}$ at
fourth order in perturbation theory.

Figure 9: (a) A contribution to the “anomalous” OTOC at second order. (b) A
contribution to the original OTOC at fourth order, which contains (a) as a
subdiagram. All propagators are bare.

We begin with

$F_{\bar{c}}\sim\Tr[\rho^{1/2}\\{\tilde{c}_{i}(it_{1}),\tilde{c}_{j}(0)\\}\rho^{1/2}\\{\tilde{c}_{i}(it_{2}),\tilde{c}_{j}(0)\\}^{\dagger}]$
(88)

which is depicted on the left side of the second equation in Fig. 2. The
tildes indicate that these are Heisenberg picture operators, which are related
to the interaction picture operators $c(\tau)$ by
$\tilde{c}(\tau)=U^{\dagger}c(\tau)U$, where

$U=\mathcal{T}\exp\quantity(-\frac{i}{\sqrt{MN}}t_{kk^{\prime}}^{\alpha}\int_{0}^{\tau}\differential{\tau^{\prime}}c_{k}^{\dagger}(\tau^{\prime})c_{k^{\prime}}(\tau^{\prime})\phi_{\alpha}(\tau^{\prime}))$
(89)

with flavor summations implied. If we expand $U^{\dagger}c(\tau)U$ to second
order in the random coupling, we find

$[U^{\dagger}c_{i}(t_{1})U]^{(2)}=\frac{-1}{MN}t_{kk^{\prime}}^{\alpha}t_{ll^{\prime}}^{\beta}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}[c^{\dagger}_{l}(it^{\prime\prime})c_{l^{\prime}}(it^{\prime\prime})\phi_{\beta}(it^{\prime\prime}),[c^{\dagger}_{k}(it^{\prime})c_{k^{\prime}}(it^{\prime})\phi_{\alpha}(it^{\prime}),c_{i}(it_{1})]]$
(90)

as detailed in [40] and [41], and as claimed in Eq. (87). Suppressing
prefactors, integrations, and flavors indices for brevity, there is a diagram
(Fig. 9(a)) corresponding to the contraction

$\displaystyle F_{\bar{c}}=\frac{\omega_{0}^{6}}{M}\times$
$\displaystyle\wick{\int\\!\\!\differential{t^{\prime}}\\!\differential{t^{\prime\prime}}\\!\differential{t^{\prime\prime\prime}}\\!\differential{t^{\prime\prime\prime\prime}}\\{[\c{2}c^{\dagger}\c{8}c\c{5}\phi(it^{\prime\prime}),[\c{1}c^{\dagger}\c{7}c\c{5}\phi(t^{\prime}),\c{1}c_{i}(it_{1})]],\c{2}c_{j}(0)\\}\\{[\c{8}c^{\dagger}\c{4}c\c{6}\phi({\beta}/{2}+it^{\prime\prime\prime\prime}),[\c{7}c^{\dagger}\c{3}c\c{6}\phi({\beta}/{2}+it^{\prime\prime\prime}),\c{3}c_{i}^{\dagger}({\beta}/{2}+it_{2})]],\c{4}c_{j}^{\dagger}({\beta}/{2})\\}]}$
$\displaystyle+\cdots,$ (91)

where we reiterate that the notation with both commutators and contractions is
explained in Appendix B. Applying the Feynman rules developed in Appendix B,
we find that (91) produces an integral with the structure

$\displaystyle F_{\bar{c}}=\frac{\omega_{0}^{6}}{M}$
$\displaystyle\int\\!\\!\differential{t^{\prime}}\\!\differential{t^{\prime\prime}}\\!\differential{t^{\prime\prime\prime}}\\!\differential{t^{\prime\prime\prime\prime}}G^{R}(t_{1}-t^{\prime})D^{R}(t^{\prime}-t^{\prime\prime})G^{A}(t^{\prime\prime})G^{r}(t^{\prime}-t^{\prime\prime\prime})G^{r}(t^{\prime\prime}-t^{\prime\prime\prime\prime})G^{A}(t_{2}-t^{\prime\prime\prime})D^{R}(t^{\prime\prime\prime}-t^{\prime\prime\prime\prime})G^{R}(t^{\prime\prime\prime\prime}).$
$\displaystyle+\cdots$ (92)

This diagram is the last term of the second equation in Fig. 2 with the shaded
rectangle removed. As a ladder diagram similar to that in Fig. 5(c), it is
straightforward to show that this diagram is of $O(1/M)$, i.e., at the same
order as other leading diagrams in the series. This contribution is included
in the integral equation (III), as can be seen via an iterative expansion.

As is indicated in the top equation of Fig. 2, the existence of $F_{\bar{c}}$
leads to an additional contribution to $F_{c}$, which at leading order is
diagrammatically shown in Fig. 9(b). Again this diagram is at $O(1/M)$ and can
be explicitly written in contractions just like Eq. (91). Writing it
explicitly in terms of contracted nested commutators will unfortunately take
too much space, but the corresponding Feynman diagram is obtained by combining
two leading diagrams for $F_{\bar{c}}$. Applying the Feynman rule developed in
Appendix B, we get

$\displaystyle F_{c}=$
$\displaystyle\cdots+\frac{\omega_{0}^{12}}{M}\int\\!\\!\differential{\bar{t}^{\prime}}\\!\differential{\bar{t}^{\prime\prime}}\\!\differential{\bar{t}^{\prime\prime\prime}}\\!\differential{\bar{t}^{\prime\prime\prime\prime}}G^{R}(t_{1}-\bar{t}^{\prime})D^{R}(\bar{t}^{\prime}-\bar{t}^{\prime\prime})G^{r}(\bar{t}^{\prime}-\bar{t}^{\prime\prime\prime})G^{r}(\bar{t}^{\prime\prime}-\bar{t}^{\prime\prime\prime\prime})G^{A}(t_{2}-\bar{t}^{\prime\prime\prime})D^{R}(\bar{t}^{\prime\prime\prime}-\bar{t}^{\prime\prime\prime\prime})$
$\displaystyle\int\\!\\!\differential{t^{\prime}}\\!\differential{t^{\prime\prime}}\\!\differential{t^{\prime\prime\prime}}\\!\differential{t^{\prime\prime\prime\prime}}G^{R}(\bar{t}^{\prime\prime\prime\prime}-t^{\prime})D^{R}(t^{\prime}-t^{\prime\prime})G^{A}(t^{\prime\prime})G^{l}(t^{\prime}-t^{\prime\prime\prime})G^{l}(t^{\prime\prime}-t^{\prime\prime\prime\prime})G^{A}(\bar{t}^{\prime\prime}-t^{\prime\prime\prime})D^{R}(t^{\prime\prime\prime}-t^{\prime\prime\prime\prime})G^{R}(t^{\prime\prime\prime\prime})$
$\displaystyle+\cdots,$ (93)

which is also obtained by iteratively expanding the integral equation (III).

## References

  * [1] A. Kitaev, A simple model of quantum holography, KITP strings seminar and Entanglement 2015 program (Feb 12, Apr 7, and May 27, 2015), https://online.kitp.ucsb.edu/online/entangled15/.
  * Chowdhury  _et al._ [2022] D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, Sachdev-Ye-Kitaev models and beyond: Window into non-Fermi liquids, Rev. Mod. Phys. 94, 035004 (2022), arXiv:2109.05037 [cond-mat.str-el] .
  * Fu  _et al._ [2017] W. Fu, D. Gaiotto, J. Maldacena, and S. Sachdev, Supersymmetric sachdev-ye-kitaev models, Phys. Rev. D 95, 026009 (2017).
  * Rosenhaus [2019] V. Rosenhaus, An introduction to the SYK model, J. Phys. A 52, 323001 (2019), arXiv:1807.03334 [hep-th] .
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the sachdev-ye-kitaev model, Phys. Rev. D 94, 106002 (2016).
  * Kitaev and Suh [2018] A. Kitaev and S. J. Suh, The soft mode in the Sachdev-Ye-Kitaev model and its gravity dual, JHEP 05, 183, arXiv:1711.08467 [hep-th] .
  * Patel  _et al._ [2022] A. A. Patel, H. Guo, I. Esterlis, and S. Sachdev, Universal theory of strange metals from spatially random interactions, arXiv e-prints , arXiv:2203.04990 (2022), arXiv:2203.04990 [cond-mat.str-el] .
  * Esterlis  _et al._ [2021] I. Esterlis, H. Guo, A. A. Patel, and S. Sachdev, Large-$n$ theory of critical fermi surfaces, Phys. Rev. B 103, 235129 (2021).
  * Aldape  _et al._ [2022] E. E. Aldape, T. Cookmeyer, A. A. Patel, and E. Altman, Solvable theory of a strange metal at the breakdown of a heavy fermi liquid, Phys. Rev. B 105, 235111 (2022).
  * Hussey  _et al._ [1998] N. E. Hussey, A. P. Mackenzie, J. R. Cooper, Y. Maeno, S. Nishizaki, and T. Fujita, Normal-state magnetoresistance of ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$, Phys. Rev. B 57, 5505 (1998).
  * Stewart [2001] G. R. Stewart, Non-fermi-liquid behavior in $d$\- and $f$-electron metals, Rev. Mod. Phys. 73, 797 (2001).
  * Keimer  _et al._ [2014] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, High Temperature Superconductivity in the Cuprates, arXiv e-prints , arXiv:1409.4673 (2014), arXiv:1409.4673 [cond-mat.supr-con] .
  * Varma  _et al._ [1989] C. M. Varma, P. B. Littlewood, S. Schmitt-Rink, E. Abrahams, and A. E. Ruckenstein, Phenomenology of the normal state of cu-o high-temperature superconductors, Phys. Rev. Lett. 63, 1996 (1989).
  * Cha  _et al._ [2020] P. Cha, N. Wentzell, O. Parcollet, A. Georges, and E.-A. Kim, Linear resistivity and sachdev-ye-kitaev (syk) spin liquid behavior in a quantum critical metal with spin-1/2 fermions, Proceedings of the National Academy of Sciences 117, 18341 (2020), https://www.pnas.org/content/117/31/18341.full.pdf .
  * Guo  _et al._ [2020] H. Guo, Y. Gu, and S. Sachdev, Linear in temperature resistivity in the limit of zero temperature from the time reparameterization soft mode, Annals Phys. 418, 168202 (2020), arXiv:2004.05182 [cond-mat.str-el] .
  * Lee [2018] S.-S. Lee, Recent developments in non-fermi liquid theory, Annual Review of Condensed Matter Physics 9, 227 (2018), https://doi.org/10.1146/annurev-conmatphys-031016-025531 .
  * Wang [2020] Y. Wang, Solvable strong-coupling quantum-dot model with a non-fermi-liquid pairing transition., Physical review letters 124 1, 017002 (2020).
  * Esterlis and Schmalian [2019] I. Esterlis and J. Schmalian, Cooper pairing of incoherent electrons: An electron-phonon version of the sachdev-ye-kitaev model, Phys. Rev. B 100, 115132 (2019).
  * Choi  _et al._ [2022] W. Choi, O. Tavakol, and Y. B. Kim, Pairing instabilities of the Yukawa-SYK models with controlled fermion incoherence, SciPost Physics 12, 151 (2022), arXiv:2110.02968 [cond-mat.str-el] .
  * Hauck  _et al._ [2020] D. Hauck, M. J. Klug, I. Esterlis, and J. Schmalian, Eliashberg equations for an electron–phonon version of the sachdev–ye–kitaev model: Pair breaking in non-fermi liquid superconductors, Annals of Physics 417, 168120 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Kim  _et al._ [2021] J. Kim, E. Altman, and X. Cao, Dirac fast scramblers, Phys. Rev. B 103, L081113 (2021).
  * Schmalian [2022] J. Schmalian, Holographic superconductivity of a critical fermi surface (2022).
  * Wang and Chubukov [2020] Y. Wang and A. V. Chubukov, Quantum Phase Transition in the Yukawa-SYK Model, Phys. Rev. Res. 2, 033084 (2020), arXiv:2005.07205 [cond-mat.str-el] .
  * Pan  _et al._ [2021] G. Pan, W. Wang, A. Davis, Y. Wang, and Z. Y. Meng, Yukawa-syk model and self-tuned quantum criticality, Phys. Rev. Research 3, 013250 (2021).
  * Wang  _et al._ [2021] W. Wang, A. Davis, G. Pan, Y. Wang, and Z. Y. Meng, Phase diagram of the spin-$\frac{1}{2}$ yukawa–sachdev-ye-kitaev model: Non-fermi liquid, insulator, and superconductor, Phys. Rev. B 103, 195108 (2021).
  * Azeyanagi  _et al._ [2018] T. Azeyanagi, F. Ferrari, and F. I. S. Massolo, Phase diagram of planar matrix quantum mechanics, tensor, and sachdev-ye-kitaev models, Phys. Rev. Lett. 120, 061602 (2018).
  * Smit  _et al._ [2021] R. L. Smit, D. Valentinis, J. Schmalian, and P. Kopietz, Quantum discontinuity fixed point and renormalization group flow of the sachdev-ye-kitaev model, Phys. Rev. Research 3, 033089 (2021).
  * Kim  _et al._ [2020] J. Kim, X. Cao, and E. Altman, Low-rank sachdev-ye-kitaev models, Phys. Rev. B 101, 125112 (2020).
  * Stanford [2016] D. Stanford, Many-body chaos at weak coupling, JHEP 10, 009, arXiv:1512.07687 [hep-th] .
  * Gu  _et al._ [2017] Y. Gu, X.-L. Qi, and D. Stanford, Local criticality, diffusion and chaos in generalized Sachdev-Ye-Kitaev models, JHEP 05, 125, arXiv:1609.07832 [hep-th] .
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Quasiclassical Method in the Theory of Superconductivity, Soviet Journal of Experimental and Theoretical Physics 28, 1200 (1969).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, JHEP 08, 106, arXiv:1503.01409 [hep-th] .
  * Anninos  _et al._ [2016] D. Anninos, T. Anous, and F. Denef, Disordered quivers and cold horizons, Journal of High Energy Physics 2016, 71 (2016).
  * Anninos  _et al._ [2015] D. Anninos, T. Anous, P. de Lange, and G. Konstantinidis, Conformal quivers and melting molecules, Journal of High Energy Physics 2015, 66 (2015).
  * Sachdev [2015] S. Sachdev, Bekenstein-Hawking Entropy and Strange Metals, Phys. Rev. X 5, 041025 (2015), arXiv:1506.05111 [hep-th] .
  * Murugan  _et al._ [2017] J. Murugan, D. Stanford, and E. Witten, More on supersymmetric and 2d analogs of the syk model, Journal of High Energy Physics 2017, 146 (2017).
  * Marcus and Vandoren [2018] E. Marcus and S. Vandoren, A new class of syk-like models with maximal chaos, Journal of High Energy Physics 2019, 1 (2018).
  * Sorokhaibam [2020] N. Sorokhaibam, Phase transition and chaos in charged SYK model, Journal of High Energy Physics 2020, 55 (2020), arXiv:1912.04326 [hep-th] .
  * Murthy  _et al._ [2021] C. Murthy, A. Pandey, I. Esterlis, and S. A. Kivelson, A stability bound on the $t$-linear resistivity of conventional metals (2021).
  * Chowdhury and Swingle [2017] D. Chowdhury and B. Swingle, Onset of many-body chaos in the $o(n)$ model, Phys. Rev. D 96, 065005 (2017).
  * Tikhanovskaya  _et al._ [2022] M. Tikhanovskaya, S. Sachdev, and A. A. Patel, Maximal quantum chaos of the critical fermi surface, Phys. Rev. Lett. 129, 060601 (2022).