# Thermodynamics and dynamics of coupled complex SYK models

Jan C. Louw<EMAIL_ADDRESS>Institute for
Theoretical Physics, Georg-August-Universität Göttingen, Friedrich-Hund-Platz
1, 37077 Göttingen, Germany  Linda M. van Manen<EMAIL_ADDRESS>Friedrich-Schiller-Universität, Institute for Theoretical Physics,
Max‐Wien‐Platz 1, 07743 Jena, Germany  Rishabh Jha  rishabh.jha@uni-
goettingen.de Institute for Theoretical Physics, Georg-August-Universität
Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany

###### Abstract

It has been known that the large-$q$ complex SYK model falls under the same
universality class as that of van der Waals (mean-field) which is also shared
by a variety of black holes. At the same time, it also saturates the
Maldacena-Shenker-Stanford (MSS) bound and is thus maximally chaotic. This
work establishes the robustness of shared universality class and quantum chaos
for SYK-like models by extending to a system of coupled large-$q$ complex SYK
models of different orders. We provide a detailed derivation of thermodynamic
(critical exponents) properties observing a phase transition and dynamic
(Lyapunov exponent) properties via the out-of-time correlator (OTOC)
calculations. Our analysis reveals that, despite the introduction of an
additional scaling parameter through interaction strength ratios, the system
undergoes a continuous phase transition at low temperatures, similar to that
of a single SYK model. The critical exponents align with the Landau-Ginzburg
(mean-field) universality class, shared with van der Waals gases and various
AdS black holes. Furthermore, we demonstrate that the coupled SYK system
remains maximally chaotic in the large-$q$ limit at low temperatures, adhering
to the Maldacena-Shenker-Stanford (MSS) bound, a feature consistent with
single large-$q$ complex SYK model. These findings open avenues for broader
inquiries into the universality and chaos in complex quantum systems by
showing that our coupled SYK system belong to the same universality class as
that of van der Waals and various AdS black holes while saturating the MSS
bound of quantum chaos.

††preprint: APS/123-QED

##  I Introduction

The Sachdev-Ye-Kitaev (SYK) model, initially conceptualized in the context of
condensed matter physics, has rapidly emerged as a pivotal tool in exploring
quantum gravity and chaos. Characterized by its non-Fermi liquid behavior and
lack of quasiparticle excitations, the SYK model exemplifies a quantum many-
body system that is exactly solvable in certain limits. The large-$q$ complex
SYK model, in particular, has been instrumental in revealing connections
between quantum chaos, black hole thermodynamics, and the holographic
principle. We define quantum chaos in terms of exponential growth in the out-
of-time-correlators (OTOCs). Studies have shown that the SYK model exhibits
maximal chaos, saturating the Maldacena-Shenker-Stanford (MSS) bound, and
parallels the thermodynamic properties of black holes, particularly in the
context of the holographic duality [1, 2, 3].

It was recently shown that the large-$q$ SYK model has a second order phase
transition, which matches the universality class of a wide range of black
holes [4]. Similar phase transitions are also found the finite-$q$ case [5,
6]. Given the analytic solvability of the large-$q$ model, one was able to
find the exact thermodynamics and find its exact holographic dual gravitation
model [3]. This duality is in terms of the thermodynamics as well as the
dynamics, namely the Lyapunov exponents, of the respective models. Here the
thermodynamics is in the grand canonical ensemble. On the SYK side, this means
that we have a chemical potential coupled to the charge density. The phase
transition results from a varying conjugate field (chemical potential)
inducing a jump in charge density. The analogy is most easily seen in a
standard $3+1$-dimensional asymptotically Anti-de-Sitter (AdS3+1) black hole
system, where the black hole charge is the conjugate field and the surface
charge density is the order parameter (see analogy 2 in [7]). Indeed, these
two parameters are mapped onto one another in the thermodynamic dictionary
between the $0+1$-dimensional SYK dot and the AdS1+1 charged black hole system
[3]. The jump in charge densities of the SYK dot is then directly related to a
jump in black hole density. In particular, since the SYK dot has no size, the
charge density (charge per flavors of fermions) jump actually corresponds to a
change in charge. On the black hole side, the charge is merely a conjugate
field, i.e., it is a set parameter and the jump happens in the surface charge
density due to the black hole transition between small and large horizon
radii.

For all cases, the coexistence line between the two phases terminates at the
critical point, e.g., $(\mu_{c},T_{c})$ for the SYK model [4] or
$(Q_{c},1/T_{c})$ for AdS3+1 [7]. These critical points correspond to the van
der Waals second order phase transition (Landau-Ginzburg mean field critical
exponents). Note the difference that the SYK model has no transition at
temperatures above $T_{c}$, while the AdS3+1 case has no transition for
temperatures below $T_{c}$. Both transitions also have a non-trivial end-point
corresponding to a first order chaotic-to-regular phase transition.

There is a caveat though. Only for the uncharged case does the smaller black
hole evaporate such that we are left with regular (non-chaotic) non-
interacting thermal radiation for a particular value of temperature [8]. On
the SYK side, however, there exists a whole range of temperatures where we
have a chaotic to regular (non-chaotic phase transition) [4] (see section V.4
below). Although the thermodynamics between the SYK model and
thermodynamically holographic dual does match, there is a short-coming of this
duality that it does not reproduce the standard Hawking-Page transition [8].
But at the same time, the SYK model has proved instrumental in exploring
connections between quantum chaos and black holes [9, 10]. The natural
question to ask is how robust are these quantum chaotic and thermodynamic
properties of SYK-like models.

In this paper, we study the robustness of the chaotic and thermodynamic
properties of a system of coupled large-$q$ complex SYK models and the
associated dynamical inequivalence to black holes in the context of Hawking-
Page transition. We do this by considering the sum of two non-commuting SYK
models of different orders. This has the effect of introducing an additional
length scale. As expected from relevance, in a renormalization context, the
lower order term dominates at low temperatures. Note that typical physical
systems will always have at least two competing terms, one associated with its
kinetic energy, the other with particle interactions. A strongly dominating
kinetic term results in a free, hence non-chaotic, system. As such, it seems
reasonable, that such a term would have a large effect on the thermodynamics.
Here, interactions are of order $q$, while we add a pseudo-kinetic term of
order $q/2$. The term “pseudo-kinetic” is motivated by various findings. For
one, such a term coupling a lattice of dots at equilibrium yields the standard
intermediate temperature resistivity related to the $q=4$, hence quadratic
hopping $q/2=2$, case [11, 12]. This is in agreement with the observation that
the large-$q$ expansions converge rather fast, leading to a large quantitative
and qualitative overlap even with the $q=4$ case as discussed at length in
[13]. This then gives us the advantage of an analytically solvable system via
a $1/q$ expansion, as opposed to the $q=4$ case which is only numerically
solvable. This is our motivation in calling the $q/2$ term in our model below
as “pseudo-kinetic” term because of this convergence and qualitative
similarity with the finite $q$ case such as $q=4$ while retaining the
analytical tractability of the system.

Since the two Hamiltonians do not commute, one might expect rather different
thermodynamics. The addition of a new scaling parameter, the ratio of
interaction strengths, adds a layer of complexity to the model. Despite this,
we find that the system retains the characteristic continuous phase transition
of single large-$q$ SYK model for all ratios in coupling strength and the
phase transition persists in the low-temperature regime. The critical point
shifts for each value of the ratio of the respective coupling strengths.
However, we find that this variant of coupled large-$q$ SYK models in fact has
the same partition function around the critical point. In other words, the
coupling of two SYK models still flows to the same model under
renormalization. The critical exponents observed in this coupled system
intriguingly align with the Landau-Ginzburg (mean-field) universality class,
which is a shared feature with van der Waals gases and a spectrum of AdS black
holes [7]. This finding underscores a potentially universal behavior
underlying these disparate physical systems.

Moving to characterize the chaos in our system of coupled SYK models, it seems
reasonable that the introduction of a new length scale in the system, namely
the ratio of coupling constants of the two SYK terms in the Hamiltonian, might
change the dynamics at lower temperatures. This would be reflected in some
novel coupling found between the $q$ and $q/2$ related kernels defining the
time evolutions of the OTOCs. Surprisingly, we find that, like the self
energies [1], the kernels are also additive for our coupled SYK system. The
chaotic nature does not change much, since the dominating SYK term just leads
to the same maximal chaos as a single SYK term in the low-temperature limit.
As such, we find that the dynamics of our coupled SYK system has a rather
similar phase diagram to that of the single SYK model. This retention of
maximal chaos, as denoted by adherence to the MSS bound, aligns with the
behavior observed in single large-$q$ complex SYK model [1]. The persistence
of maximal chaos in the coupled system is not just a continuation of the SYK
model’s established traits but also highlights the robust nature of chaos in
these quantum systems. This work establishes this robustness for a coupled SYK
system not just for the case of chaos but also for the shared universality
with the Landau-Ginzburg mean field critical exponents.

##  II Model

We consider the following coupled complex large-$q$ SYK system Hamiltonian in
equilibrium:

${\cal H}=J_{q}{\cal H}_{q}+J_{q/2}{\cal H}_{q/2}$ (1)

where ${\cal H}_{\kappa q}$ ($\kappa\in\\{\frac{1}{2},1\\}$) is large-$q$
complex SYK model Hamiltonian given by

${\cal H}_{\kappa
q}=\sum\limits_{\begin{subarray}{c}\\{\bm{\mu}\\}_{1}^{\kappa q/2}\\\
\\{\bm{\nu}\\}_{1}^{\kappa
q/2}\end{subarray}}^{N}X^{\bm{\mu}}_{\bm{\nu}}c^{{\dagger}}_{\mu_{1}}\cdots
c^{{\dagger}}_{\mu_{\kappa q/2}}c_{\nu_{\kappa
q/2}}^{\vphantom{{\dagger}}}\cdots c_{\nu_{1}}^{\vphantom{{\dagger}}}+\text{
h.c. }$ (2)

summing over $\\{\bm{\mu}\\}_{1}^{q/2}\equiv 1\leq\mu_{1}<\cdots<\mu_{\kappa
q/2}\leq{N}$. Here $X^{\bm{\mu}}_{\bm{\nu}}$ is a random matrix whose
components are derived from a Gaussian ensemble with zero mean and variance
given by

$\overline{|X|^{2}}=\frac{(\kappa q)^{-2}((\kappa q/2)!)^{2}}{(N/2)^{\kappa
q-1}}$ (3)

