# Near-Equilibrium Approach to Transport in Complex Sachdev-Ye-Kitaev Models

Cristian Zanoci<EMAIL_ADDRESS>Department of Physics, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts
02139, USA  Brian Swingle<EMAIL_ADDRESS>Department of Physics,
Brandeis University, Waltham, Massachusetts 02453, USA

###### Abstract

We study the non-equilibrium dynamics of a one-dimensional complex Sachdev-Ye-
Kitaev chain by directly solving for the steady state Green’s functions in
terms of small perturbations around their equilibrium values. The model
exhibits strange metal behavior without quasiparticles and features diffusive
propagation of both energy and charge. We explore the thermoelectric transport
properties of this system by imposing uniform temperature and chemical
potential gradients. We then expand the conserved charges and their associated
currents to leading order in these gradients, which we can compute numerically
and analytically for different parameter regimes. This allows us to extract
the full temperature and chemical potential dependence of the transport
coefficients. In particular, we uncover that the diffusivity matrix takes on a
simple form in various limits and leads to simplified Einstein relations. At
low temperatures, we also recover a previously known result for the Wiedemann-
Franz ratio. Furthermore, we establish a relationship between diffusion and
quantum chaos by showing that the diffusivity eigenvalues are upper bounded by
the chaos propagation rate at all temperatures. Our work showcases an
important example of an analytically tractable calculation of transport
properties in a strongly interacting quantum system and reveals a more general
purpose method for addressing strongly coupled transport.

##  I Introduction

The study of quantum systems out of equilibrium can shed light on many
outstanding questions regarding thermalization, transport, and quantum many-
body chaos in condensed matter theories. In addition to these conceptual
problems, there are many practical questions motivated by recent experimental
advances in ultracold atoms and solid state systems, which present new avenues
for investigating the non-equilibrium dynamics of many-body systems. In
particular, quantum transport has garnered a lot of attention recently in an
attempt to uncover new features of the dynamical processes governing the
behavior of strongly interacting systems out of equilibrium. Despite numerous
efforts, practical calculations of transport coefficients in quantum many-body
systems remain challenging from both a theoretical and technical standpoint
[1], especially at low temperatures.

One-dimensional models have emerged as prototypical examples for studying
transport phenomena, partly due to their computational tractability. They
usually consist of interacting particles or spins on a lattice that are driven
away from equilibrium by certain external biases. The system then relaxes to a
steady state dictated by its microscopic dynamics, which typically involves
the transport of conserved quantities according to local conservation laws [2,
3]. These conserved charges and their associated currents are the quantities
of interest.

A common implementation of this idea involves connecting the system to baths
that drive it towards a desired steady state, where many transport properties
are easily available [1, 4, 5]. However, reaching this non-equilibrium steady
state (NESS) in the hydrodynamic limit can be practically challenging [6],
since most numerical techniques are usually limited to small systems and short
evolution times. If we could bypass simulating the open-system non-equilibrium
dynamics entirely and instead access the emergent NESS directly, we would be
able to immediately find all the transport properties of the system.

For a general class of models, we have previously shown that the local Green’s
functions in NESS are only slightly perturbed from their equilibrium values in
the case of weak driving [7]. This allowed us to find these non-equilibrium
corrections explicitly in terms of the equilibrium Green’s functions, without
having to solve for the open-system dynamics. Our method is equivalent to a
first order expansion in the local gradients, and thus falls under the
umbrella of linear response theory. In this approximation, the temperature and
chemical potential differences across the system are assumed to be small
compared to their average values. Conveniently, most experimental setups
studying transport in many-body systems also operate in the linear-response
regime.

The class of models in question consists of lattices built from the Sachdev-
Ye-Kitaev (SYK) model [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
describing fermions with random all-to-all $q$-body interactions. In this
paper, we will focus specifically on the complex fermion version of SYK [10,
11, 8, 9, 12, 20, 21, 22, 18, 23, 24, 19], which has an additional conserved
global $U(1)$ charge. This model displays a multitude of remarkable
properties, ranging from an emergent approximate conformal symmetry at low
temperatures [13, 12, 14] to maximal many-body chaos [25]. In fact, the SYK
model is holographically dual to extremal charged black holes with AdS2
horizons [26, 27, 12, 13, 14, 28, 29, 15, 30], and has a residual entropy
directly connected to the Bekenstein-Hawking entropy of these black holes [26,
27, 12]. The model and its many variations [31, 32, 33, 34, 35, 36, 37] belong
to a class of systems realizing holographic quantum matter without
quasiparticles, an thus represent a valuable platform for studying non-Fermi
liquid behavior [10, 11, 12, 19]. Given the interesting physical properties of
the SYK family of models, several experimental implementations have been
recently proposed [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].

The non-equilibrium dynamics of SYK models has been previously studied through
various quench protocols [50, 51, 52, 53, 54] or through couplings to external
baths [55, 56, 7, 57, 58, 59, 60] and Lindblad operators [61, 62]. In
particular, several questions pertaining to transport and chaos in higher-
dimensional lattices of coupled SYK clusters have been addressed [31, 21, 33,
32, 34, 35, 37, 36, 7]. These include many indicative properties of strange
metals, such as diffusive propagation of energy [31, 21, 33] and resistivity
that scales linearly with temperature [33, 36]. Moreover, it was shown that
the same time-reparametrization field is responsible for the propagation of
both low-energy modes and quantum chaos [31, 21], thus leading to a connection
between energy diffusion and the butterfly velocity [63, 64, 65, 66, 67, 68,
69, 70, 71, 72, 73].

Nonetheless, the problem of characterizing transport for arbitrary model
parameters remains mostly unsolved. Many of the previous approaches relied on
the large $q$ limit or the low-temperature Schwarzian effective action to
describe the energy and charge fluctuations [31, 68, 32, 70]. These methods
have a limited range of applicability and often do not lead to explicit
solutions for the transport coefficients. In this paper, we propose a more
general approach based on the expansion of the SYK Green’s functions in the
near-equilibrium regime, in the presence of constant temperature and chemical
potential gradients. This allows us to compute the charge and energy currents
throughout an SYK chain, and hence determine the associated diffusivities and
conductivities. This method has the immediate advantage of delivering
numerical results for any values of temperature, chemical potential, and $q$.
Additionally, we obtain closed-form solutions in the limits of large $q$ and
$q=2$.

This work represents a natural extension of our previous analysis of energy
diffusion in Majorana SYK models [7]. Since complex fermions feature both
charge and energy conservation, we were able to fully characterize the
combined thermoelectric response of a strange metal and expose some of its
most fascinating aspects. First, we studied the interplay of transport with
various thermodynamic quantities and the physics of phase transitions. Second,
we showed that the diffusivity matrix takes on a particular form at both high
and low temperatures, as well as in the large $q$ limit. Third, we verified
that the Wiedemann-Franz ratio approaches a known constant at zero temperature
[21]. Last, but not least, we related the eigenvalues of the diffusivity
matrix to an upper bound set by chaos $D_{\pm}\leq v_{B}^{2}/\lambda_{L}$ at
all temperatures and saturated in the conformal limit [31, 21, 70, 7]. Our
results provide concrete values for the transport coefficients that can be
measured in the aforementioned experiments. But most importantly, they suggest
a promising path towards more general studies of out-of-equilibrium phenomena
in strongly interacting systems in which one directly accesses non-equilibrium
steady states of interest.

The rest of the paper is structured as follows. In Sec. II we introduce our
one-dimensional SYK model. Subsequently, in Sec. III we describe in detail our
approach to studying the equilibrium, non-equilibrium, and chaotic properties
of this model. In Sec. IV we review the phase diagram of the complex SYK model
and present our main results for the transport coefficients as a function of
temperature and chemical potential. We also discuss a chaos bound on
diffusivities in that section. We then provide a brief discussion of our
findings and comment on possible extensions in Sec. V. The details of our
calculations are available in the Appendix.

##  II Model

The building block of our model is a complex SYK cluster [10, 11, 8, 9, 12,
20, 21, 22, 18, 23, 24, 19] with random all-to-all $q-$body interactions among
$N$ fermions in $(0+1)$ dimensions. In order to study transport in this model,
we generalize it to an infinitely-long one-dimensional chain (see Fig. 1),
where each site $x$ is an SYK cluster characterized by the Hamiltonian

$H_{0}^{x}=\sum_{\\{i\\},\\{j\\}}J_{i_{1}\ldots i_{\frac{q}{2}}j_{1}\ldots
j_{\frac{q}{2}}}^{(0)}(c_{i_{1}}^{x})^{\dagger}\cdots(c_{i_{\frac{q}{2}}}^{x})^{\dagger}c_{j_{1}}^{x}\cdots
c_{j_{\frac{q}{2}}}^{x},$ (1)

where $q$ is an even integer and $\\{i\\}$ denotes the restricted sum over
$1\leq i_{1}<\cdots<i_{\frac{q}{2}}\leq N$. The neighboring sites of the chain
interact via a similar Hamiltonian

$H_{1}^{x,x+1}=\sum_{\\{i\\},\\{j\\}}J_{i_{1}\ldots
j_{\frac{q}{2}}}^{(1)}(c_{i_{1}}^{x})^{\dagger}\cdots(c_{i_{\frac{q}{2}}}^{x})^{\dagger}c_{j_{1}}^{x+1}\cdots
c_{j_{\frac{q}{2}}}^{x+1}+\mathrm{h.c.}.$ (2)

The fermions obey the standard anti-commutation relations
$\\{(c_{i}^{x})^{\dagger},c_{j}^{x^{\prime}}\\}=\delta_{ij}\delta_{xx^{\prime}}$.
Note that we choose the interaction term to consist of the same number of
fermion operators from each site and one can consider adding more general
interactions that may change the transport properties of the model [21]. The
SYK couplings are complex, independent Gaussian random variables with zero
mean obeying

$\displaystyle J_{i_{1}\ldots i_{\frac{q}{2}}j_{1}\ldots
j_{\frac{q}{2}}}^{(0,1)}$ $\displaystyle=(J_{j_{1}\ldots
j_{\frac{q}{2}}i_{1}\ldots i_{\frac{q}{2}}}^{(0,1)})^{*},$ (3)
$\displaystyle\langle(J_{i_{1}\ldots i_{\frac{q}{2}}j_{1}\ldots
j_{\frac{q}{2}}}^{(0,1)})^{2}\rangle$
$\displaystyle=\frac{J_{0,1}^{2}(q/2)!(q/2-1)!}{N^{q-1}}.$ (4)

The numerical coefficients are chosen to cancel additional factors in the path
integral and the powers of $N$ ensure the correct scaling of extensive
thermodynamic variables.

For future convenience, we write the total system Hamiltonian in terms of bond
operators $H^{x,x+1}$ acting on consecutive sites $(x,x+1)$

$H=\sum_{x}H^{x,x+1}=\sum_{x}\left(\frac{1}{2}H_{0}^{x}+\frac{1}{2}H_{0}^{x+1}+H_{1}^{x,x+1}\right).$
(5)

This Hamiltonian is invariant under the particle-hole symmetry and has a
globally conserved $U(1)$ charge density $Q=\sum_{x}Q_{x}$, where the local
charge density $Q_{x}\in(-1/2,1/2)$ is defined as in Ref. [21]

$Q_{x}=\frac{1}{N}\sum_{i}\langle(c_{i}^{x})^{\dagger}c_{i}^{x}\rangle-\frac{1}{2}.$
(6)

The only other conserved quantity is the energy and it is precisely the
interplay between the transport properties of these two conserved charges that
we aim to study.

We should mention that similar higher-dimensional SYK models have been
previously studied in the context of transport [31, 68, 21, 32, 33, 34, 35,
37, 36, 7], quantum chaos [74, 59, 75, 76, 69], and quantum phase transitions
[77, 78, 32, 79, 80]. Most importantly, a generalization of our setup to
arbitrary graphs coupled to thermal reservoirs should be straightforward [7].

##  III Methods

In this section, we derive the equations governing the equilibrium and non-
equilibrium dynamics of our model, review the definitions of various
thermodynamic, transport, and chaos-related quantities that we report later in
the paper, and show how these observables simplify in the limit of small and
large $q$.

###  III.1 Equilibrium

We begin with the equilibrium description of our model, which is most easily
done in imaginary time. The SYK chain maintains all the exactly solvable
properties of a single SYK cluster in the large-$N$ limit [12, 21, 33, 18]. We
are interested in the grand canonical partition function $Z=\Tr
e^{-\beta(H-\mu Q)}$, where $\beta=1/T$ is the inverse temperature and $\mu$
is the chemical potential that fixes the value of $Q$. In equilibrium, these
parameters are constant (site-independent) throughout the system. Due to the
self-averaging property of this model at large $N$ [12], it is sufficient to
consider the replica-diagonal partition function
$Z=\int[\mathop{}\\!\mathrm{d}G_{x}][\mathop{}\\!\mathrm{d}\Sigma_{x}]e^{-S}$,
for which the Euclidean effective action, after integrating out the fermions,
becomes

$\displaystyle S=\sum_{x}(S_{x,x+1}+S_{x}),$ (7) $\displaystyle
S_{x,x+1}=-\frac{2NJ_{1}^{2}}{q}\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}(-G_{x}(\tau_{1},\tau_{2})G_{x+1}(\tau_{2},\tau_{1}))^{\frac{q}{2}},$
(8)
$S_{x}=-N\log\det\big((\partial_{\tau}-\mu)\delta(\tau_{1}-\tau_{2})+\Sigma_{x}(\tau_{1},\tau_{2})\big{missing})-N\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}\Big{(}\Sigma_{x}(\tau_{1},\tau_{2})G_{x}(\tau_{2},\tau_{1})+\frac{J_{0}^{2}}{q}(-G_{x}(\tau_{1},\tau_{2})G_{x}(\tau_{2},\tau_{1}))^{\frac{q}{2}}\Big{)}.$
(9)

For each site $x$, we defined the Euclidean time-ordered fermion two-point
function

$G_{x}(\tau_{1},\tau_{2})=-\frac{1}{N}\sum_{i=1}^{N}\langle
Tc_{i}^{x}(\tau_{1})c_{i}^{x}(\tau_{2})^{\dagger}\rangle,$ (10)

and the fermion self-energy $\Sigma_{x}(\tau_{1},\tau_{2})$ as the associated
Lagrange multiplier. In the large-$N$ limit, the saddle point of this
effective action produces the Schwinger-Dyson (SD) equations of motion

$\displaystyle G_{x}(i\omega_{n})$
$\displaystyle=\frac{1}{i\omega_{n}+\mu-\Sigma_{x}(i\omega_{n})},$ (11)
$\displaystyle\Sigma_{x}(\tau)$
$\displaystyle=(-1)^{\frac{q}{2}-1}G_{x}(\tau)^{\frac{q}{2}}\Big{(}J_{0}^{2}G_{x}(-\tau)^{\frac{q}{2}-1}$
$\displaystyle+J_{1}^{2}G_{x-1}(-\tau)^{\frac{q}{2}-1}+J_{1}^{2}G_{x+1}(-\tau)^{\frac{q}{2}-1}\Big{)},$
(12)

where $\omega_{n}=(2n+1)\pi/\beta$ is a Matsubara frequency and we assumed
time-translation invariance $\tau=\tau_{1}-\tau_{2}$ in equilibrium. At half
filling ($\mu=0$), the Green’s function and all the other quantities derived
from it are identical to those of a Majorana SYK model. The only difference is
that complex fermions have twice as many degrees of freedom, leading to a
trivial doubling of all the extensive observables.

As we have previously shown [7], for a uniform chain in equilibrium, the
Green’s functions take on the site-independent value $G(\tau)$ and the
effective on-site coupling become $J=\sqrt{J_{0}^{2}+2J_{1}^{2}}$. In other
words, the SD equations for an interacting SYK cluster have the exact same
form as those of an isolated $(0+1)$ dimensional SYK model with coupling $J$

$\begin{split}G(i\omega_{n})&=\frac{1}{i\omega_{n}+\mu-\Sigma(i\omega_{n})},\\\
\Sigma(\tau)&=J^{2}G(\tau)^{\frac{q}{2}}(-G(-\tau))^{\frac{q}{2}-1}.\end{split}$
(13)

This system of SD equations can be solved numerically using the method
described in Appendix A.

###  III.2 Thermodynamics