Our main focus will be the thermodynamics of this model, defined by the
partition function ${\cal Z}=e^{-\beta N\Omega}$. The thermodynamics of this
model is, in fact, equivalent 111The equivalence can be seen by considering
the Lagrangian in Section II.B in [15], for a uniformly coupled lattice, i.e.,
site independent coupling $D_{i},J_{i}\to D,J$. The remaining terms in the
Lagrangian are the Green’s functions associated with the $i^{\text{th}}$ dot
${\cal G}_{i}(t_{1},t_{2})$. Since the lattice is at equilibrium, there is no
heat or charge flow, meaning that each dot has equal temperature and charge.
This leads to the Green’s functions also being independent of $i$. Although
the results in [15] are valid only when the $r$-body hopping on nearest
lattice sites is of the order $r={\cal O}(q^{0})$, the formalism and
methodology developed there is valid for any value of $r$, including $r=q/2$.
to a uniformly coupled lattice of SYK dots with $q/2$-body hopping with
$q$-body inter-dot interactions as considered in [15]. Here the grand
potential (per lattice site) $\Omega$ may, in fact, be found by studying the
Green’s functions which are defined as follows:

$\mathcal{G}\left(t_{1},t_{2}\right)\equiv\frac{-1}{N}\sum_{\mu=1}^{N}\left\langle\mathcal{T}_{\mathcal{C}}c_{\mu}\left(t_{1}\right)c_{\mu}^{\dagger}\left(t_{2}\right)\right\rangle$
(4)

where $\mathcal{T}_{\mathcal{C}}$ is the time ordering operator on the Keldysh
contour [16]. This Green’s function may be separated into the greater
($t_{1}>t_{2}$) and lesser ($t_{1}>t_{2}$) Green’s functions, which depict the
order of the operators. We use the ansatz for the Green’s functions of
large-$q$ complex SYK model [17, 4, 15]

${\cal G}^{\gtrless}(t_{1},t_{2})=({\cal Q}\mp
1/2)e^{g^{\gtrless}(t_{1},t_{2})/q}$ (5)

with the boundary conditions

$g^{\gtrless}(t,t)=0.$ (6)

Before proceeding further, we introduce some notations introduced earlier in
the literature [17, 4, 15] that will serve us later in our analysis. Our main
focus of interest will be in the so-called “symmetric” and “asymmetric” little
$g$ defined as follows:

$g^{\pm}(t_{1},t_{2})\equiv\frac{g^{>}(t_{1},t_{2})\pm
g^{<}(t_{2},t_{1})}{2}.$ (7)

Notice the reverse time argument in $g^{<}$ on the right-hand side where we
recall that
$g^{\gtrless}\left(t_{1},t_{2}\right)^{*}=g^{\gtrless}\left(t_{2},t_{1}\right)$.
We introduce the effective coupling coefficients for both the terms in the
Hamiltonian in eq. (1)

${\cal J}_{\kappa q}\equiv(1-4{\cal Q}^{2})^{\kappa q/4-1/2}J_{\kappa q}$ (8)

where again $\kappa\in\\{\frac{1}{2},1\\}$.

We defer the derivation of differential equations for $g^{\pm}$ to appendix A.
Recall that we are only interested in the equilibrium situation, implying that
the Green’s functions are only dependent on time differences, $t_{12}\equiv
t_{1}-t_{2}$, i.e., $g^{\pm}(t_{1},t_{2})\equiv g^{\pm}(t_{12})$. Using
appendix A, we get as solution for the Green’s function $g^{-}$ as follows:

$g^{-}(t_{1},t_{2})=g^{-}(t_{1}-t_{2})=g^{-}(t_{12}=2{\cal
Q}{\imath}\alpha\,t_{12}.$ (9)

Here, using eq. (89), we see that $\alpha$ is related to the derivative of the
symmetric part

$\dot{g}^{+}(0)={\imath}\alpha$ (10)

This symmetric part $g^{+}(t_{12})$ is the solution to a second order
differential equation. Since we only have to focus on time differences, focus
on the Green’s function with a single time argument which we label $t$

$\ddot{g}^{+}(t)=-2{\cal J}_{q}^{2}e^{g_{+}(t)}-2{\cal
J}_{q/2}^{2}e^{g_{+}(t)/2}$ (11)

Fortunately, this differential equation can be solved in closed form and the
solution is given by

$e^{g^{+}(t)/2}=\frac{1}{\left(\beta\mathcal{J}_{q/2}\right)^{2}}\frac{(\pi
v)^{2}}{1+\sqrt{A^{2}+1}\cos(\pi v(1/2-\imath t/\beta))}$ (12)

with the dimensionless parameter

$A\equiv\frac{\pi
v\beta\mathcal{J}_{q}}{\left(\beta\mathcal{J}_{q/2}\right)^{2}}$ (13)

We impose the boundary condition for $g$, namely $g^{+}(0)=0$ to get the
closure relation for $v$ as follows:

$\pi
v=\sqrt{\left(\beta\mathcal{J}_{q}\right)^{2}+\left(\frac{\left(\beta\mathcal{J}_{q/2}\right)^{2}}{\pi
v}\right)^{2}}\cos(\pi v/2)+\frac{\left(\beta\mathcal{J}_{q/2}\right)^{2}}{\pi
v}$ (14)

Using the first equation in eq. (89), we see that at equal times, we have

$\dot{g}^{+}(0)={\imath}\alpha$ (15)

where $\alpha$ is a constant in equilibrium whose details can be found in
appendix A. For convenience, we reproduce the definition of $\alpha$ here from
eq. (86). $\alpha$ is related to the expectation value of energy per particle
given by

$(1-{\cal Q}^{2})\alpha(t)=\epsilon_{q/2}(t)/2+\epsilon_{q}(t)$ (16)

where in the large-$q$ limit, we have ${\cal Q}\sim{\cal
O}(\frac{1}{\sqrt{q}})$ and $\epsilon_{\kappa q}(t)\equiv q^{2}\frac{J_{\kappa
q}\langle{\cal H}_{\kappa q}\rangle}{N}$ where $\kappa\in\\{\frac{1}{2},1\\}$.

This allows us to explicitly calculate $\alpha$ using the solution for $g^{+}$
in eq. (12) to get

$\beta\alpha=-2\sqrt{\frac{\left(\beta{\cal
J}_{q/2}\right)^{4}}{\pi^{2}v^{2}}+\left(\beta{\cal J}_{q}\right)^{2}}\sin(\pi
v/2)$ (17)

where we have made use of the closure relation eq. (14).

##  III Equation of State, Grand Potential and Phase Diagram

We start by calculating the equation of state that is true for all
temperature, and then identifying the region of interest where we can observe
a phase transition. Following the logic in [17, 4], we impose the Kubo-Martin-
Schwinger (KMS) relation [18] on the full Green’s function of the system that
leads us to the equation of state for our system given by

$\beta\mu=\ln\left[\frac{1+2{\cal Q}}{1-2{\cal Q}}\right]-2{\cal
Q}\beta\alpha/q$ (18)

We immediately notice that if $\beta\alpha={\cal O}(q^{0})$ and ${\cal
Q}\sim{\cal O}(q^{-1/2})$, the second term will be suppressed in large-$q$
limit. Therefore, at intermediate and large temperatures, we are left with the
equation of state of a free system only. But if we rescale the temperature as
${\tilde{\beta}}=\beta/q$ such that ${\tilde{\beta}}\sim{\cal O}(q^{0})$, then
we allow for the second term in the equation of state to contribute as well.
Given this re-scaling, we are accordingly interested in the low-
temperregime,regime as we are considering the large-$q$ limit.

Note that the effective interactions are exponentially suppressed as ${\cal
J}_{\kappa q}\sim e^{-\kappa q{\cal Q}^{2}}$ (using eq. (8)) for non-rescaled
charge densities. As such, our first step considers charge densities

${\cal Q}=\frac{\tilde{{\cal Q}}}{q^{1/2}},\qquad\tilde{{\cal Q}}={\cal
O}(q^{0})$ (19)

As noted above, we also focus on the low-temperature regime

$T={\tilde{T}}/q={\cal O}(q^{-1}),\quad\beta=q\tilde{\beta}$ (20)

and a rescaled chemical potential

$\mu=\tilde{\mu}q^{-3/2}$ (21)

where there is a competition between the interaction terms and the free system
contribution. To better understand the equation of state in the low-
temperature regime, we need to first find $v$ using eq. (14) in that limit.

###  III.1 Finding v in low-temperature limit

We must find the solution $v$ to eq. (14) in the low-temperature regime. In
the case where we set $J_{q/2}=0$, we would be left purely with $v\sim 1$.
This changes here, however. Let us define

$U\equiv\tilde{\beta}{\cal J}_{q},\qquad K\equiv\tilde{\beta}{\cal J}_{q/2}$
(22)

leaving

$A=\frac{\pi vU}{qK^{2}}={\cal O}(q^{-1}).$ (23)

With this, the eq. (14) takes the form

$\pi v=\sqrt{\frac{U^{2}}{q^{2}}+\left(\frac{K^{2}}{\pi
v}\right)^{2}}q^{2}\cos(\pi v/2)+\frac{q^{2}K^{2}}{\pi v},$ (24)

where we can interpret $U$ as an interaction term while $K$ as a pseudo-
kinetic term. We take both of the dimensionless couplings to be of order
$q^{0}$. From the above we note how the pseudo-kinetic term, corresponding to
$q/2$-body interactions, is the relevant term, since it enters in at a higher
order of $q$. If we however first set $K=0$, then this would yield the
standard closure relation [1] $\pi v=Uq\cos(\pi v/2)$ solved by $v=1+{\cal
O}(1/q)$. Since this special case has already been treated [4, 3] extensively,
we focus on the case of non-zero $K$ at large-$q$

$\pi v=\left[\frac{K^{2}}{\pi v}+\frac{U^{2}\pi v}{2q^{2}K^{2}}+{\cal
O}(1/q^{4})\right]q^{2}\cos(\pi v/2)+\frac{q^{2}K^{2}}{\pi v},$ (25)

with the solution

$v=2-\frac{4}{qK}\sqrt{2+(U/K)^{2}}+{\cal O}(q^{-2}).$ (26)

We realize that for our model, $v\in[0,2]$. We know that in the $K=0$ case
that the low-temperature limit would yield $v\to 1$ [1]. We will return to
this doubling of $v$ in (26) from the $K=0$ case when we discuss Lyapunov
exponents in Section V.

###  III.2 Equation of state

Since we have $\alpha$ in the equation of state (eq. (18)), we make use of eq.
(26) together with eq. (17) to get $\alpha$ in the low-temperature limit as
follows:

$\beta\alpha/q=-2\sqrt{\frac{q^{2}K^{4}}{\pi^{2}v^{2}}+U^{2}}\sin(\pi
v/2)\sim-2\sqrt{2K^{2}+U^{2}}.$ (27)

Then using the scaling of different quantities mentioned in eqs. (19), (20)
and (21), we are left with the equation of state in the low-temperature limit

$\tilde{\beta}\tilde{\mu}=4\tilde{{\cal Q}}\Big{[}1+\sqrt{2K({\cal
Q})^{2}+U({\cal Q})^{2}}\Big{]}$ (28)

which notably reduces to the standard large-$q$ complex SYK model for $K=0$
[4] as it should.

We now consider the equation of state where we first substitute from the
respective definitions (namely, eqs. (8) and (22) in large-$q$ scaled
temperature limit)

$U={\tilde{\beta}}J_{q}\exp\left(-{\tilde{\mathcal{Q}}}^{2}\right),K={\tilde{\beta}}J_{q/2}\exp\left(-\frac{{\tilde{\mathcal{Q}}}^{2}}{2}\right),{\tilde{\beta}}=\frac{1}{{\tilde{T}}}$
(29)

to get

$\frac{{\tilde{\mu}}}{J_{q}}=4\tilde{{\cal
Q}}\left[\sqrt{2\left(\frac{J_{q/2}}{J_{q}}\right)^{2}e^{-{\tilde{\mathcal{Q}}}^{2}}+e^{-2{\tilde{\mathcal{Q}}}^{2}}}+{\tilde{T}}/J_{q}\right]$
(30)

We introduce a new dimensionless parameter of the system that acts as a new
scale compared to a single large-$q$ complex SYK model, namely

$\gamma\equiv\frac{J_{q/2}}{J_{q}}$ (31)

which makes our equation of state as

$\frac{{\tilde{\mu}}}{J_{q}}=4\tilde{{\cal
Q}}\left[\sqrt{2\gamma^{2}e^{-{\tilde{\mathcal{Q}}}^{2}}+e^{-2{\tilde{\mathcal{Q}}}^{2}}}+{\tilde{T}}/J_{q}\right]$
(32)

Let us now consider this function plotted in Fig. 1 for various temperatures,
having set $J_{q}=1$. We clearly observe a phase transition in the (scaled)
low-temperature regime where there exists a critical temperature $T_{c}$ below
which there are three solutions while above which there is only one. This is
very reminiscent of the van der Waals phase transition and indeed, as we will
see in section IV, the critical exponents do belong to the same universality
class as that of van der Waals for all values of $\gamma=J_{q/2}/J_{q}$.

Figure 1: Plots of the equation of state (eq. (32)) for various temperatures
for $\gamma^{2}=0.28$ and $J_{q}=1$. The blue and orange curves are both one-
to-one functions, each corresponding to the supercritical phase at $T\geq
T_{c}$. The three colors in the bottom curve correspond to $3$ different
phases. Pink for the low density gaseous phase, green for the unstable phase,
and blue for the high density liquid phase.

###  III.3 Grand potential

The grand potential is given by $\Omega=E-\mu\mathcal{Q}-T\mathcal{S}$ where
${\cal S}$ is the entropy density. We already have all the required relations
except the entropy density.

The total energy $E=q^{2}\langle{\cal H}\rangle$ can be calculated via the
effective action by considering functional derivatives. This process yields
integrals over the self-energy, which we evaluate in appendix B to yield the
expression in eq. (107) which we reproduce here for convenience

$q\beta
E=-\frac{4K^{2}}{U}\coth^{-1}\left[\sqrt{2\frac{K^{2}}{U^{2}}+1}\right]-2U\sqrt{2\frac{K^{2}}{U^{2}}+1}.$
(33)

Motivated by the energy contributions, we consider an ansatz for ${\cal S}$
where the leading order is the free entropy of the system (without any
interaction) and seek the leading order in $1/q$ correction. To find the
leading order correction to the free entropy contribution, we use the
Maxwell’s relation

$-\left(\frac{{\partial}{\cal S}}{{\partial}{\cal
Q}}\right)_{T,J_{q},J_{q/2}}=\left(\frac{{\partial}\mu}{{\partial}T}\right)_{{\cal
Q},J_{q},J_{q/2}}$ (34)

where we already know $\mu$ using the equation of state. We find that the
${\cal O}(1/q)$ contribution of entropy density is at the same order of $q$ as
the energy term and the explicit expression for entropy density is given by

${\cal S}=\ln 2-2\tilde{Q}^{2}/q+{\cal O}(q^{-2})$ (35)

We proceed to calculate the dimensionless contribution $q\beta\Omega$ which is
of order ${\cal O}(q^{0})$. We label this as $q\beta{\tilde{\Omega}}$ where we
subtract off the logarithmic contribution to the entropy, namely
$\tilde{\beta}{\tilde{\Omega}}=q\beta\Omega+q\ln 2$, thereby leaving us with
after plugging the entropy density (recall $\tilde{{\cal
Q}}\equiv\sqrt{q}{\cal Q}={\cal O}(q^{0})$ in eq. (19))

$q\beta{\tilde{\Omega}}=q\beta E-q\beta\mu{\cal Q}+2\tilde{{\cal Q}}^{2}$ (36)

Next we use the equation of state at low-temperature given in eq. (28)
combined with the scaling for $\beta$ and $\mu$ in eqs. (20) and (21),
respectively. The resulting expression for $q\beta\mu{\cal Q}$ is of the order
${\cal O}(q^{0})$ and is given by

$q\beta\mu{\cal Q}=4{\tilde{\mathcal{Q}}}^{2}[1+U\sqrt{1+2(K/U)^{2}}]$ (37)

We already have the expression for $q\beta E$ in eq. (33) (or eq. (107) as
derived in appendix B) where we see that $q\beta E\sim{\cal O}(q^{0})$.
Therefore, plugging the equation of state and the energy contributions, we
finally get for the grand potential

$\tilde{\beta}{\tilde{\Omega}}=-2(1+2{\tilde{\mathcal{Q}}}^{2})U\sqrt{1+2(K/U)^{2}}-\frac{4K^{2}}{U}\coth^{-1}\left[\sqrt{[2(K/U)^{2}+1]}\right]-2\tilde{{\cal
Q}}^{2}$ (38)

where $q\beta{\tilde{\Omega}}\sim{\cal O}(q^{0})$.

Next we wish to the explicit dependence of ${\tilde{\Omega}}$ on
${\tilde{\mathcal{Q}}}$ and ${\tilde{T}}$. We use eq. (29) and recall from eq.
(31) that $\gamma\equiv J_{q/2}/J_{q}$ to get as the final expression for the
grand potential

$\tilde{\Omega}/J_{q}=-2(1+2\tilde{\mathcal{Q}}^{2})e^{-\tilde{\mathcal{Q}}^{2}}\sqrt{1+2\gamma^{2}e^{\tilde{\mathcal{Q}}^{2}}}-4\gamma^{2}\coth^{-1}\left[\sqrt{2\gamma^{2}e^{\tilde{\mathcal{Q}}^{2}}+1}\right]-2\tilde{\mathcal{Q}}^{2}\tilde{T}/J_{q}$
(39)

where in the limit $\gamma\rightarrow 0$, we recover the grand potential for a
single large-$q$ complex SYK model (eq. (6) of [4]) as it should.

##  IV Phase diagram

Having obtained the grand potential at low temperature, we are in a position
to study the phase transition as depicted in Fig. 1. There exists a critical
temperature ${\tilde{T}_{c}}$ that corresponds to the point of inflection, or
in other words, when both the first derivative and the second derivative of
${\tilde{\mu}}$ with respect to ${\tilde{\mathcal{Q}}}$ vanish.

To find ${\tilde{T}_{c}}$, we start with finding the critical charge density
${\tilde{\mathcal{Q}}_{c}}$ which is obtained by imposing that the second
derivative of ${\tilde{\mu}}$, given in eq. (32), with respect to ${\cal Q}$
vanish, namely ${\tilde{\mu}}^{(2)}({\tilde{\mathcal{Q}}_{c}})=0$. This
equation is satisfied trivially for ${\tilde{\mathcal{Q}}_{c}}=0$ or when

$2\gamma^{4}e^{2{\tilde{\mathcal{Q}}_{c}}^{2}}\left({\tilde{\mathcal{Q}}_{c}}^{2}-3\right)+3\gamma^{2}e^{{\tilde{\mathcal{Q}}_{c}}^{2}}\left(2{\tilde{\mathcal{Q}}_{c}}^{2}-3\right)+2{\tilde{\mathcal{Q}}_{c}}^{2}-3=0$
(40)

Since we cannot find a simple closed form solution for
${\tilde{\mathcal{Q}}_{c}}$, we solve for $\gamma$ as a function of
${\tilde{\mathcal{Q}}_{c}}$. The positive solution is given by

$\gamma^{2}({\tilde{\mathcal{Q}}_{c}})=e^{-{\tilde{\mathcal{Q}}_{c}}^{2}}\frac{-6{\tilde{\mathcal{Q}}_{c}}^{2}+9-\sqrt{20{\tilde{\mathcal{Q}}_{c}}^{4}-36{\tilde{\mathcal{Q}}_{c}}^{2}+9}}{4\left({\tilde{\mathcal{Q}}_{c}}^{2}-3\right)}$
(41)

We note that $\gamma^{2}({\tilde{\mathcal{Q}}_{c}})\in[0,\infty)$ and
therefore

${\tilde{\mathcal{Q}}_{c}}\in\left[\sqrt{\frac{3}{2}},\sqrt{3}\right].$ (42)

Indeed for the case of $\gamma^{2}({\tilde{\mathcal{Q}}_{c}})=0$, we get
${\tilde{\mathcal{Q}}_{c}}=\sqrt{3/2}$ which matches with the case of single
large-$q$ complex SYK model [4] as expected.

The critical temperature is found by solving
$\frac{\partial{\tilde{\mu}}}{\partial{\tilde{\mathcal{Q}}_{c}}}=0$, which
yields

$\frac{{\tilde{T}}_{c}}{J_{q}}=e^{-{\tilde{\mathcal{Q}}_{c}}^{2}}\frac{2{\tilde{\mathcal{Q}}_{c}}^{2}-1+2e^{{\tilde{\mathcal{Q}}_{c}}^{2}}\gamma^{2}({\tilde{\mathcal{Q}}_{c}})({\tilde{\mathcal{Q}}_{c}}^{2}-1)}{\sqrt{1+2e^{{\tilde{\mathcal{Q}}_{c}}^{2}}\gamma^{2}({\tilde{\mathcal{Q}}_{c}}^{2})}}$
(43)

With this, we can now plot the parts of the EOS and grand potential
corresponding to the various phases as given in Fig. 2.

Figure 2: Plots of the equation of state (eq. (32)) and grand potential (eq.
(39)) for $\gamma^{2}=0.28$, at temperature $T=0.6T_{c}$, and $J_{q}=1$.

###  IV.1 Universality Class

Finally, we go near the critical point on the phase diagram and expand both
the chemical potential (eq. (32)) and grand potential (eq. (39)) around the
critical point, including
${\tilde{\mu}_{c}}\equiv{\tilde{\mu}}({\tilde{\mathcal{Q}}_{c}})|_{T=T_{c}}$,
to get the critical exponents. To do so, we first define the reduced shifted
set of variables as follows:

$m\equiv\frac{1}{a_{c}}\left[\frac{{\tilde{\mu}}}{{\tilde{\mu}_{c}}}-1\right],\;\,\rho\equiv\frac{{\tilde{\mathcal{Q}}}}{{\tilde{\mathcal{Q}}_{c}}}-1,\;\,{\tilde{T}}\equiv{\tilde{T}_{c}}\left[1+\frac{t}{b_{c}}\right]$
(44)

such that the critical point corresponds to $m=\rho=t=0$. We have further
defined

$\displaystyle a_{c}$
$\displaystyle\equiv\frac{1+4{\tilde{\mathcal{Q}}_{c}}^{2}(3-5{\tilde{\mathcal{Q}}_{c}}^{2}/3)+3(2{\tilde{\mathcal{Q}}_{c}}^{2}-3)\sqrt{(2{\tilde{\mathcal{Q}}_{c}}^{2}/3-1)(10{\tilde{\mathcal{Q}}_{c}}^{2}/3-1)}}{4}$
(45) $\displaystyle b_{c}$
$\displaystyle\equiv\frac{1-4{\tilde{\mathcal{Q}}_{c}}^{2}+24(2{\tilde{\mathcal{Q}}_{c}}^{2}/3)^{2}-\sqrt{(2{\tilde{\mathcal{Q}}_{c}}^{2}/3-1)(10{\tilde{\mathcal{Q}}_{c}}^{2}/3-1)}}{16(2{\tilde{\mathcal{Q}}_{c}}^{2}/3)^{2}a({\tilde{\mathcal{Q}}_{c}})}.$

We have also rescaled the reduced shifted variables such that the overlap with
the solved case of a single SYK model becomes more apparent. This is the
origin of the constants
$a({\tilde{\mathcal{Q}}_{c}}),b({\tilde{\mathcal{Q}}_{c}})$. Such a rescaling
and shifting also enters into the grand potential

$f\equiv\frac{\tilde{\Omega}-\tilde{\Omega}_{c}}{{\tilde{\mu}_{c}}{\tilde{\mathcal{Q}}_{c}}a_{c}}+m-t/3$
(46)

We will focus on the expansion around the critical order parameter,
corresponding to $\rho=0$. Using the explicit expressions in eq. (39) for the
grand potential and eq. (32) for the equation of state, we expand the reduced
grand potential $f$ as well as the reduced equation of state $m$ to yield

$\displaystyle f$ $\displaystyle=-\rho^{2}\frac{t+(3\rho/2)^{2}}{3}+{\cal
O}(\rho^{5})$ (47) $\displaystyle m$
$\displaystyle=2t/3+\rho(2t/3+\rho^{2})+{\cal O}(\rho^{4})$

respectively. At this point, we note that we are left with the same equations
as a single large-$q$ complex SYK model [4]. As such, we also have to consider
field mixing to give the correct scaling function and critical exponents. By
this, we mean that the correct ordering field $\mathscr{h}$ is not the
chemical potential, but rather also gains a contribution from the temperature
field $\mathscr{h}\equiv m-2t/3$. We then solve for $\rho$ in terms of
$\mathscr{h}$ and $t$ in various regimes of small $\mathscr{h}$ and $t$ and
use this to express the reduced grand potential. For instance, around
$\mathscr{h}=0$ for small values of $\mathscr{h}$ (but not necessarily small
$t$) we get

$\displaystyle f(t,\mathscr{h})$
$\displaystyle=-\frac{|t|^{2}}{6}-|\mathscr{h}||2t/3|^{1/2}-\frac{3\mathscr{h}^{2}}{2}|t|^{-1}+\mathcal{O}\left(\mathscr{h}^{3}|t|^{-5/2}\right)$
(48) $\displaystyle\rho(t,\mathscr{h})$
$\displaystyle=\operatorname{sgn}(\mathscr{h})|2t/3|^{1/2}+2\mathscr{h}|2t/3|^{-1}+\mathcal{O}\left(\mathscr{h}^{2}t^{-5/2}\right)$

respectively, where note that we are considering $\mathscr{h}$ to be small but
not necessarily $t$ such that the term inside the big-${\cal O}$ notation is
small. For small $t$ around $t=0$ (but not necessarily small $\mathscr{h}$),
we also have

$\displaystyle f(t,\mathscr{h})$
$\displaystyle=-\frac{3}{4}|\mathscr{h}|^{4/3}\left[1+\mathcal{O}\left(t\mathscr{h}^{-2/3}\right)\right]$
(49) $\displaystyle\rho(t,\mathscr{h})$
$\displaystyle=-\partial_{\hbar}f(t,\mathscr{h})$

where $f(t,\mathscr{h})$ used in the derivative is the one that has been
expanded for small $t$ in eq. (49). Then we have the critical exponents
$\alpha$ calculated through specific heat
$C_{\mathscr{h}}\propto-\partial_{t}^{2}f(t,0)\propto|t|^{-\alpha}$ (using $f$
in eq. (48)), $\beta$ through order parameter $\rho(t,0)\propto|t|^{\beta}$
(using $\rho$ in eq. (48)), $\gamma$ through susceptibility
$\left.\chi_{\mathscr{h}}\propto\partial_{\mathscr{h}}^{2}f(t,\mathscr{h})\right|_{\mathscr{h}=0}\propto|t|^{-\gamma}$
(using $f$ in eq. (48)) and $\delta$ through
$\rho(0,\mathscr{h})\propto\mathscr{h}^{1/\delta}$ (using $\rho$ in eq. (49)).
The critical exponents are given in table 1. These are mean-field critical
exponents, and therefore our model belongs to the same universality class as
van der Waals.

Table 1: Critical exponents $\alpha$ | $\beta$ | $\gamma$ | $\delta$
---|---|---|---
0 | $\frac{1}{2}$ | 1 | 3

Surprisingly, even though the expressions for the equation of state (eq. (32))
and the grand potential (eq. (39)) are different from the case of a single
large-$q$ complex SYK model due to the presence of another length scale
$\gamma$, we get the exact same scaling of both of them around the critical
point as in the case of a single large-$q$ SYK model. In other words, all the
above critical point expansions of the reduced equation of state $m$, the
reduced charge density or equivalently the order parameter $\rho$ and the
reduced grand potential $f$ are the same as that for a single large-$q$
complex SYK model. Therefore, our system of coupled large-$q$ SYK models
belong to the same (Landau-Ginzburg mean-field) universality class as that of
a single large-$q$ complex SYK model, which in turn belongs to the van der
Waals (mean-field) universality class.

##  V Dynamics and chaos

In this section we will compute the Lyapunov exponent for our coupled system
of complex large-$q$ SYK models and show that the model saturates the
Maldacena-Shenker-Stanford (MSS) bound at low temperature [19], thereby being
maximally chaotic to order $\order{1/q}$ in the low temperature limit.

###  V.1 4-point function

In order to obtain the Lyapunov exponent, we first analyze the four point
function for the coupled system. In the large-$N$ limit, the correlator for
the complex fermions $c$ is given by [1, 20]

$\displaystyle\frac{1}{N^{2}}\sum_{j,k}^{N}\langle T$ $\displaystyle
c_{j}(\tau_{1})c_{j}^{\dagger}(\tau_{2})c_{k}^{\dagger}(\tau_{3})c_{k}(\tau_{4})\rangle=$
(50) $\displaystyle{\cal G}(\tau_{12}){\cal
G}(\tau_{34})+\frac{1}{N}\mathcal{F}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})+\dots$