Once we obtain the solutions of the equilibrium SD equations, we can compute
any thermodynamic variable. Our goal here is twofold. First, we would like to
derive expressions for the common thermodynamic functions, such as entropy and
heat capacity, that would help us identify a phase transition in the complex
SYK model and assess its impact on the transport coefficients and chaos [81,
82, 83, 53, 23, 84]. Second, we need to find the susceptibility matrix, which
relates diffusivities and conductivities [85, 86, 87, 63], as described in
more detail in Sec. III.4.

Since our system is homogeneous, we can focus on the thermodynamic properties
of a single SYK node with interaction strength $J$. In what follows, all the
extensive quantities are replaced by their densities per particle (i.e.
divided by $N$). As is the case for most thermodynamic problems, our starting
point is the grand canonical potential $\Omega$. In the large-$N$ limit,
$\Omega$ is approximated by evaluating the action on the saddle-point
equations of motion

$\begin{split}\Omega&=-T\log{Z}=T\Bigg{[}\sum_{n}\log\left(\frac{G(i\omega_{n})}{G_{0}(i\omega_{n})}\right)\\\
&-\frac{q-1}{q}\sum_{n}\Sigma(i\omega_{n})G(i\omega_{n})-\log\left(2\cosh(\mu/2T)\right)\Bigg{]},\end{split}$
(14)

where $G_{0}(i\omega_{n})=(i\omega_{n}+\mu)^{-1}$ is the free fermion Green’s
function. Here we have regularized the logarithm by adding and subtracting the
free fermion contribution [21, 33, 18], and evaluated the last term using the
Matsubara frequency summation [85, 86, 87]. The free energy is given by a
Legendre transform $F=\Omega+\mu Q$. Recall that the charge and chemical
potential can be obtain from their respective ensembles at fixed temperature

$\begin{split}Q&=-\left(\frac{\partial\Omega}{\partial\mu}\right)_{T}=\frac{1}{2}\left(G(0^{+})-G(\beta^{-})\right),\\\
\mu&=\left(\frac{\partial F}{\partial
Q}\right)_{T}=-\partial_{\tau}G(0^{+})-\partial_{\tau}G(\beta^{-}),\end{split}$
(15)

where the second equalities in terms of Green’s function are derived in Ref.
[12, 18]. Note that our definition of the grand canonical potential in Eq.
(14) gives us exactly the charge density introduced in Eq. (6).

The entropy is computed as the first derivative of the potential, using the
standard thermodynamic identities

$S=-\left(\frac{\partial\Omega}{\partial T}\right)_{\mu}=-\left(\frac{\partial
F}{\partial T}\right)_{Q}.$ (16)

A striking feature of the SYK model is its non-zero residual entropy $S_{0}$
in the limit of zero temperature, which is not due to an exponentially large
ground state degeneracy, but rather because of the exponentially small level
spacing all the way down to the ground state [11, 8, 9, 12]. We will use this
entropy to distinguish between an SYK-like ground state and a trivial one in
Sec. IV.1, and will also relate it to the thermopower in Appendix D.

The static susceptibility matrix $\chi$ relates the change in macroscopic
observables due to the variation of the underlying microscopic quantities

$\begin{pmatrix}\nabla Q\\\ \nabla E-\mu\nabla
Q\end{pmatrix}=\begin{pmatrix}\chi_{11}&\chi_{12}\\\
\chi_{21}&\chi_{22}\end{pmatrix}\begin{pmatrix}\nabla\mu\\\ \nabla
T\end{pmatrix},$ (17)

and can be cast in terms of the second derivative of the grand potential [85,
86, 87]

$\chi=\begin{pmatrix}-\left(\frac{\partial^{2}\Omega}{\partial\mu^{2}}\right)_{T}&-\left(\frac{\partial^{2}\Omega}{\partial\mu\partial
T}\right)_{\mu,T}\\\ -T\left(\frac{\partial^{2}\Omega}{\partial
T\partial\mu}\right)_{T,\mu}&-T\left(\frac{\partial^{2}\Omega}{\partial
T^{2}}\right)_{\mu}\end{pmatrix}.$ (18)

By virtue of equality of mixed partial derivatives, the off-diagonal elements
are always related by $\chi_{21}=T\chi_{12}$. The diagonal elements can be
interpreted as the charge compressibility

$K\equiv\chi_{11}=-\left(\frac{\partial^{2}\Omega}{\partial\mu^{2}}\right)_{T}=\left(\frac{\partial
Q}{\partial\mu}\right)_{T},$ (19)

and heat capacity at constant chemical potential

$C_{\mu}\equiv\chi_{22}=-T\left(\frac{\partial^{2}\Omega}{\partial
T^{2}}\right)_{\mu}=T\left(\frac{\partial S}{\partial T}\right)_{\mu}.$ (20)

The heat capacity at fixed charge can be related to the other entries in the
susceptibility matrix via the thermodynamic identity [63]

$C_{Q}=T\left(\frac{\partial S}{\partial
T}\right)_{Q}=C_{\mu}-\frac{T\chi_{12}^{2}}{K}.$ (21)

Finally, the linear-in-$T$ coefficient of the specific heat is simply defined
as $\gamma=C_{Q}/T$.

Both $K$ and $\gamma$ play an important role in transport. They appear as the
coefficients in the low-temperature Schwarzian effective action used to
describe charge and energy fluctuations [21]. We will also show that in this
conformal limit, the ratio of energy to charge diffusivities is governed by
$K/\gamma$.

###  III.3 Non-equilibrium

Although the Euclidean time formulation works well for thermodynamics, it is
not suitable for non-equilibrium dynamics, due to the problems arising from
analytic continuation to zero frequency. Therefore, the non-equilibrium
evolution of a quantum many-body system is better described in real-time using
the Schwinger-Keldysh formalism [88, 89]. Following the derivation in Refs.
[33, 58], we can write down the terms in a Lorentzian effective action after
integrating out the fermions, just as we did in imaginary time

$S_{x,x+1}=\frac{2iNJ_{1}^{2}}{q}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}G_{x}(t_{2},t_{1})^{\frac{q}{2}}G_{x+1}(t_{1},t_{2})^{\frac{q}{2}},$
(22)
$S_{x}=-iN\log\det\left(\partial_{t}\delta_{\mathcal{C}}(t_{1},t_{2})+i\Sigma_{x}(t_{1},t_{2})\right)+iN\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}\left(\frac{J_{0}^{2}}{q}G_{x}(t_{2},t_{1})^{\frac{q}{2}}G_{x}(t_{1},t_{2})^{\frac{q}{2}}-\Sigma_{x}(t_{1},t_{2})G_{x}(t_{2},t_{1})\right),$
(23)

where $\mathcal{C}$ denotes the closed-time Keldysh contour consisting of a
positive and a negative branch [88, 89]. Recall that in the Schwinger-Keldysh
formalism, the contour-ordered Green’s functions $G_{x}(t_{1},t_{2})$ are
actually $2\times 2$ matrices, where each entry corresponds to a placement of
the time arguments on either branch of the contour. The off-diagonal entries
of this matrix are the well-known greater and lesser Green’s functions

$\displaystyle G_{x}^{>}(t_{1},t_{2})\equiv G_{x}(t_{1}^{-},t_{2}^{+})$
$\displaystyle=-\frac{i}{N}\sum_{i=1}^{N}\langle
c_{i}(t_{1}^{-})c_{i}^{\dagger}(t_{2}^{+})\rangle,$ (24) $\displaystyle
G_{x}^{<}(t_{1},t_{2})\equiv G_{x}(t_{1}^{+},t_{2}^{-})$
$\displaystyle=\frac{i}{N}\sum_{i=1}^{N}\langle
c_{i}^{\dagger}(t_{2}^{-})c_{i}(t_{1}^{+})\rangle,$ (25)

where $t_{i}^{\pm}$ live on the positive and negative branch respectively. The
contour-ordered Green’s functions above are related to the more conventional
retarded, advanced, and Keldysh Green’s functions via a Keldysh rotation [88,
89]

$\displaystyle G_{x}^{R}(t_{1},t_{2})$
$\displaystyle=\Theta(t_{1}-t_{2})\big{(}G_{x}^{>}(t_{1},t_{2})-G_{x}^{<}(t_{1},t_{2})\big{)},$
(26) $\displaystyle G_{x}^{A}(t_{1},t_{2})$
$\displaystyle=\Theta(t_{2}-t_{1})\big{(}G_{x}^{<}(t_{1},t_{2})-G_{x}^{>}(t_{1},t_{2})\big{)},$
(27) $\displaystyle G_{x}^{K}(t_{1},t_{2})$
$\displaystyle=G_{x}^{>}(t_{1},t_{2})+G_{x}^{<}(t_{1},t_{2}),$ (28)

The corresponding self-energies are defined in a similar manner.

For the purposes of our analysis, we will only consider states in thermal
equilibrium or steady states weakly perturbed from equilibrium [7]. In both
cases, the fermion Green’s functions become time-translation invariant and
satisfy the identity [90]

$G_{x}^{\gtrless}(t_{1},t_{2})=G_{x}^{\gtrless}(t=t_{1}-t_{2})=-G_{x}^{\gtrless}(-t)^{*}.$
(29)

Furthermore, their values at $t=0$ are related to the local charge and
chemical potential via

$\begin{split}Q_{x}&=-iG_{x}^{>}(0)+\frac{1}{2}=-iG_{x}^{<}(0)-\frac{1}{2},\\\
\mu_{x}&=\partial_{t}G_{x}^{>}(0)-\partial_{t}G_{x}^{<}(0),\end{split}$ (30)

which follow from an analytic continuation of Eq. (15) to real time.

To obtain the Schwinger-Dyson equations governing the real-time dynamics of
the system, we can look for large-$N$ saddle point solutions of the Lorentzian
action [33, 58]

$\begin{split}G_{x}^{R}(\omega)&=\frac{1}{i\omega-\Sigma_{x}^{R}(\omega)},\\\
\Sigma_{x}^{\gtrless}(t)&=G_{x}^{\gtrless}(t)^{\frac{q}{2}}\Big{(}J_{0}^{2}G_{x}^{\lessgtr}(-t)^{\frac{q}{2}-1}\\\
&+J_{1}^{2}G_{x-1}^{\lessgtr}(-t)^{\frac{q}{2}-1}+J_{1}^{2}G_{x+1}^{\lessgtr}(-t)^{\frac{q}{2}-1}\Big{)}.\end{split}$
(31)

Note that these equations are only valid for the time-translation invariant
case. For more general non-equilibrium setups, one has to derive the full
Kadanoff-Baym equations and solve them numerically [7].

We emphasize that the real-time action does not involve a chemical potential,
since $\mu_{x}$ is a property of the state, rather than the Hamiltonian. Often
times this issue is addressed by explicitly adding a mass term
$-\mu_{x}\sum_{i=1}^{N}(c_{i}^{x})^{\dagger}c_{i}^{x}$ to the Hamiltonian [12,
33, 60, 58, 57]. This works in imaginary time, where it is simply equivalent
to working in the grand canonical ensemble. However, the chemical potential
and mass term are not equivalent in real time [83], with the Green’s functions
being typically off by a factor of $e^{i\mu_{x}t}$, which can lead to
incorrect dynamics and transport properties.

Similarly, Eq. (31) does not have an explicit dependence on chemical potential
or temperature, in contrast to its Euclidean-time counterpart in Eq. (12).
Hence there are infinitely many distinct saddle point solutions of the real-
time SD equation to which we can converge, each corresponding to different
$\mu_{x}$ and $\beta_{x}$ . To circumvent this problem, we use the
fluctuation-dissipation theorem (FDT) [85, 86, 87] to set the values of these
parameters

$\frac{iG_{x}^{K}(\omega)}{A_{x}(\omega)}=\tanh(\frac{\beta_{x}(\omega-\mu_{x})}{2}),$
(32)

where $A_{x}(\omega)=-2\imaginary G_{x}^{R}(\omega)$ is the spectral function
and $G_{x}^{K}(\omega)$ is the Fourier transform of the Keldysh Green’s
function. Fixing the local temperature and chemical potential in this way is
applicable to both the equilibrium and near-thermal steady states under
consideration [7]. The FDT, together with Eq. (31), form a closed set of
equations that can be solved iteratively to find a unique solution (see
Appendix A).

So far we have assumed that each cluster has a well defined chemical potential
$\mu_{x}$ and inverse temperature $\beta_{x}$. This is indeed true for a
uniform chain in equilibrium with $\mu_{x}=\mu$ and $\beta_{x}=\beta$.
Similarly to the imaginary-time version, the real-time site-independent
solution $G^{\gtrless}(t)$ is the same as that of a single SYK node with an
effective coupling $J=\sqrt{J_{0}^{2}+2J_{1}^{2}}$. We have also shown that in
the presence of a small bias throughout the chain, the NESS Green’s functions
at late times are only slightly perturbed from their equilibrium values [7],
as long as we are still in the linear response regime. This bias can be
introduced by either directly coupling the system to baths at different
chemical potentials and temperatures [59, 56, 55, 57, 58, 7, 60], or by
introducing an effective coupling to the environment through Lindblad
operators [61, 62]. As we will discuss in detail in the next section, to
extract the thermoelectric transport coefficients, it is enough to impose a
uniform chemical potential or temperature gradient along the chain. Thus one
can define local parameters that are ever so slightly perturbed from their
equilibrium values

$\begin{split}\mu_{x}&=\mu+x\nabla\mu,\\\
\beta_{x}&=\beta+x\nabla\beta,\end{split}$ (33)

with $|\nabla\mu|\ll\mu$ and $|\nabla\beta|\ll\beta$. This allows us to write
the near-equilibrium Green’s functions in terms of an extra site-dependent
correction

$G_{x}^{\gtrless}(t)=G^{\gtrless}(t)+xF_{\mu,\beta}^{\gtrless}(t),$ (34)

where $|F_{\mu,\beta}^{\gtrless}(t)|\ll|G^{\gtrless}(t)|$ are the non-
equilibrium contributions proportional to the gradients

$\begin{split}F_{\mu}^{\gtrless}(t)&=\frac{\mathop{}\\!\mathrm{d}G^{\gtrless}(t)}{\mathop{}\\!\mathrm{d}\mu}\nabla\mu,\\\
F_{\beta}^{\gtrless}(t)&=\frac{\mathop{}\\!\mathrm{d}G^{\gtrless}(t)}{\mathop{}\\!\mathrm{d}\beta}\nabla\beta.\end{split}$
(35)

The subscripts refer to whether the perturbation is due a chemical potential
or a temperature gradient. To first order, these contributions can be summed
to characterize the system’s response to any mixed thermoelectric bias. Eq.
(34) is analogous to a gradient expansion in hydrodynamics. We see that to
access the non-equilibrium transport physics in the linear response regime, it
is sufficient to solve the SD equations in real time at equilibrium.

Figure 1: Schematic depiction of the infinite one-dimensional SYK chain in
non-equilibrium. Each site contains $N$ fermions with intra-cluster coupling
$J_{0}$ and inter-cluster coupling $J_{1}$. The system is subject to uniform
biases $\nabla\mu$ and $\nabla\beta$. In the steady state, homogeneous charge
and energy currents $j^{Q,E}$ flow through the chain. The local Green’s
functions are perturbed from their equilibrium values $G(t)$ by the linear
response contributions $F(t)$.

###  III.4 Transport

Our model has charge and energy as the only two conserved quantities. These
are expressed in terms of the local on-site charge density $Q_{x}$ and the on-
bond energy density $E_{x,x+1}=\langle
H^{x,x+1}\rangle=(E_{0}^{x}+E_{0}^{x+1})/2+E_{1}^{x,x+1}$ introduced in Sec.
II. The charge density can be computed from the real-time Green’s functions
using Eq. (30), while the energy density has contributions from

$\displaystyle E_{0}^{x}$
$\displaystyle=-2i\frac{J_{0}^{2}}{q}\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{1}\left(G_{x}^{>}(t,t_{1})^{\frac{q}{2}}G_{x}^{<}(t_{1},t)^{\frac{q}{2}}-\mathrm{h.c.}\right),$
$\displaystyle E_{1}^{x,x+1}$
$\displaystyle=-2i\frac{J_{1}^{2}}{q}\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{1}\Big{(}G_{x}^{>}(t,t_{1})^{\frac{q}{2}}G_{x+1}^{<}(t_{1},t)^{\frac{q}{2}}$
$\displaystyle+G_{x}^{<}(t_{1},t)^{\frac{q}{2}}G_{x+1}^{>}(t,t_{1})^{\frac{q}{2}}-\mathrm{h.c.}\Big{)}.$
(36)

These conserved quantities have an associated charge current density
$j_{x}^{Q}=\mathop{}\\!\mathrm{d}Q_{x}/\mathop{}\\!\mathrm{d}t$ and energy
current density
$j_{x}^{E}=\mathop{}\\!\mathrm{d}E_{x,x+1}/\mathop{}\\!\mathrm{d}t$
respectively. The formulas for the currents flowing across a site $x$ can be
derived by combining the continuity equation and Heisenberg’s equation of
motion [2, 3], resulting in

$\displaystyle j_{x}^{Q}$ $\displaystyle=i[Q_{x},H^{x,x+1}],$ (37)
$\displaystyle j_{x}^{E}$ $\displaystyle=i[H^{x-1,x},H^{x,x+1}].$ (38)

Computing the expectation value of these commutators in the Schwinger-Keldysh
formalism is more involved and we provide a derivation in Appendix B. Our
general formulas for the currents are given by Eqs. (90, 92-94).

In equilibrium, there are no currents flowing through the system. To observe a
finite current, we have to introduce a small bias, accomplished, for instance,
by connecting the chain to reservoirs at its two ends [7]. In the long-time
limit, when the system reaches its steady state, the currents become uniform
throughout the chain $j^{Q,E}\equiv\langle j_{x}^{Q,E}\rangle$, as shown in
Fig. 1. In the linear response regime, the gradients are also small and
constant as in Eq. (33). Therefore, we can use Eq. (34) to write all the
quantities of interest in terms of the equilibrium Green’s functions and to
first order in non-equilibrium corrections $F_{\mu,\beta}^{\gtrless}(t)$. For
example, the charge gradient becomes

$\nabla Q=-iF_{\mu,\beta}^{>}(0),$ (39)

while the energy gradient is given by

$\nabla
E=2J^{2}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[\left(-G^{>}(t)G^{<}(t)^{*}\right)^{\frac{q}{2}}\left(\frac{F_{\mu,\beta}^{>}(t)}{G^{>}(t)}+\frac{F_{\mu,\beta}^{<}(t)^{*}}{G^{<}(t)^{*}}\right)\right].$
(40)

Similarly, the charge current in Eq. (90) becomes

$j^{Q}=\frac{qJ^{2}}{2}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Re}\left[\left(-G^{>}(t)G^{<}(t)^{*}\right)^{\frac{q}{2}}\left(\frac{F_{\mu,\beta}^{>}(t)}{G^{>}(t)}-\frac{F_{\mu,\beta}^{<}(t)^{*}}{G^{<}(t)^{*}}\right)\right],$
(41)

and the energy current in Eqs. (92-94) simplifies to

$j^{E}=\frac{1}{2}J_{1}^{2}J^{2}\real(j_{++}+j_{+-}),$ (42) $\displaystyle
j_{++}$
$\displaystyle=q\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(G^{>}(t)G^{<}(t)^{*}G^{>}(t^{\prime})G^{<}(t^{\prime})^{*}\right)^{\frac{q}{2}-1}$
$\displaystyle\cdot\Big{(}G^{<}(t^{\prime}-t)^{*}\left(G^{<}(t)^{*}F_{\mu,\beta}^{>}(t^{\prime})-G^{>}(t^{\prime})F_{\mu,\beta}^{<}(t)^{*}\right)-G^{>}(t^{\prime}-t)\left(G^{<}(t^{\prime})^{*}F_{\mu,\beta}^{>}(t)-G^{>}(t)F_{\mu,\beta}^{<}(t^{\prime})^{*}\right)\Big{)},$
(43) $\displaystyle j_{+-}$
$\displaystyle=-q\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(G^{>}(t)^{*}G^{<}(t)G^{>}(t^{\prime})G^{<}(t^{\prime})^{*}\right)^{\frac{q}{2}-1}$
$\displaystyle\cdot\Big{(}G^{<}(t^{\prime}-t)^{*}\left(G^{>}(t)^{*}F_{\mu,\beta}^{>}(t^{\prime})-G^{>}(t^{\prime})F_{\mu,\beta}^{>}(t)^{*}\right)-G^{>}(t^{\prime}-t)\left(G^{<}(t^{\prime})^{*}F_{\mu,\beta}^{<}(t)-G^{<}(t)F_{\mu,\beta}^{<}(t^{\prime})^{*}\right)\Big{)}.$
(44)

Within linear response, the currents are related to the conjugate gradients
via the conductivity matrix $L$

$\begin{pmatrix}j^{Q}\\\ j^{E}-\mu
j^{Q}\end{pmatrix}=-\begin{pmatrix}\sigma&\alpha\\\ \alpha
T&\bar{\kappa}\end{pmatrix}\begin{pmatrix}\nabla\mu\\\ \nabla T\end{pmatrix},$
(45)

where $\sigma\equiv L_{11}$ is the electrical conductivity, $\alpha\equiv
L_{12}$ is the thermoelectric conductivity, and $\kappa\equiv
L_{22}-L_{12}L_{21}/L_{11}=\bar{\kappa}-\alpha^{2}T/\sigma$ is the thermal
conductivity [85, 86, 87, 63]. The off-diagonal elements are constrained by
the Onsager reciprocal relation $L_{21}=TL_{12}$. The quantity $j^{E}-\mu
j^{Q}$ is referred to as the heat current [63]. Eq. (45) contains three
unknown transport coefficients. To solve it, we will consider two different
setups (see Fig. 1): one with $\nabla\mu=\text{const}$ and $\nabla T=0$ (or
equivalently $\nabla\beta=0$), and the other with $\nabla\mu=0$ and $\nabla
T=\text{const}$. This will give us a system of equations from which we can
easily derive $\sigma$, $\alpha$, and $\bar{\kappa}$ (or $\kappa$). We will
refer to the non-equilibrium contribution in each scenario as
$F_{\mu}^{\gtrless}(t)$ and $F_{\beta}^{\gtrless}(t)$ respectively (see Eq.
(35)).

The SYK model is known to exhibit diffusive transport [21]. The hydrodynamic
relations defining the diffusivity matrix $D$ are given by

$\begin{pmatrix}j^{Q}\\\ j^{E}-\mu
j^{Q}\end{pmatrix}=-\begin{pmatrix}D_{11}&D_{12}\\\
D_{21}&D_{22}\end{pmatrix}\begin{pmatrix}\nabla Q\\\ \nabla E-\mu\nabla
Q\end{pmatrix},$ (46)

where $D_{11}$ is the charge diffusion constant, $D_{22}$ is the thermal (not
energy!) diffusion constant, and the off-diagonal elements describe mixed
transport [63]. The diffusivity matrix can be diagonalized, with eigenvalues
$D_{\pm}$ describing the coupled diffusion of charge and heat. It is these
modes that govern the dynamics of charge and energy fluctuations in the system
and are thus more physically relevant than the individual entries in $D$ [63,
21]. In particular, in order for the fluctuations to decay, we must have that
$D_{\pm}\geq 0$, while the matrix elements $D_{ij}$ can be negative.

Finally, combining Eqs. (17, 45, 46) yields the generalized Einstein relation
$L=D\chi$. In the absence of coupling between the charge and energy carriers
(e.g. at $\mu=0$), we have $\chi_{12}=\alpha=0$ and recover the standard
Einstein relations for charge $\sigma=D_{11}K$ and energy
$\kappa=D_{22}C_{\mu}$ transport [63]. However, more generally, one has the
coupled relations given by the full matrix equation.

###  III.5 Chaos

The non-Fermi liquid phase described by the SYK model is known to be highly
chaotic [13, 14] and even saturates a bound on chaos at low temperatures [25].
In such maximally chaotic theories, energy dynamics and diffusion are
fundamentally related to chaos [72, 73]. Moreover, the thermal diffusion
constant of SYK models in the conformal limit is directly controlled by the
butterfly velocity [31, 21, 7, 70], thus realizing a conjectured bound on
diffusion in incoherent metals [63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73].
In this section, we review the many-body chaos properties of the SYK model and
will later show that chaos provides an upper bound on diffusivity in SYK
chains at any temperature and chemical potential. We will mostly follow our
analysis of the Majorana SYK model [7].

We begin by introducing the out-of-time-order correlation function (OTOC),
which has been widely used as a measure of chaos in quantum systems [91, 92,
93, 25, 14, 94]. The regularized OTOC in real time is defined as

$C(x,t_{1},t_{2})=\frac{1}{N^{2}}\sum_{i,j=1}^{N}\Tr[yc_{j}^{x}(t_{1})^{\dagger}yc_{i}^{0}(0)^{\dagger}yc_{j}^{x}(t_{2})yc_{i}^{0}(0)],$
(47)

where $y=e^{-\beta H/4}/Z^{1/4}$ evenly spaces the fermionic fields along the
thermal circle [37]. To leading order, the OTOC can be written as

$C(x,t_{1},t_{2})=\mathcal{F}_{d}-\frac{\mathcal{F}(x,t_{1},t_{2})}{N},$ (48)

where $\mathcal{F}_{d}$ is just a constant corresponding to the disconnected
correlator and $\mathcal{F}$ is the first order contribution stemming from the
contraction of $c_{i}^{0}$ with $c_{j}^{x}$ [95]. For a chaotic system with a
large hierarchy of timescales between thermalization and scrambling, we expect
$\mathcal{F}$ to scale exponentially as $e^{\lambda_{L}t}$, where
$t_{1}=t_{2}=t$ is in the Lyapunov regime $\beta\lesssim t\lesssim\beta\ln N$
[95]. The Lyapunov exponent $\lambda_{L}$ determines the rate of growth of an
operator under Heisenberg evolution and serves as a quantum mechanical measure
for information scrambling in phase space [93]. For an isolated SYK cluster,
$\mathcal{F}(t_{1},t_{2})$ is determined by summing over a set of ladder
diagrams [14, 77, 75, 22, 96, 37, 83], resulting in the self-consistency
equation

$\mathcal{F}(t_{1},t_{2})=\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}t_{3}\mathop{}\\!\mathrm{d}t_{4}K^{R}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}(t_{3},t_{4}),$
(49)

where $K_{R}$ is the retarded kernel

$\begin{split}K^{R}(t_{1},t_{2},t_{3},t_{4})&=(q-1)J^{2}G^{R}(t_{1}-t_{3})G^{A}(t_{4}-t_{2})G^{W}(t_{3}-t_{4})^{\frac{q}{2}-1}G^{W}(t_{4}-t_{3})^{\frac{q}{2}-1}\\\
&=(q-1)J^{2}G^{R}(t_{1}-t_{3})G^{R}(t_{2}-t_{4})^{*}|G^{W}(t_{3}-t_{4})|^{q-2}.\end{split}$
(50)

Here $G^{W}$ is the Wightman Green’s function and we used the symmetry
properties $G^{A}(-t)=G^{R}(t)^{*}$ and $G^{W}(-t)=G^{W}(t)^{*}$ to simplify
the expression. For fermionic systems, the Wightman propagator is related to
the spectral function in frequency space via [37]

$G^{W}(\omega)=\frac{A(\omega)}{2\cosh(\beta\omega/2)}.$ (51)

To determine the Lyapunov exponent, we follow the prescription in Ref. [97],
which works for both Majorana and complex SYK. We define a variant of the
kernel with a parameter $\alpha<0$

$K_{\alpha}^{R}(t,t^{\prime})=\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}se^{\alpha
s}K^{R}\left(s+\frac{t}{2},s-\frac{t}{2},\frac{t^{\prime}}{2},-\frac{t^{\prime}}{2}\right).$
(52)

This operator can be cast in matrix form, with its largest eigenvalue
$k_{R}(\alpha)$ depending on $\alpha$. The Lyapunov exponent is then
determined by the equation $k_{R}(-\lambda_{L})=1$. This condition is
equivalent to $\mathcal{F}$ being an eigenvector of the kernel $K^{R}$ with
eigenvalue one. The Lyapunov exponent of a $(0+1)$-d SYK model is known to
saturate the bound $\lambda_{L}\leq 2\pi/\beta$ at low temperatures [13, 14,
25].

For spatially extended systems, such as our one-dimensional chain, the
operators can also grow in space. Chaos propagation in a translation-invariant
system is described by the Fourier transform of the momentum-space OTOC

$\mathcal{F}(x,t)\sim\int_{-\infty}^{\infty}\frac{\mathop{}\\!\mathrm{d}p}{2\pi}\frac{e^{\lambda_{L}(p)t+ipx}}{\cos(\lambda_{L}(p)\beta/4)},$
(53)

where $\lambda_{L}(p)$ is the momentum-dependent Lyapunov exponent [37]. In
the hydrodynamic limit, this integral can be evaluated using a saddle point
approximation. Depending on the parameters of our model, the integral can
either pick up a contribution solely from the saddle point $p_{s}$, or from
both the saddle point and the momentum-space pole $p_{1}$, both of which are
located on the imaginary axis $p_{s,1}=i|p_{s,1}|$ [97, 37]. In either case,
the result can be written as

$\mathcal{F}(x,t)\sim e^{\lambda_{L}(p_{s,1})(t-|x|/v_{B})},$ (54)

where the butterfly velocity $v_{B}$ is defined as

$v_{B}=\begin{cases}\frac{\lambda_{L}(p_{s})}{|p_{s}|}&\mbox{if
}|p_{s}|<|p_{1}|\\\ \frac{2\pi}{\beta|p_{1}|}&\mbox{if
}|p_{s}|>|p_{1}|\end{cases}$ (55)

Physically, the butterfly velocity $v_{B}$ defines a light-cone that bounds
the speed of operator growth in space [92]. It can also be viewed as a
temperature-dependent extension of the Lieb-Robinson velocity [98].

To compute the butterfly velocity, it is enough to find the momenta $p_{s,1}$
according to Ref. [97]. For a uniform SYK chain, the retarded kernel in
momentum space factorizes $K^{R}(p)=s(p)K^{R}$, where
$s(p)=1+\frac{qJ_{1}^{2}}{2(q-1)J^{2}}(\cos(p)-1)$ is the spatial kernel [31]
and $K^{R}$ is the kernel for a single cluster defined in Eq. (50). Therefore
the eigenvalues of the kernel also factorize
$k_{R}(p,\alpha)=s(p)k_{R}(\alpha)$. The momentum-dependent Lyapunov exponent
can be obtained by solving the equation $k_{R}(p,-\lambda_{L}(p))=1$. At
$p=0$, we recover our previous formula for the Lyapunov exponent of a single
cluster. Generally, this equation has to be solved numerically by repeatedly
diagonalizing the kernel in Eq. (52) and using the bisection method, although
closed-form solutions are available in some limits (see Sec. III.7). Once we
have the entire function $\lambda_{L}(p)$, the location of the saddle can then
be found by solving
$\lambda_{L}(p_{s})=p_{s}\lambda_{L}^{{}^{\prime}}(p_{s})$, while $p_{1}$ is
the momentum at which the Lyapunov exponent attains its maximum value
$\lambda_{L}(p_{1})=2\pi/\beta$.

We can now use the newly introduced measures of chaos to define a
characteristic chaos diffusivity $v_{B}^{2}/\lambda_{L}$, which is known to be
closely related to the thermal diffusion constant in strange metals [21, 31,
7, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]. In fact, we will show that for
all systems under consideration, the chaos diffusivity provides an upper bound
$D_{\pm}\leq v_{B}^{2}/\lambda_{L}$, just as in the Majorana case [7].

###  III.6 $q=2$ limit

Although the Green’s functions generally do not have a closed-form
representation, they simplify significantly in the limits of small and large
$q$, which we discuss next. We start with the special case of $q=2$
corresponding to free fermions, where the system has a quasiparticle
description [21] and the Hamiltonian becomes integrable and non-chaotic [99,
100].

The SD equations are quadratic and can be solved exactly [50]. The spectral
function in equilibrium is given by

$A(\omega)=\frac{2}{J}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}}\quad\mathrm{for}\
|\omega|<2J,$ (56)

and the Green’s functions can be obtained from the fluctuation-dissipation
theorem

$G^{\gtrless}(\omega)=\mp\frac{iA(\omega)}{1+e^{\mp\beta(\omega-\mu)}},$ (57)