where the indices $j,k$ are the fermionic sites, and
$\mathcal{F}=\sum_{n}\mathcal{F}_{n}$ is the unregularized out-of-time-order
correlator (OTOC). The first term is a disconnected diagram, and we are
interested in the $\order{\frac{1}{N}}$ diagrams. Here, the disorder average
is non-vanishing if the number of indices for the two couplings
$X_{\mu_{1}\dots\mu_{\kappa q/2}}^{\nu_{1}\dots\nu_{\kappa q/2}}$ (where
$\kappa=\\{\frac{1}{2},1\\}$) are the same. As a result, the average over the
product of the $q$ interaction coupling with the $q/2$ interaction is always
zero. E.g. for $q=4$, the disorder average
$\expectationvalue{X_{\mu_{1}\mu_{2}}^{\nu_{1}\nu_{2}}\;X_{\mu_{1}}^{\nu_{1}}}=0$.

Unlike the single large-$q$ complex SYK model, where a single
$\mathcal{F}_{n}$ in the sum denotes a single diagram [20], in our model each
$\mathcal{F}_{n}$ represents a sum of ladder diagrams as depicted in Fig. (3).
This is after resumming the perturbative terms to get the ladder diagrams in
terms of dressed propagators [21].

$\tau_{3}$$\tau_{2}$$\tau_{1}$$\tau_{4}$$\mathcal{F}_{0}$+$\tau_{3}$$\tau_{2}$$\tau_{1}$$\tau_{4}$$J_{q}$$J_{q}$+$\tau_{3}$$\tau_{2}$$\tau_{1}$$\tau_{4}$$J_{q/2}$$J_{q/2}$+$\mathcal{F}_{1}$$J_{q}$$J_{q}$$J_{q}$$J_{q}$+$J_{q/2}$$J_{q/2}$$J_{q/2}$$J_{q/2}$+$J_{q}$$J_{q}$$J_{q/2}$$J_{q/2}$+$J_{q}$$J_{q}$$J_{q/2}$$J_{q/2}$$\mathcal{F}_{2}$
Figure 3: The diagrams related to the $\order{\frac{1}{N}}$ part of the
correlator for our system of coupled complex large-$q$ SYK models. The diagram
is drawn for $q=8$, and also consist of diagram with
$\tau_{3}\leftrightarrow\tau_{4}$, which have been omitted here. There are
$q-2$ lines connecting the two horizontal rails out of which Half of them
($q/2-1$ lines) run in one direction while the remaining half ($q/2-1$ lines)
run in the opposite.

According to standard techniques [2, 20, 1], each diagram $\mathcal{F}_{n+1}$
is generated by applying a kernel $K$ to the previous diagram
$\mathcal{F}_{n}$, i.e.,

$\mathcal{F}_{n+1}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})=\int d\tau
d\tau^{\prime}K(\tau_{1},\tau_{2},\tau,\tau^{\prime})\mathcal{F}_{n}(\tau,\tau^{\prime},\tau_{3},\tau_{4})$
(51)

or $\mathcal{F}_{n+1}=K\odot\mathcal{F}_{n}$ in short. We can recursively
repeat the process and assuming asymptotic exponential growth of the OTOC. The
zero rung diagrams are also negligible at large times, which gives an
eigenvalue problem of type $\mathcal{F}=K\odot F$ (see below for more
details). The kernel for the coupled system at equilibrium is given by