followed by an inverse Fourier transform

$G^{\gtrless}(t)=\mp\frac{i}{\pi
J}\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{e^{-i\omega
t}}{1+e^{\mp\beta(\omega-\mu)}}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}}.$
(58)

Finally, we can derive the non-equilibrium contributions from Eq. (35). It is
straightforward to check that $F_{\mu,\beta}^{<}(t)=F_{\mu,\beta}^{>}(t)$ and

$\displaystyle\real\left[G^{<}(t)-G^{>}(t)\right]$ $\displaystyle=0,$ (59)
$\displaystyle\imaginary\left[G^{<}(t)-G^{>}(t)\right]$
$\displaystyle=\frac{B_{1}(2Jt)}{Jt},$ (60)

where $B_{1}$ is the Bessel function of the first kind. With this in mind, we
arrive at the following simplified formulas for the energy gradient and
currents

$\displaystyle\nabla
E=2J\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\frac{B_{1}(2Jt)}{t}\real
F_{\mu,\beta}^{>}(t),$ (61) $\displaystyle j^{Q}$
$\displaystyle=-\frac{J_{1}^{2}}{J}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\frac{B_{1}(2Jt)}{t}\imaginary
F_{\mu,\beta}^{>}(t),$ (62) $\displaystyle j^{E}$
$\displaystyle=-J_{1}^{2}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\frac{B_{1}(2Jt)}{t}\frac{B_{1}(2J(t^{\prime}-t))}{t^{\prime}-t}\real
F_{\mu,\beta}^{>}(t^{\prime}).$ (63)

The charge gradient is still given by Eq. (39). The equations above contain
all the necessary information to calculate the conductivities and
diffusivities numerically at arbitrary $\mu$ and $\beta$ using the non-
equilibrium setups described in Sec. III.4. Moreover, we were able to compute
closed-form results for these transport coefficients in the limit of zero and
infinite temperatures, as described in Appendix C.

###  III.7 Large $q$ limit

We now turn to the opposite limit of large $q$, where an analytic
approximation for the Green’s function at all temperatures is available [14,
21, 96, 101]. To leading order in $1/q$, the Green’s functions for an SYK
model in equilibrium can be expanded as

$\displaystyle G^{>}(t)$
$\displaystyle=-\frac{i}{e^{\beta\mu}+1}\left(1+\frac{g(t)}{q}+\cdots\right)\approx-\frac{ie^{g(t)/q}}{e^{\beta\mu}+1},$
(64) $\displaystyle G^{<}(t)$
$\displaystyle=\frac{i}{e^{-\beta\mu}+1}\left(1+\frac{g(t)^{*}}{q}+\cdots\right)\approx\frac{ie^{g(t)^{*}/q}}{e^{-\beta\mu}+1},$
(65)

where “…” denotes higher order terms, $g(t)$ is a function of order one
satisfying $g(t)^{*}=g(-t)$ and $g(0)=0$ [96]. With this ansatz, the SD
equations are equivalent to a differential equation for $g(t)$

$-\frac{\partial^{2}g(t)}{\partial
t^{2}}=iq\left(\Sigma^{>}(t)+\Sigma^{<}(t)^{*}\right)=2\mathcal{J}^{2}e^{g(t)},$
(66)

where $\mathcal{J}=J\sqrt{q2^{1-q}\cosh^{2-q}(\beta\mu/2)}$ is the effective
coupling. Notice that the original theory has two independent scales $\beta J$
and $\beta\mu$, while the new differential equation only depends on the
combined $\beta\mathcal{J}$ [96]. This holds even after including higher-order
terms in $1/q$ and seems to be an artefact of this expansion [101]. The large
$q$ limit is well defined only when we adjust the original coupling $J$ such
that the re-scaled interaction $\mathcal{J}$ is kept finite as $q\to\infty$.
This implies that $J$ has to be a function of $\beta\mu$, which makes the
comparison with the numerical results at finite $q$ and constant $J$ a bit
more complicated. A direct comparison is possible for $\mu=0$, where we
recover some of our findings for Majorana fermions [7].

The solution to Eq. (66) is of the form found in Ref. [14]

$e^{g(t)}=\frac{\cos^{2}(\pi v/2)}{\cosh^{2}\left(\frac{\pi
v}{\beta}\left(\frac{i\beta}{2}+t\right)\right)},$ (67)

where $v$ satisfies

$\beta\mathcal{J}=\frac{\pi v}{\cos(\pi v/2)}.$ (68)

This gives us the full dependence of the Green’s functions on $\mu$ and
$\beta$. We find it convenient to write the derivatives
$F_{\mu,\beta}^{\gtrless}$ as follows

$\displaystyle\frac{F_{\mu,\beta}^{>}(t)}{G^{>}(t)}$
$\displaystyle=\frac{f_{\mu,\beta}(t)}{q}-\frac{C_{\mu,\beta}}{1+e^{-\beta\mu}},$
(69) $\displaystyle\frac{F_{\mu,\beta}^{<}(t)}{G^{<}(t)}$
$\displaystyle=\frac{f_{\mu,\beta}(t)^{*}}{q}+\frac{C_{\mu,\beta}}{1+e^{\beta\mu}},$
(70)

where $C_{\mu}=\beta\nabla\mu$, $C_{\beta}=\mu\nabla\beta$, and
$f_{\mu,\beta}(t)$ are non-equilibrium contributions to $g(t)$ in the presence
of small gradients (see Ref. [7] for the Majorana case)

$\begin{split}f_{\mu}(t)&=\frac{\mathop{}\\!\mathrm{d}g(t)}{\mathop{}\\!\mathrm{d}\mu}\nabla\mu,\\\
f_{\beta}(t)&=\frac{\mathop{}\\!\mathrm{d}g(t)}{\mathop{}\\!\mathrm{d}\beta}\nabla\beta.\end{split}$
(71)

Our non-equilibrium observables simplify drastically in terms of these
functions

$\displaystyle\nabla Q$
$\displaystyle=\frac{C_{\mu,\beta}}{4\cosh^{2}(\beta\mu/2)},$ (72)
$\displaystyle\nabla E$
$\displaystyle=\frac{\mathcal{J}^{2}}{q^{2}\cosh^{2}(\beta\mu/2)}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[e^{g(t)}\left(2f_{\mu,\beta}(t)-qC_{\mu,\beta}\tanh(\beta\mu/2)\right)\right],$
(73) $\displaystyle j^{Q}$
$\displaystyle=-\frac{\mathcal{J}_{1}^{2}C_{\mu,\beta}}{4\cosh^{2}(\beta\mu/2)}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Re}\left[e^{g(t)}\right]=-\frac{\mathcal{J}_{1}^{2}C_{\mu,\beta}\cos(\pi
v/2)}{4\mathcal{J}\cosh^{2}(\beta\mu/2)},$ (74) $\displaystyle j^{E}$
$\displaystyle=\frac{1}{2}\mathcal{J}_{1}^{2}\mathcal{J}^{2}\real(j_{++}+j_{+-}),$
(75) $\displaystyle j_{++}$
$\displaystyle=\frac{i}{q^{2}\cosh^{2}(\beta\mu/2)}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}e^{g(t)+g(t^{\prime})}\left(f_{\mu,\beta}(t^{\prime})-f_{\mu,\beta}(t)-qC_{\mu,\beta}\tanh(\beta\mu/2)\right),$
(76) $\displaystyle j_{+-}$
$\displaystyle=\frac{i}{q^{2}\cosh^{2}(\beta\mu/2)}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}e^{g(t)^{*}+g(t^{\prime})}\left(f_{\mu,\beta}(t^{\prime})-f_{\mu,\beta}(t)^{*}\right).$
(77)

Notice that the charge gradient and current do not have any dependence on
$f_{\mu,\beta}(t)$ and we can already obtain closed-form expressions for them.
On the other hand, the energy gradient and current require an explicit
calculation of $f_{\mu,\beta}(t)$ for different biases, which we defer to
Appendix C. Nevertheless, we managed to compute the diffusivity and
conductivity matrices analytically for arbitrary $\mu$ and $\beta$ in the
large $q$ limit and will present our results in the next section.

Finally, we comment on the chaos characteristics in this approximation. The
Lyapunov exponent has been previously computed in the large $q$ limit [96,
97]. Since the same derivation applies for both Majorana and complex fermions,
one finds, in our notation

$k_{R}(\alpha)=\frac{8(\pi v)^{2}}{\alpha(\alpha-2\pi v/\beta)\beta^{2}}.$
(78)

Hence the momentum-dependent Lyapunov exponent is given by

$\begin{split}\lambda_{L}(p)&=\frac{\pi
v}{\beta}\left(\sqrt{1+8s(p)}-1\right)\\\ &=\frac{\pi
v}{\beta}\left(\sqrt{9+4\frac{J_{1}^{2}}{J^{2}}\left(\cos(p)-1\right)}-1\right),\end{split}$
(79)

from which the butterfly velocity can be found numerically by solving the
appropriate equations from Sec. III.5. For a single SYK cluster we recover the
well-known answer $\lambda_{L}(0)=2\pi v/\beta$ [96, 97]. We immediately see
that in the low-temperature limit $v\to 1$ and the system is maximally chaotic
[13, 14, 25]. Moreover, in this limit, the butterfly velocity approaches
$v_{B}=\pi\mathcal{J}_{1}/\sqrt{3\beta\mathcal{J}}$ and the thermal diffusion
constant saturates the chaos bound
$D_{+}=D_{22}=v_{B}^{2}/\lambda_{L}=\pi\mathcal{J}_{1}^{2}/6\mathcal{J}$ (see
Eq. (86)).

##  IV Results

Figure 2: Phase diagram of the complex SYK model for $q=4-8$, in units of
$J=\sqrt{3\cdot 2^{q-1}/q}$. The lines correspond to first order phase
transitions between a high-entropy SYK-like phase and a low-entropy harmonic
oscillator-like phase. They start at $(T=0,\mu_{*})$ and end at the critical
point $(T_{c},\mu_{c})$. At high temperatures, the system enters a
supercritical phase.

We report our results on the thermodynamic, transport, and chaos properties of
the complex SYK model in the following sections. We show that there are two
distinct phases in equilibrium, each leading to very different scalings of our
observables. We then study the dependence of the diffusivity and conductivity
matrices on chemical potential and temperature, in relation to the
aforementioned phases. Lastly, we investigate a bound on diffusion imposed by
the chaotic dynamics of the system.

In order to emphasize that our methods are applicable to a range of
parameters, we display results for different interaction orders $q$. To this
extent, we fix the rescaled couplings
$\mathcal{J}_{0,1}|_{\mu=0}=J_{0,1}\sqrt{q2^{1-q}}=1$, which sets the results
for different $q$ on equal footing and allows for a direct comparison to
previously reported values for Majorana fermions [7]. It also keeps $J$
independent of both $\mu$ and $\beta$. Furthermore, we only focus on the
regime with $\mu\geq 0$, since the sign of $\mu$ can be changed by simply
swapping the roles of the creation and annihilation operators in our model.

###  IV.1 Equilibrium phase diagram

We begin by investigating the phase diagram of the SYK model at finite
chemical potential and temperature. In this model, a first order phase
transition arises as a result of the competition between a high-entropy SYK-
like phase and a low-entropy harmonic oscillator-like phase, and it has been
extensively studied in the literature [81, 82, 84, 53, 23].

Figure 3: Phase transition characteristics of the SYK model for $q=4$. Grand
canonical potential (a), charge (b), entropy (c), and charge compressibility
(d) as a function of temperature. Signatures of a first order phase transition
are seen for $0.52\leq\mu\leq 0.83$. The thermodynamic derivatives $Q$ and $S$
experience a jump where the potential $\Omega$ changes slope, while the
compressibility $K$ shows a discontinuity.

The SD equations can have two distinct solutions depending on the values of
$\mu$ and $\beta$. At small $\mu$, the behavior is similar to the Majorana SYK
case. The model has a trivial perturbative expansion around the maximally
mixed state at high temperatures, and a non-trivial conformal regime with an
emergent approximate time-reparametriation symmetry at low temperatures [14].
The latter regime also features a finite zero-temperature entropy and maximal
chaos, reminiscent of nearly extremal black holes. Therefore, we label this
region as the high-entropy or SYK-like phase [81, 82].

On the other hand, in the limit of large $\mu$, the model behaves like a set
of weakly coupled harmonic oscillators and the ground state is given by the
unique Fock vacuum all the way to zero temperature [81, 82]. Hence the system
is non-chaotic, has a vanishingly small entropy at low temperatures and an
exponentially decaying Euclidean two-point function $G(\tau)\sim
e^{-\mu\tau}$. We will refer to this as the low-entropy or harmonic
oscillator-like phase [81, 82].

Figure 4: Temperature dependence of diffusivity $D$ for $q=2$ and multiple
values of $\mu$. Individual entries of the diffusivity matrix (a-d) and its
eigenvalues (e-f) approach finite values at zero temperature in accordance
with Eq. (127).

The two solutions are separated by a finite first order phase transition line,
which starts at $(T=0,\mu_{*})$ and culminates at a critical point
$(T_{c},\mu_{c})$ with asymmetric $q$-dependent critical exponents [81, 82,
84]. For $q=4$, we have $\mu_{*}\approx 0.52$, $T_{c}\approx 0.16$, and
$\mu_{c}\approx 0.83$. At the critical point, the two solutions are identical
and the transition becomes second order. For $T>T_{c}$, the SD equations have
only one solution, corresponding to a high-temperature perturbative regime,
and the system is in a supercritical phase [81, 82]. The high- and low-entropy
phases can be smoothly connected by going around the critical point, which
emphasizes that there is no sharp distinction between them. We summarize these
findings in the phase diagram of Fig. 2. Note that all of our units are
rescaled by a trivial factor of $J$ compared the diagrams in Refs. [81, 82].

On the transition line, the values of $\Omega$ for the two solutions are equal
and the two phases can coexist (see Fig. 3(a)). Upon crossing the line, the
Green’s function jumps from one solution to the other, causing a discontinuity
in the first order derivatives of the potential. This is illustrated in Fig.
3(b-c), where the charge and entropy show clear signs of a first order phase
transition for $\mu_{*}\leq\mu\leq\mu_{c}$. Consequently, second derivatives
experience a singularity at the transition point, as exemplified by the charge
compressibility $K$ in Fig. 3(d). Notice that the $T=0$ state has a finite
entropy and compressibility below $\mu_{*}$, while above $\mu_{*}$ it has
maximal charge and zero entropy and compressibility. This is consistent with
our previous description of the two phases.

The same qualitative behavior is observed for all values of $q\geq 4$, with
the transition line shrinking rapidly as $q$ increases (see Fig. 2). We expect
this transition to disappear completely in the infinite-$q$ limit, as can be
seen explicitly from the thermodynamic potential in Eq. (129). The low-entropy
solution becomes favorable when the second term switches sign from negative to
positive, which never happens at finite temperatures because
$\tan\left(\frac{\pi v}{2}\right)>\frac{\pi v}{4}$ for all $v\in(0,1)$.
Analogously, there is no phase transition in the case of $q=2$ either. The
grand canonical potential in Eq. (96) and its derivatives are smooth,
continuous functions, and the Green’s function always converges to its free-
fermion value.

Figure 5: Temperature dependence of conductivity matrix $L$ for $q=2$ and
several values of $\mu$. Electrical conductivity (a) saturates to a constant
at zero temperature, while thermoelectric (b) and thermal (c) conductivities
approach zero as $1/\beta$. The $\mu$-dependence at low temperatures is
quadratic, in agreement with Eq. (128). (d) The Wiedemann-Franz ratio
converges to its free-fermion value of $\pi^{2}/3$.

###  IV.2 Near-equilibrium transport

The presence of a phase transition has important consequences for both the
transport coefficients and the Lyapunov exponents discussed next. In the high-
entropy phase, we find that transport is diffusive and the Lyapunov exponent
is non-vanishing. Since the high- and low-entropy phases are smoothly
connected by going around the critical point via the super-critical phase, we
expect that the dynamics is diffusive and chaotic throughout the phase
diagram. However, we do observe extreme changes in the diffusivities and
Lyapunov exponents in the vicinity of the phase transition line. Moreover,
while these properties are expected to be non-vanishing, they can be very
small and quite difficult to ascertain numerically. Therefore, when
applicable, we will restrict our analysis to the SYK-like phase, where our
quantities of interest are more straightforward to obtain. Lastly, for all the
parameter regimes considered below, we checked numerically using the same
open-system setup as in Ref. [7], that the NESS solutions of the full
Kadanoff-Baym equations in the presence of weak driving indeed take on the
form in Eq. (34). Thus our ansatz is justified.

Figure 6: Temperature dependence of diffusivity $D$ for $q=4$ and multiple
values of $\mu$. The matrix elements $D_{ij}$ (a-d) and eigenvalues $D_{\pm}$
(e-f) are shown in the SYK-like phase. The eigenvalues reach a constant at low
temperatures and its variation with $\mu$ for different $q$ is shown in the
inset. The diffusivity decreases for both larger $\mu$ and $q$.

We now proceed with our results for the simplest free-fermion case of $q=2$.
We numerically compute all the integrals in Sec. III.6 and extract the
transport coefficients. Their values are plotted in Fig. 4 and Fig. 5 as a
function of inverse temperature and for different $\mu\in[0,2J]$. In the
limits of zero and infinite temperature, we were able to find the linear
response functions analytically and obtained exact solutions for both $D$ and
$L$ in Appendix C. It is easy to check that they agree with our numerical
results in Fig. 4 and Fig. 5 in the corresponding limits. We will elaborate
below on the specific structure of the diffusivity and conductivity matrices
in these limits.

Next, we discuss our results for $q\geq 4$. Since all these cases are very
similar, we focus on $q=4$ in the main panels of Fig. 6 and Fig. 7, with the
understanding that the same conclusions hold for larger $q$. At small $\mu$,
we recover the same behavior as in the Majorana case [7]. For $\mu=0.7$ and
$\mu=0.83$, we encounter the phase transition within our range of
temperatures, and the transport coefficients drop close to zero abruptly. At
large $\mu$, we avoid the phase transition completely and directly enter the
low-entropy phase. In this case, the conductivities and $D_{\pm}$ smoothly
decrease as we lower the temperature. Notice that $D_{11}$ can become negative
as we approach the low-entropy phase, which seems troubling at first. However,
recall that only the eigenvalues $D_{\pm}$ are required to be positive to
ensure the decay of charge and energy fluctuations, which we verify to be the
case in Fig. 6(e-f).

Figure 7: Temperature dependence of conductivity matrix $L$ for $q=4$ and
several values of $\mu$. Electrical (a), thermoelectric (b), and thermal (c)
conductivities are shown in the SYK-like phase. The insets display their low-
temperature asymptotic behavior as a function of $\mu$. (d) The Wiedemann-
Franz ratio converges to $\pi^{2}/12$ for $q=4$. The inset confirms that the
scaling generalizes to $\kappa\beta/\sigma\to 4\pi^{2}/3q^{2}$ for other
values of $q$.

There are a lot of features that emerge from the structure of the diffusivity
matrix for both $q=2$ (Fig. 4) and $q\geq 4$ (Fig. 6). At high temperature,
$D_{12}=0$ and $D_{21}\sim\mu$. The eigenvalues $D_{\pm}$ are identical to the
diagonal entries $D_{11,22}$ and approach finite $\mu$-independent values. At
low temperature and away from the phase transition, the other off-diagonal
entry vanishes $D_{21}=0$ and the eigenvalues $D_{+}=D_{22}$ and
$D_{-}=D_{11}$ converge to $\beta$-independent values. The $\mu$-dependence of
these numbers for different $q$ is shown in the inset of Fig. 6(e-f) and in
Eq. (127). We see that the diffusivity decreases with both $\mu$ and $q$.
These constraints on the diffusivity matrix and the generalized Einstein
relations are enough to conclude that

$\begin{split}\sigma&=D_{11}K,\\\ \kappa&=D_{22}\gamma/\beta,\end{split}$ (80)

at both high and low temperatures. These non-trivial relations are checked
explicitly for $q=2$ in Appendix C. The same dependence among transport
coefficients was found for holographic theories and the SYK chain in the
conformal limit [21]. There it was attributed to the interplay between the
global $U(1)$ charge and the emergent PSL$(2,{\rm I\\!R})$ symmetry. It is
interesting that here we see the same structure also emerge at infinite
temperature.

The conductivity matrix can be examined in the same way (see Fig. 5 and Fig.
7). At high temperature, all the conductivities are zero, while at low
temperature, $\sigma$ is finite and $\kappa$ decays as $1/\beta$. In fact, as
$T\to 0$ we observe a linear-in-T resistivity $\sigma^{-1}$, above a
background residual resistivity $\sigma_{0}^{-1}$, according to the prediction
in [36, 24]

$\frac{1}{\sigma}=\frac{1}{\sigma_{0}}\left(1+4\alpha_{G}\frac{T}{J}\right),$
(81)

where $\alpha_{G}$ is a known numerical constant [14]. A linear fit to our
data yields $\alpha_{G}\approx 0.194$ for $\mu=0$ and $q=4$, which is very
close to the literature value $\alpha_{G}\approx 0.187$ (e.g. Fig. 9 in [14]).
This linear-in-$T$ resistivity is a common feature of many non-Fermi liquid
models [11, 33, 34, 35, 36, 19]. The $\mu$-dependence of the conductivities at
low temperature is available in the inset of Fig. 7 and in Eq. (128). We
notice that $\alpha$ scales linearly with $\mu$, while $\sigma$ and
$\bar{\kappa}$ have a dependence that is closer to quadratic. We can also
combine these results to show that our model has a non-vanishing thermopower
all the way to zero temperature, as discussed further in Appendix D

These observations about the structure of the conductivity matrix lead us to
believe that the Wiedemann-Franz ratio $\kappa\beta/\sigma$ approaches a
constant at zero temperature. Indeed we find numerically in Fig. 7(d) for
$q=4$ and analytically from Eq. (128) for $q=2$ that

$\lim_{\beta\to\infty}\frac{\kappa\beta}{\sigma}=\frac{4\pi^{2}}{3q^{2}},$
(82)

in agreement with the results of Ref. [21]. This also holds for other values
of $q$, as long as we are still in the SYK-like phase, as shown in the inset
of Fig. 7(d). The slight deviations at larger values of $\mu$ are caused by
our inability to numerically reach low enough temperatures without crossing
the phase transition. At zero temperature, we can combine the two results
above to find that the ratio of diffusivities obeys [21]

$\frac{D+}{D_{-}}=\frac{D_{22}}{D_{11}}=\frac{4\pi^{2}}{3q^{2}}\frac{K}{\gamma}.$
(83)

Finally, we are ready to present our findings in the large $q$ limit,
following the derivation in Sec. III.7 and Appendix C. We find that
$D_{12}=0$, which together with the Einstein relations, is enough to conclude
that Eq. (80) holds for all values of $\mu$ and $\beta$. Moreover, we can
combine the conductivities in Eq. (144) and Eq. (147) to arrive at the
Wiedemann-Franz ratio

$\frac{\kappa\beta}{\sigma}=\frac{4\pi^{2}v^{2}}{3q^{2}}.$ (84)

At zero temperature, $v\to 1$ and we recover the results in Eq. (82) and Eq.
(83). Therefore, all the previously found features of transport at finite $q$
are also applicable to the infinite $q$ regime. In addition, this expansion
provides compact solutions for all the transport coefficients over the entire
parameter range (see Appendix C).

Figure 8: Comparison between finite and infinite $q$ diffusivities and
conductivities at $\mu=0$. (a) The diffusivity eigenvalues $D_{\pm}$ (dashed
and dotted lines) are upper bounded by the chaos propagation rate
$v_{B}^{2}/\lambda_{L}$ (solid lines) for all $q$. (b) The electrical
conductivity $\sigma$ (solid lines) and thermal conductivity $\kappa$ (dashed
lines) show little variation with $q$.

In order to make a fair comparison to the finite $q$ results, we restrict
ourselves to the case $\mu=0$, for reasons explained in Appendix C. The charge
and energy diffusion modes decouple and we are left with a diagonal
diffusivity matrix

$\displaystyle D_{-}$
$\displaystyle=D_{11}=\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}\cos(\frac{\pi
v}{2}),$ (85) $\displaystyle D_{+}$
$\displaystyle=D_{22}=\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\left(\frac{\pi
v}{2}\sin(\frac{\pi v}{2})+\cos(\frac{\pi v}{2})\right).$ (86)

The temperature dependence enters the expressions implicitly through $v$ (see
Eq. (68)). We plot these results in Fig. 8(a), together with the finite $q$
values obtained numerically by following the prescription in Sec. III.4. All
the curves obey the same pattern and the agreement with the $q\to\infty$
result clearly improves with increasing $q$. The energy diffusion constant
$D_{22}$ agrees with our previous answer for Majorana fermions [7], as
expected for $\mu=0$. Similarly, the electrical and thermal conductivities are
given by

$\displaystyle\sigma$
$\displaystyle=\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}^{2}}\frac{\pi v}{4},$
(87) $\displaystyle\kappa$ $\displaystyle=\frac{\mathcal{J}_{1}^{2}}{\beta
q^{2}\mathcal{J}^{2}}\frac{(\pi v)^{3}}{3}.$ (88)

These are shown in Fig. 8(b) next to their finite $q$ counterparts. The
agreement is quite good even for moderate values of $q$.

###  IV.3 Chaos

In this section, we further explore the connections between transport and
many-body chaos. We numerically diagonalize the kernel introduced in Sec.
III.5 and extract the Lyapunov exponent $\lambda_{L}$ and the butterfly
velocity $v_{B}$. Both of these quantities only weakly depend on the chemical
potential $\mu$ in the SYK-like phase. In particular, we checked that in the
limit of infinite temperature and finite charge density, both $\lambda_{L}$
and $v_{B}$ saturate the bounds proposed in Ref. [69]. Upon approaching the
phase transition, they decay exponentially with $\mu$ and tend to zero in the
low-entropy phase [96, 83]. This is not surprising, since the conserved $U(1)$
charge constrains the phase-space dynamics of the system. A large chemical
potential eventually renders the system integrable, as manifested by a
transition to the harmonic oscillator-like phase, for which very weak chaotic
behavior is expected.

The ratio $v_{B}^{2}/\lambda_{L}$ exhibits a similar behavior. In Fig. 9 we
plot its temperature dependence in the SYK-like phase for $q=4$ and compare it
to the diffusivity eigenvalues $D_{\pm}$, since these are more physically
relevant than the diagonal entries of $D$. Our results indicate that
$D_{\pm}\leq v_{B}^{2}/\lambda_{L}$ at all temperatures, suggesting that chaos
upper bounds diffusion. This inequality generalizes the bound we previously
found for energy diffusion in Majorana SYK chains [7]. We observe that in the
limit of zero temperature and for $\mu<\mu_{*}$, the thermal diffusivity
saturates the chaos bound $D_{22}=D_{+}=v_{B}^{2}/\lambda_{L}$. This
remarkable result is a consequence of the fact that the same
reparameterization degrees of freedom are responsible both for thermal
diffusion and the OTOC chaos dynamics [14, 31]. The charge diffusivity, on the
other hand, is not easily related to chaos in this model [21]. In Fig. 8(a) we
verify that the same results hold for other values of $q$, as well as in the
large $q$ limit.

The fact that the SYK chain reaches this equality in the conformal limit has
been previously shown for both Majorana and complex fermions [31, 21].
However, our method for calculating the diffusivities at arbitrary $\mu$ and
$\beta$ allows us to confirm the inequality $D_{\pm}\leq
v_{B}^{2}/\lambda_{L}$ beyond the conformal or large $q$ limits [70]. A
similar bound has been found for other families of models as well [68, 102,
69]. We should mention that this inequality is by no means universal, since
there are examples of theories where it holds in the opposite direction [64,
65, 66]. A more rigorous upper bound on diffusivity can be written in the form
of $D\leq v^{2}\tau_{\text{eq}}$ [67, 103, 104, 71], where $v\sim v_{B}$ is
the operator growth velocity and $\tau_{\text{eq}}$ is the local equilibration
timescale, which can be much larger than the Lyapunov timescale
$1/\lambda_{L}$ [67].

Figure 9: The diffusivity eigenvalues $D_{\pm}$ (dashed and dotted lines) and
the chaos bound $v_{B}^{2}/\lambda_{L}$ (solid lines) for $q=4$ and different
chemical potentials. The bound is saturated in the conformal limit. All three
quantities drop to zero outside the SYK-like phase.

##  V Discussion

This work has described the thermodynamic, transport, and chaos properties of
an SYK chain with general $q$-body interactions. Our main result is a detailed
analysis of the near-equilibrium response of local Green’s functions to small
external biases. More specifically, we expanded the Green’s function of each
SYK cluster to first order in the non-equilibrium corrections
$F_{\mu,\beta}^{\gtrless}$ due to constant chemical potential and temperature
gradients. We were then able to express all the conserved charges and their
associated currents in terms of these functions. The calculations were carried
out analytically for $q=2$ and $q\to\infty$, and numerically for all other
values of $q$ using the solutions of the Schwinger-Dyson equation described in
Appendix A. This allowed us to fully characterize the mixed thermoelectric
response of the model in terms of its diffusivity, conductivity, and
susceptibility matrices. Moreover, we showed that the eigenvalues of the
diffusivity matrix satisfy the inequality $D_{\pm}\leq v_{B}^{2}/\lambda_{L}$
at all temperatures, with equality achieved for $D_{+}$ in the conformal
limit. This result generalizes our previous bound on energy diffusion in the
case of Majorana fermions [7] and establishes a connection between transport
and chaos in the SYK model.

Our analysis has revealed new features in the structure of the transport
coefficients. In particular, we showed that one of the off-diagonal entries of
the diffusivity matrix approaches zero at both high and low temperatures, as
well as in the large $q$ limit. Together with the Einstein relations, this
results in a simplified expression for the conductivities in Eq. (80), which
was previously established only for SYK and holographic models in the
conformal limit [21]. Additionally, we showed that the Wiedemann-Franz ratio
approaches the finite value $4\pi^{2}/3q^{2}$ at zero temperature, in
agreement with Ref. [21]. For $q=2$ we recover the universal Fermi liquid
prediction in the form of the Lorenz number $\pi^{2}/3$. We should emphasize
that for $q\geq 4$, this result is not universal and depends on the specific
choice of interaction between clusters [21].

Although our methods are valid for arbitrary values of $\beta$ and $\mu$, we
have to be careful when interpreting our results close to the phase transition
between the high- and low-entropy phases and in the low-entropy phase.
Specifically, various transport quantities experience a sudden drop or
divergence when crossing the transition. Moreover, since the low-entropy phase
has very small values of the diffusivity and Lyapunov exponent, we cannot
always robustly study the low temperature limit of our observables past the
phase transition. Hence, in some instances, we have to restrict ourselves to
small values of $\mu$, where the SYK-like phase extends all the way to zero
temperature. Given the weakly coupled nature of the low-entropy phase, other
analytical methods may be useful in that regime, if the physics is of
interest.

Our work paves the way for further analytical and numerical studies of linear
response in quantum many-body systems. In this paper, we focused on the zero-
frequency response of a uniform one-dimensional chain, but generalizations
should be straightforward. For example, the frequency dependence of the
transport coefficients can be extracted by imposing a time-dependent
oscillatory bias [7, 52]. Our methods are also suited for other higher-
dimensional non-Fermi liquid models built from SYK clusters [35, 34], or more
general theories with tractable local Green’s functions. The same ideas can in
principle be applied to study transport in more conventional spin systems [1,
4, 5, 6], where the NESS is approximated as a tensor network, although the
details of this calculation are more complicated.

Despite its success predicting the main features of transport, linear response
theory has some limitations as a probe of non-equilibrium dynamics in quantum
systems. It would be interesting to investigate the effect of strong driving
on our SYK system, where non-linear effects, such as Joule heating, play an
important role. To capture the physics beyond linear response, one would have
to solve the full Kadanoff-Baym for the system out-of-equilibrium [7]. Probing
non-linear transport and out-of-equilibrium phase transitions are both
interesting future directions.