$K(\tau_{1},\tau_{2};\tau,\tau^{\prime})\equiv-\sum_{\kappa=\\{\frac{1}{2},1\\}}\frac{2^{\kappa
q-1}J_{\kappa q}^{2}}{\kappa q}(\kappa q-1)\;{\cal G}(\tau_{1}-\tau){\cal
G}(\tau_{2}-\tau^{\prime}){\cal G}^{>}(\tau-\tau^{\prime})^{\kappa q/2-1}{\cal
G}^{<}(\tau^{\prime}-\tau)^{\kappa q/2-1}$ (52)

and the generation of the ladder diagrams with the kernel is illustrated in
Fig. (4). A notable result is that our kernel does not merely add a rung to
the ladder diagram. It adds a sum of diagrams, where each diagram is extended
with one rung.

$J_{q}$$J_{q}$+$J_{q/2}$$J_{q/2}$=$J_{q}$$J_{q}$+$J_{q/2}$$J_{q/2}$$\times$
Figure 4: In the large-$N$ limit, the diagrams of order $\order{\frac{1}{N}}$
in the out-of-time-correlator (OTOC) for the coupled SYK system can be
constructed with a kernel $K$. The OTOC is then written as
$\mathcal{F}=\mathcal{F}_{0}+K\mathcal{F}$, where
$\mathcal{F}=\sum\limits_{n=0}^{\infty}\mathcal{F}_{n}$, where
$\mathcal{F}_{0}$ is defined in Fig. 3 and becomes suppressed at large time
due to exponential growth of the OTOC. This figure shows the process
$\mathcal{F}_{1}=K\mathcal{F}_{0}$ for $q=8$, where the term in parentheses is
the kernel.

###  V.2 Chaos regime

We wish to obtain the exponential growth rate of the OTOC at late times. For
this, we utilize the standard approach [2, 19] by considering the following
regularized out-of-time-order correlator (OTOC) in real time as

$\mathcal{F}(t_{1},t_{2},t_{3},t_{4})=\Tr\left[y\;c_{j}(t_{1})\,y\,c_{k}^{\dagger}(t_{3})\,y\,c_{j}^{\dagger}(t_{2})\,y\,c_{k}(t_{4})\right],$
(53)

having defined $y\equiv\rho(\beta)^{1/4}$, where $\rho(\beta)$ is the thermal
density matrix. In this approach, the fermions are separated by a quarter of
the thermal circle, and a real-time separation between the pairs of fermions.
More on the specifics of this approach can be found in appendix C. We notice
that the kernel generates all the infinite diagrams, with the exception of the
free propagator, i.e.,

$\sum\limits_{n=1}^{\infty}\mathcal{F}_{n}=K\odot\sum\limits_{n=0}^{\infty}\mathcal{F}_{n}.$
(54)

Adding $\mathcal{F}_{0}$ on both sides and assuming asymptotic exponential
growth rate of OTOC, the zero rung diagram $\mathcal{F}_{0}$ gets suppressed
in the limit $t_{1},t_{2}\to\infty$. We get an eigenvalue equation of the type
[20, 1] $\mathcal{F}=K\odot\mathcal{F}$, namely

$\mathcal{F}(t_{1},t_{2},t_{3},t_{4})=\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$\infty$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t$\\!
}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$\infty$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t^{\prime}$\\!
}K_{R}(t_{1},t_{2};t,t^{\prime})\mathcal{F}(t,t^{\prime},t_{3},t_{4}).$ (55)

The exponential growth of the OTOC is determined by the real-time part of
$\mathcal{F}$ , and the diagrams describing the OTOC are generated by the
retarded kernel. Given the time translational invariance for our equilibrium
setup, we have (with the notation from section II, namely $t_{ij}\equiv
t_{i}-t_{j}$)

$K_{R}(t_{1},t_{2};t_{3},t_{4})=\sum_{\kappa=\\{\frac{1}{2},1\\}}\frac{2^{\kappa
q-1}J_{\kappa q}^{2}}{\kappa q}(\kappa q-1){\cal G}_{R}(t_{13}){\cal
G}_{R}(-t_{24}){\cal G}_{lr}^{>}(t_{34})^{\kappa q/2-1}{\cal
G}_{lr}^{<}(-t_{34})^{\kappa q/2-1},$ (56)

where ${\cal G}_{R}$ and ${\cal G}_{lr}$ are the retarded propagator and the
Wightman correlator respectively, defined as

$\displaystyle{\cal G}_{R}(t)$ $\displaystyle\equiv[{\cal G}^{>}(t)-{\cal
G}^{<}(t)]\;\theta(t)\stackrel{{\scriptstyle
q\rightarrow\infty}}{{=}}\theta(t),$ (57) $\displaystyle{\cal G}^{>}_{lr}(t)$
$\displaystyle\equiv{\cal G}^{>}(\imath t+\beta/2),$ $\displaystyle{\cal
G}^{<}_{lr}(t)$ $\displaystyle\equiv{\cal G}^{<}(\imath t-\beta/2).$

Note the relation in frequency domain between the Matsubara Green’s function
and real time retarded Green’s function, namely ${\cal
G}({\imath}\omega_{n}\to\omega+{\imath}0^{+})={\cal G}_{R}(\omega)$ [22].
Compared with the kernel in eq. (52), we see that the overall minus sign
vanishes. This is due to the transformation $\tau\to it$ and $\tau^{\prime}\to
it^{\prime}$, leading to $d\tau d\tau^{\prime}=-dtdt^{\prime}$ in the integral
of type of eq. (55), where the minus sign is absorbed in the retarded kernel.
For a detailed discussion on the retarded kernel and the chaos exponents in
real time, see [23] (in particular, section 8).

The explicit form of the Wightman function can be found with the Green’s
function in eqs. (5) and (12). Combined with the effective coupling (8), we
have (recall $\cosh(-t)=\cosh(t)$)

$K_{R}(t_{1},t_{2};t_{3},t_{4})=\theta(t_{13})\theta(t_{24})B(t_{34}),$ (58)

where

$\displaystyle B(t_{34})\equiv$ $\displaystyle\frac{2(\pi
v)^{2}A^{2}\;\theta(t_{13})\theta(t_{24})}{\beta^{2}\left[1+\sqrt{A^{2}+1}\cosh(\frac{\pi
vt_{34}}{\beta})\right]^{2}}$ (59) $\displaystyle+\frac{(\pi
v)^{2}\;\theta(t_{13})\theta(t_{24})}{\beta^{2}\left[1+\sqrt{A^{2}+1}\cosh(\frac{\pi
vt_{34}}{\beta})\right]}.$

Recall the definition of $A$ from eq. (13), namely $A\equiv\frac{\pi
v\beta\mathcal{J}_{q}}{\left(\beta\mathcal{J}_{q/2}\right)^{2}}$ $=\frac{\pi
vU}{qK^{2}}$ where we used definitions in eq. (22) and the scaling for $\beta$
in eq. (20). It is important to notice that even if the kernel applies
multiple diagrams, the total sum of infinite diagrams in $\mathcal{F}$ does
not change with applying the retarded kernel. Furthermore, the zero rung
ladder diagram gets suppressed at larger times. This implies that in the
coupled SYK model, $\mathcal{F}$ is also an eigenfunction of the integral
transform that can be further simplified by taking the two derivatives
$\partial_{t_{1}}\partial_{t_{2}}$ on both sides of the integral equation (eq.
(55)) to get

$\partial_{t_{1}}\partial_{t_{2}}\mathcal{F}(t_{1},t_{2},t_{3},t_{4})=\int_{-\infty}^{\infty}dt\int_{-\infty}^{\infty}dt^{\prime}\;\delta(t_{1}-t)\delta(t_{2}-t^{\prime})B(t_{34})\;\mathcal{F}(t,t^{\prime},t_{3},t_{4}).$
(60)

The delta functions enforces $t=t_{1}$ and $t^{\prime}=t_{2}$. Due to the
theta function and our interest in large times $t_{1}$ and $t_{2}$, we can put
$t_{3}=t_{4}=0$ without loss of generality. For brevity, we write
$\mathcal{F}(t_{1},t_{2},t_{3}=0,t_{4}=0)=\mathcal{F}(t_{1},t_{2})$.

###  V.3 The Lyapunov exponent

To solve this integral in real time, assuming a chaotic behavior at asymptotic
times, we make the standard ansatz [2, 1, 20]

$\mathcal{F}(t_{1},t_{2})=e^{\lambda_{L}(t_{1}+t_{2})/2}f(t_{12}),$ (61)

where $\lambda_{L}$ is the growth exponent, known as the Lyapunov exponent.
With the ansatz (61), and assuming the function $\mathcal{F}(t_{1},t_{2})$
decays at the boundaries $t_{1},t_{2}\to-\infty$, this reduces to a second
order differential equation

$\left[\partial_{t_{1}}\partial_{t_{2}}-B(t_{12})\right]f(t_{12})=-\frac{\lambda_{L}^{2}}{4}f(t_{12}).$
(62)

In fact, it is a “Schrödinger equation” with a potential

$V(x)\equiv-2B(t_{12})/(\pi vT)^{2},\qquad x\equiv\pi vTt_{12}/2.$ (63)

From the chain rule we then have

$\partial_{t_{1}}\partial_{t_{2}}=\frac{\partial
x}{{\partial}t_{1}}\frac{\partial
x}{{\partial}t_{2}}{\partial}_{x}^{2}=-\frac{(\pi
vT)^{2}}{2}{\partial}_{x}^{2}/2$

which leaves us with

$\left[-{\partial}_{x}^{2}/2+V(x)\right]f(t_{12})=-\frac{\lambda_{L}^{2}}{2(\pi
vT)^{2}}f(t_{12}).$ (64)

While there typically exists an entire spectrum of Lyapunov exponents, the
systems chaos is defined by the maximum $\lambda_{L}$. This corresponds to the
minimal negative eigenvalue, which would correspond to the ground state energy
$E_{0}$. As such, we have a ground state wave function $\psi_{0}(x)\equiv
f(t_{12})$ and energy

$E_{0}\equiv-\frac{\lambda_{L}^{2}}{2[\pi vT]^{2}}.$ (65)

This leaves us with the more familiar form of the Schrödinger equation

$\left[-\frac{{\partial}_{x}^{2}}{2}+V(x)\right]\psi_{0}(x)=E_{0}\psi_{0}(x).$
(66)

####  V.3.1 Perturbation theory

Let us now consider the above Schrödinger type equation (66) and in particular
the potential, which may be written as

$\displaystyle V(x)=$
$\displaystyle\frac{1}{(\sqrt{A^{2}+1}-1)/2+\sqrt{A^{2}+1}\,/V_{0}(x)}$ (67)
$\displaystyle-\frac{A^{2}}{[(\sqrt{A^{2}+1}-1)/2+\sqrt{A^{2}+1}\,/V_{0}(x)]^{2}}$