The SYK chain discussed in this paper is a solvable theoretical model
displaying some of the major properties of a non-Fermi liquid [19]. However,
its non-Fermi liquid behavior is yet to be observed experimentally. In recent
years, multiple experimental realizations [38, 39, 40, 41, 42, 43, 44, 45] and
quantum simulations [46, 47, 48, 49] of SYK have been proposed. These include
ultracold atom experiments [40, 41], Majorana modes at the interface of a
topological insulator and superconductor [42], semiconductor wires coupled
through a disordered quantum dot [43], superconducting circuits [44], and
graphene flakes [45]. The latter configuration, based on the zeroth Landau
level in graphene flakes with irregular boundaries subject to strong magnetic
fields, is especially well suited for probing mesoscopic transport in the
complex SYK model [60, 105]. One could use this setup to look for signatures
of a linear-in-$T$ resistivity at low temperatures according to Eq. (81). It
has also been suggested that measurements of the thermopower can serve as an
indicator of the non-vanishing residual entropy at low temperatures [105].
This opens up the possibility of directly comparing our theoretical
predictions with actual experimental data.

###### Acknowledgements.

We are grateful to Nikolay Gnezdilov, Aavishkar Patel, Nilakash Sorokhaibam,
and Maria Tikhanovskaya for helpful discussions related to the numerical
solutions of the SYK model. C.Z. acknowledges financial support from the
Harvard-MIT Center for Ultracold Atoms through NSF Grant No. PHY-1734011. The
work of B.S. is supported in part by AFOSR grant FA9550-19-1-0360.

## References

  * Bertini  _et al._ [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Rev. Mod. Phys. 93, 025003 (2021).
  * Zotos  _et al._ [1997] X. Zotos, F. Naef, and P. Prelovsek, Phys. Rev. B 55, 11029 (1997).
  * Kapustin and Spodyneiko [2021] A. Kapustin and L. Spodyneiko, Phys. Rev. B 104, 035150 (2021).
  * Weimer  _et al._ [2021] H. Weimer, A. Kshetrimayum, and R. Orús, Rev. Mod. Phys. 93, 015008 (2021).
  * Landi  _et al._ [2021] G. T. Landi, D. Poletti, and G. Schaller, Non-equilibrium boundary driven quantum systems: models, methods and properties (2021), arXiv:2104.14350 .
  * Zanoci and Swingle [2021] C. Zanoci and B. Swingle, Phys. Rev. B 103, 115148 (2021).
  * Zanoci and Swingle [2022] C. Zanoci and B. Swingle, Phys. Rev. Research 4, 023001 (2022).
  * Georges  _et al._ [2000] A. Georges, O. Parcollet, and S. Sachdev, Phys. Rev. Lett. 85, 840 (2000).
  * Georges  _et al._ [2001] A. Georges, O. Parcollet, and S. Sachdev, Phys. Rev. B 63, 134406 (2001).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Phys. Rev. Lett. 70, 3339 (1993).
  * Parcollet and Georges [1999] O. Parcollet and A. Georges, Phys. Rev. B 59, 5341 (1999).
  * Sachdev [2015] S. Sachdev, Phys. Rev. X 5, 041025 (2015).
  * Kitaev [2015] A. Kitaev, in  _Proceedings of KITP: A Simple Model of Quantum Holography_ , Entanglement in Strongly-Correlated Quantum Matter (Kavli Institute of Theoretical Physics, 2015).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Phys. Rev. D 94, 106002 (2016).
  * Kitaev and Suh [2018] A. Kitaev and S. J. Suh, J. High Energy Phys. 2018 (5), 183.
  * Sarosi [2018] G. Sarosi, Proc. Sci. Modave2017, 001 (2018).
  * Rosenhaus [2019] V. Rosenhaus, J. Phys. A Math. Theor. 52, 323001 (2019).
  * Gu  _et al._ [2020] Y. Gu, A. Kitaev, S. Sachdev, and G. Tarnopolsky, J. High Energy Phys. 2020 (2), 157.
  * Chowdhury  _et al._ [2021] D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, Sachdev-ye-kitaev models and beyond: A window into non-fermi liquids (2021), arXiv:2109.05037 .
  * Fu and Sachdev [2016] W. Fu and S. Sachdev, Phys. Rev. B 94, 035135 (2016).
  * Davison  _et al._ [2017] R. A. Davison, W. Fu, A. Georges, Y. Gu, K. Jensen, and S. Sachdev, Phys. Rev. B 95, 155131 (2017).
  * Bulycheva [2017] K. Bulycheva, J. High Energy Phys. 2017 (12), 69.
  * Tikhanovskaya  _et al._ [2021a] M. Tikhanovskaya, H. Guo, S. Sachdev, and G. Tarnopolsky, Phys. Rev. B 103, 075141 (2021a).
  * Tikhanovskaya  _et al._ [2021b] M. Tikhanovskaya, H. Guo, S. Sachdev, and G. Tarnopolsky, Phys. Rev. B 103, 075142 (2021b).
  * Maldacena  _et al._ [2016a] J. Maldacena, S. H. Shenker, and D. Stanford, J. High Energy Phys. 2016 (8), 106.
  * Sachdev [2010a] S. Sachdev, Phys. Rev. Lett. 105, 151602 (2010a).
  * Sachdev [2010b] S. Sachdev, J. Stat. Mech. Theory Exp. 2010, P11022 (2010b).
  * Almheiri and Polchinski [2015] A. Almheiri and J. Polchinski, J. High Energy Phys. 2015 (11), 14.
  * Maldacena  _et al._ [2016b] J. Maldacena, D. Stanford, and Z. Yang, Prog. Theor. Exp. Phys. 2016 (2016b).
  * Engelsöy  _et al._ [2016] J. Engelsöy, T. G. Mertens, and H. Verlinde, J. High Energy Phys. 2016 (7), 139.
  * Gu  _et al._ [2017a] Y. Gu, X.-L. Qi, and D. Stanford, J. High Energy Phys. 2017 (5), 125.
  * Jian  _et al._ [2017] C.-M. Jian, Z. Bi, and C. Xu, Phys. Rev. B 96, 115122 (2017).
  * Song  _et al._ [2017] X.-Y. Song, C.-M. Jian, and L. Balents, Phys. Rev. Lett. 119, 216601 (2017).
  * Patel  _et al._ [2018] A. A. Patel, J. McGreevy, D. P. Arovas, and S. Sachdev, Phys. Rev. X 8, 021049 (2018).
  * Chowdhury  _et al._ [2018] D. Chowdhury, Y. Werman, E. Berg, and T. Senthil, Phys. Rev. X 8, 031024 (2018).
  * Guo  _et al._ [2020] H. Guo, Y. Gu, and S. Sachdev, Ann. Phys. 418, 168202 (2020).
  * Guo  _et al._ [2019] H. Guo, Y. Gu, and S. Sachdev, Phys. Rev. B 100, 045140 (2019).
  * Franz and Rozali [2018] M. Franz and M. Rozali, Nat. Rev. Mater. 3, 491 (2018).
  * Rahmani and Franz [2019] A. Rahmani and M. Franz, Rep. Prog. Phys. 82, 084501 (2019).
  * Danshita  _et al._ [2017] I. Danshita, M. Hanada, and M. Tezuka, Prog. Theor. Exp. Phys. 2017 (2017).
  * Wei and Sedrakyan [2021] C. Wei and T. A. Sedrakyan, Phys. Rev. A 103, 013323 (2021).
  * Pikulin and Franz [2017] D. I. Pikulin and M. Franz, Phys. Rev. X 7, 031006 (2017).
  * Chew  _et al._ [2017] A. Chew, A. Essin, and J. Alicea, Phys. Rev. B 96, 121119 (2017).
  * Yang  _et al._ [2018] F. Yang, L. Henriet, A. Soret, and K. Le Hur, Phys. Rev. B 98, 035431 (2018).
  * Chen  _et al._ [2018] A. Chen, R. Ilan, F. de Juan, D. I. Pikulin, and M. Franz, Phys. Rev. Lett. 121, 036403 (2018).
  * García-Álvarez  _et al._ [2017] L. García-Álvarez, I. L. Egusquiza, L. Lamata, A. del Campo, J. Sonner, and E. Solano, Phys. Rev. Lett. 119, 040501 (2017).
  * Luo  _et al._ [2019] Z. Luo, Y.-Z. You, J. Li, C.-M. Jian, D. Lu, C. Xu, B. Zeng, and R. Laflamme, NPJ Quantum Inf. 5, 53 (2019).
  * Babbush  _et al._ [2019] R. Babbush, D. W. Berry, and H. Neven, Phys. Rev. A 99, 040301 (2019).
  * Behrends and Béri [2022] J. Behrends and B. Béri, Phys. Rev. Lett. 128, 106805 (2022).
  * Eberlein  _et al._ [2017] A. Eberlein, V. Kasper, S. Sachdev, and J. Steinberg, Phys. Rev. B 96, 205123 (2017).
  * Bhattacharya  _et al._ [2019] R. Bhattacharya, D. P. Jatkar, and N. Sorokhaibam, J. High Energy Phys. 2019 (7), 66.
  * Kuhlenkamp and Knap [2020] C. Kuhlenkamp and M. Knap, Phys. Rev. Lett. 124, 106401 (2020).
  * Samui and Sorokhaibam [2021] T. Samui and N. Sorokhaibam, J. High Energy Phys. 2021 (4).
  * Louw and Kehrein [2022] J. C. Louw and S. Kehrein, Phys. Rev. B 105, 075117 (2022).
  * Almheiri  _et al._ [2019] A. Almheiri, A. Milekhin, and B. Swingle, Universal constraints on energy flow and syk thermalization (2019), arXiv:1912.04912 .
  * Zhang [2019] P. Zhang, Phys. Rev. B 100, 245104 (2019).
  * Cheipesh  _et al._ [2021] Y. Cheipesh, A. I. Pavlov, V. Ohanesjan, K. Schalm, and N. V. Gnezdilov, Phys. Rev. B 104, 115134 (2021).
  * Haldar  _et al._ [2020] A. Haldar, P. Haldar, S. Bera, I. Mandal, and S. Banerjee, Phys. Rev. Research 2, 013307 (2020).
  * Chen  _et al._ [2017] Y. Chen, H. Zhai, and P. Zhang, J. High Energy Phys. 2017 (7), 150.
  * Can  _et al._ [2019] O. Can, E. M. Nica, and M. Franz, Phys. Rev. B 99, 045419 (2019).
  * Sá  _et al._ [2021] L. Sá, P. Ribeiro, and T. Prosen, Lindbladian dissipation of strongly-correlated quantum matter (2021), arXiv:2112.12109 .
  * Kulkarni  _et al._ [2021] A. Kulkarni, T. Numasawa, and S. Ryu, Syk lindbladian (2021), arXiv:2112.13489 .
  * Hartnoll [2014] S. A. Hartnoll, Nat. Phys. 11, 54 (2014).
  * Blake [2016a] M. Blake, Phys. Rev. Lett. 117, 091601 (2016a).
  * Blake [2016b] M. Blake, Phys. Rev. D 94, 086014 (2016b).
  * Blake  _et al._ [2017] M. Blake, R. A. Davison, and S. Sachdev, Phys. Rev. D 96, 106008 (2017).
  * Hartman  _et al._ [2017] T. Hartman, S. A. Hartnoll, and R. Mahajan, Phys. Rev. Lett. 119, 141601 (2017).
  * Gu  _et al._ [2017b] Y. Gu, A. Lucas, and X.-L. Qi, SciPost Phys. 2, 018 (2017b).
  * Chen  _et al._ [2020] X. Chen, Y. Gu, and A. Lucas, SciPost Phys. 9, 71 (2020).
  * Choi  _et al._ [2021] C. Choi, M. Mezei, and G. Sárosi, J. High Energy Phys. 2021 (2), 207.
  * Hartnoll and Mackenzie [2021] S. A. Hartnoll and A. P. Mackenzie, Planckian dissipation in metals (2021), arXiv:2107.07802 .
  * Blake  _et al._ [2018] M. Blake, H. Lee, and H. Liu, J. High Energy Phys. 2018 (10), 127.
  * Blake and Liu [2021] M. Blake and H. Liu, J. High Energy Phys. 2021 (5), 229.
  * Gu  _et al._ [2017c] Y. Gu, A. Lucas, and X.-L. Qi, J. High Energy Phys. 2017 (9), 120.
  * Zhang [2017] P. Zhang, Phys. Rev. B 96, 205138 (2017).
  * Bentsen  _et al._ [2019] G. Bentsen, Y. Gu, and A. Lucas, Proc. Natl. Acad. Sci. USA 116, 6689 (2019).
  * Banerjee and Altman [2017] S. Banerjee and E. Altman, Phys. Rev. B 95, 134302 (2017).
  * Haldar  _et al._ [2018] A. Haldar, S. Banerjee, and V. B. Shenoy, Phys. Rev. B 97, 241106 (2018).
  * Jian and Yao [2017] S.-K. Jian and H. Yao, Phys. Rev. Lett. 119, 206602 (2017).
  * Cai  _et al._ [2018] W. Cai, X.-H. Ge, and G.-H. Yang, J. High Energy Phys. 2018 (1), 76.
  * Azeyanagi  _et al._ [2018] T. Azeyanagi, F. Ferrari, and F. I. S. Massolo, Phys. Rev. Lett. 120, 061602 (2018).
  * Ferrari and Schaposnik Massolo [2019] F. Ferrari and F. I. Schaposnik Massolo, Phys. Rev. D 100, 026007 (2019).
  * Sorokhaibam [2020] N. Sorokhaibam, J. High Energy Phys. 2020 (7), 55.
  * Cao  _et al._ [2021] S. Cao, Y.-C. Rui, and X.-H. Ge, Thermodynamic phase structure of complex sachdev-ye-kitaev model and charged black hole in deformed jt gravity (2021), arXiv:2103.16270 .
  * de Groot and Mazur [1984] S. de Groot and P. Mazur,  _Non-equilibrium Thermodynamics_, Dover Books on Physics (Dover Publications, New York, 1984).
  * Kubo  _et al._ [1985] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical Physics II: Nonequilibrium Statistical Mechanics_, Springer series in solid-state sciences (Springer-Verlag, New York, 1985).
  * Forster [1975] D. Forster,  _Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions_ (W. A. Benjamin, Advanced Book Program, 1975).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, Cambridge, England, 2011).
  * 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, 2013).
  * Babadi  _et al._ [2015] M. Babadi, E. Demler, and M. Knap, Phys. Rev. X 5, 041005 (2015).
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP 28, 1200 (1969).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, J. High Energy Phys. 2014 (3), 67.
  * Shenker and Stanford [2015] S. H. Shenker and D. Stanford, J. High Energy Phys. 2015 (5), 132.
  * Kobrin  _et al._ [2021] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, Phys. Rev. Lett. 126, 030602 (2021).
  * Mezei and Sárosi [2020] M. Mezei and G. Sárosi, J. High Energy Phys. 2020 (1), 186.
  * Bhattacharya  _et al._ [2017] R. Bhattacharya, S. Chakrabarti, D. P. Jatkar, and A. Kundu, J. High Energy Phys. 2017 (11), 180.
  * Gu and Kitaev [2019] Y. Gu and A. Kitaev, J. High Energy Phys. 2019 (2), 75.
  * Roberts and Swingle [2016] D. A. Roberts and B. Swingle, Phys. Rev. Lett. 117, 091602 (2016).
  * García-García  _et al._ [2018] A. M. García-García, B. Loureiro, A. Romero-Bermúdez, and M. Tezuka, Phys. Rev. Lett. 120, 241603 (2018).
  * Haque and McClarty [2019] M. Haque and P. A. McClarty, Phys. Rev. B 100, 115122 (2019).
  * Tarnopolsky [2019] G. Tarnopolsky, Phys. Rev. D 99, 026010 (2019).
  * Lucas and Steinberg [2016] A. Lucas and J. Steinberg, J. High Energy Phys. 2016 (10), 143.
  * Lucas [2021] A. Lucas, Constraints on hydrodynamics from many-body quantum chaos (2021), arXiv:1710.01005 .
  * Han and Hartnoll [2018] X. Han and S. A. Hartnoll, Phys. Rev. Lett. 121, 170601 (2018).
  * Kruchkov  _et al._ [2020] A. Kruchkov, A. A. Patel, P. Kim, and S. Sachdev, Phys. Rev. B 101, 205148 (2020).
  * Brian and Liang [2018] S. Brian and F. Liang, Sci. Adv. 4, eaat2621 (2018).
  * Kozii  _et al._ [2019] V. Kozii, B. Skinner, and L. Fu, Phys. Rev. B 99, 155123 (2019).
  * Han  _et al._ [2020] F. Han, N. Andrejevic, T. Nguyen, V. Kozii, Q. T. Nguyen, T. Hogan, Z. Ding, R. Pablo-Pedro, S. Parjan, B. Skinner, A. Alatas, E. Alp, S. Chi, J. Fernandez-Baca, S. Huang, L. Fu, and M. Li, Nat. Commun. 11, 6167 (2020).

##  Appendix A Numerical solutions of the SD equations

As mentioned in the main text, the Schwinger-Dyson equations in both real and
imaginary time are solved numerically for the case of a single isolated SYK
cluster in equilibrium. Our approach in Euclidean time is almost identical to
the original one for Majorana fermions [14]:

  1. 1.

Initialize $G(i\omega_{n})$ with the free-fermion propagator
$(i\omega_{n}+\mu)^{-1}$ and compute its inverse Fourier transform $G(\tau)$.

  2. 2.

Calculate $\Sigma(\tau)$ using the second line in Eq. (13) and Fourier
transform it to $\Sigma(i\omega_{n})$.

  3. 3.

Compute a new Green’s function $\tilde{G}(i\omega_{n})$ from the first line in
Eq. (13) and get its inverse Fourier transform $\tilde{G}(\tau)$.

  4. 4.

Perform a weighted update
$G(\tau)\leftarrow(1-\alpha)G(\tau)+\alpha\tilde{G}(\tau)$ with $\alpha=0.3$.

  5. 5.

Repeat steps 2-4 until the iterative procedure converges
$\max_{\tau}|G(\tau)-\tilde{G}(\tau)|<\epsilon$ with $\epsilon=10^{-12}$.

The imaginary time domain is discretized as $\tau_{n}=nd\tau$, where
$d\tau=\beta/M$, $-M\leq n\leq M$, and $M=2^{20}$. Similarly, the Matsubara
frequencies are given by $\omega_{n}=(2n+1)d\omega$, with $d\omega=\pi/\beta$
and $-M\leq n\leq M-1$. All the thermodynamic properties are then derived from
the grand canonical potential in Eq. (14).

Our algorithm for obtaining the real-time Green’s functions is an extension of
the method used in Ref. [50] for Majorana fermions:

  1. 1.

Initialize $G^{\gtrless}(t)$ with the $q=2$ result in Eq. (58).

  2. 2.

Calculate $\Sigma^{\gtrless}(t)$ using the second line in Eq. (31) (for a
single cluster $G_{x}^{\gtrless}=G_{x\pm 1}^{\gtrless}=G^{\gtrless}$).
Evaluate the retarded self-energy
$\Sigma^{R}(t)=\Theta(t)(\Sigma^{>}(t)-\Sigma^{<}(t))$ and its Fourier
transform $\Sigma^{R}(\omega)$.

  3. 3.

Compute $G^{R}(\omega)$ from the first line in Eq. (13) and find the spectral
function $A(\omega)=-2\imaginary G^{R}(\omega)$.

  4. 4.

Determine a new Green’s function $\tilde{G}^{\gtrless}(\omega)$ from the FDT
in Eq. (57) and get its inverse Fourier transform $\tilde{G}^{\gtrless}(t)$.

  5. 5.

Perform a weighted update
$G^{\gtrless}(t)\leftarrow(1-\alpha)G^{\gtrless}(t)+\alpha\tilde{G}^{\gtrless}(t)$
with $\alpha=0.7$.

  6. 6.

Repeat steps 2-5 until the iterative procedure converges
$\max_{t}|G^{\gtrless}(t)-\tilde{G}^{\gtrless}(t)|<\epsilon$ with
$\epsilon=10^{-5}$.

The real time domain is discretized as $t_{n}=ndt$, where $dt=0.05$, $-M\leq
n\leq M$, and $M=10^{4}$. Similarly, Fourier transform frequencies are given
by $\omega_{n}=nd\omega$, with $d\omega=2\pi/Mdt$ and $-M\leq n\leq M$. At
large values of $\mu$ and $\beta$, the above procedure can experience
convergence issues. In order to mitigate this problem, we perform an
additional annealing step, where we start at a high temperature and gradually
lower it while re-running the algorithm with the Green’s functions initialized
to the previously converged values from a higher temperature run.

##  Appendix B Charge and energy currents

In this section, we present a derivation of the formulas for the charge and
energy currents introduced in Sec. III.4. The non-equilibrium expectation
value of any operator can be computed in the Keldysh formalism using the
generating functional [88, 89]. For a one-dimensional chain, the charge
current at site $x$ is given by

$\begin{split}j_{x}^{Q}&=i[Q_{x},H^{x,x+1}]=i[Q_{x},H_{1}^{x,x+1}]\\\
&=i\sum_{\\{i\\},\\{j\\}}\left(J_{i_{1}\ldots
j_{\frac{q}{2}}}^{(1)}(c_{i_{1}}^{x})^{\dagger}\cdots(c_{i_{\frac{q}{2}}}^{x})^{\dagger}c_{j_{1}}^{x+1}\cdots
c_{j_{\frac{q}{2}}}^{x+1}-\mathrm{h.c.}\right),\end{split}$ (89)

and its expectation value is similar to the on-bond energy

$\begin{split}\langle
j_{x}^{Q}\rangle&=-\frac{J_{1}^{2}}{2}\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{1}\Big{(}G_{x}^{>}(t,t_{1})^{\frac{q}{2}}G_{x+1}^{<}(t_{1},t)^{\frac{q}{2}}\\\
&-G_{x}^{<}(t_{1},t)^{\frac{q}{2}}G_{x+1}^{>}(t,t_{1})^{\frac{q}{2}}+\mathrm{h.c.}\Big{)}.\end{split}$
(90)

This result is equivalent to the tunneling current formula derived in Ref.
[57] for $q=2$.

In a similar fashion, the energy current at site $x$ is

$\begin{split}j_{x}^{E}&=i[H^{x-1,x},H^{x,x+1}]\\\
&=i\left(\frac{1}{2}[H_{0}^{x},H_{1}^{x,x+1}-H_{1}^{x-1,x}]+[H_{1}^{x-1,x},H_{1}^{x,x+1}]\right).\end{split}$
(91)

The calculation of its expectation value mirrors the one for the energy
current in the Majorana case [7], except that now we have extra contributions
form the Hermitian conjugate terms in the interaction Hamiltonian

$\langle
j_{x}^{E}\rangle=\frac{1}{2}J_{1}^{2}J^{2}\real(j_{++}^{x-1,x,x+1}+j_{+-}^{x-1,x,x+1}),$
(92)

$\displaystyle j_{++}^{x-1,x,x+1}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{<}(t_{1},t)^{\frac{q}{2}}G_{x}^{<}(t_{2},t)^{\frac{q}{2}-1}G_{x}^{>}(t,t_{1})^{\frac{q}{2}-1}G_{x}^{>}(t_{2},t_{1})G_{x+1}^{>}(t,t_{2})^{\frac{q}{2}}$
$\displaystyle-\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{>}(t,t_{1})^{\frac{q}{2}}G_{x}^{>}(t,t_{2})^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t)^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t_{2})G_{x+1}^{<}(t_{2},t)^{\frac{q}{2}}$
$\displaystyle+\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{<}(t_{2},t)^{\frac{q}{2}}G_{x}^{<}(t_{1},t)^{\frac{q}{2}-1}G_{x}^{>}(t,t_{2})^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t_{2})G_{x+1}^{>}(t,t_{1})^{\frac{q}{2}}$
$\displaystyle-\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{>}(t,t_{2})^{\frac{q}{2}}G_{x}^{>}(t,t_{1})^{\frac{q}{2}-1}G_{x}^{<}(t_{2},t)^{\frac{q}{2}-1}G_{x}^{>}(t_{2},t_{1})G_{x+1}^{<}(t_{1},t)^{\frac{q}{2}},$
(93) $\displaystyle j_{+-}^{x-1,x,x+1}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{<}(t_{1},t)^{\frac{q}{2}}G_{x}^{>}(t_{2},t)^{\frac{q}{2}-1}G_{x}^{>}(t,t_{1})^{\frac{q}{2}-1}G_{x}^{>}(t_{2},t_{1})G_{x+1}^{<}(t,t_{2})^{\frac{q}{2}}$
$\displaystyle-\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{>}(t,t_{1})^{\frac{q}{2}}G_{x}^{<}(t,t_{2})^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t)^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t_{2})G_{x+1}^{>}(t_{2},t)^{\frac{q}{2}}$
$\displaystyle+\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{>}(t_{2},t)^{\frac{q}{2}}G_{x}^{<}(t_{1},t)^{\frac{q}{2}-1}G_{x}^{<}(t,t_{2})^{\frac{q}{2}-1}G_{x}^{<}(t_{1},t_{2})G_{x+1}^{>}(t,t_{1})^{\frac{q}{2}}$
$\displaystyle-\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{x-1}^{<}(t,t_{2})^{\frac{q}{2}}G_{x}^{>}(t,t_{1})^{\frac{q}{2}-1}G_{x}^{>}(t_{2},t)^{\frac{q}{2}-1}G_{x}^{>}(t_{2},t_{1})G_{x+1}^{<}(t_{1},t)^{\frac{q}{2}}.$
(94)

These expressions simplify greatly for a uniform chain in the near-equilibrium
linear response regime, as shown in Sec. III.4.

##  Appendix C Exact calculations of the transport coefficients

In Sec. III.4 we introduced a simple non-equilibrium correction
$F_{\mu,\beta}^{\gtrless}(t)$ to the Green’s functions in the linear response
regime. We showed that the conserved quantities and their currents can be
expanded to first order in this function. Next, we will consider special cases
where it is possible to analytically compute the equilibrium Green’s
functions, and hence also $F_{\mu,\beta}^{\gtrless}(t)$. In particular, we
will provide detailed derivations for the susceptibility, diffusivity, and
conductivity matrices in the limit of zero and infinite temperature for $q=2$,
as well as in the $q\to\infty$ limit at arbitrary temperature. This will be a
continuation of our discussion of these limits in Sec. III.6 and Sec. III.7.

###  C.1 $q=2$ limit

The SYK Hamiltonian for $q=2$ is equivalent to a random hopping model. The
partition function can be computed directly from the free-fermion picture
[14]. After fixing the reference energy level to match our convention for the
charge in Eq. (6), we have

$Z=\prod_{|\omega|<2J}e^{-\beta\mu/2}\left(1+e^{-\beta(\omega-\mu)}\right)$
(95)

Note that unlike Majoranas, complex fermions are not paired up when performing
the product over all modes. The thermodynamic potential becomes

$\Omega=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega
g(\omega)\left(\log(1+e^{-\beta(\omega-\mu)})-\frac{\beta\mu}{2}\right),$ (96)

where we introduced the normalized density of states
$g(\omega)=A(\omega)/2\pi$. The entries of the susceptibility matrix can be
found by taking second derivatives of the potential

$\displaystyle\chi_{11}$
$\displaystyle=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{\beta}{4\pi
J\cosh^{2}(\beta(\omega-\mu)/2)}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}},$
(97) $\displaystyle\chi_{12}$
$\displaystyle=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{\beta^{2}(\omega-\mu)}{4\pi
J\cosh^{2}(\beta(\omega-\mu)/2)}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}},$
(98) $\displaystyle\chi_{22}$
$\displaystyle=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{\beta^{2}(\omega-\mu)^{2}}{4\pi
J\cosh^{2}(\beta(\omega-\mu)/2)}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}}.$
(99)

Now switching over to transport, from Eq. (58) we can deduce

$\displaystyle F_{\beta}^{>}(t)$
$\displaystyle=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{i(\mu-\omega)e^{-i\omega
t}\nabla\beta}{4\pi
J\cosh^{2}(\beta(\omega-\mu)/2)}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}},$
(100) $\displaystyle F_{\mu}^{>}(t)$
$\displaystyle=\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{i\beta
e^{-i\omega t}\nabla\mu}{4\pi
J\cosh^{2}(\beta(\omega-\mu)/2)}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}}.$
(101)

All these integrals can always be performed numerically, but in the special
case of high or low temperatures, we can evaluate them analytically.

####  C.1.1 Infinite temperature limit

First consider $\beta\to 0$ and expand everything to leading order in $\beta$.
For instance, the susceptibility matrix becomes

$\chi=\begin{pmatrix}\frac{\beta}{4}&-\frac{\beta^{2}\mu}{4}\\\
-\frac{\beta\mu}{4}&\frac{\beta^{2}(J^{2}+\mu^{2})}{4}\end{pmatrix},$ (102)

and hence $\gamma=\beta^{3}J^{2}/4$. In the case when we bias our chain with a
temperature gradient, we can approximate

$F_{\beta}^{>}(t)=-\frac{B_{2}(2Jt)}{2t}\nabla\beta+i\frac{\mu
B_{1}(2Jt)}{4Jt}\nabla\beta,$ (103)

where $B_{1,2}$ are Bessel functions of the first kind. Plugging this into
Eqs. (39, 61-63), we find

$\displaystyle\nabla Q$ $\displaystyle=\frac{\mu}{4}\nabla\beta,$ (104)
$\displaystyle\nabla E$ $\displaystyle=-\frac{J^{2}}{4}\nabla\beta,$ (105)
$\displaystyle j^{Q}$ $\displaystyle=-\frac{2\mu J_{1}^{2}}{3\pi
J}\nabla\beta,$ (106) $\displaystyle j^{E}$
$\displaystyle=\frac{8J_{1}^{2}J}{15\pi}\nabla\beta.$ (107)

Similarly, in the presence of a chemical potential gradient we can write

$F_{\mu}^{>}(t)=\frac{\beta^{3}\mu B_{2}(2Jt)}{4t}\nabla\mu+i\frac{\beta
B_{1}(2Jt)}{4Jt}\nabla\mu,$ (108)

and therefore deduce a new set of observables

$\displaystyle\nabla Q$ $\displaystyle=\frac{\beta}{4}\nabla\mu,$ (109)
$\displaystyle\nabla E$ $\displaystyle=\frac{\beta^{3}\mu J^{2}}{8}\nabla\mu,$
(110) $\displaystyle j^{Q}$ $\displaystyle=-\frac{2\beta J_{1}^{2}}{3\pi
J}\nabla\mu,$ (111) $\displaystyle j^{E}$ $\displaystyle=-\frac{4\beta^{3}\mu
J_{1}^{2}J}{15\pi}\nabla\mu.$ (112)

Together these form a set of four equations each for the diffusivity and
conductivity matrices. The solutions are given by

$\displaystyle D$ $\displaystyle=\begin{pmatrix}\frac{8J_{1}^{2}}{3\pi J}&0\\\
-\frac{8\mu J_{1}^{2}}{15\pi J}&\frac{32J_{1}^{2}}{15\pi J}\end{pmatrix},$
(113) $\displaystyle L$ $\displaystyle=\begin{pmatrix}\frac{2\beta
J_{1}^{2}}{3\pi J}&-\frac{2\beta^{2}\mu J_{1}^{2}}{3\pi J}\\\ -\frac{2\beta\mu
J_{1}^{2}}{3\pi
J}&\frac{8\beta^{2}J_{1}^{2}J}{15\pi}+\frac{2\beta^{2}\mu^{2}J_{1}^{2}}{3\pi
J}\end{pmatrix},$ (114)

in agreement with the results in Fig. 4 and Fig. 5. Note that
$D_{21}=\mu(D_{22}-D_{11})$. In addition, we can check that the generalized
Einstein relation $L=D\chi$ is indeed satisfied and that the thermal
conductivity takes the form

$\kappa=\frac{8\beta^{2}J_{1}^{2}J}{15\pi}=\frac{D_{22}\gamma}{\beta}.$ (115)

####  C.1.2 Zero temperature limit

Now take the opposite limit of $\beta\to\infty$ and expand in $1/\beta$. To
leading order, the susceptibility is given by