where $A\equiv\frac{\pi
v\beta\mathcal{J}_{q}}{\left(\beta\mathcal{J}_{q/2}\right)^{2}}$ $=\frac{\pi
vU}{qK^{2}}$. We now consider this for small $A$, with corrections
$\psi_{0}=\psi_{0}^{(0)}+A^{2}\psi_{0}^{(1)}+\ldots$ and
$E_{0}=E_{0}^{(0)}+A^{2}E_{0}^{(1)}+\ldots$ obtained via perturbation theory.

For $A=0$, we are left with the Pöschl-Teller potential
$V_{0}(x)=-\sech^{2}(x)$ [24] and the ground state is

$\psi_{0}^{(0)}(x)=\sech(x)/\sqrt{2}.$ (68)

with corresponding eigenvalue is $E_{0}^{(0)}=-1/2$. The leading order
perturbed potential is

$V(x)=V_{0}(x)-\frac{V_{0}(x)}{4}(2+33V_{0}(x))A^{2}+{\cal O}(A^{4})$ (69)

As such, first order perturbation theory yields the energy correction
$E_{0}=-1/2+A^{2}E^{(1)}_{0}+{\cal O}(A^{4})$, with

$E^{(1)}_{0}=-\bra{\psi^{(0)}_{0}}\frac{\hat{V}_{0}}{4}(2+33\hat{V}_{0})\ket{\psi^{(0)}_{0}}.$
(70)

Inserting the identity $\int dx|x\rangle\langle x|=1$ yields

$\displaystyle E^{(1)}_{0}$
$\displaystyle=\int_{\mbox{\scriptsize$$}}^{\mbox{\scriptsize$$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}x$\\!
}\frac{V_{0}^{2}(x)}{8}(2+33V_{0}(x))$ (71)

This integrates to $E^{(1)}_{0}=-61/15$, hence the energy is lowered. Since we
have the Lyapunov exponent (65)

$\lambda_{L}=\sqrt{2|E_{0}|}\pi vT,$ (72)

one might assume that this leads to an increase to the exponents, hence
allowing one to surpass the MSS bound. Note however that this correction is
actually at quadratic order since

$A=\frac{\pi vU}{qK^{2}}$ (73)

whereas there is another correction to $v$ in eq. (26) which reduces the
Lyapunov exponent at linear order, namely

$v=2-\frac{4}{qK}\sqrt{2+(U/K)^{2}}+{\cal O}(q^{-2})$ (74)

As such the MSS bound $\lambda_{L}\to 2\pi T$ remains in tact in the low $T$
limit. We can now return to the topic of $v$ being twice that of the $q$-body
SYK model. The point is that $v$ is not necessarily a direct measure of chaos.
While we have $\lambda_{L}=\sqrt{2|E_{0}|}\pi vT$ remains true, if we had
instead first set $\gamma=0$, we would have found a potential

$V(x)\sim-\frac{1}{[1/2+1/V_{0}(x)]^{2}}=V_{0}(2x)/4.$

This then leads to a ground state energy of $E_{0}=2$ and $v\to 1$ in the low-
temperature regime [1]. Said another way, while $v$ might saturate to either
$1$ or $2$, in either of these cases, the Lyapunov exponent will saturate to
the MSS bound.

###  V.4 Relation to the phase transitions at “very” low temperature

Given that the MSS bound is always saturated over the phase transition
corresponding to the temperature regime $T={\cal O}(q^{-1})$ (20), we conclude
that both liquid and gaseous phases are maximally chaotic. In other words, we
have a maximally chaotic to maximally chaotic phase transition similar to
those found in small to large black hole phase transitions [3, 7].

This, however, changes in the rescaled regime of very low temperature

$T=\bar{T}q^{-2}={\cal O}(q^{-2}),$ (75)

which is the regime one must consider to carefully take the zero temperature
limit. In this regime, much like the standard $q$-body SYK model, the gaseous
and liquid densities diverge [4]. In fact, the liquid phase, the blue line in
Fig. 1, is at $Q_{\ell}={\cal O}(q^{0})$. Now recall that the effective SYK
interactions have a Gaussian suppression in charge density. As such, the large
charge density of the liquid eventually fully suppresses the SYK interactions,
leaving a free non-interacting liquid

$\tilde{\Omega}_{\ell}=\frac{\bar{T}}{2}\ln[1-4{\cal
Q}_{\ell}^{2}],\quad\bar{\mu}_{\ell}=2\bar{T}\tanh^{-1}(2{\cal Q}_{\ell}),$
(76)

where we have defined $\bar{\mu}\equiv\mu q^{-2}$.

With this in mind, let us consider the mathematical details around the phase
transition in this rescaled regime. The full EOS (18) may be written as

$\bar{\beta}\bar{\mu}=\ln\left[\frac{1+2{\cal Q}}{1-2{\cal Q}}\right]-2{\cal
Q}q\bar{\beta}\alpha.$ (77)

For the case where $\gamma=0$, we know that, at such temperatures and chemical
potentials, the phase transition is from $Q={\cal O}(1/q)\to{\cal O}(q^{0})$
[4]. Let us consider the EOS and the grand potential here for such charge
densities. As for the $\gamma=0$ case, for charge densities $Q={\cal
O}(q^{0})$, the Lyapunov exponent is exponentially suppressed, $v\sim
e^{-qQ^{2}}\to 0$, hence the system is free. Here we set $J_{q}=1$, since we
will not focus on the case where $J_{q}=0$.

For $Q=\bar{Q}q^{-1}$, we find the gaseous phase still remains in the
maximally chaotic phase. This leaves us with the gaseous part of the EOS

$\bar{\mu}_{g}\sim 4\bar{{\cal Q}}_{g}\sqrt{2\gamma^{2}+1}$ (78)

and the grand potential (39)

$\tilde{\Omega}_{g}=-2\sqrt{1+2\gamma^{2}}-4\gamma^{2}\coth^{-1}\sqrt{2\gamma^{2}+1}.$
(79)

The phase transition will be where (79) and (78) are equal to (76). Equating
the two expressions, we find that

${\cal
Q}_{\ell}=\frac{\sqrt{1-e^{2\bar{\beta}\tilde{\Omega}_{g}}}}{2}\xrightarrow[\bar{T}\to
0]{}1/2$

and

$\bar{{\cal Q}}_{g}=\bar{T}\frac{\tanh^{-1}(2{\cal
Q}_{\ell})}{2\sqrt{2\gamma^{2}+1}}\xrightarrow[\bar{T}\to
0]{}-\frac{\tilde{\Omega}_{g}}{2\sqrt{1+2\gamma^{2}}}$

We may now see explicitly how, for the liquid phase $Q_{\ell}={\cal
O}(q^{0})$, in this regime, both $K$ and $U$ are smaller that
$q\bar{T}e^{-qQ^{2}}$. The exponential suppression dominates over any
polynomial growth in $q$, implying that the closure relation (24)

$\pi v=\sqrt{\frac{U^{2}}{q^{2}}+\left(\frac{K^{2}}{\pi
v}\right)^{2}}q^{2}\cos(\pi v/2)+\frac{q^{2}K^{2}}{\pi v},$ (80)

reduces to $\pi v\sim 0$. As such, we have also confirmed that the free liquid
state is indeed regular (non-chaotic), with a Lyapunov exponent of
$\lambda_{L}=0$. As such, we are left with a maximally chaotic to regular
phase transition, which is similar to that of the standard Hawking-Page
transition from a large black hole to non-interacting radiation [8] but as
stated in the introduction, instead of having a particular value of
temperature like in Hawking-Page, we have an entire region of temperatures in
the very low temperature regime in our coupled SYK case.

##  VI Conclusion and Outlook

Motivated from the study of a single large-$q$ complex SYK model, we
considered a system of coupled SYK models in this work and studied its
thermodynamic and dynamic properties. Despite having a new scale in the
system, namely the ratio of interaction strengths of the two SYK terms in the
Hamiltonian, we found that there exists a continuous phase transition at low
temperatures just like a single SYK and the associated critical exponents
correspond to the same universality class as Landau-Ginzburg (mean-field)
exponents, also shared by van der Waals gas. Coincidentally, various AdS black
holes also belong to the same universality class [7, 25]. We will comment more
on this below. Finally, we calculated the Lyapunov exponent, and we again
found that the coupled SYK system also saturates the Maldacena-Shenker-
Stanford (MSS) bound in large-$q$ limit at low temperature, thereby making our
model maximally chaotic, surprisingly matching the behavior of a single SYK
model.

The coupled SYK system’s adherence to the established Landau-Ginzburg
universality class and its demonstration of maximal chaos are not merely
mathematical curiosities. They provide a critical lens through which the
fundamental nature of quantum chaos and gravity can be examined. This
universality across different models and configurations raises the question of
how one can obtain deviating results. Given that we flow to the exact same
grand potential scaling relations near a critical point for the combination of
two SYK models, the natural question to ask is to what extent this result can
be generalized. Does it hold for any system of coupled large-$q$ type SYK
models with the Hamiltonian given by

${\cal H}=\sum_{\kappa>0}J_{\kappa}{\cal H}_{\kappa q}$ (81)

or does it extend beyond this, or is it just a coincidence that we found for
our coupled case? Given that the kinetic SYK model is integrable, one would
imagine this to change upon adding a quadratic term to the above Hamiltonian.
This term would dominate in the low-temperature regime, which would at least
mean that the model is unlikely to have a maximally chaotic to maximally
chaotic second order phase transition, as observed in this work. But without
any kinetic term, the future work can be to find if at large-$q$ and low
temperature limits, all coupled SYK systems as in eq. (81) belong to the same
universality class as Landau-Ginzburg and if they are also maximally chaotic.

Looking ahead, our study naturally leads to the exploration of holographic
mappings of these coupled systems to $1+1$-dimensional gravity, analogous to
what has been achieved for a single SYK model [3]. As we have found that the
scaling relations near the critical point are the same as a single SYK model,
essentially what we have to do is to match our grand potential of the coupled
SYK system to that of a single SYK model. Then we use the dictionary developed
in [3] to the deformed JT gravity to find the necessary modifications on the
gravity side. We expect to find a phase transition on the gravity side between
a small and a large black hole (for example, see [7]) that corresponds to the
chaotic-to-chaotic continuous phase transition as obtained for our model. We
leave this for the future work.

## Acknowledgments

JCL and RJ would like to thank Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation))-217133147/SFB 1073 (Project No. B03) for supporting this
work. LMvM gratefully acknowledges the support of the Volkswagen foundation,
and by the Deutsche Forschungsgemeinschaft (DFG) under Grant No 406116891
within the Research Training Group RTG 2522/1.

## References

  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D 94, 106002 (2016).
  * [2] A. Kitaev, A simple model of quantum holography, Talks given at “KITP: Entanglement in Strongly-Correlated Quantum Matter”,(2015) (Part 1, Part 2).
  * Louw  _et al._ [2023] J. C. Louw, S. Cao, and X.-H. Ge, Matching partition functions of deformed Jackiw-Teitelboim gravity and the complex SYK model, Phys. Rev. D 108, 10.1103/PhysRevD.108.086014 (2023).
  * Louw and Kehrein [2023] J. C. Louw and S. Kehrein, Shared universality of charged black holes and the complex large-$q$ Sachdev-Ye-Kitaev model, Phys. Rev. B 107, 075132 (2023).
  * Azeyanagi  _et al._ [2018] T. Azeyanagi, F. Ferrari, and F. I. Schaposnik Massolo, Phase Diagram of Planar Matrix Quantum Mechanics, Tensor, and Sachdev-Ye-Kitaev Models, Phys. Rev. Lett. 120, 061602 (2018).
  * Ferrari and Schaposnik Massolo [2019] F. Ferrari and F. I. Schaposnik Massolo, Phases of melonic quantum mechanics, Phys. Rev. D 100, 026007 (2019).
  * Kubizňák and Mann [2012] D. Kubizňák and R. B. Mann, P - V criticality of charged AdS black holes, J. High Energy Phys. 2012 (7), 1.
  * Hawking and Page [1982] S. W. Hawking and D. N. Page, Thermodynamics of black holes in anti-de Sitter space, Commun. Math. Phys. 87, 577 (1982).
  * Cotler  _et al._ [2017] 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, J. High Energy Phys. 2017 (5), 1.
  * 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).
  * Patel  _et al._ [2018] A. A. Patel, M. J. Lawler, and E.-A. Kim, Coherent Superconductivity with a Large Gap Ratio from Incoherent Metals, Phys. Rev. Lett. 121, 187001 (2018).
  * Cha  _et al._ [2020] P. Cha, A. A. Patel, E. Gull, and E.-A. Kim, Slope invariant $T$-linear resistivity from local self-energy, Phys. Rev. Res. 2, 033434 (2020).
  * Louw [2023] J. Louw,  _Analytic studies on the SYK models: from instantaneous thermalization to holographic duality_ , Ph.D. thesis, Goettingen University (2023).
  * Note [1] The equivalence can be seen by considering the Lagrangian in Section II.B in [15], for a uniformly coupled lattice, i.e., site independent coupling $D_{i},J_{i}\to D,J$. The remaining terms in the Lagrangian are the Green’s functions associated with the $i^{\text{th}}$ dot ${\cal G}_{i}(t_{1},t_{2})$. Since the lattice is at equilibrium, there is no heat or charge flow, meaning that each dot has equal temperature and charge. This leads to the Green’s functions also being independent of $i$. Although the results in [15] are valid only when the $r$-body hopping on nearest lattice sites is of the order $r={\cal O}(q^{0})$, the formalism and methodology developed there is valid for any value of $r$, including $r=q/2$.
  * Jha and Louw [2023] R. Jha and J. C. Louw, Dynamics and charge fluctuations in large-$q$ Sachdev-Ye-Kitaev lattices, Phys. Rev. B 107, 235114 (2023).
  * Kamenev and Levchenko [2009] A. Kamenev and A. Levchenko, Keldysh technique and non-linear $\sigma$-model: basic principles and applications, Adv. Phys. (2009).
  * Louw and Kehrein [2022] J. C. Louw and S. Kehrein, Thermalization of many many-body interacting Sachdev-Ye-Kitaev models, Phys. Rev. B 105, 075117 (2022).
  * Stefanucci and van Leeuwen [2013] G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, Cambridge, England, UK, 2013).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Standford, A bound on chaos, J. High Energy Phys. (106).
  * Bhattacharya  _et al._ [2017] R. Bhattacharya, S. Chakrabarti, D. P. Jatkar, and A. Kundu, SYK model, chaos and conserved charge, J. High Energy Phys. 2017 (11), 1.
  * Klebanov and Tarnopolsky [2017] I. R. Klebanov and G. Tarnopolsky, Uncolored random tensors, melon diagrams, and the Sachdev-Ye-Kitaev models, Phys. Rev. D 95, 046004 (2017).
  * Jishi [2013] R. A. Jishi,  _Feynman Diagram Techniques in Condensed Matter Physics_ (Cambridge University Press, Cambridge, England, UK, 2013).
  * Murugan  _et al._ [2017] J. Murugan, D. Stanford, and E. Witten, More on supersymmetric and 2d analogs of the SYK model, J. High Energy Phys. 2017 (8), 1.
  * Pöschl and Teller [1933] G. Pöschl and E. Teller, Bemerkungen zur Quantenmechanik des anharmonischen Oszillators, Z. Phys. 83, 10.1007/BF01331132 (1933).
  * Majhi and Samanta [2017] B. R. Majhi and S. Samanta, P-V criticality of AdS black holes in a general framework, Phys. Lett. B 773, 203 (2017).
  * Standford [2016] D. Standford, Many-body chaos at weak coupling, J. High Energy Phys. (9).

##  Appendix A Derivation of differential equation for Green’s function

For our model in eq. (1), the Dyson equation is ${\cal G}^{-1}={\cal
G}_{0}^{-1}-\Sigma$, where ${\cal G}_{0}$ is the free Green’s functions.
Further the self-energy $\Sigma$ here is just the sum of the individual self-
energies $\Sigma=\Sigma_{q}+\Sigma_{q/2}$ where $\Sigma_{\kappa q}$
($\kappa=\\{1/2,1\\}$) is defined below [1].

We start with the definition of effective coupling coefficients in eq. (8).
This allows us to define another quantity

${\cal L}^{>}(t_{1},t_{2})\equiv\sum_{\kappa\in\\{\frac{1}{2},1\\}}2{\cal
J}_{\kappa}^{2}e^{\kappa g_{+}(t_{1},t_{2})},{\cal L}^{<}(t_{1},t_{2})={\cal
L}^{>}(t_{1},t_{2})^{*}$ (82)

where we reproduce the definition of the effective couplings for the
Hamiltonian in eq. (1) as

${\cal J}_{\kappa q}\equiv(1-4{\cal Q}^{2})^{\kappa q/4-1/2}J_{\kappa q}$ (83)

This contributes in the self-energy [17] of both terms in the Hamiltonian as

$\Sigma_{\kappa q}^{\gtrless}(t_{1},t_{2})=\frac{1}{q}{\cal
L}^{\gtrless}(t_{1},t_{2}){\cal G}^{\gtrless}(t_{1},t_{2})$ (84)

We are interested in obtaining the differential equations for $g^{\pm}$ whose
solutions will solve our model. We start with the differential equations for
$g^{\gtrless}$ (plugging eq. (20) into eq. (19) in [17]) where the Kadanoff-
Baym equations reduce to the following:

$\partial_{t_{1}}g^{\gtrless}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}dt_{3}{\cal
L}^{\gtrless}(t_{1},t_{3})+\underbrace{2{\cal G}^{\lessgtr}(0,0)}_{2{\cal
Q}\pm 1}{\imath}\alpha(t_{1}).$ (85)

Here $\alpha$ is related to the expectation value of energy per particle given
by

$(1-{\cal Q}^{2})\alpha(t)=\epsilon_{q/2}(t)/2+\epsilon_{q}(t)$ (86)

where in the large-$q$ limit, we have ${\cal Q}={\cal O}(q^{-1/2})$ and
$\epsilon_{\kappa q}(t)$ is defined as

$\epsilon_{\kappa q}(t)\equiv q^{2}\frac{J_{\kappa q}\langle{\cal H}_{\kappa
q}\rangle}{N}$ (87)

where $\kappa\in\\{\frac{1}{2},1\\}$.

We take the differential equation for $g^{<}(t_{1},t_{2})$ and take its
complex conjugate to get the differential equation for time reversed
$g^{<}(t_{2},t_{1})$ given by (recall ${\cal L}^{<}(t_{1},t_{2})={\cal
L}^{>}(t_{1},t_{2})^{*}$)

${\partial}_{t_{1}}g^{<}(t_{2},t_{1})=\int_{\mbox{\scriptsize$t_{1}$}}^{\mbox{\scriptsize$t_{2}$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t_{3}$\\!
}{\cal L}^{>}(t_{1},t_{3})-\left(2{\cal Q}-1\right){\imath}\alpha(t_{1})$ (88)

We add and subtract this to the differential equation for $g^{>}(t_{1},t_{2})$
to get

$\displaystyle{\partial}_{t_{1}}g^{+}(t_{1},t_{2})$
$\displaystyle=\int_{t_{1}}^{t_{2}}dt_{3}{\cal
L}^{>}(t_{1},t_{3})+{\imath}\alpha(t_{1})$ (89)
$\displaystyle{\partial}_{t_{1}}g^{-}(t_{1},t_{2})$ $\displaystyle=2{\cal
Q}{\imath}\alpha(t_{1})$

We recall that we are only interested in the equilibrium situation, implying
that all expectation values are constant, hence $\alpha(t)$ is constant as
well as the Green’s functions are only dependent on time differences
$g^{\pm}(t_{1},t_{2})\equiv g^{\pm}(t_{1}-t_{2})$. Therefore, we solve the
above differential equation for $g^{-}$ to get

$g^{-}(t_{1},t_{2})=2{\cal Q}{\imath}\alpha\cdot(t_{1}-t_{2})$ (90)

and obtain a second order differential equation for $g^{+}$ given by

${\partial}_{t_{1}}{\partial}_{t_{2}}g^{+}(t_{1}-t_{2})={\cal
L}^{>}(t_{1},t_{2})$ (91)

We finally plug the definition of ${\cal L}^{>}$ from eq. (82) above to get
the differential equation for $g^{+}$ (where $t_{1}-t_{2}\equiv t$)

$\ddot{g}^{+}(t)=-2{\cal J}_{q}^{2}e^{g_{+}(t)}-2{\cal
J}_{q/2}^{2}e^{g_{+}(t)/2}$ (92)

Thus we have solved for $g^{-}$ (eq. (90)) and obtained the differential
equation for $g^{+}$ whose solution is presented in the main text in eq. (12),
thereby solving the system.

##  Appendix B Energy contributions

The energy contributions in eq. (86) can also be derived by the integral
expression for $\alpha(t)$ [15]. Asymptotically in the large-$q$ limit using
eq. (5) and the scaling of ${\cal Q}$ mentioned in eq. (19), we get to leading
order in $1/q$ (see eqs. (19) to (20) in [17])

$\alpha(t_{1})={\imath}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$t_{1}$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t_{3}$\\!
}\left(q\Sigma^{>}(t_{1}-t_{3})-q\Sigma^{<}(t_{1}-t_{3})\right)$ (93)

where using eq. (84), we have