$\chi=\begin{pmatrix}\frac{1}{\pi
J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}&-\frac{\pi\mu}{12J^{3}\beta\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}\\\
-\frac{\pi\mu}{12J^{3}\beta^{2}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}&\frac{\pi}{3\beta
J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\end{pmatrix},$ (116)

and $\gamma=\frac{\pi}{3J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}$. The non-
equilibrium contribution from a temperature gradient is

$F_{\beta}^{>}(t)=\frac{\pi e^{-i\mu
t}}{3J\beta^{3}}\Bigg{(}\frac{i\mu}{4J^{2}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}-t\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\Bigg{)}\nabla\beta,$
(117)

leading to

$\displaystyle\nabla Q$ $\displaystyle=\frac{\pi\mu}{12(\beta
J)^{3}}\frac{1}{\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}\nabla\beta,$ (118)
$\displaystyle\nabla E$
$\displaystyle=-\frac{\pi}{3\beta^{3}J}\frac{1-2\left(\frac{\mu}{2J}\right)^{2}}{\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}\nabla\beta,$
(119) $\displaystyle j^{Q}$ $\displaystyle=-\frac{\pi\mu
J_{1}^{2}}{6\beta^{3}J^{4}}\nabla\beta,$ (120) $\displaystyle j^{E}$
$\displaystyle=\frac{\pi
J_{1}^{2}}{3\beta^{3}J^{2}}\left(1-3\left(\frac{\mu}{2J}\right)^{2}\right)\nabla\beta.$
(121)

Analogously, the contribution due to a chemical potential gradient can be
expanded as

$F_{\mu}^{>}(t)=i\frac{e^{-i\mu t}}{\pi
J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\nabla\mu,$ (122)

and hence our final observables are

$\displaystyle\nabla Q$ $\displaystyle=\frac{1}{\pi
J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\nabla\mu,$ (123)
$\displaystyle\nabla E$ $\displaystyle=\frac{\mu}{\pi
J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\nabla\mu,$ (124) $\displaystyle
j^{Q}$ $\displaystyle=-\frac{J_{1}^{2}}{\pi
J^{2}}\left(1-\left(\frac{\mu}{2J}\right)^{2}\right)\nabla\mu,$ (125)
$\displaystyle j^{E}$ $\displaystyle=-\frac{\mu J_{1}^{2}}{\pi
J^{2}}\left(1-\left(\frac{\mu}{2J}\right)^{2}\right)\nabla\mu.$ (126)

Solving for the transport coefficients, we conclude that

$\displaystyle D$
$\displaystyle=\begin{pmatrix}\frac{J_{1}^{2}}{J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}&-\frac{\mu
J_{1}^{2}}{4J^{3}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}\\\
-\frac{\pi^{2}\mu
J_{1}^{2}}{12\beta^{2}J^{3}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}}&\frac{J_{1}^{2}}{J}\sqrt{1-\left(\frac{\mu}{2J}\right)^{2}}\end{pmatrix},$
(127) $\displaystyle L$ $\displaystyle=\begin{pmatrix}\frac{J_{1}^{2}}{\pi
J^{2}}\left(1-\left(\frac{\mu}{2J}\right)^{2}\right)&-\frac{\pi\mu
J_{1}^{2}}{6\beta J^{4}}\\\ -\frac{\pi\mu
J_{1}^{2}}{6\beta^{2}J^{4}}&\frac{\pi J_{1}^{2}}{3\beta
J^{2}}\left(1-\left(\frac{\mu}{2J}\right)^{2}\right)\end{pmatrix},$ (128)

which agrees with our numerics in Fig. 4 and Fig. 5. Again, we find that
$D_{21}=\mu(D_{22}-D_{11})=0$ to first order in $1/\beta$. The non-zero
contribution above is a higher-order correction necessary to satisfy the
Einstein relation. Eq. (80) holds as expected with
$\kappa=\bar{\kappa}=D_{22}\gamma/\beta$. Our answer agrees with the free-
fermion calculation in Ref. [33] that found $\sigma=1/\pi$ and $\kappa=\pi
T/3$ in units of $J\approx J_{1}=1$ and $\mu=0$, thus providing an independent
consistency check of our methods.

###  C.2 Large $q$ limit

We now return to the large $q$ analysis of Sec. III.7 and discuss the
thermodynamic properties of our model in this limit. The grand canonical
potential has been previously derived in Ref. [21]

$\begin{split}\Omega&=-\frac{1}{\beta}\log\left(2\cosh(\beta\mu/2)\right)\\\
&-\frac{2\pi v}{\beta q^{2}\cosh^{2}(\beta\mu/2)}\left(\tan(\frac{\pi
v}{2})-\frac{\pi v}{4}\right).\end{split}$ (129)

To leading order in $1/q$, the charge becomes

$Q=\frac{1}{2}\tanh(\frac{\beta\mu}{2})+\mathcal{O}\left(\frac{1}{q^{2}}\right).$
(130)

It follows that at low temperatures, the chemical potential should scale as
$\mu\sim T$ to maintain a constant charge. This stems from a failure of the
infinite $q$ and infinite $\beta$ limits to commute, which is an inherent
shortcoming of this expansion [96].

Taking a second derivative, we find the susceptibilities

$\displaystyle\chi_{11}$
$\displaystyle=\frac{\beta}{4\cosh^{2}(\beta\mu/2)}+\mathcal{O}\left(\frac{1}{q^{2}}\right),$
(131) $\displaystyle\chi_{12}$
$\displaystyle=-\frac{\beta^{2}\mu}{4\cosh^{2}(\beta\mu/2)}+\mathcal{O}\left(\frac{1}{q^{2}}\right),$
(132) $\displaystyle\chi_{22}$
$\displaystyle=\frac{(\beta\mu)^{2}}{4\cosh^{2}(\beta\mu/2)}\left(1+\left(\frac{2\pi
v}{q\beta\mu}\right)^{2}\frac{1}{\frac{\pi v}{2}\tan(\frac{\pi
v}{2})+1}\right),$ (133)

where the second term in $\chi_{22}$ is necessary to obtain the leading order
contribution to

$\gamma=\frac{\beta}{q^{2}\cosh^{2}(\beta\mu/2)}\frac{(\pi v)^{2}}{\frac{\pi
v}{2}\tan\left(\frac{\pi v}{2}\right)+1}.$ (134)

The prefactors can also be written in terms of the charge
$\cosh^{2}(\beta\mu/2)=(1-4Q^{2})^{-1}$. Note that in the zero-temperature
limit, at fixed charge, the compressibility $\chi_{11}$ diverges with $\beta$,
which is unphysical. The correct behavior can be recovered by keeping the next
order term in the large $q$ expansion and taking the temperature to zero first
[21]. This is another example where the order in which we take the limits
matters.

In the case of transport, we first imagine maintaining a constant chemical
potential gradient across the chain held at a fixed temperature. The non-
equilibrium contribution defined in Eq. (71) can be computed via the chain
rule

$\begin{split}f_{\mu}(t)&=-(q-2)\frac{\beta^{2}\mathcal{J}}{4}\tanh(\frac{\beta\mu}{2})\frac{\cos(\frac{\pi
v}{2})}{\frac{\pi v}{2}\tan(\frac{\pi v}{2})+1}\\\
&\cdot\left(\left(1-\frac{2it}{\beta}\right)\tan(\frac{\pi v}{2}-\frac{i\pi
vt}{\beta})-\tan(\frac{\pi v}{2})\right).\end{split}$ (135)

The energy gradient is obtained by direct integration

$\nabla
E=\frac{\beta\mathcal{J}\tanh(\frac{\beta\mu}{2})}{2q\cosh^{2}(\frac{\beta\mu}{2})}\Bigg{(}\frac{\frac{\pi
v}{2}\cos(\frac{\pi v}{2})}{\frac{\pi v}{2}\tan(\frac{\pi
v}{2})+1}+\sin\left(\frac{\pi v}{2}\right)\Bigg{)}\nabla\mu.$ (136)

However, the energy current is a bit more subtle. It turns out that the
overall contribution from the terms proportional $f_{\mu}$ is zero, so we only
have to consider the last term in Eq. (76). We use the identity

$\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\,\mathrm{Im}\left[e^{g(t)+g(t^{\prime})}\right]=-\frac{\sin(\pi
v)}{2\mathcal{J}^{2}},$ (137)

and finally arrive at

$j^{E}=-\frac{\beta\mathcal{J}_{1}^{2}\tanh(\frac{\beta\mu}{2})}{4q\cosh^{2}(\frac{\beta\mu}{2})}\sin(\pi
v)\nabla\mu.$ (138)

The other two observables related to charge transport are given by Eqs. (72,
74).

Next, we will consider a slightly simpler setup where we impose both a
temperature and chemical potential gradient, but maintain a constant charge
$\nabla Q=0$. From Eq. (130) this is equivalent to holding $\beta\mu$ constant
along the chain and setting $\beta\nabla\mu=-\mu\nabla\beta$. In linear
response, this corresponds to an additive contribution from both $f_{\mu}$ and
$f_{\beta}$, which we denote by

$f_{Q}(t)=\frac{(1+i\mathcal{J}t\sin(\frac{\pi v}{2}))\tan(\frac{\pi
v}{2}-\frac{i\pi vt}{\beta})-\tan(\frac{\pi v}{2})}{\frac{\beta}{\pi
v}\left(\frac{\pi v}{2}\tan(\frac{\pi v}{2})+1\right)}\nabla\beta$ (139)

This is exactly the answer we found for a Majorana chain [7], up to a constant
prefactor. Therefore, we can follow the same calculations to find

$\displaystyle\nabla E$
$\displaystyle=-\frac{\mathcal{J}}{q^{2}\beta\cosh^{2}(\beta\mu/2)}\frac{\pi
v\cos(\frac{\pi v}{2})}{\frac{\pi v}{2}\tan(\frac{\pi v}{2})+1}\nabla\beta,$
(140) $\displaystyle j^{E}$
$\displaystyle=\frac{\mathcal{J}_{1}^{2}}{q^{2}\beta\cosh^{2}(\beta\mu/2)}\frac{\pi
v}{3}\cos^{2}\left(\frac{\pi v}{2}\right)\nabla\beta.$ (141)

Note that the charge current vanishes as expected, since
$j^{Q}\sim(C_{\mu}+C_{\beta})=0$ for this setup.

Finally, combining all the results into a system of equations, we deduce the
diagonal diffusivity entries

$\displaystyle D_{11}$
$\displaystyle=\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}\cos(\frac{\pi v}{2}),$
(142) $\displaystyle D_{22}$
$\displaystyle=\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\left(\frac{\pi
v}{2}\sin(\frac{\pi v}{2})+\cos(\frac{\pi v}{2})\right),$ (143)

and the off-diagonal values $D_{12}=0$ and $D_{21}=\mu(D_{22}-D_{11})$.
Similarly, the conductivities are

$\displaystyle\sigma$
$\displaystyle=\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}^{2}}\frac{\pi
v}{4\cosh^{2}(\beta\mu/2)},$ (144) $\displaystyle\alpha$
$\displaystyle=-\beta\mu\sigma=-\beta\mu\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}^{2}}\frac{\pi
v}{4\cosh^{2}(\beta\mu/2)},$ (145) $\displaystyle\bar{\kappa}$
$\displaystyle=\beta\mu^{2}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}^{2}}\frac{\pi
v}{4\cosh^{2}(\beta\mu/2)}\left(1+\frac{4}{3}\left(\frac{\pi
v}{q\beta\mu}\right)^{2}\right),$ (146)

and we can check the dependence in Eq. (80)

$\kappa=\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}^{2}}\frac{(\pi v)^{3}}{3\beta
q^{2}\cosh^{2}(\beta\mu/2)}=\frac{D_{22}\gamma}{\beta}.$ (147)

The large $q$ approximation has given us remarkably simple closed-form answers
for all the transport coefficients. The generalized Einstein relation now can
be checked explicitly. Notice that the conductivities, just like the
susceptibilites, are suppressed by a factor of $\cosh^{2}(\beta\mu/2)$. This
again suggests that we should scale our parameters to maintain a finite value
of $\beta\mu$. In our numerics at finite $q$, we maintain both $J$ and $\mu$
constant while sweeping a wide range of temperatures. Therefore, in order to
make the results of the large $q$ approximation consistent in this regime, we
will restrict ourselves to the case $\mu=0$ when comparing to fixed $q$
results.

##  Appendix D Thermopower

Figure 10: Temperature dependence of the thermopower (solid lines) for the SYK
model with $q=4$. At low temperatures in the SYK-like phase, $\Theta$
approaches the conformal answer $2\pi\mathcal{E}$ (dashed lines). The inset
shows the thermoelectric figure of merit $ZT$ for the same parameters.

If a temperature gradient is applied to across a material with free charge
carriers, a potential gradient will arise as a result of the carriers’ motion
from hot to cold areas. The magnitude of this thermoelectric effect is
characterized by the thermopower $\Theta$, also known as the Seebeck
coefficient. The thermopower is defined as the ratio of the induced potential
gradient $\nabla\mu$ to the applied temperature gradient $\nabla T$ after the
system has reached a steady state with no charge current [86, 85, 87]. Eq.
(45) implies

$\Theta=-\frac{\nabla\mu}{\nabla T}=\frac{\alpha}{\sigma}.$ (148)

Materials with high thermopower are very important for building efficient
thermoelectric generators and coolers. A useful metric for quantifying the
effectiveness of a thermoelectric material for practical applications is the
dimensionless thermoelectric figure of merit

$ZT=\frac{\sigma\Theta^{2}T}{\kappa}=\frac{\alpha^{2}T}{\sigma\kappa}=\frac{\bar{\kappa}}{\kappa}-1.$
(149)

For conventional metals and insulators, $ZT$ is at most of order one [106].
Recently, it has been shown that the thermopower of Dirac and Weyl semimetals
in an external magnetic field grows linearly with the magnitude of the field
and can reach extremely high values [106, 107, 108]. This makes it possible to
have $ZT\gg 1$ in these materials.

We now calculate the thermopower and thermoelectric figure of merit for the
complex SYK model using our previous results for the conductivity matrix. We
showcase our findings for $q=4$ in Fig. 10. The thermopower reaches a non-zero
constant at low temperatures and increases with $\mu$ as we approach the phase
transition. The thermoelectric figure of merit is of order one and has a
similar dependence on $\beta$ and $\mu$. This is consistent with a previous
analysis of thermopower in SYK models and holographic theories [21, 105],
where the authors showed that in the low temperature limit

$\Theta=2\pi\mathcal{E},$ (150)

where $\mathcal{E}$ is the particle-hole asymmetry of the fermionic spectral
function [9, 18, 21], which can be expressed in terms of charge via a
Luttinger-Ward identity

$\displaystyle e^{2\pi\mathcal{E}}$
$\displaystyle=\frac{\sin(\pi/q+\theta)}{\sin(\pi/q-\theta)},$ (151)
$\displaystyle Q$
$\displaystyle=-\frac{\theta}{\pi}-\left(\frac{1}{2}-\frac{1}{q}\right)\frac{\sin(2\theta)}{\sin(2\pi/q)}.$
(152)

Note that with our conventions, both $\mathcal{E}$ and $\theta$ are negative.
By combining this result with the Wiedemann-Franz ratio in Eq. (82), we find
that the figure of merit should scale as $ZT=3q^{2}\mathcal{E}^{2}$ at low
temperatures. The asymmetry diverges as the charge becomes larger when leaving
the SYK-like phase, which explains the scaling in Fig. 10.

The non-vanishing thermopower in the limit of zero temperature is a remarkable
feature of the SYK model. It can be related to the existence of a finite zero-
temperature entropy $S_{0}$ via an exact Kelvin formula [9, 21]. Recently, it
has been suggested that thermopower measurements can serve as a direct probe
of the residual entropy $S_{0}$ [105]. This low-temperature behavior is in
stark contrast to that of a Fermi liquid, whose thermopower vanishes linearly
with $T$. In fact, we can see this explicitly for the $q=2$ SYK, where Eq.
(128) implies $\Theta\sim T$ and $ZT\sim T^{2}$.

The thermoelectric parameters can be computed exactly in the large $q$ limit.
We find the simple expressions $\Theta=-\beta\mu$ and
$ZT=\frac{3q^{2}(\beta\mu)^{2}}{4(\pi v)^{2}}$. The thermopower is directly
related to $Q$ through Eq. (130) and is constant at fixed charge. At high
temperatures, $\pi v\approx\beta\mathcal{J}$ and the figure of merit remains
finite with $ZT\sim q^{2}\mu^{2}$. This precisely matches our results in the
inset of Fig. 10. Given this scaling of $ZT$, it seems possible to further
increase the thermodynamic figure of merit by going to higher $q$ SYK models.