$q\Sigma_{\kappa q}^{>}(t_{1},t_{2})\sim-{\cal L}^{>}(t_{1},t_{2})/2$ (94)

and ${\cal L}^{<}(t_{1},t_{2})={\cal L}^{>}(t_{1},t_{2})^{*}$, leaving

$\alpha(0)={\operatorname{Im}}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$0$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t$\\!
}{\cal L}^{>}(-t)$ (95)

Note from eq. (86) that ${\cal Q}^{2}\sim{\cal O}(1/q)$ and therefore to
leading order in $1/q$, we have $\alpha$ equal to the weighted energy sum
(which is also a constant as we are considering equilibrium situation) given
by

$\alpha(0)=\sum_{\kappa=\frac{1}{2},1}\kappa\epsilon_{\kappa q}(0)$ (96)

Therefore, we can isolate the individual energy integrals as

$\kappa\epsilon_{\kappa
q}(0)={\operatorname{Im}}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$0$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t$\\!
}2{\cal J}^{2}_{\kappa q}e^{\kappa g^{+}(-t)}$ (97)

We first evaluate $\epsilon_{q/2}$ where we substitute eq. (12) for
$e^{g^{+}(-t)}$ to get

$\displaystyle\beta\epsilon_{q/2}(t_{1})$
$\displaystyle={\operatorname{Im}}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$0$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}t$\\!
}4\beta^{2}T{\cal J}_{q/2}^{2}e^{g^{+}(-t)/2}$ (98)
$\displaystyle={\operatorname{Im}}\int_{\mbox{\scriptsize$-\infty$}}^{\mbox{\scriptsize$0$}}\\!\\!\\!\mbox{\small$\,\mathrm{d}x$\\!
}\frac{4\pi v}{1+B\cos(\pi v/2+{\imath}x))}$

where we have let $x=\pi vTt$ and $B=\sqrt{A^{2}+1}$. Recall that $A$ is given
by eq. (23)

$A=\frac{\pi vU}{qK^{2}}.$ (99)

The integral is solved to get

$-\frac{8{\imath}\pi
v}{A}\tanh^{-1}\left[\left(\sqrt{A^{-2}+1}-A^{-1}\right)\tan(\pi
v/4+{\imath}x/2)\right]$ (100)

where $x\rightarrow-\infty$ is a purely real contribution thereby dropping out
(since we are interested in the imaginary part of the integral only) and at
$x=0$, we are left with

$-\frac{8{\imath}\pi
v}{A}\tanh^{-1}\left[\left(\sqrt{A^{-2}+1}-A^{-1}\right)\tan(\pi v/4)\right]$
(101)

Therefore we have the final expression as

$\beta\epsilon_{q/2}(t_{1})=-\frac{8\pi
v}{A}\tanh^{-1}\left[\left(\sqrt{A^{-2}+1}-A^{-1}\right)\tan(\pi v/4)\right]$
(102)

In the limiting case of ${\cal J}_{q}\rightarrow 0$, we are left with
$A\rightarrow 0$, thereby leading to

$\beta\epsilon_{q/2}(t_{1})=-4\pi v\tan(\pi v/4)$ (103)

which is the correct energy density for a single large-$q/2$ complex SYK
model. Furthermore, considering the (scaled) low-temperature case

$T={\tilde{T}}/q={\cal O}(q^{-1}),\quad\beta=q\tilde{\beta},$ (104)

we find $A={\cal O}(q^{-1})$ implying that the obtained result is still the
correct leading order result.

Next we proceed to calculate $\epsilon_{q}$. We know that to the leading order
in $1/q$, we have $\alpha=\epsilon_{q}+\frac{1}{2}\epsilon_{q/2}$. Here we
already know $\alpha$ from eq. (27) which is at low-temperature and
$\epsilon_{q/2}$ from eq. (102). Therefore we get

$\displaystyle E$ $\displaystyle=\langle{\cal
H}/N\rangle=\epsilon_{q/2}/q^{2}+\epsilon_{q}/q^{2}$ (105)
$\displaystyle=\epsilon_{q/2}/q^{2}+(\alpha-\epsilon_{q/2}/2)/q^{2}$

which can be re-arranged as follows:

$q\beta E=\beta\epsilon_{q/2}/(2q)+\beta\alpha/q$ (106)

where using the scaling for $\beta$ in eq. (104), we immediately see that
$\alpha$ and $q\beta E$ are both at order ${\cal O}(q^{0})$. Then we
substitute eqs. (27) and (102) to finally get the full expression for the
energy in the (scaled) low-temperature limit as

$q\beta
E=-\frac{4K^{2}}{U}\coth^{-1}\left[\sqrt{2\frac{K^{2}}{U^{2}}+1}\right]-2U\sqrt{2\frac{K^{2}}{U^{2}}+1}$
(107)

##  Appendix C Chaos theory of the SYK model

In chaos theory, we wish to analyze the evolution of the systems operator
$W(t)$ after a small perturbation $V(0)$ in the initial conditions. The
correlation function describing this evolution is the out-of-time-order
correlator (OTOC), generally written as

$c(t)=\expectationvalue{[W(t),V(0)]^{2}}_{\beta},$ (108)

where $\expectationvalue{\dots}_{\beta}$ is the thermal expectation value,
which can be written as $\expectationvalue{{\cal
O}}_{\beta}=\Tr[\rho(\beta){\cal O}]=\frac{1}{Z}\Tr[e^{\beta H}{\cal O}]$. It
has been noticed that the correlator grows exponentially with some growth rate
$\lambda_{L}$ [2],

$c(t)\sim e^{\lambda_{L}t},$ (109)

and the goal is to find this grow rate. The exponential growth of the OTOC
depends on real time, and thus we wish to compute the correlator in real time.

A common and convenient _choice_ of writing the OTOC in SYK models is [19]

$\mathcal{F}(t_{1},t_{2})=\Tr\left[y\;W(t)\;y\;V(0)\;y\;W(t)\;y\;V(0)\right],\quad
y\equiv\rho(\beta)^{1/4}$ (110)

where $\rho(\beta)$ is the thermal density matrix. This choice is to make
computation easier, and will not affect the exponent [19, 26]. For simplicity,
we have illustrated the OTOC with a Keldysh-Schwinger contour in figure (5).

$\beta$$W(t)$$W(t)$$V(0)$$V(0)$Re(t)$\beta$$W(t)$$W(t)$$V(0)$$V(0)$ Figure 5:
On the left is the Keldysh-Schwinger contour for the out-of-time-correlator
(OTOC) for the chaotic system. The perturbation $V(0)$ and generic operator
$W(t)$ are separated by a large real time and a quarter on the thermal circle,
as is shown in the picture on the right.

Here, the operators are separated by a quarter on the thermal circle (for
details, see [1]), and additionally a real-time separation between the pair of
operators. We are interested in regimes where $t\rightarrow\infty$, such that
the zero rung diagrams, can be neglected. The paths in real time are described
by the retarded propagator defined as (57), and the Wightman Greens function
(57) describes a pair of operators separated by $\beta/2$ on the thermal
circle. If we evaluate the Greens function (12) with $t\rightarrow
it+\beta/2$, we see that the function becomes a function of real time. Hence,
this operation can be seen as an analytic continuation to real time.


* [33] M. V. Lawson.  Inverse semigroups.  World Scientific Publishing Co., Inc., River Edge, NJ, 1998.  The theory of partial symmetries.
  * [34] G. Lehrer and R. Zhang.  The second fundamental theorem of invariant theory for the orthogonal group.  Ann. of Math. (2), 176(3):2031–2054, 2012.
  * [35] S. Lipscomb.  Symmetric inverse semigroups, volume 46 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 1996.
  * [36] V. Maltcev and V. Mazorchuk.  Presentation of the singular part of the Brauer monoid.  Math. Bohem., 132(3):297–323, 2007.
  * [37] P. Martin.  Temperley-Lieb algebras for nonplanar statistical mechanics—the partition algebra construction.  J. Knot Theory Ramifications, 3(1):51–82, 1994.
  * [38] P. Martin.  On diagram categories, representation theory and statistical mechanics.  In Noncommutative rings, group rings, diagram algebras and their applications, volume 456 of Contemp. Math., pages 99–136. Amer. Math. Soc., Providence, RI, 2008.
  * [39] P. Martin and V. Mazorchuk.  Partitioned binary relations.  Math. Scand., 113(1):30–52, 2013.
  * [40] V. Mazorchuk and S. Srivastava.  Jucys–Murphy elements and Grothendieck groups for generalized rook monoids.  J. Comb. Algebra, 6(1):185–222, 2022.
  * [41] E. H. Moore.  Concerning the abstract groups of order $k!$ and $\frac{1}{2}k!$ holohedrically isomorphic with the symmetric and the alternating substitution-groups on $k$ letters.  Proc. Lond. Math. Soc., 28:357–366, 1896/97.
  * [42] W. D. Munn.  Matrix representations of inverse semigroups.  Proc. London Math. Soc. (3), 14:165–181, 1964.
  * [43] R. J. Plemmons and M. T. West.  On the semigroup of binary relations.  Pacific J. Math., 35:743–753, 1970.
  * [44] L. M. Popova.  Defining relations in some semigroups of partial transformations of a finite set (in Russian).  Uchenye Zap. Leningrad Gos. Ped. Inst., 218:191–212, 1961.
  * [45] L. M. Popova.  Defining relations of a semigroup of partial endomorphisms of a finite linearly ordered set (in Russian).  Leningrad. Gos. Ped. Inst. Učen. Zap., 238:78–88, 1962.
  * [46] E. F. Robertson, N. Ruškuc, and M. R. Thomson.  On finite generation and other finiteness conditions for wreath products of semigroups.  Comm. Algebra, 30(8):3851–3873, 2002.
  * [47] E. F. Robertson, N. Ruškuc, and M. R. Thomson.  Finite generation and presentability of wreath products of monoids.  J. Algebra, 266(2):382–392, 2003.
  * [48] E. F. Robertson, N. Ruškuc, and J. Wiegold.  Generators and relations of direct products of semigroups.  Trans. Amer. Math. Soc., 350(7):2665–2685, 1998.
  * [49] L. Solomon.  Representations of the rook monoid.  J. Algebra, 256(2):309–342, 2002.
  * [50] B. Steinberg.  Möbius functions and semigroup representation theory. II. Character formulas and multiplicities.  Adv. Math., 217(4):1521–1557, 2008.
  * [51] B. Steinberg.  Representation theory of finite monoids.  Universitext. Springer, Cham, 2016.
  * [52] H. N. V. Temperley and E. H. Lieb.  Relations between the “percolation” and “colouring” problem and other graph-theoretical problems associated with regular planar lattices: some exact results for the “percolation” problem.  Proc. Roy. Soc. London Ser. A, 322(1549):251–280, 1971.