# Energy Transport in Sachdev-Ye-Kitaev Networks Coupled to Thermal Baths

Cristian Zanoci<EMAIL_ADDRESS>Department of Physics, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts
02139, USA  Brian Swingle<EMAIL_ADDRESS>Condensed Matter Theory Center and
Joint Quantum Institute, Department of Physics, University of Maryland,
College Park, Maryland 20742, USA  Department of Physics, Brandeis University,
Waltham, Massachusetts 02453, USA

###### Abstract

We develop a general framework for studying the equilibrium and non-
equilibrium properties of arbitrary networks of Sachdev-Ye-Kitaev clusters
coupled to thermal baths. We proceed to apply this technique to the problem of
energy transport, which is known to be diffusive due to the strange metal
behavior of these models. We use the external baths to impose a temperature
gradient in the system and study the emerging non-equilibrium steady state
using the Schwinger-Keldysh formalism. We consider two different
configurations for the baths, implementing either a boundary or bulk driving,
and show that the latter leads to a significantly faster convergence to the
steady state. This setup allows us to compute both the temperature and
frequency dependence of the diffusion constant. At low temperatures, our
results agree perfectly with the previously known values for diffusivity in
the conformal limit. We also establish a relationship between energy transport
and quantum chaos by showing that the diffusion constant is upper bounded by
the chaos propagation rate at all temperatures. Moreover, we find a simple
analytical form for the non-equilibrium Green’s functions in the linear
response regime and use it to derive exact closed-form expressions for the
diffusion constant in various limits. We mostly focus on uniform one-
dimensional chains, but we also discuss higher-dimensional generalizations.

##  I Introduction

The subject of non-equilibrium dynamics spans across multiple branches of
physics, ranging from condensed matter to quantum gravity. Recently, it has
drawn a lot of attention in an attempt to answer fundamental questions
regarding thermalization, quantum many-body chaos, and transport in strongly
interacting systems. The latter is particularly useful for unveiling the
properties of quantum matter and exploring the dynamical processes governing
the behavior of quantum systems out of equilibrium. Despite recent efforts,
practical calculations of transport coefficients in quantum many-body systems
remain challenging from both a conceptual and technical standpoint [1]. On one
hand, a general hydrodynamic description of transport is possible at a
macroscopic level, but relating the parameters of this theory to the
microscopic parameters of the underlying model is difficult except in special
cases, for example, at weak coupling. On the other hand, numerical methods for
simulating non-equilibrium dynamics can directly address a wide class of
microscopic models, but their applicability is typically limited to small
systems and short time scales. Therefore, it would be very useful to develop a
framework that allows us to efficiently extract the transport properties of
small systems, while simultaneously offering insights about the non-
equilibrium dynamics in the hydrodynamic limit.

In the context of one-dimensional models, open-system dynamics has shown
promising results for bridging the gap between macroscopic effective theories
and microscopic models, especially when combined with tensor network
techniques [1, 2, 3]. Within this framework, the system is coupled to external
baths that drive the system towards a desired steady state, from which various
transport properties are readily available. However, for a generic non-
integrable model, these simulations can suffer from entanglement growth and
very slow convergence at low temperatures [4]. Moreover, very little is known
about the structure of the emergent non-equilibrium states. Our goal is to
study a solvable open-system model and show how one can mitigate the
aforementioned convergence problems, while also developing an analytical
understanding of the non-equilibrium steady state (NESS) and elucidating the
relationship between transport and other many-body phenomena, such as quantum
chaos.

Looking more broadly, we expect the results about the structure of steady
states and convergence rates in the solvable model to generalize to a variety
of more physical models. Combined with tensor network methods, these lessons
may significantly improve our ability to compute transport properties from
microscopic theories. There are also a growing number of principled numerical
techniques [5, 6, 7, 8, 9, 10], besides the open-system approach, which might
also benefit from additional information about the structure of current
carrying states.

Our solvable model is built from the Sachdev-Ye-Kitaev (SYK) model [11, 12,
13, 14, 15, 16, 17, 18, 19, 20], which in recent years has emerged as a
paradigmatic example of a strongly coupled, yet exactly solvable quantum many-
body system. The model describes $N$ fermions with random all-to-all $q$-body
interactions and displays a multitude of remarkable properties. The
Hamiltonian belongs to a class of systems realizing holographic quantum matter
without quasiparticle excitations [14]. In the limit of large $N$ and low
temperatures, the system has an emergent approximate time reparameterization
conformal symmetry [11, 17]. The SYK model also exhibits many-body chaos [11,
17] and saturates a universal bound on the relevant Lyapunov exponent [21], a
feature shared by black holes in Einstein gravity. In fact, the SYK model is
holographically dual to gravitational theories of black holes with near-
horizon AdS2 geometry [22, 23, 24, 25, 11, 26, 27, 18, 28]. On the gravity
side, this connection can be used to study unsolved questions related to
holography and black holes, such as black hole evaporation [29, 30, 31, 28].
On the condensed matter side, the SYK model presents a valuable platform for
studying non-Fermi liquid behavior [12, 13, 14] and thermalization [32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42], among other phenomena.

Subsequent works [43, 44, 45, 46, 47, 48] extended the SYK model to higher
spatial dimensions by coupling individual SYK clusters to their neighbors with
similar SYK interactions, thus building lattices that exhibit strange metal
behavior. The generalized models remain exactly solvable via a saddle point
expansion and retain many features of the original SYK, including local
criticality and maximal chaos [43, 44]. The spatial locality of the models
allows for the investigation of transport and chaos in these systems. This
revealed many indicative properties of a strange metal, such as diffusive
propagation of energy [43, 44, 45] and resistivity that scales linearly with
temperature [45]. Moreover, it was shown that the same time-reparametrization
field is responsible for the propagation of both low-energy modes and quantum
chaos in this system [43], thus leading to a connection between the diffusion
constant and the butterfly velocity [49, 50, 51, 52, 53, 54, 55].

Alongside these theoretical developments, several experimental proposals for
the SYK model have been put forward [56, 57], including realizations in
ultracold atoms [58], Majorana modes at the interface of a topological
insulator and superconductor [59], semiconductor wires coupled through a
disordered quantum dot [60], and graphene flakes with irregular boundaries
[61]. These experiments open up the possibility of directly studying the
transport properties of the model.

In this paper, we study energy transport in SYK networks driven out of
equilibrium by coupling to external baths. The system-reservoir interaction in
SYK models has been previously discussed in the context of thermalization [40,
39, 41], and only recently in the context of transport [62]. Most of the
previous approaches to transport used the explicit form of the fermion four-
point function and relied on the low-temperature or large $q$ limits to
compute the diffusion constant [43, 55, 48]. One advantage of coupling the
system to baths is that it allows us to study transport for arbitrary values
of $q$ and inverse temperature $\beta$. We consider two configurations with
baths attached either on the boundary or throughout the bulk, and show that
the latter leads to faster convergence. We study the resulting non-equilibrium
dynamics using the Schwinger-Keldysh formalism [63, 64] and derive the
associated Kadanoff-Baym equations governing the approach to a steady state.
Once we reach the NESS, we compute the current and energy gradient across the
system to determine the diffusion coefficient.

We rely on this setup to investigate both DC and AC transport. For the former,
the baths are used to impose a constant temperature imbalance on the system,
while for the latter, we oscillate the temperature of the baths to create a
time-dependent drive. For the DC case, we find that the diffusivity increases
with $\beta$ until it saturates at a constant value matching exactly the
conformal answer [43]. We also show that the diffusion constant is related to
chaos via the inequality $D\leq v_{B}^{2}/\lambda_{L}$, which holds at all
temperatures and becomes an equality in the conformal limit. This behavior is
consistent with studies of related systems [55, 65, 53, 54]. For the AC case,
we compute the frequency dependence of the diffusivity and find that it
decreases exponentially at high frequencies.

Additionally, our results indicate that the Green’s functions in NESS are only
slightly perturbed from their equilibrium values, allowing us to re-write the
equations of motion in terms of these new non-equilibrium corrections without
referencing the baths. In the case of DC transport, we were able to find this
correction explicitly in terms of the equilibrium Green’s function and derived
analytic formulas for the diffusion constant in the limit of large $q$, as
well as for $q=2$. In the case of AC transport, we found that the NESS
contribution has a more complicated form which can only be estimated
numerically. Our results computed via this formalism agree well with the ones
obtained from coupling the system to baths and represent one of the first
attempts at characterizing the NESS analytically.

The rest of the paper is structured as follows. In Sec. II we define our
higher-dimensional generalization of the SYK model coupled to baths. We
describe in detail our approach to non-equilibrium dynamics and show how to
compute the different quantities related to transport. In Sec. III we briefly
review the equilibrium properties of our system. Next, in Sec. IV we present
our results for DC transport on different lattices, while focusing
specifically on the temperature dependence of the diffusion coefficients and
the structure of the emerging NESS. Similarly, Sec. V is devoted to frequency-
dependent transport. Finally, we give a short discussion and outlook in Sec.
VI.

##  II Setup

We begin by introducing our family of models, the equations governing their
equilibrium and non-equilibrium dynamics, the different approximations which
significantly simplify these equations, and the general setup and relevant
quantities used to study the transport properties of these models.

###  II.1 SYK models coupled to thermal baths

The SYK model [11, 12, 13, 15, 16, 14, 17, 18, 19, 20] is a strongly
interacting fermion model in $(0+1)$ dimensions. We consider a generalization
of this model on arbitrary graphs $\mathcal{G}=(V,E)$. Each vertex $u\in V$ is
an SYK cluster (quantum dot) of $N$ Majorana fermions $\psi_{i}^{u}$ with
random all-to-all $q_{S}-$body interactions given by the Hamiltonian

$H_{0}^{u}=i^{q_{S}/2}\sum_{\\{i\\}}J_{i_{1}\ldots
i_{q_{S}}}^{(0)}\psi_{i_{1}}^{u}\ldots\psi_{i_{q_{S}}}^{u},$ (1)

where $\\{i\\}$ denotes the restricted sum over $1\leq
i_{1}<\cdots<i_{q_{S}}\leq N$. The vertices of the graph connected by an edge
$(u,v)\in E$ are coupled to each other via the interaction Hamiltonian

$H_{1}^{uv}=i^{q_{I}/2}\sum_{\\{i\\},\\{j\\}}J_{i_{1}\ldots
i_{\frac{q_{I}}{2}}j_{1}\ldots
j_{\frac{q_{I}}{2}}}^{(1)}\psi_{i_{1}}^{u}\ldots\psi_{i_{\frac{q_{I}}{2}}}^{u}\psi_{j_{1}}^{v}\ldots\psi_{j_{\frac{q_{I}}{2}}}^{v}.$
(2)

The Majorana fermions obey the standard anti-commutation relations
$\\{\psi_{i}^{u},\psi_{j}^{v}\\}=\delta_{ij}\delta_{uv}$. Note that it is not
essential that the interaction term contains the same number of fermions from
both sites and a generalization to more generic interactions should be
straightforward. Similar higher-dimensional SYK models have been previously
studied in the context of transport [43, 55, 44, 45, 46, 47, 66, 67], quantum
chaos [68, 69, 70], and quantum phase transitions [71, 72, 48, 73, 74].

In order to study the non-equilibrium properties of our system, we couple a
subset of vertices to thermal baths [69, 40, 39, 62, 41]. Each bath is modeled
as an SYK cluster of $M$ Majorana fermions $\chi_{i}^{u}$ with all-to-all
$q_{B}-$body interactions given by the Hamiltonian

$H_{B}^{u}=i^{q_{B}/2}\sum_{\\{i\\}}J_{i_{1}\ldots
i_{q_{B}}}^{(B)}\chi_{i_{1}}^{u}\ldots\chi_{i_{q_{B}}}^{u}.$ (3)

The system-bath interaction is of the same form as the inter-cluster coupling
on the graph

$H_{SB}^{u}=i^{q_{I}/2}\sum_{\\{i\\},\\{j\\}}V_{i_{1}\ldots
i_{\frac{q_{I}}{2}}j_{1}\ldots
j_{\frac{q_{I}}{2}}}^{(u)}\psi_{i_{1}}^{u}\ldots\psi_{i_{\frac{q_{I}}{2}}}^{u}\chi_{j_{1}}^{u}\ldots\chi_{j_{\frac{q_{I}}{2}}}^{u}.$
(4)

Since not all vertices are necessarily coupled to a bath, we use $V^{(u)}\neq
0$ to indicate the presence of a bath and set $V^{(u)}=0$ otherwise. The baths
are coupled to the system at time $t=0$.

All the SYK couplings are independent Gaussian random variables with zero mean
and variances

$\displaystyle\langle(J_{i_{1}\ldots i_{q_{S}}}^{(0)})^{2}\rangle$
$\displaystyle=\frac{J_{0}^{2}(q_{S}-1)!}{N^{q_{S}-1}},$ (5)
$\displaystyle\langle(J_{i_{1}\ldots i_{q_{I}/2}j_{1}\ldots
j_{q_{I}/2}}^{(1)})^{2}\rangle$
$\displaystyle=\frac{J_{1}^{2}(q_{I}/2)!(q_{I}/2-1)!}{N^{q_{I}-1}},$ (6)
$\displaystyle\langle(J_{i_{1}\ldots i_{q_{B}}}^{(B)})^{2}\rangle$
$\displaystyle=\frac{J_{B}^{2}(q_{B}-1)!}{M^{q_{B}-1}},$ (7)
$\displaystyle\langle(V_{i_{1}\ldots i_{q_{I}/2}j_{1}\ldots
j_{q_{I}/2}}^{(u)})^{2}\rangle$
$\displaystyle=\frac{V_{u}^{2}(q_{I}/2)!(q_{I}/2-1)!}{M^{q_{I}/2}N^{q_{I}/2-1}}.$
(8)

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, such as the energy.

The total Hamiltonian is simply the sum of all the terms

$H=\sum_{u\in V}\big{(}H_{0}^{u}+H_{B}^{u}+H_{SB}^{u}\big{)}+\sum_{(u,v)\in
E}H_{1}^{uv}.$ (9)

Our model represents a unified framework for studying both equilibrium and
non-equilibrium properties of arbitrary SYK lattices coupled to baths. We can
reproduce many of the previously studied SYK configurations [43, 44, 45, 46,
47, 66, 67, 68, 70, 71, 72, 48, 73, 74, 69, 40, 39, 62, 41, 42] by adjusting
the topology of the graph, the type of intra- and inter-cluster interactions,
or the nature of the baths in our model.

###  II.2 Equilibrium

The generalized SYK model maintains all the exactly solvable properties of the
original model in the large-$N$ limit [17]. Strictly speaking, in the presence
of quenched disorder, we would have to introduce replicas in the path
integral. However, in the large-$N$ limit, our model self-averages and the
interaction between replicas is suppressed [17, 43]. Therefore, it is
sufficient to consider the replica-diagonal partition function [17, 43], for
which the Euclidean effective action can be written as

$\displaystyle S=\sum_{(u,v)\in E}S_{uv}+\sum_{u\in V}S_{u},$ (10)
$\displaystyle
S_{uv}=-\frac{NJ_{1}^{2}}{q_{I}}\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}G_{u}^{S}(\tau_{1},\tau_{2})^{q_{I}/2}G_{v}^{S}(\tau_{1},\tau_{2})^{q_{I}/2},$
(11)
$\begin{split}S_{u}&=-\frac{N}{2}\log\det\left(\partial_{\tau}-\Sigma_{u}^{S}(\tau_{1},\tau_{2})\right)-\frac{M}{2}\log\det\left(\partial_{\tau}-\Sigma_{u}^{B}(\tau_{1},\tau_{2})\right)-\frac{NV_{u}^{2}}{q_{I}}\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}G_{u}^{S}(\tau_{1},\tau_{2})^{q_{I}/2}G_{u}^{B}(\tau_{1},\tau_{2})^{q_{I}/2}\\\
&+\frac{N}{2}\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}\left(\Sigma_{u}^{S}(\tau_{1},\tau_{2})G_{u}^{S}(\tau_{1},\tau_{2})-\frac{J_{0}^{2}}{q_{S}}G_{u}^{S}(\tau_{1},\tau_{2})^{q_{S}}\right)+\frac{M}{2}\int\mathop{}\\!\mathrm{d}\tau_{1}\mathop{}\\!\mathrm{d}\tau_{2}\left(\Sigma_{u}^{B}(\tau_{1},\tau_{2})G_{u}^{B}(\tau_{1},\tau_{2})-\frac{J_{B}^{2}}{q_{B}}G_{u}^{B}(\tau_{1},\tau_{2})^{q_{B}}\right).\end{split}$
(12)

For each vertex $u$, we introduced the Euclidean time-ordered fermion two-
point functions

$\displaystyle G_{u}^{S}(\tau_{1},\tau_{2})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\langle
T\psi_{i}(\tau_{1})\psi_{i}(\tau_{2})\rangle,$ (13) $\displaystyle
G_{u}^{B}(\tau_{1},\tau_{2})$ $\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\langle
T\chi_{i}(\tau_{1})\chi_{i}(\tau_{2})\rangle,$ (14)

and the fermion self-energies $\Sigma_{u}^{S,B}$ as the associated Lagrange
multipliers. In the large-$N$ limit, the saddle point of this effective action
produces the Schwinger-Dyson (SD) equations of motion

$\displaystyle 1$
$\displaystyle=(\partial_{\tau}-\Sigma_{u}^{S,B}(\tau))G_{u}^{S,B}(\tau),$
(15) $\displaystyle\Sigma_{u}^{S}(\tau)$
$\displaystyle=J_{0}^{2}G_{u}^{S}(\tau)^{q_{S}-1}+V_{u}^{2}G_{u}^{S}(\tau)^{q_{I}/2-1}G_{u}^{B}(\tau)^{q_{I}/2}$
$\displaystyle+J_{1}^{2}\sum_{v\in
V}A_{uv}G_{u}^{S}(\tau)^{q_{I}/2-1}G_{v}^{S}(\tau)^{q_{I}/2},$ (16)
$\displaystyle\Sigma_{u}^{B}(\tau)$
$\displaystyle=J_{B}^{2}G_{u}^{B}(\tau)^{q_{B}-1}+\frac{N}{M}V_{u}^{2}G_{u}^{S}(\tau)^{q_{I}/2}G_{u}^{B}(\tau)^{q_{I}/2-1},$
(17)

where $A_{uv}$ is the adjacency matrix of the graph and we assumed time-
translation symmetry $\tau=\tau_{1}-\tau_{2}$ in equilibrium.

If the bath is much larger than the system $M\gg N$, we can neglect the back-
reaction on the bath and drop the second term in Eq. (17). Hence, even though
the smaller system is affected by its coupling to the bath through the term
proportional to $V_{u}^{2}$, the larger bath can be approximated as decoupled
[39, 40]. As a result, the bath’s Green’s function always matches its
equilibrium value, even when coupled to the system. We will operate under this
assumption in the remainder of the paper.

###  II.3 Non-equilibrium

The real-time evolution of a quantum many-body system coupled to baths can be
described using the Schwinger-Keldysh formalism [63, 64]. Following the
derivation in Refs. [37, 38], we can write down the terms in the Lorentzian
action

$S_{uv}=\frac{iNJ_{1}^{2}}{q_{I}}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}G_{u}^{S}(t_{2},t_{1})^{q_{I}/2}G_{v}^{S}(t_{1},t_{2})^{q_{I}/2},$
(18) $\displaystyle S_{u}$
$\displaystyle=-\frac{iN}{2}\log\det\left(\partial_{t}\delta_{\mathcal{C}}(t_{1},t_{2})+i\Sigma_{u}^{S}(t_{1},t_{2})\right)+\frac{iN}{2}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}\left(\frac{J_{0}^{2}}{q_{S}}G_{u}^{S}(t_{2},t_{1})^{q_{S}/2}G_{u}^{S}(t_{1},t_{2})^{q_{S}/2}-\Sigma_{u}^{S}(t_{1},t_{2})G_{u}^{S}(t_{2},t_{1})\right)$
$\displaystyle-\frac{iM}{2}\log\det\left(\partial_{t}\delta_{\mathcal{C}}(t_{1},t_{2})+i\Sigma_{u}^{B}(t_{1},t_{2})\right)+\frac{iM}{2}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}\left(\frac{J_{B}^{2}}{q_{B}}G_{u}^{B}(t_{2},t_{1})^{q_{B}/2}G_{u}^{B}(t_{1},t_{2})^{q_{B}/2}-\Sigma_{u}^{B}(t_{1},t_{2})G_{u}^{B}(t_{2},t_{1})\right)$
$\displaystyle+\frac{iNV_{u}^{2}}{q_{I}}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}G_{u}^{S}(t_{2},t_{1})^{q_{I}/2}G_{u}^{B}(t_{1},t_{2})^{q_{I}/2},$
(19)

where $\mathcal{C}=\mathcal{C}^{+}\cup\mathcal{C}^{-}$ is the closed-time
Keldysh contour with the positive (forward) branch $\mathcal{C}^{+}$ going
from $-\infty$ to $+\infty$ and the negative (backward) branch
$\mathcal{C}^{-}$ going from $+\infty$ to $-\infty$. Recall that in the
Schwinger-Keldysh formalism, the contour-ordered Green’s functions are
actually $2\times 2$ matrices

$G_{u}^{S}(t_{1},t_{2})=\begin{pmatrix}G_{u}^{S}(t_{1}^{+},t_{2}^{+})&G_{u}^{S}(t_{1}^{+},t_{2}^{-})\\\
G_{u}^{S}(t_{1}^{-},t_{2}^{+})&G_{u}^{S}(t_{1}^{-},t_{2}^{-})\end{pmatrix},$
(20)

where $t_{i}^{+}$ lives on the positive branch and $t_{i}^{-}$ lives on the
negative branch. We will be mainly interested in the greater and lesser
Green’s functions

$\displaystyle G_{u}^{>}(t_{1},t_{2})\equiv G_{u}^{S}(t_{1}^{-},t_{2}^{+})$
$\displaystyle=-\frac{i}{N}\sum_{i=1}^{N}\langle\psi_{i}(t_{1}^{-})\psi_{i}(t_{2}^{+})\rangle,$
(21) $\displaystyle G_{u}^{<}(t_{1},t_{2})\equiv
G_{u}^{S}(t_{1}^{+},t_{2}^{-})$
$\displaystyle=-\frac{i}{N}\sum_{i=1}^{N}\langle\psi_{i}(t_{1}^{+})\psi_{i}(t_{2}^{-})\rangle,$
(22)

through which we can define the other two components

$\displaystyle G_{u}^{S}(t_{1}^{+},t_{2}^{+})$
$\displaystyle=\theta(t_{1}-t_{2})G_{u}^{>}(t_{1},t_{2})+\theta(t_{2}-t_{1})G_{u}^{<}(t_{1},t_{2}),$
$\displaystyle G_{u}^{S}(t_{1}^{-},t_{2}^{-})$
$\displaystyle=\theta(t_{1}-t_{2})G_{u}^{<}(t_{1},t_{2})+\theta(t_{2}-t_{1})G_{u}^{>}(t_{1},t_{2}).$
(23)

In practice, it is often more convenient to work with the retarded, advanced,
and Keldysh Green’s functions

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

which are related to the contour-ordered Green’s functions above via a Keldysh
rotation [63, 64]. The corresponding self-energies are defined in a similar
manner. For Majorana fermions in thermal equilibrium, the greater and lesser
Green’s functions are related via

$G_{u}^{>}(t_{1},t_{2})=-G_{u}^{<}(t_{2},t_{1})=-G_{u}^{>}(t_{2},t_{1})^{*},$
(27)

which holds true even out-of-equilibrium as long as the time evolution starts
from thermal equilibrium [75]. Therefore, all other two-point functions can be
computed from $G_{u}^{>}(t_{1},t_{2})$ and we can focus solely on studying
this component of the Green’s function. Furthermore, due to the Majorana anti-
commutation relations, the greater Green’s function at equal time reduces to
$G_{u}^{>}(t,t)=-i/2$.

To obtain the Schwinger-Dyson equations in real time, we can use the large-$N$
saddle point equations for the Lorentzian action [37, 38], yielding

$\displaystyle\Sigma_{u}^{S,B}(t_{1},t_{2})$
$\displaystyle=\partial_{t}\delta_{\mathcal{C}}(t_{1},t_{2})-(G_{u}^{S,B}(t_{1},t_{2}))^{-1},$
(28) $\displaystyle\Sigma_{u}^{>}(t_{1},t_{2})$
$\displaystyle=-i^{q_{S}}J_{0}^{2}G_{u}^{>}(t_{1},t_{2})^{q_{S}-1}$
$\displaystyle-i^{q_{I}}V_{u}^{2}G_{u}^{>}(t_{1},t_{2})^{q_{I}/2-1}G_{u,B}^{>}(t_{1},t_{2})^{q_{I}/2}$
(29) $\displaystyle-i^{q_{I}}J_{1}^{2}\sum_{v\in
V}A_{uv}G_{u}^{>}(t_{1},t_{2})^{q_{I}/2-1}G_{v}^{>}(t_{1},t_{2})^{q_{I}/2},$
$\displaystyle\Sigma_{u,B}^{>}(t_{1},t_{2})$
$\displaystyle=-i^{q_{B}}J_{B}^{2}G_{u,B}^{>}(t_{1},t_{2})^{q_{B}-1}.$ (30)

Furthermore, one can apply convolutions on either side of the first SD
equation, use the Langreth rules and integration by parts [63, 64, 37, 38], to
obtain the Kadanoff-Baym (KB) equations for the system’s Green’s functions

$\displaystyle i\partial_{t_{1}}G_{u}^{>}(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}t_{3}\left(\Sigma_{u}^{R}(t_{1},t_{3})G_{u}^{>}(t_{3},t_{2})+\Sigma_{u}^{>}(t_{1},t_{3})G_{u}^{A}(t_{3},t_{2})\right),$
(31) $\displaystyle-i\partial_{t_{2}}G_{u}^{>}(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}t_{3}\left(G_{u}^{R}(t_{1},t_{3})\Sigma_{u}^{>}(t_{3},t_{2})+G_{u}^{>}(t_{1},t_{3})\Sigma_{u}^{A}(t_{3},t_{2})\right).$
(32)

In order to make the causal structure of the KB equations more explicit and
for future use in our numerical calculations, we use the definitions of
advanced and retarded propagators to write them as

$\begin{split}i\partial_{t_{1}}G_{u}^{>}(t_{1},t_{2})&=\int_{-\infty}^{t_{1}}\mathop{}\\!\mathrm{d}t_{3}\left(\Sigma_{u}^{>}(t_{1},t_{3})+\Sigma_{u}^{>}(t_{3},t_{1})\right)G_{u}^{>}(t_{3},t_{2})-\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{3}\Sigma_{u}^{>}(t_{1},t_{3})\left(G_{u}^{>}(t_{3},t_{2})+G_{u}^{>}(t_{2},t_{3})\right),\\\
-i\partial_{t_{2}}G_{u}^{>}(t_{1},t_{2})&=\int_{-\infty}^{t_{1}}\mathop{}\\!\mathrm{d}t_{3}\left(G_{u}^{>}(t_{1},t_{3})+G_{u}^{>}(t_{3},t_{1})\right)\Sigma_{u}^{>}(t_{3},t_{2})-\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{3}G_{u}^{>}(t_{1},t_{3})\left(\Sigma_{u}^{>}(t_{3},t_{2})+\Sigma_{u}^{>}(t_{2},t_{3})\right).\end{split}$
(33)

Thus the evolution of $G_{u}^{>}(t_{1},t_{2})$ only depends on
$G_{u}^{>}(t,t^{\prime})$ evaluated at earlier times $t<t_{1}$ and
$t^{\prime}<t_{2}$, making the causal structure evident.

These equations are applicable to any non-equilibrium situations. They
describe the generic time evolution of a Green’s function, with the only
information about the physical system encoded through the definition of self-
energy in Eq. (29). In our case, we couple the system to baths at different
temperatures and observe the steady states that arise at late times.

The KB equations are generally not tractable analytically and one typically
has to resort to numerical integration. Our approach is similar to the one
used in Refs. [37, 66, 38, 39, 40, 42, 41]. We begin by solving the SD
equations self-consistently in Lorentzian time for each decoupled cluster [37,
66]. This gives us the initial equilibrium conditions for
$G_{u}^{>}(t_{1},t_{2})$ at all times $t_{1},t_{2}<0$. Then, at
$t_{1}=t_{2}=0$, we turn on the bath couplings $V_{u}$ and solve the KB
equations on a discrete grid with spacing $\mathop{}\\!\mathrm{d}t=0.05$ (in
units of inverse coupling) in the $(t_{1},t_{2})$ plane. The grid size is
usually set to $1000L\times 1000L$, where $L$ is the size of the system. In
the case of frequency-dependent driving, the time step is further decreased
such that $\omega\mathop{}\\!\mathrm{d}t\ll 1$ and the sub-sampling effects
are negligible. The equations are solved using a predictor-corrector
integration scheme [76, 41]. Since the Green’s functions typically decay
exponentially away from the diagonal $t_{1}=t_{2}$, the calculations can be
sped up by restricting our attention to a strip $|t_{1}-t_{2}|\lesssim
10\beta$ and setting all the Green’s functions to zero outside that strip
[39]. Furthermore, Eq. (27) reduces the computation time in half because we
only have to solve the KB equations for $G_{u}^{>}(t_{1},t_{2})$ with
$t_{1}>t_{2}$.

Although the KB equations do not have a closed-form solution, they simplify
significantly in the conformal and large $q$ limits, which we discuss next.

###  II.4 Conformal limit

At low temperatures $\beta J_{0}\gg 1$ and large timescales $\tau J_{0}\gg 1$,
the system develops an emergent conformal symmetry [17]. In Euclidean time,
the Green’s function for an isolated SYK cluster can be found explicitly

$G_{u}^{S}(\tau)=b\left(\frac{\pi}{\beta\sin(\pi\tau/\beta)}\right)^{2/q_{S}}\,\mathrm{sgn}(\tau),$
(34)

where $b$ is given by the equation

$b^{q_{S}}=\frac{1}{\pi
J_{0}^{2}}\left(\frac{1}{2}-\frac{1}{q_{S}}\right)\tan(\pi/q_{S}).$ (35)

Many properties of the SYK system, such as the fermion four-point function and
the chaos exponent, can be computed exactly in this limit [17]. We will use
the conformal answer when studying the equilibrium properties of our system in
Sec. III.

###  II.5 Large $q$ limit

In order to obtain an analytic approximation for the Green’s function at both
small and large energies, we take the large $q$ limit [17, 77]. For the sake
of simplicity, we set $q_{S}=q_{I}=q_{B}=q$ throughout this section. To
leading order in $1/q$, the system’s Green’s function and self-energy are
given by

$\displaystyle G_{u}^{>}(t_{1},t_{2})$
$\displaystyle=-\frac{i}{2}\left(1+\frac{g_{u}(t_{1},t_{2})}{q}+\cdots\right)\approx-\frac{i}{2}e^{\frac{g_{u}(t_{1},t_{2})}{q}},$
(36) $\displaystyle\Sigma_{u}^{>}(t_{1},t_{2})$
$\displaystyle=-\frac{i}{q}\Big{(}\mathcal{J}_{0}^{2}e^{g_{u}(t_{1},t_{2})}+\mathcal{V}_{u}^{2}e^{\frac{1}{2}(g_{u}(t_{1},t_{2})+g_{u}^{B}(t_{1},t_{2}))}$
$\displaystyle+\mathcal{J}_{1}^{2}e^{g_{u}(t_{1},t_{2})/2}\sum_{v}A_{uv}e^{g_{v}(t_{1},t_{2})/2}\Big{)},$
(37)

where “…” denotes higher order terms in $1/q$, $g_{u}(t_{1},t_{2})$ is a
function of order one satisfying $g_{u}(t,t)=0$, and
$\mathcal{J}_{0,1}=J_{0,1}\sqrt{q2^{1-q}}$,
$\mathcal{V}_{u}=V_{u}\sqrt{q2^{1-q}}$ are re-scaled couplings. The large $q$
limit is well defined only when we adjust the original couplings
($J_{0},J_{1},V_{u}$) such that the re-scaled couplings
($\mathcal{J}_{0},\mathcal{J}_{1},\mathcal{V}_{u}$) are kept fixed as
$q\to\infty$.

We can use these expressions to write the KB equations in terms of the new
functions $g_{u}(t_{1},t_{2})$ [37]. To leading order in $1/q$, we have

$\displaystyle\partial_{t_{1}}g_{u}(t_{1},t_{2})$
$\displaystyle=2iq\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{3}\Sigma_{u}^{>}(t_{1},t_{3})$
$\displaystyle-
iq\int_{-\infty}^{t_{1}}\mathop{}\\!\mathrm{d}t_{3}\left(\Sigma_{u}^{>}(t_{1},t_{3})+\Sigma_{u}^{>}(t_{3},t_{1})\right),$
(38) $\displaystyle\partial_{t_{2}}g_{u}(t_{1},t_{2})$
$\displaystyle=2iq\int_{-\infty}^{t_{1}}\mathop{}\\!\mathrm{d}t_{3}\Sigma_{u}^{>}(t_{3},t_{2})$
$\displaystyle-
iq\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{3}\left(\Sigma_{u}^{>}(t_{3},t_{2})+\Sigma_{u}^{>}(t_{2},t_{3})\right).$
(39)

Notice that in the first equation, $t_{2}$ only appears in the integration
limit on the right-hand side (and similarly for $t_{1}$ in the second
equation). This is a feature of the large $q$ expansion and it does not hold
more generally. Therefore, the second time derivative of either equation takes
on a simple form

$\frac{\partial^{2}g_{u}(t_{1},t_{2})}{\partial t_{1}\partial
t_{2}}=2iq\Sigma_{u}^{>}(t_{1},t_{2}).$ (40)

If we additionally assume time translation invariance, then we have

$\begin{split}-\frac{\partial^{2}g_{u}(t)}{\partial
t^{2}}&=2\Big{(}\mathcal{J}_{0}^{2}e^{g_{u}(t)}+\mathcal{V}_{u}^{2}e^{\frac{1}{2}(g_{u}(t)+g_{u}^{B}(t))}\\\
&+\mathcal{J}_{1}^{2}e^{g_{u}(t)/2}\sum_{v}A_{uv}e^{g_{v}(t)/2}\Big{)}.\end{split}$
(41)

This results in a system of coupled ODEs that can be solved numerically. For a
single isolated cluster, the solution can be obtained from the imaginary time
formula derived in Ref. [17] by analytic continuation

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

where $v$ satisfies

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

As shown in Ref. [37], given the symmetries of the KB equation and equilibrium
initial conditions (see Eq. (27)), the large $q$ solution always obeys

$g_{u}(-t)=g_{u}(t)^{*},$ (44)

which is indeed true for the equilibrium solution $g^{(0)}(t)$.

###  II.6 Observables

Figure 1: Schematic diagram of the non-equilibrium setups under study. (a) The
boundary-driven SYK chain is connected at both ends to thermal baths at
inverse temperatures $\beta_{L}$ and $\beta_{R}$. Each system site contains
$N$ Majorana fermions and the hopping between neighboring sites is set by
$J_{1}$, while the baths consist of $M$ Majorana fermions coupled to the chain
with interaction $V$. In NESS, a homogeneous current $j$ flows through the
bulk. (b) Same as (a), but with baths attached at every site throughout the
chain. The temperature of the baths varies linearly from left to right. (c)
Two-dimensional SYK lattice connected to baths at its two vertical edges,
resulting in a net uniform horizontal current $j_{x}$.

Once we have access to the greater Green’s functions, we can easily compute
any observable in our system. We are interested in the transport of conserved
quantities and in the case of the SYK Hamiltonian only the energy is
conserved. It is convenient to decompose the system Hamiltonian in terms of
bond operators $H_{S}=\sum_{(u,v)\in E}H^{uv}=\sum_{(u,v)\in
E}(H_{0}^{u}/d_{u}+H_{0}^{v}/d_{v}+H_{1}^{uv})$ acting on sites $(u,v)$, where
$d_{u}=\sum_{u^{\prime}}A_{uu^{\prime}}$ is the degree of vertex $u$. We then
define separately the on-site energy per particle of each cluster

$E_{0}^{u}=-i^{q_{S}+1}\frac{J_{0}^{2}}{q_{S}}\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{1}\left(G_{u}^{>}(t,t_{1})^{q_{S}}-G_{u}^{>}(t_{1},t)^{q_{S}}\right),$
(45)

and similarly the interaction energy between clusters

$\displaystyle E_{1}^{uv}$
$\displaystyle=-i^{q_{I}+1}\frac{2J_{1}^{2}}{q_{I}}\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{1}\Big{(}G_{u}^{>}(t,t_{1})^{q_{I}/2}G_{v}^{>}(t,t_{1})^{q_{I}/2}$
$\displaystyle-
G_{u}^{>}(t_{1},t)^{q_{I}/2}G_{v}^{>}(t_{1},t)^{q_{I}/2}\Big{)}.$ (46)

The on-bond energy then becomes $E_{uv}\equiv\langle
H^{uv}\rangle=E_{0}^{u}/d_{u}+E_{0}^{v}/d_{v}+E_{1}^{uv}$. The formula for the
associated local energy current $j_{u}$ per particle flowing from $u$ to $v$
can be derived by combining the continuity equation at site $u$ with
Heisenberg’s equation of motion [78]

$j_{u}=i\sum_{u^{\prime}}A_{u^{\prime}u}[H^{u^{\prime}u},H^{uv}].$ (47)

Computing the expectation value of these commutators in the Schwinger-Keldysh
formalism is more complicated. We provide a detailed derivation in Appendix A
and our general formula for the current is given by Eq. (87).

In this paper, we will mainly be interested in regular lattices (one- and two-
dimensional) with baths acting on the boundaries (see Fig. 1). The right and
left baths are held at different inverse temperatures
$\beta_{R,L}=\beta\pm\delta\beta$, where $\beta$ is the average bath
temperature and $\delta\beta$ is a small temperature imbalance driving the
system out of equilibrium. We typically choose $\delta\beta=0.1\beta$ to make
sure we are in the linear response regime. In the long-time limit, when the
system reaches its NESS, the energy current becomes uniform in the bulk
$j\equiv\langle j_{u}\rangle$. For the particular case of a uniform one-
dimensional chain (Fig. 1(a)), the current reads

$j=\frac{i^{q_{S}+q_{I}}}{2}J_{1}^{2}J^{2}\real(j_{++}^{x-1,x,x+1}+j_{+-}^{x-1,x,x+1}),$
(48)

where we labeled the sites by their horizontal coordinate $x$. Similarly, the
energy gradient is simply the difference in energy on two consecutive bonds

$\nabla E=E_{x,x+1}-E_{x-1,x}.$ (49)

For a uniform setup with $q_{S}=q_{I}$, our model is a diffusive metal [43].
We expect that transport in such a system is governed by Fourier’s law
$j=-D\nabla E=-D\Delta E/L$ with a temperature-dependent diffusion constant
$D$. Here $\Delta E$ denotes the energy difference across the system and $L$
is the linear size of the system. More generally, a system can exhibit non-
diffusive transport where the current scales with system size as $j\sim
1/L^{\gamma}$ with $\gamma\neq 1$ [1]. We explicitly verify in our numerical
simulations that the transport is indeed diffusive at all temperatures, except
for the case $q_{S}<q_{I}$, when the system becomes an insulator at low
temperatures (see Sec. IV.5).

In order to accurately describe the temperature dependence of the diffusion
constants, we need to introduce an effective local temperature. This
temperature can in principle be very different from the average bath
temperature at which we drive the system [4, 79]. At sufficiently long times
after turning on the bath couplings, the system reaches a nearly-thermal state
where the Green’s functions become time-translation invariant. This allows us
to use the fluctuation-dissipation theorem (FDT) [63, 64, 37] to define an
effective temperature $\beta_{u}^{-1}$ similarly to the thermal equilibrium
case

$\frac{iG_{u}^{K}(\omega)}{A_{u}(\omega)}=\tanh(\frac{\beta_{u}\omega}{2}),$
(50)

where $A_{u}(\omega)=-2\imaginary G_{u}^{R}(\omega)$ is the spectral function
and $G_{u}^{K}(\omega)$ is the Fourier transform of the Keldysh Green’s
function. Note that the FDT holds only for low frequencies, since high
frequencies are affected by the size of the discretization step
$\mathop{}\\!\mathrm{d}t$. One can extract the local temperature by fitting
Eq. (50) over a small frequency interval [37, 38]. Alternatively, we can
access the low-frequency limit directly by taking the slope of the ratio of
Green’s functions at $\omega=0$ [40, 62]

$\beta_{u}=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\omega}\left(\frac{2iG_{u}^{K}(\omega)}{A_{u}(\omega)}\right)\bigg{\rvert}_{\omega=0}.$
(51)

This simple relation unambiguously defines the local temperature of an SYK
cluster, in contrast with conventional spin systems, where defining a non-
equilibrium temperature is more challenging [4]. Since the system is many-body
chaotic, we expect it to reach a nearly-thermal state relatively fast, with a
local temperature closely matching the bath temperature $\beta$ [39, 40].

Finally, we would like to study how fast we approach the non-equilibrium
steady state. This would allow us to compare the effectiveness of different
driving setups. The convergence time to NESS can be defined for any of our
observables, such as the current or the diffusion constant. Suppose we measure
an observable $z_{t}$ at each step $t$ during our time evolution after
coupling the system to the baths. To assess the convergence of a sequence of
measurements $\\{z_{t}\\}_{t=0}^{T}$, we define the convergence time in terms
of the error relative to the asymptotic value
$z_{\infty}=\lim\limits_{T\to\infty}z_{T}$ of a given observable

$\epsilon(t)=\left|1-\frac{z_{t}}{z_{\infty}}\right|.$ (52)

Note that this metric is invariant under rescaling, but not under a
translation by a constant of the data. At late times, we expect the
convergence to be exponentially fast $\epsilon(t)\sim e^{-\Gamma t}$, where
$\Gamma$ is the convergence rate [39]. We can also define a convergence time
$t_{\mathrm{NESS}}$ as the minimum time such that $\epsilon(t)\leq\eta$ for
all $t\geq t_{\mathrm{NESS}}$, where $\eta$ is a set threshold. Notice that
the definition of $t_{\mathrm{NESS}}$ does not assume exponential convergence
to NESS and is therefore more general. By comparing convergence times, we can
quantitatively assess whether a certain bath configuration is more efficient
at inducing the non-equilibrium steady state. In Sec. IV.2 we will show that
this is indeed the case for bulk driving.

##  III Equilibrium

We begin by studying the equilibrium setup, where all the bath and system
clusters are at the same inverse temperature $\beta$. This will help us
understand the changes that a local Green’s function undergoes by simply
connecting to its neighbors and the baths, before even introducing non-
equilibrium effects caused by temperature imbalance. Since there is no
dynamics, we look for solutions of the SD equations in Euclidean time. We will
focus on the case $q_{S}=q_{I}=q_{B}=q$, but the more general setup was
thoroughly discussed in Refs. [39, 40]. This case corresponds to a marginal
interaction, such that the effective $q$ does not change, while the effective
on-site interaction $J_{0}$ gets renormalized to $J_{u}$. In other words, the
SD equation for an interacting SYK cluster has the exact same form as that of
an isolated $(0+1)$ dimensional SYK model with coupling $J_{u}$

$\Sigma_{u}^{S}(\tau)=J_{u}^{2}G_{u}^{S}(\tau)^{q-1}.$ (53)

We will consider the conformal limit for simplicity, but the same qualitative
results can be obtained by solving the SD equations numerically using the
approach described in Ref. [17] with minor adjustments to the self-energy
formula as given by Eq. (16). If $q$ is the same throughout the system, then
the time-dependence cancels completely in the SD equations in the conformal
limit. From Eq. (35), it is easy to see that the Green’s function depends on
the coupling as $G_{u}^{S}\sim J_{u}^{-2/q}$ and $\Sigma_{u}^{S}\sim
J_{u}^{2/q}$. Therefore, Eq. (16) simplifies to

$\frac{J_{0}^{2}}{J_{u}^{2}}+\frac{J_{1}^{2}}{J_{u}}\sum_{v}\frac{A_{uv}}{J_{v}}+\frac{V_{u}^{2}}{J_{B}J_{u}}=1,$
(54)

for each site $u$. This is a system of quadratic equations for the effective
coupling $J_{u}$, whose general solution can always be found numerically.

Figure 2: Equilibrium properties of an SYK chain of length $L=20$. The
effective on-site coupling $J_{x}$ is shown as a function of system site in
the case of boundary driving (dash-dotted line), strong bulk driving (solid
line), and weak bulk driving (dotted line). The asymptotic values for the
coupling deep in the bulk are shown with dashed lines. We set $J_{0}=J_{1}=1$.

We further specialize to the case of a one-dimensional chain of $L$ SYK
clusters, labeled by their lattice position $x$, with nearest-neighbor
couplings and open boundary conditions (i.e. the only non-zero entries of the
adjacency matrix are $A_{x,x-1}=A_{x,x+1}=1$). We then investigate
configurations in which the baths are either coupled at the two ends of the
chain (boundary driving, see Fig. 1(a)) or at every site of the chain (bulk
driving, Fig. 1(b)). The system-bath coupling always has the same value $V$.
This setup is appropriate for studying the transport properties of a one-
dimensional SYK system [43, 55, 44, 45, 46, 47, 66, 67] and can be easily
generalized to higher dimensions (Fig. 1(c)).

Our results for the effective coupling $J_{x}$ obtained by solving the system
of equations (54) is shown in Fig. 2. In the case of boundary driving, we see
a small edge effect due to the presence of the baths, but otherwise the
effective coupling approaches its site-independent value
$J=\sqrt{J_{0}^{2}+2J_{1}^{2}}$ deep in the bulk. Note that this effective
coupling is slightly different from the one derived in Ref. [43] due to a
distinct choice of normalization for the variance of the coupling constant
$J_{1}$. When baths are present on every site, their interaction with the
system will strongly renormalize the coupling throughout the chain. For
example, if we choose $J_{B}=J$, then the effective coupling away from the
boundaries converges to $J=\sqrt{J_{0}^{2}+2J_{1}^{2}+V^{2}}$, which is very
different from the boundary-driven case. In non-equilibrium, this bulk driving
would significantly alter the transport coefficients.

We consider two approaches for mitigating the unwanted effects of bulk
driving, while preserving some of its appealing features, such as fast NESS
convergence rates. First, we could adjust the couplings
$\tilde{V}=\tilde{J}_{0}=J_{0}/\sqrt{2}$ such that the effective on-site
interaction in the bulk becomes the same as in the case of boundary driving
$J=\sqrt{J_{0}^{2}+2J_{1}^{2}}$. We refer to this as strong bulk driving. The
issue with this solution is that it slightly changes the properties of the
underlying system by reducing $J_{0}$. Alternatively, we could take the limit
of weak bulk driving $V\ll J_{0},J_{1}$, which leaves the effective coupling
mostly unchanged (see Fig. 2). However, the rate at which energy is exchanged
with the baths will also decrease, leading to slightly longer convergence
times. The weak driving scheme serves as a middle ground between the strong
and boundary drivings. In the next section, we will further explore the
advantages of each setup in the context of non-equilibrium dynamics.

##  IV DC Transport

We now proceed with our results for the non-equilibrium case in the presence
of time-independent (DC) driving. Following the procedure described in Sec.
II, we couple the system to baths at different temperatures and measure the
energy gradients, local temperatures, and currents that arise as a result. We
first consider the case of boundary driving and investigate the properties of
the emerging NESS. We then show that bulk driving leads to approximately the
same diffusion constants, but with much faster convergence times. Based on
these findings, we postulate a simple ansatz for the non-equilibrium Green’s
function and show that it perfectly captures the transport properties of the
SYK system. Using this ansatz, we derive a closed-form solution for the
temperature-dependent diffusion coefficient in the large $q$ limit. Throughout
these sections, we focus on one-dimensional chains and set
$q_{S}=q_{I}=q_{B}=q$ for simplicity. We discuss the case with $q_{S}\neq
q_{I}$ separately in Sec. IV.5 and conclude our analysis of DC transport with
a generalization to higher-dimensional lattices. To emphasize that our methods
are applicable to a wide range of parameters, we display results for various
values of $\beta$ and $q$. We also fix $\mathcal{J}_{0}=\mathcal{J}_{1}=1$ in
order to treat the results for different $q$ on equal footing.

###  IV.1 Boundary-driven SYK chain

In the case of boundary driving, we introduce a small temperature imbalance
$\delta\beta$ at the two ends of the chain. Fig. 3 showcases our findings for
$q=4$, but similar results hold for all values of $q$. We should mention that
although the SYK model with $q=2$ corresponds to free fermions, is non-
integrable, and does not satisfy ETH [35, 36, 38], it still slowly thermalizes
when coupled to external baths and conducts energy in the same way as its
analogs with higher $q$. In Fig. 3(a) we plot the on-bond energy $E_{x,x+1}$,
normalized by its equilibrium value $E_{0}$, as a function of lattice site
$x$. The energy profiles are linear in the bulk, with small boundary effects
that become more prominent at low temperatures. Moreover, we find that local
Green’s functions obey the FDT with inverse temperatures that interpolate
linearly between $\beta-\delta\beta$ and $\beta+\delta\beta$ at the two ends,
as shown in Fig. 3(b). The temperature in the middle of the chain precisely
agrees with the average bath temperature. This is different from chaotic spin
models, where one usually finds a much higher local temperature in the bulk
[4]. We also find that the current scales linearly with inverse system size at
all temperatures, as displayed in Fig. 3(c). Notice that the scaling exponent
$\gamma$ is close to one, which confirms that our model is a diffusive metal.

Figure 3: DC NESS transport properties of the boundary-driven SYK chain with
$q=4$ at different bath temperatures $\beta$. (a) Spatial energy profiles,
rescaled to the equilibrium energy $E_{0}$ at temperature $\beta$, for a chain
of length $L=50$. (b) Same as (a), but for the local inverse temperature. Both
energy and temperature profiles are linear in the bulk. (c) Scaled energy
current $j/\Delta E$ as a function of system size. Symbols denote numerical
values and lines represent fits to the scaling $j/\Delta E=-D/L^{\gamma}$. The
values of $\gamma$ suggest diffusive transport according to Fourier’s law. (d)
Derivative of the equilibrium Green’s function
$\mathop{}\\!\mathrm{d}G_{0}(t)/\mathop{}\\!\mathrm{d}\beta$ (dashed lines)
showing agreement with the non-equilibrium contribution $F(t)$ (solid lines)
to the Green’s functions in Eq. (56), as extracted numerically from NESS. Both
are scaled to unit norm. Figure 4: Temperature dependence of the DC diffusion
constant for different values of $q$. Filled circles represent numerical
values for the boundary-driven SYK chain in both panels. (a) Empty circles
(shaded diamonds) correspond to the case of weak (strong) bulk driving. Solid
lines represent the diffusion constant computed using the ansatz
$F(t)=dG_{0}(t)/d\beta$, while the dashed lines indicate the conformal limit
derived in Ref. [43]. (b) Dashed lines show the diffusion constants in the
large $q$ limit.

In addition, we observe a very fast convergence of the diffusion constant with
system size, which is likely due to the all-to-all nature of the couplings in
the SYK model. Therefore, when computing the temperature dependence of the
diffusion coefficient, we can restrict our attention to moderately-sized
systems with $L=20$ sites. Our results for $D(\beta)$ at various $q$ can be
found in Fig. 4(a). The transport coefficient increases with $\beta$ and
approaches a $q$-dependent constant both in the limit of zero and infinite
temperature. The answer for the diffusion constant of an SYK chain in the low-
temperature conformal limit has been previously derived in Ref. [43] by
studying the leading order contributions to the fermionic four-point
functions. Their derivation can be generalized to the case of arbitrary $q$ in
a straightforward way and we obtain

$D=\frac{q}{q-1}\frac{\pi}{2\alpha_{K}}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}},$
(55)

where $\mathcal{J}=J\sqrt{q2^{1-q}}$ is the rescaled effective coupling
introduced in Sec. III and $\alpha_{K}$ is a numerical constant whose value is
close to $3$ (e.g. Fig. 9 in [17]). We recover the original expression in Ref.
[43] by setting $q=4$. The conformal answer in Eq. (55) is indicated by dashed
lines in Fig. 4(a). We observe a striking agreement with our numerical results
at large $\beta$.

Most importantly, the authors also showed that in the conformal limit, the
energy diffusion constant $D$ is related to the butterfly velocity $v_{B}$ via
the simple relation $D=\beta v_{B}^{2}/2\pi$, which realises a conjectured
bound on diffusion in incoherent metals [49, 50, 51, 52, 53, 54]. We further
explore the connections between transport and many-body chaos in Appendix B,
where we show that the energy diffusion is upper bounded by chaos $D\leq
v_{B}^{2}/\lambda_{L}$, with $\lambda_{L}$ denoting the Lyapunov exponent. The
same conclusion has been reached for the energy diffusion constant in
inhomogeneous SYK chains [55] and the charge diffusion constant in a
holographic model [65].

###  IV.2 Bulk-driven SYK chain

Our non-equilibrium driving scheme with baths on the edges can be used in
principle to extract the diffusion constant of arbitrary models at any
temperature. However, in practice, the times for which we have to numerically
evolve the KB equations until convergence can be quite long, especially at low
temperatures and for the $q=2$ case. The main computational cost comes from
the fact that the influence of the bath has to propagate from the boundary all
the way to the middle of the system. Therefore, it seems natural to search for
an alternative setup which circumvents this problem.

Motivated by the linearity of the temperature profile along the chain, we
propose to couple each site $x$ to a bath at inverse temperature $\beta_{x}$
which linearly interpolates between $\beta_{L}$ and $\beta_{R}$, as shown in
Fig. 1(b). We have already studied the equilibrium properties of the system in
the presence of this bulk driving in Sec. III. We identified two scenarios in
which the effective coupling $\mathcal{J}$ in the bulk remains mostly
unchanged. In the case of strong bulk driving, we set
$\mathcal{V}=\mathcal{J}_{0}=1/\sqrt{2}$ and in the case of weak bulk driving,
we choose $\mathcal{V}=0.2$ and $\mathcal{J}_{0}=1$. We expect that the bulk
driving does not significantly alter the properties of our system even in non-
equilibrium and its diffusion coefficient remains unchanged as long as
$\mathcal{J}$ and $\mathcal{J}_{1}$ stay the same in the bulk. Indeed, in the
low-temperature limit we can confirm this directly using Eq. (55). For the
more general case, we solve the KB equations in the presence of bulk driving
and plot our results in Fig. 4(a). We see that the diffusion constants for
both strong and weak driving agree well with the results obtained via boundary
driving, with the weak driving showing a slightly better match at low
temperatures.

Figure 5: Convergence time of the non-equilibrium dynamics for boundary (solid
circles) and bulk (empty circles) driving at various temperatures. The
boundary-driven systems exhibit a power law scaling $t_{\mathrm{NESS}}\sim
L^{z}$ with system size (solid lines), while both weak and strong bulk driving
lead to nearly-constant convergence times (dashed lines). The convergence
threshold was set to $\eta=10^{-3}$.

However, the main appeal of using these alternative driving schemes is their
fast convergence times. We compare the times required to reach NESS in Fig. 5.
For the setup with baths at the boundaries, we find a clear power law
dependence of $t_{\mathrm{NESS}}$ on system size. Notice that the convergence
time scales superlinearly with $L$ and the exponent is somewhat higher at low
temperatures. On the other hand, $t_{\mathrm{NESS}}$ is approximately constant
in the case of weak driving and the advantage over boundary driving is already
noticeable around $L\approx 20$. Furthermore, strong bulk driving leads to an
even lower size- and temperature-independent convergence time. These scalings
are roughly the same for all $q$. Hence, there is a clear computational
advantage in adding baths to the bulk with little to no impact on the
transport properties of the system. We expect that a similar speed-up can be
achieved in other spin models.

It is important to emphasize that using bulk driving to exactly match the
results of boundary driving was only possible due to the simple, predictable
nature of the temperature profile across the system. Moreover, we relied on
the fact that the baths are similar to the system SYK clusters (i.e. same $q$,
similar $\beta$ etc.). This may no longer be the case in other, more complex
non-equilibrium scenarios. In fact, we will briefly discuss this in the
context of AC transport in Sec. V.3.

###  IV.3 Non-equilibrium ansatz

Our analysis of the steady-state observables in the previous sections suggests
a relatively simple ansatz for the local non-equilibrium Green’s function that
emerges at late times. In equilibrium, the Green’s function deep in the bulk
takes on a site-independent and time translation-invariant form $G^{(0)}(t)$,
which is the solution to the SD equation in Lorentzian time with an effective
coupling $J$ and inverse temperature $\beta$. In the presence of a small
uniform bias, we expect the non-equilibrium Green’s function to have an extra
site-dependent correction proportional to the temperature gradient.
Additionally, we expect the Green’s function to become time translation-
invariant at late times, which we indeed observe in our simulations.
Therefore, we look for a NESS solution to the KB equations of the form

$G_{x}^{>}(t)=G^{(0)}(t)+xF(t),$ (56)

where $|F(t)|\ll|G^{(0)}(t)|$ is the non-equilibrium contribution. Similarly,
we can expand the self-energy as

$\Sigma_{x}^{>}(t)=-i^{q}J^{2}G^{(0)}(t)^{q-2}\left(G^{(0)}(t)+x(q-1)F(t)\right).$
(57)

These expressions already lead to a much simpler version of the KB equations
in (33) for the non-equilibrium contribution $F(t)$. However, we can find an
explicit solution for $F(t)$ by exploiting the fact that the inverse
temperature gradient $\nabla\beta$ is constant in the bulk (see Fig. 3(b)).
Assuming that the effective coupling $J$ also remains constant, we can perform
a first-order expansion in the small gradient

$G_{x}^{>}(t)\approx
G^{(0)}(t)+x\frac{\mathop{}\\!\mathrm{d}G^{(0)}(t)}{\mathop{}\\!\mathrm{d}\beta}\nabla\beta,$
(58)

This is akin to a gradient expansion in hydrodynamics. By comparing this to
Eq. (56), without loss of generality, we can identify
$F(t)\equiv\mathop{}\\!\mathrm{d}G^{(0)}(t)/\mathop{}\\!\mathrm{d}\beta$.
Notice that the overall magnitude of $F(t)$ does not matter, since to first
order, both the energy gradient and the current will be proportional to $F$.
This is consistent with the expectation that the exact value of the
temperature bias $\delta\beta$ should not affect the transport properties, as
long as we are within the linear response regime. A comparison between $F(t)$
extracted numerically as the difference between two consecutive Green’s
functions in NESS and the derivative
$\mathop{}\\!\mathrm{d}G^{(0)}(t)/\mathop{}\\!\mathrm{d}\beta$ is presented in
Fig. 3(d). As we argued above, the two functions coincide.

Furthermore, one can easily verify using Eq. (56) and Eq. (57) that if
$G^{(0)}(t)$ is a solution to the KB equations in (33), then so is
$G_{x}^{>}(t)$ in Eq. (58), to first order in the gradient. We emphasize that
$F(t)=\mathop{}\\!\mathrm{d}G^{(0)}(t)/\mathop{}\\!\mathrm{d}\beta$ is a
solution to the KB equations, but it is not necessarily the unique solution.
One could in principle devise more complicated non-equilibrium setups and
initial conditions where Eq. (56) holds, but the non-equilibrium correction
$F(t)$ has a more complicated form. This is why the numerical comparison in
Fig. 3(d) is crucial.

Now we can use this ansatz to obtain a numerical solution for the diffusion
constant at arbitrary temperature. To first order in $F$, the energy gradient
in Eq. (49) becomes

$\nabla
E=-2i^{q+2}J^{2}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[G^{(0)}(t)^{q-1}F(t)\right]=\frac{\mathop{}\\!\mathrm{d}E_{0}}{\mathop{}\\!\mathrm{d}\beta}.$
(59)

where $E_{0}$ is the equilibrium energy computed from $G^{(0)}(t)$ and we have
used the identities in Eq. (27) to write the integral only in terms of
functions evaluated at $t\geq 0$. Similarly, the current in Eq. (48) becomes

$j=\frac{1}{2}J_{1}^{2}J^{2}\real(j_{++}+j_{+-}),$ (60)

$\begin{split}j_{++}+j_{+-}&=-2iq\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}G^{(0)}(t^{\prime})^{q-2}G^{(0)}(t^{\prime}-t)\left(G^{(0)}(t^{\prime})\,\mathrm{Im}\left[F(t)G^{(0)}(t)^{q-2}\right]-F(t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t)^{q-1}\right]\right)\\\
&=-\frac{2iq}{q-1}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}G^{(0)}(t^{\prime})^{2q-2}G^{(0)}(t^{\prime}-t)\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\beta}\left(\frac{\,\mathrm{Im}\left[G^{(0)}(t)^{q-1}\right]}{G^{(0)}(t^{\prime})^{q-1}}\right).\end{split}$
(61)

Note that both of these quantities are independent of $x$ and we can compute
the diffusion constant $D=-j/\nabla E$ solely in terms of the equilibrium
function $G^{(0)}(t)$.

In general, the Green’s function does not have a closed-form representation,
except for the case of free fermions $q=2$, which we discuss in Appendix C.
Therefore, we have to numerically solve the SD equation to get $G^{(0)}(t)$
[37, 66], then perform the integrals above to obtain $D$. The final result is
shown as the solid lines in Fig. 4(a) and it agrees perfectly with diffusion
constants extracted from solving the full KB equations in the presence of
baths (solid circles in Fig. 4(a)).

Our non-equilibrium ansatz significantly reduces the complexity of computing
the transport properties of the SYK chain compared to the boundary- or bulk-
driven setups. Moreover, under certain circumstances, such as the large $q$
limit discussed below, we can use this ansatz to find simple expressions for
the diffusion constant.

###  IV.4 Large $q$ limit

We now apply our construction in the large $q$ limit introduced in Sec. II.5.
For a 1D chain, Eq. (41) away from the boundary can be written as

$-\frac{\partial^{2}g_{x}(t)}{\partial
t^{2}}=2\mathcal{J}_{0}^{2}e^{g_{x}(t)}+2\mathcal{J}_{1}^{2}e^{\frac{g_{x}(t)}{2}}\left(e^{\frac{g_{x-1}(t)}{2}}+e^{\frac{g_{x+1}(t)}{2}}\right).$
(62)

In equilibrium, $g^{(0)}(t)$ derived in Eq. (42) is a uniform solution.
However, in the presence of a small temperature gradient, we look for a
general solution of the form

$g_{x}(t)=g^{(0)}(t)+xf(t),$ (63)

where $|f(t)|\ll|g^{(0)}(t)|$ is the non-equilibrium correction. Just as
before, we identify
$f(t)\equiv\mathop{}\\!\mathrm{d}g^{(0)}(t)/\mathop{}\\!\mathrm{d}\beta$. In
terms of this new function we have $F(t)=G^{(0)}(t)f(t)/q$. Since $f$ is
proportional to the imposed temperature gradient, we can make it arbitrarily
small and only keep terms linear in $f$. The KB equation becomes

$\frac{\partial^{2}f(t)}{\partial t^{2}}=-2\mathcal{J}^{2}e^{g^{(0)}(t)}f(t),$
(64)

where $\mathcal{J}^{2}=\mathcal{J}_{0}^{2}+2\mathcal{J}_{1}^{2}$ and we
canceled the equilibrium terms. We supplement this ODE with the initial
condition $f(0)=0$, imposed by $g_{x}(0)=g^{(0)}(0)=0$. However, since the
equation is of second order, we have to also specify $f^{\prime}(0)$, on which
we have (almost) no restrictions. This is a consequence of the aforementioned
freedom in choosing the overall scale of $f$. Note that $f(t)$ also has to
obey the identity in Eq. (44).

Given our guess for $f(t)$ and Eq. (42), we find the following non-equilibrium
contribution

$f(t)=-\tan(\frac{\pi v}{2})+\Big{(}1+i\mathcal{J}t\sin(\frac{\pi
v}{2})\Big{)}\tan(\frac{\pi v}{2}-\frac{i\pi vt}{\beta})$ (65)

It is easy to check that this is indeed a solution to Eq. (64) and that it
satisfies all the conditions mentioned above. In the large $q$ limit, the
energy gradient reads

$\begin{split}\nabla
E&=\frac{\mathcal{J}^{2}}{q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[e^{g^{(0)}(t)}f(t)\right]\\\
&=-\frac{1}{2q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[\frac{\partial^{2}f(t)}{\partial
t^{2}}\right]\\\
&=\frac{1}{2q^{2}}\,\mathrm{Im}\left[f^{\prime}(0)-f^{\prime}(\infty)\right]\\\
&=-\frac{\mathcal{J}}{2q^{2}}\cos(\frac{\pi v}{2}),\end{split}$ (66)

where we used Eq. (64) in the second line. The formula for the current is a
bit more involved

$\begin{split}j&=\frac{\mathcal{J}_{1}^{2}\mathcal{J}^{2}}{2q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\Bigg{[}\,\mathrm{Re}\left[e^{\left(1-\frac{1}{q}\right)g^{(0)}(t)}f(t)\right]\,\mathrm{Im}\left[e^{\left(1-\frac{1}{q}\right)g^{(0)}(t^{\prime})+g^{(0)}(t^{\prime}-t)/q}\right]\\\
&-\,\mathrm{Re}\left[e^{\left(1-\frac{1}{q}\right)g^{(0)}(t)}\right]\,\mathrm{Im}\left[e^{\left(1-\frac{1}{q}\right)g^{(0)}(t^{\prime})+g^{(0)}(t^{\prime}-t)/q}f(t^{\prime})\right]\Bigg{]}.\end{split}$
(67)

Unfortunately, this double integral does not have a closed-form solution for
arbitrary $q$ and $\beta$. However, we are able to approximate it in both the
$q\to\infty$ limit with arbitrary $\beta$ and in the $\beta\to 0$ and
$\beta\to\infty$ limits with arbitrary $q$. The details of this computation
are provided in Appendix C. For example, the $q\to\infty$ limit yields

$D=\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\left(\frac{\pi v}{2}\sin(\frac{\pi
v}{2})+\cos(\frac{\pi v}{2})\right).$ (68)

The temperature dependence enters this expression implicitly through $v$ (see
Eq. (43)). Note that this equation exactly matches the one derived in Ref.
[54] from the energy density two-point function, thus providing an independent
consistency check for our NESS calculation.

Our non-equilibrium ansatz led to a remarkably simple formula for the
diffusion constant in the limit of infinite $q$. We plot this result in Fig.
4(b). In the same figure, we also show the diffusion constants obtained by
numerically integrating Eq. (67). The large $q$ expansion curves follow the
same trend as the previously obtained exact results. The agreement between the
two is significantly better at higher temperatures and larger $q$. The
diffusion constant clearly decreases with $q$ and we expect it to eventually
approach the $q=\infty$ result (black line in Fig. 4(b)).

###  IV.5 Models with $q_{S}\neq q_{I}$

In all the examples previously studied in this paper, we found that our
systems behaved like diffusive metals with a finite energy diffusion constant
at both zero and infinite temperatures. However, this is not the case for all
SYK chains. In this section, we will discuss a family of models with
$q_{S}\neq q_{I}$, which were first introduced in the context of metal to
insulator transitions [48]. First, we consider a model with $q_{I}=2q_{S}=2q$,
where $q\geq 2$. This model is equivalent to isolated SYK clusters with
$q_{S}=q$ in the IR limit, since the inter-cluster coupling becomes
irrelevant. Hence, the system should become an insulator at low temperatures
[48]. To confirm this prediction, we compute the diffusion constant as a
function of temperature using the boundary-driven setup for $q=4$. The results
are shown in Fig. 6(a). We can see that the diffusion constant approaches zero
in the low-temperature limit. In fact, we can further observe that $D$
decreases quadratically with temperature, as depicted in the inset. This
agrees with the large $q$ calculation of Ref. [48], which predicts a decay
$D\sim\beta^{-2}$.

Figure 6: Temperature dependence of the diffusion constant for other SYK
systems. (a) One-dimensional chain with $q_{S}>q_{I}$ (red) and $q_{S}<q_{I}$
(blue). For the latter, the inset shows a quadratic fit $D\sim\beta^{-2}$ in
the conformal limit. (b) Two-dimensional uniform square lattice. Filled
circles represent numerical values obtained from the KB equations, while empty
circles denote values inferred from one-dimensional data. Dashed lines
indicate the conformal limit in Eq. (70).

Next, we consider the opposite model where $q_{S}=2q_{I}=2q$. In this model,
we expect the on-site coupling $\mathcal{J}_{0}$ to become irrelevant in the
IR, where the low-temperature transport should be governed by the intra-
cluster coupling $\mathcal{J}_{1}$. The system is expected to behave like a
diffusive metal, which we confirm numerically in Fig. 6(a) for the case $q=4$.
The diffusion constant approaches a constant both at the high and low
temperatures. According to the large $q$ analysis of Ref. [48], the diffusion
constant should converge to

$D=\frac{\pi\mathcal{J}_{1}}{3\sqrt{1+\mathcal{J}_{0}^{2}/8\mathcal{J}_{1}^{2}}}\approx
0.99,$ (69)

for $\mathcal{J}_{0}=\mathcal{J}_{1}=1$. Although our results seem to converge
to a somewhat smaller value, this discrepancy is on par with the one in our
other large $q$ predictions in Fig. 4(b).

###  IV.6 Higher-dimensional generalizations

So far, we have focused on one-dimensional SYK chains, but our analysis can be
easily generalized to higher-dimensional lattices. For example, we can
consider a two-dimensional $L_{x}\times L_{y}$ square lattice with baths
attached along the vertical boundaries, as shown in Fig. 1(c). Let the lattice
dimensions be $L_{x}=15$ and $L_{y}=4$, with periodic boundary conditions in
the $y$-direction. The boundary driving will impose a current $j_{x}$ in the
$x$-direction, but there will be no net current in the $y$-direction due to
symmetry. We time-evolve the KB equations in the usual way and use Eq. (89) to
compute the current by summing contributions from all four neighbors of a
given site.

Our results for the square lattice are given in Fig. 6(b). We observe the same
temperature dependence as in the one-dimensional case. The asymptotic value at
low temperatures can be compared again to the conformal limit answer in Ref.
[43]. For a translation invariant lattice, we expect the pole determining the
diffusion constant to get equal contributions from both the $p_{x}^{2}$ and
$p_{y}^{2}$ momenta, thus doubling the value in Eq. (55)

$D=\frac{q}{q-1}\frac{\pi}{\alpha_{K}}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}},$
(70)

where $\mathcal{J}=\sqrt{\mathcal{J}_{0}^{2}+4\mathcal{J}_{1}^{2}}$ is now the
effective coupling in 2D. This conformal answer is depicted by dashed lines in
Fig. 6(b) and it agrees well with our data.

Given the simplicity of the low-temperature answer above, one could attempt to
estimate the two-dimensional diffusion constant at arbitrary temperatures by
appropriately rescaling the one-dimensional data

$D_{\mathrm{2D}}=2\sqrt{\frac{\mathcal{J}_{0}^{2}+2\mathcal{J}_{1}^{2}}{\mathcal{J}_{0}^{2}+4\mathcal{J}_{1}^{2}}}D_{\mathrm{1D}},$
(71)

where we took into account the difference in effective coupling. This guess is
represented by empty circles in Fig. 6(b) and it is in surprisingly good
agreement with the actual 2D data, even though the diffusion constant at high
temperatures has a more complicated dependence on $\mathcal{J}$. This result
suggests that diffusion in higher dimensional SYK lattices is very similar to
its one-dimensional counterpart.

##  V AC Transport

We turn our attention to the frequency dependence of the diffusion coefficient
$D(\omega)$. This quantity is fundamental, as it describes the response of the
system to an external time-dependent bias with angular frequency $\omega$. The
frequency dependence is typically computed via the Kubo formula [80], which
relates transport coefficients to current autocorrelation functions. However,
our framework is more suitable for directly measuring the current that emerges
as a result of periodic driving.

In order to induce a time-dependent energy gradient in our system, we must
couple it to baths at inverse temperatures
$\beta_{R,L}=\beta\pm\delta\beta\cos(\omega t)$. Since it is non-trivial to
treat a time-dependent temperature directly within the Schwinger-Keldysh
formalism, we choose instead to oscillate the on-site coupling of the baths
$J_{B}\pm\delta J_{B}\cos(\omega t)$, where $\delta J_{B}\ll J_{B}$. Note that
the KB equations remain virtually unchanged in the presence of a time-
dependent coupling [37, 38, 42], which only changes the local energy scale for
the baths and leads to an oscillating temperature. In the limit of $\omega\to
0$ we exactly recover the temperature imbalance $\beta\pm\delta\beta$ imposed
in the DC case.

Figure 7: AC NESS transport properties of a boundary-driven SYK chain of
length $L=11$ with $q=4$ at $\beta=1$ and $\omega=5$. Time-dependent (a)
energy, (b) local inverse temperature, and (c) current for consecutive sites
$x$, showing a uniform amplitude decay $c=0.13$ and phase shift $\delta=1.74$.
The means were subtracted for easier visualization. (d) The non-equilibrium
contributions $F_{1,2}(t)$ to the Green’s functions extracted directly from
the boundary-driven NESS (solid lines) and computed using the ansatz in
Appendix D (dashed lines).

We proceed in a manner similar to the DC case. First, we attach the baths at
the boundary of our SYK chain and study the resulting NESS. We find that all
observables oscillate at frequency $\omega$ with an amplitude that becomes
more attenuated the further we measure into the bulk. Based on these insights,
we propose a new ansatz for the non-equilibrium Green’s function and solve the
KB equations numerically using this functional form. We show that the ansatz
describes well the AC transport in SYK models. Last, we briefly comment on the
setup in the presence of bulk driving. Just as in the case of DC transport, we
set $q_{S}=q_{I}=q_{B}=q$ and $\mathcal{J}_{0}=\mathcal{J}_{1}=1$.

###  V.1 Boundary-driven SYK chain

Figure 8: Frequency dependence of the AC diffusion constant for (a) $q=4$ and
(b) $q=6$ at various temperatures. Filled circles represent numerical values
for the boundary-driven SYK chain. Empty circles correspond to the diffusion
constant computed using the ansatz in Eq. (72). The insets show an exponential
scaling $D\sim e^{-\nu\omega}$ at high frequencies.

In the boundary-driven case, the inverse temperatures of the two baths at the
ends oscillate out of phase by a half-period, thus always creating an energy
gradient and current, whose direction changes in time. A similar configuration
has been used to study magnetization transport in a quantum spin chain subject
to an oscillating magnetic field bias at the boundaries [81]. Our results for
a sample chain driven at $\omega=5$ are shown in Fig. 7 and similar outcomes
are seen for all regimes of parameters. All local observables experience an
oscillatory time-dependence. In Fig. 7(a-c) we plot the oscillations of the
bond energy, inverse temperature, and current around their equilibrium values
for three consecutive sites in the bulk. An interesting feature is that the
amplitude of these oscillations is small and decays at a steady rate $c$ as we
approach the middle of the chain from either end. Moreover, we observe a
constant phase shift $\delta$ from site to site. This indicates an attenuation
in the bulk of the wave modes created at the boundary.

Notice that the current and energy gradient oscillate approximately in phase
at every site. Therefore, we can use the ratio of their amplitudes, which is
independent of the dampening $c$, to define a diffusion constant
$D(\omega)=|j(\omega)/\nabla E(\omega)|$ at each site. Even though the fast
decay prevents us from numerically studying large systems, we checked that the
diffusion constant converges rapidly with system size and is uniform
throughout the chain. We plot our results for $D(\omega)$ at different
temperatures in Fig. 8. Note that in the zero-frequency limit, we recover the
DC values of the diffusion constants from Fig. 4(a). Our simulations are
limited by the time discretization step at large frequencies, where we have to
keep $\mathop{}\\!\mathrm{d}t$ small to avoid sub-sampling errors, while also
evolving the KB equations long enough to achieve convergence. Nonetheless,
even in the limited frequency range attainable numerically we observe a clear
exponential decay $D\sim e^{-\nu\omega}$ at $\omega\gg J,\beta^{-1}$. The
exponent $\nu$ is small and decreases with $\beta$. Moreover, we notice that
$\nu$ approaches a $q$-independent value at low temperatures, suggesting a
possible universal scaling in the IR.

The exponential decay of the diffusion constant at large frequencies can be
explained using recent developments in the theory of periodically-driven many-
body systems [82, 83, 84, 85]. Intuitively, the energy exchange between the
system and the bath (or its neighbors) becomes inefficient at high
frequencies, since it requires significant changes in the many-body state in a
very short time [42]. The amount of energy outputted by the bath during
periodic driving is on the order of $\omega$, while the energy that can be
absorbed by a local degree of freedom (fermion) is of order $J$. At high
frequencies, many fermions (roughly $\omega/J$) would have to work together to
absorb this energy over one period. However, such a many-body effect is
exponentially suppressed for a Hamiltonian with few-body interactions [85].
Therefore, the rate of energy exchange between clusters becomes exponentially
small, which in turn leads to an exponentially small diffusion coefficient.

###  V.2 Non-equilibrium ansatz

Similarly to the DC case, we use our observations about the NESS to devise an
ansatz for the local non-equilibrium Green’s function under periodic driving.
For convenience, we perform a change of variables from $(t_{1},t_{2})$ to
$t=t_{1}-t_{2}$ and $T=(t_{1}+t_{2})/2$. At late times, the Green’s functions
are no longer time-translation invariant and depend on $T$ in addition to $t$.
However, we find that $T$ only enters the oscillatory component of the Green’s
function. The most general ansatz consistent with these observations and the
symmetries of the Green’s function (see Eq. (27)) is given by

$\begin{split}G_{x}^{>}(t,T)&=G^{(0)}(t)+c^{x}F(t)\cos(\omega
T-x\delta+\phi(t))\\\ &\equiv G^{(0)}(t)+c^{x}\left(F_{1}(t)\cos(\omega
T-x\delta)+F_{2}(t)\sin(\omega T-x\delta)\right),\end{split}$ (72)

where we allowed for a time-dependent phase $\phi(t)$ in addition to the
relative phase $\delta$ between sites. The parameter $c$ controls the ratio of
amplitudes between neighboring sites. We find it easier to work directly with
the sine and cosine components $F_{1,2}(t)$, rather than the phase $\phi(t)$.
In the linear response regime, we have $|F_{1,2}(t)|\ll|G^{(0)}(t)|$.

We can use this guess to write a set of KB equations for the non-equilibrium
functions $F_{1,2}(t)$. However, the additional factors stemming from the AC
driving prevent us from finding an explicit solution in terms of the
equilibrium Green’s functions. Nevertheless, we can solve the system of
equations iteratively to find both $F_{1,2}(t)$ and $(c,\delta)$, as detailed
in Appendix D. A comparison between $F_{1,2}(t)$ extracted directly from the
NESS Green’s functions and those computed numerically in Appendix D is shown
in Fig. 7(d). We see a very good agreement between the two. Furthermore, we
can use the ansatz to compute $D(\omega)$ at various temperatures. We observe
remarkable agreement with the previous data in Fig. 8. The slight deviations
at high frequencies are due to the extremely small amplitudes of the current
and energy gradient in the case of boundary driving.

###  V.3 Bulk-driven SYK chain

Most of the difficulties related to the computation of $D(\omega)$ in the
previous two sections arise from the dampening of the oscillations in the
bulk. In particular, this has led us to introduce two additional parameters,
$c$ and $\delta$, in our non-equilibrium ansatz, which further complicated the
calculations. Therefore, it seems natural to consider adding baths in the
bulk. Indeed, this would lead to more uniform observables across the chain,
corresponding to the limit $c\to 1$ and $\delta\to 0$. Note that a small, but
finite $\delta$ is still required to have non-zero current and can be imposed
by staggering the baths. However, a new problem emerges at large frequencies.
In this regime, there will be a significant amplitude decay and phase shift
between the oscillations in the bath and those in the system at every site.
Since the KB equations in the bulk now involve the bath’s Green’s function, we
would have to take into account both the NESS contributions from the system
$F_{S}(t)$, as well as those from the bath $F_{B}(t)$ due to the time-
dependent temperature. The two functions would be related by a new set of
parameters $(c^{\prime},\delta^{\prime})$ that need to be determined
numerically. Therefore, the attenuation problem persists even in the case of
bulk driving and we do not get a clear advantage in this case. It is worth
noting that in the case of time-independent bulk driving, the bath at fixed
temperature did not have a NESS correction $F(t)$ to its Green’s function, and
hence its influence could be fully captured through renormalizing the
effective coupling $J$.

##  VI Discussion

In this paper, we discussed a framework for studying the non-equilibrium
properties of a generalized SYK model on arbitrary graphs coupled to thermal
baths. We then used this formalism to investigate energy transport in SYK
lattices with various $q$-body interactions. Our main focus was on one-
dimensional chains, although we also discussed generalizations to higher
dimensions. We showed that the transport coefficients can be equivalently
computed either by applying a small bias at the boundaries of the system, or
by attaching baths to every site in the bulk. Using this setup, we verified
that all the models feature diffusive transport, and numerically computed the
temperature dependence of the diffusion constants. At low temperatures, we
were able to show that the diffusivity approaches a constant value which
exactly matches the conformal limit prediction [43]. We also showed that
energy diffusion is upper bounded by chaos at all temperatures $D\leq
v_{B}^{2}/\lambda_{L}$, with equality holding in the conformal regime. It
remains to be seen whether this inequality holds for other generalizations of
the SYK model.

Our analysis of the non-equilibrium Green’s functions that emerge during
transport has revealed that they are only weakly perturbed from their
equilibrium values. Moreover, we were able to identify the exact functional
form of this perturbation for both frequency-independent and frequency-
dependent transport. In the DC case, we were able to derive a series of
analytical results in the large $q$ expansion of the SYK model. In particular,
we obtained a closed-form expression for the diffusion constant at all
temperatures in the $q\to\infty$ limit. In the AC case, we managed to solve
for the non-equilibrium contributions numerically and showed that the
diffusion constant decays exponentially at high frequencies. In both cases,
our non-equilibrium ansatz allowed for a very efficient computation of the
transport properties of an SYK chain without referencing the specific setup
driving the system out of equilibrium. This result is not surprising, since we
know that the transport coefficients should not depend on the details of the
drive, as long as we are in the linear response regime. Note that although we
focused on linear-response-like scenarios in this work, nothing in the general
formalism requires this, and it would be interesting to investigate the non-
linear response of such SYK networks due to strong driving.

The study of non-equilibrium SYK chains uncovered important insights about the
structure of NESS, which can be further applied to transport in more
conventional spin systems [1, 2, 3, 4]. First, we showed that adding baths in
the bulk dramatically improved the convergence time of our simulations. This
could be a key ingredient necessary to reach lower temperatures in certain
spin systems, where convergence time becomes the limiting factor [4]. Second,
we argued that the NESS can be described locally using the non-equilibrium
contribution $F$, which acts as a small perturbation on top of the equilibrium
solution. A similar approach could be applicable to spin systems where the
state is represented as a tensor network [4, 1]. One would proceed by
expanding the master equation to first order in $F$ and then solving for the
non-equilibrium contribution, in the spirit of Appendix D. Finally, it would
be interesting to study AC transport in spin chains using a periodic drive,
similar to the one described in Sec. V, and compare the results to those
obtained via the Kubo formula [81, 80, 1].

One natural extension of our analysis involves thermoelectric transport in
complex SYK models [14, 44, 45, 66, 62]. This variant of the SYK Hamiltonian
is written in terms of complex fermions and features a conserved fermion
number in addition to energy [14, 44]. Diffusive charge and energy transport
has been observed in the strongly correlated metals built from complex SYK
clusters [44, 45]. They also have other interesting properties, such as a
linear in temperature resistivity reminiscent of high-$T_{c}$ cuprates [45].
It would be instructive to study mixed thermal and electrical transport in
these models by explicitly coupling them to baths in various ways.

Given the interesting physical properties of the SYK model and its extensions,
multiple experimental realizations [58, 59, 60, 61] and quantum simulations
[86, 87] of SYK have been proposed. However, there are a few limitations that
experiments must overcome. First, all-to-all interactions and the large $N$
limit, which are crucial for exactly solving the SYK model, are hard to
achieve in real materials. Second, depending on the Hamiltonian under study,
one has to suppress the free-fermion hopping terms, while enhancing the
higher-order $q$-body random interactions between the fermions [56, 57]. To
circumvent these issues, experimental setups in solid-state systems that rely
on approximate symmetries to forbid two-body interactions, but allow higher-
order couplings, have been proposed. For example, implementations with
Majorana modes on the surface of a topological insulator [59] and
semiconductor quantum wires coupled to a quantum dot [60] have been discussed.
Additionally, graphene flakes with irregular boundaries in the presence of
strong magnetic fields have been suggested as realizations of the complex SYK
model [61, 42]. These experiments provide a promising path towards studying
out-of-equilibrium dynamics in SYK models. Actual measurements of the
diffusion constants could potentially be performed on these devices in the
near future.

###### Acknowledgements.

We would like to thank Richard Davison and Julia Steinberg for valuable
discussions. We are also grateful to Mark Mezei for pointing out a correction
to our chaos bound in an earlier version of the manuscript. 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 the AFOSR
under grant number 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).
  * 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).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Leviatan  _et al._ [2017] E. Leviatan, F. Pollmann, J. H. Bardarson, D. A. Huse, and E. Altman, Quantum thermalization dynamics with matrix-product states (2017), arXiv:1702.08894 .
  * White  _et al._ [2018] C. D. White, M. Zaletel, R. S. K. Mong, and G. Refael, Phys. Rev. B 97, 035127 (2018).
  * Rakovszky  _et al._ [2020] T. Rakovszky, C. W. von Keyserlingk, and F. Pollmann, Dissipation-assisted operator evolution method for capturing hydrodynamic transport (2020), arXiv:2004.05177 .
  * Polkovnikov [2010] A. Polkovnikov, Ann. Phys. 325, 1790 (2010).
  * Schachenmayer  _et al._ [2015] J. Schachenmayer, A. Pikovski, and A. M. Rey, Phys. Rev. X 5, 011022 (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).
  * 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).
  * 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).
  * 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).
  * Maldacena  _et al._ [2016a] J. Maldacena, S. H. Shenker, and D. Stanford, J. High Energy Phys. 2016 (8), 106.
  * Jackiw [1985] R. Jackiw, Nucl. Phys. B 252, 343 (1985).
  * Teitelboim [1983] C. Teitelboim, Phys. Lett. B 126, 41 (1983).
  * 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.
  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Page [1993] D. N. Page, Phys. Rev. Lett. 71, 3743 (1993).
  * Rocha [2010] J. V. Rocha, J. Phys. Conf. Ser. 222, 012005 (2010).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Sonner and Vielma [2017] J. Sonner and M. Vielma, J. High Energy Phys. 2017 (11), 149.
  * 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).
  * 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.
  * 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).
  * Haldar  _et al._ [2020] A. Haldar, P. Haldar, S. Bera, I. Mandal, and S. Banerjee, Phys. Rev. Research 2, 013307 (2020).
  * Kuhlenkamp and Knap [2020] C. Kuhlenkamp and M. Knap, Phys. Rev. Lett. 124, 106401 (2020).
  * Gu  _et al._ [2017a] Y. Gu, X.-L. Qi, and D. Stanford, J. High Energy Phys. 2017 (5), 125.
  * Davison  _et al._ [2017] R. A. Davison, W. Fu, A. Georges, Y. Gu, K. Jensen, and S. Sachdev, Phys. Rev. B 95, 155131 (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).
  * Jian  _et al._ [2017] C.-M. Jian, Z. Bi, and C. Xu, Phys. Rev. B 96, 115122 (2017).
  * 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).
  * Choi  _et al._ [2021] C. Choi, M. Mezei, and G. Sárosi, J. High Energy Phys. 2021 (2), 207.
  * Gu  _et al._ [2017b] Y. Gu, A. Lucas, and X.-L. Qi, SciPost Phys. 2, 018 (2017b).
  * 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).
  * 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).
  * Chen  _et al._ [2018] A. Chen, R. Ilan, F. de Juan, D. I. Pikulin, and M. Franz, Phys. Rev. Lett. 121, 036403 (2018).
  * Cheipesh  _et al._ [2021] Y. Cheipesh, A. I. Pavlov, V. Ohanesjan, K. Schalm, and N. V. Gnezdilov, Phys. Rev. B 104, 115134 (2021).
  * 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).
  * Lucas and Steinberg [2016] A. Lucas and J. Steinberg, J. High Energy Phys. 2016 (10), 143.
  * Guo  _et al._ [2019] H. Guo, Y. Gu, and S. Sachdev, Phys. Rev. B 100, 045140 (2019).
  * Can  _et al._ [2019] O. Can, E. M. Nica, and M. Franz, Phys. Rev. B 99, 045419 (2019).
  * Gu  _et al._ [2017c] Y. Gu, A. Lucas, and X.-L. Qi, J. High Energy Phys. 2017 (9), 120.
  * Chen  _et al._ [2017] Y. Chen, H. Zhai, and P. Zhang, J. High Energy Phys. 2017 (7), 150.
  * 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.
  * Babadi  _et al._ [2015] M. Babadi, E. Demler, and M. Knap, Phys. Rev. X 5, 041005 (2015).
  * Butcher [2016] J. Butcher,  _Numerical Methods for Ordinary Differential Equations_ (Wiley, New York, 2016).
  * Tarnopolsky [2019] G. Tarnopolsky, Phys. Rev. D 99, 026010 (2019).
  * Zotos  _et al._ [1997] X. Zotos, F. Naef, and P. Prelovsek, Phys. Rev. B 55, 11029 (1997).
  * Zanoci and Swingle [2016] C. Zanoci and B. G. Swingle, Entanglement and thermalization in open fermion systems (2016), arXiv:1612.04840 .
  * Kubo [1957] R. Kubo, J. Phys. Soc. Japan 12, 570 (1957).
  * Trauzettel  _et al._ [2008] B. Trauzettel, P. Simon, and D. Loss, Phys. Rev. Lett. 101, 017202 (2008).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. m. c. Huveneers, Phys. Rev. Lett. 115, 256803 (2015).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. m. c. Huveneers, Phys. Rev. B 95, 014112 (2017).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Phys. Rev. Lett. 116, 120401 (2016).
  * Kuwahara  _et al._ [2016] T. Kuwahara, T. Mori, and K. Saito, Ann. Phys. 367, 96 (2016).
  * 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).
  * 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.
  * Gu and Kitaev [2019] Y. Gu and A. Kitaev, J. High Energy Phys. 2019 (2), 75.

##  Appendix A Energy current

The current flowing from $u$ to $v$ across the edge $(u,v)$ is given by

$j_{u}=i\sum_{u^{\prime}}A_{u^{\prime}u}[H^{u^{\prime}u},H^{uv}],$ (73)

where $H^{uv}=H_{0}^{u}/d_{u}+H_{0}^{v}/d_{v}+H_{1}^{uv}$ is the on-bond
Hamiltonian. The commutator simplifies to

$\begin{split}j_{u}&=i\sum_{u^{\prime}\neq
v}A_{u^{\prime}u}\Big{(}\frac{1}{d_{u}}\left([H_{0}^{u},H_{1}^{uv}]-[H_{0}^{u},H_{1}^{uu^{\prime}}]\right)\\\
&+[H_{1}^{u^{\prime}u},H_{1}^{uv}]\Big{)}.\end{split}$ (74)

One can easily check using the anti-commutation relations that the commutator
between a product of $m$ Majorana fermions and a product of $n$ Majorana
fermions is proportional to $1-(-1)^{mn-l}$, where $l$ is the number of
identical fermions in common between the two terms. Hence the commutator
vanishes unless $l$ is odd. Furthermore, in the large $N$ limit, the leading
order contribution comes from the commutator where the terms have a single
fermion in common, with other contributions suppressed by factors of $1/N$.

To compute the non-equilibrium expectation value of any operator $O$ we employ
the generating functional [63]

$\langle O(t)\rangle=\lim\limits_{\eta\to 0}\frac{i}{2}\frac{\delta
Z[\eta]}{\delta\eta(t)}=\lim\limits_{\eta\to
0}\frac{i}{2}\frac{\delta\ln(Z[\eta])}{\delta\eta(t)}=\lim\limits_{\eta\to
0}\frac{i}{2}\frac{\delta iS_{\eta}}{\delta\eta(t)},$ (75)

where $Z[\eta]$ is the generating functional with an additional source term
$\eta(t)O(t)$. We also used the fact that $Z[\eta\to 0]=1$ in the Keldysh
formalism [63], and denoted by $S_{\eta}$ the part of the action that depends
on the source term. The factor of $1/2$ accounts for the fact that $t$ can
belong to either the positive $\mathcal{C}^{+}$ or negative $\mathcal{C}^{-}$
branch of the contour. Without loss of generality, we will assume that $t$
lives on $\mathcal{C}^{+}$ and ignore the factor of $1/2$.

Let us focus on the first commutator in Eq. (74). In order to find $S_{\eta}$,
we have to perform the disorder averaging over couplings in the path integral,
which results in a standard Gaussian integral [19]. Its contribution to the
action is given by

$iS_{\eta}=-\frac{i^{q_{S}+q_{I}}}{2}NJ_{0}^{2}J_{1}^{2}\int_{\mathcal{C}}\mathop{}\\!\mathrm{d}t_{1}\mathop{}\\!\mathrm{d}t_{2}G_{u}(t_{1},t_{2})G_{u}(t_{1},t)^{q_{S}-1}G_{u}(t_{2},t)^{\frac{q_{I}}{2}-1}G_{v}(t_{2},t)^{\frac{q_{I}}{2}}\eta(t).$
(76)

Note that due to the causal structure, we only have to consider times
$t_{1},t_{2}\leq t$. However, they can belong to either of the two contour
branches, resulting in eight different orderings, which we group into four
terms as follows:

$\displaystyle j_{++}^{uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{<}(t_{1},t_{2})\left(G_{u}^{<}(t_{1},t)G_{u}^{<}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{u}^{<}(t_{1},t)^{q_{S}-\frac{q_{I}}{2}}G_{v}^{<}(t_{2},t)^{\frac{q_{I}}{2}}-G_{v}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{u}^{<}(t_{2},t)^{q_{S}-\frac{q_{I}}{2}}\right),$
(77) $\displaystyle j_{+-}^{uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{<}(t_{1},t_{2})\left(G_{u}^{<}(t_{1},t)G_{u}^{>}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{v}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{u}^{>}(t_{2},t)^{q_{S}-\frac{q_{I}}{2}}-G_{u}^{<}(t_{1},t)^{q_{S}-\frac{q_{I}}{2}}G_{v}^{>}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(78) $\displaystyle j_{-+}^{uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{>}(t_{1},t_{2})\left(G_{u}^{>}(t_{1},t)G_{u}^{<}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{v}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{u}^{<}(t_{2},t)^{q_{S}-\frac{q_{I}}{2}}-G_{u}^{>}(t_{1},t)^{q_{S}-\frac{q_{I}}{2}}G_{v}^{<}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(79) $\displaystyle j_{--}^{uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{>}(t_{1},t_{2})\left(G_{u}^{>}(t_{1},t)G_{u}^{>}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{u}^{>}(t_{1},t)^{q_{S}-\frac{q_{I}}{2}}G_{v}^{>}(t_{2},t)^{\frac{q_{I}}{2}}-G_{v}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{u}^{>}(t_{2},t)^{q_{S}-\frac{q_{I}}{2}}\right).$
(80)

Notice that we used the identities in Eq. (23) to write the answers
exclusively in terms of greater or lesser Green’s functions. Additionally,
using Eq. (27), we obtain that $j_{++}^{uv}=(j_{--}^{uv})^{*}$ and
$j_{+-}^{uv}=(j_{-+}^{uv})^{*}$. The expectation value of the commutator is
given by the sum of the four terms above, with the appropriate proportionality
constant

$\langle[H_{0}^{u},H_{1}^{uv}]\rangle=-\frac{i^{q_{S}+q_{I}+1}}{2}J_{0}^{2}J_{1}^{2}(j_{++}^{uv}+j_{+-}^{uv}+j_{-+}^{uv}+j_{--}^{uv})=-i^{q_{S}+q_{I}+1}J_{0}^{2}J_{1}^{2}\real(j_{++}^{uv}+j_{+-}^{uv}).$
(81)

A similar expression can be obtained for $[H_{1}^{u^{\prime}u},H_{1}^{uv}]$ if
we define

$\displaystyle j_{++}^{u^{\prime}uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{<}(t_{1},t_{2})\left(G_{u}^{<}(t_{1},t)G_{u}^{<}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{u^{\prime}}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{v}^{<}(t_{2},t)^{\frac{q_{I}}{2}}-G_{v}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{u^{\prime}}^{<}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(82) $\displaystyle j_{+-}^{u^{\prime}uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{<}(t_{1},t_{2})\left(G_{u}^{<}(t_{1},t)G_{u}^{>}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{v}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{u^{\prime}}^{>}(t_{2},t)^{\frac{q_{I}}{2}}-G_{u^{\prime}}^{<}(t_{1},t)^{\frac{q_{I}}{2}}G_{v}^{>}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(83) $\displaystyle j_{-+}^{u^{\prime}uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{>}(t_{1},t_{2})\left(G_{u}^{>}(t_{1},t)G_{u}^{<}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{v}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{u^{\prime}}^{<}(t_{2},t)^{\frac{q_{I}}{2}}-G_{u^{\prime}}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{v}^{<}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(84) $\displaystyle j_{--}^{u^{\prime}uv}$
$\displaystyle=\int_{-\infty}^{t}\mathop{}\\!\mathrm{d}t_{2}\int_{-\infty}^{t_{2}}\mathop{}\\!\mathrm{d}t_{1}G_{u}^{>}(t_{1},t_{2})\left(G_{u}^{>}(t_{1},t)G_{u}^{>}(t_{2},t)\right)^{\frac{q_{I}}{2}-1}\left(G_{u^{\prime}}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{v}^{>}(t_{2},t)^{\frac{q_{I}}{2}}-G_{v}^{>}(t_{1},t)^{\frac{q_{I}}{2}}G_{u^{\prime}}^{>}(t_{2},t)^{\frac{q_{I}}{2}}\right),$
(85)

where only the index of the terms in brackets changed from $u$ to $u^{\prime}$
and $q_{S}$ was changed to $q_{I}$. The expectation value of the commutator is
now given by

$\langle[H_{1}^{u^{\prime}u},H_{1}^{uv}]\rangle=-\frac{i}{2}J_{1}^{4}(j_{++}^{u^{\prime}uv}+j_{+-}^{u^{\prime}uv}+j_{-+}^{u^{\prime}uv}+j_{--}^{u^{\prime}uv})=-iJ_{1}^{4}\real(j_{++}^{u^{\prime}uv}+j_{+-}^{u^{\prime}uv})$
(86)

Putting everything together using Eq. (74), we arrive at the following
expression for the current:

$\langle j_{u}\rangle=J_{1}^{2}\sum_{u^{\prime}\neq
v}A_{u^{\prime}u}\left(i^{q_{S}+q_{I}}\frac{J_{0}^{2}}{d_{u}}\real(j_{++}^{uv}+j_{+-}^{uv})-i^{q_{S}+q_{I}}\frac{J_{0}^{2}}{d_{u}}\real(j_{++}^{uu^{\prime}}+j_{+-}^{uu^{\prime}})+J_{1}^{2}\real(j_{++}^{u^{\prime}uv}+j_{+-}^{u^{\prime}uv})\right).$
(87)

It is important to mention that in non-equilibrium situations where the
Green’s functions are only weakly perturb from their equilibrium values
$G_{u}^{>}=G^{(0)}+F_{u}$ with $|F_{u}|\ll|G^{(0)}|$ and $q_{S}=q_{I}$, there
are significant simplifications to the formula above

$\displaystyle j_{+\pm}^{u^{\prime}uv}$
$\displaystyle=j_{+\pm}^{uv}-j_{+\pm}^{uu^{\prime}},$ (88)
$\displaystyle\langle j_{u}\rangle$
$\displaystyle=\frac{J_{1}^{2}(J_{0}^{2}+d_{u}J_{1}^{2})}{d_{u}}\sum_{u^{\prime}\neq
v}A_{u^{\prime}u}\real(j_{++}^{u^{\prime}uv}+j_{+-}^{u^{\prime}uv}).$ (89)

We verified numerically that this condition applies to all the non-equilibrium
setups studied in this paper and therefore use Eq. (89) throughout, except for
the cases when $q_{S}\neq q_{I}$, where the more general Eq. (87) applies.

##  Appendix B Chaos bound on diffusion

In this section, we review the many-body chaos properties of the SYK model.
Chaos is fundamentally related to energy fluctuations and diffusion, since the
local energy characterizes the rate of change of the quantum phase, and phase
decoherence leads to chaos [44]. 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 [88, 89, 90, 21, 17, 91]. Following [17, 43], we define the
regularized OTOC in real time

$C(x,t_{1},t_{2})=\frac{1}{N^{2}}\sum_{i,j=1}^{N}\Tr[y\psi_{j}^{x}(t_{1})y\psi_{i}^{0}(0)y\psi_{j}^{x}(t_{2})y\psi_{i}^{0}(0)],$
(90)

where $y=e^{-\beta H/4}/Z^{1/4}$ evenly spaces the fermionic fields along the
thermal circle. The OTOC measures how sensitive the system is to an initial
local perturbation created by the operator $\psi_{i}^{0}(0)$, thus
characterizing a quantum analog for the butterfly effect [89]. The leading
order contribution to the OTOC comes from a time independent constant equal to
the disconnected correlator $\mathcal{F}_{d}$ [92]. The next order
contribution comes from contracting $\psi_{i}$ with $\psi_{j}$ and is of order
$1/N$ at early times

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

For chaotic systems with a large hierarchy between thermalization and
scrambling, we expect $\mathcal{F}$ to grow exponentially as
$e^{\lambda_{L}t}$, where $t_{1}=t_{2}=t$ is in the Lyapunov regime
$\beta\lesssim t\lesssim\beta\ln N$ [92]. Here $\lambda_{L}$ is the Lyapunov
exponent determining the scrambling rate [90]. For a single SYK cluster (with
no spatial dependence), $\mathcal{F}(t_{1},t_{2})$ is determined by summing
over a set of ladder diagrams [17], leading to 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}),$
(92)

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^{R}(t_{2}-t_{4})G^{W}(t_{3}-t_{4})^{q-2}.\end{split}$
(93)

The functions $G^{R}$ and $G^{W}$ are the retarded and Wightman Green’s
functions. The Wightman propagator is related to the spectral function in
frequency space [66]

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

To determine the Lyapunov exponent, we follow the prescription introduced in
Ref. [93]. 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).$
(95)

We can view this operator as a matrix with its largest eigenvalue denoted by
$k_{R}(\alpha)$. Then the Lyapunov exponent is determined by the equation
$k_{R}(-\lambda_{L})=1$. This is equivalent to solving for $\mathcal{F}$ as an
eigenvector of the kernel $K^{R}$ with eigenvalue one. The $(0+1)$-d SYK model
is known to saturate a bound on the Lyapunov exponent at low temperatures
$\lambda_{L}=2\pi/\beta$ [17, 21].

The one-dimensional SYK chain allows us to study the chaos dynamics in space.
To characterize the spatial propagation in a translation-invariant system, it
is convenient to first introduce the Fourier transform
$\mathcal{F}(p,\omega)$. As derived in Ref. [66] using the ladder identity,
the OTOC has a pole in both frequency and momentum space

$\mathcal{F}(p,\omega)\sim\frac{1}{\cos(\lambda_{L}(p)\beta/4)}\frac{1}{\omega-i\lambda_{L}(p)},$
(96)

where $\lambda_{L}(p)$ is the momentum-dependent Lyapunov exponent. By
performing an inverse Fourier transform back to real space, we find [66]

$\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)}.$
(97)

At large distances and times, 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}|$ [93, 66]. If
$|p_{s}|<|p_{1}|$, the OTOC receives a contribution only from the saddle point

$\mathcal{F}(x,t)\sim
e^{\lambda_{L}(p_{s})t+ip_{s}|x|}=e^{\lambda_{L}(p_{s})(t-|x|/v_{s})},$ (98)

where $v_{s}\equiv\lambda_{L}(p_{s})/|p_{s}|$. Conversely, if
$|p_{s}|>|p_{1}|$, the OTOC receives a dominant contribution from the pole,
resulting in a wave-front that propagates with a maximal chaos rate
$\lambda_{L}(p_{1})=2\pi/\beta$

$\mathcal{F}(x,t)\sim
e^{\lambda_{L}(p_{1})t+ip_{1}|x|}=e^{\frac{2\pi}{\beta}(t-|x|/v_{1})},$ (99)

where $v_{1}=2\pi/\beta|p_{1}|$. We can now define the butterfly velocity [93,
66, 54]

$v_{B}=\begin{cases}v_{s}&\mbox{if }|p_{s}|<|p_{1}|\\\ v_{1}&\mbox{if
}|p_{s}|>|p_{1}|\end{cases}$ (100)

Physically, $v_{B}$ represents the growth rate of the region where operators
have large anti-commutators with the initial $\psi_{i}^{0}(0)$. We can use the
butterfly velocity to define two chaos diffusion constants, $\beta
v_{B}^{2}/2\pi$ and $v_{B}^{2}/\lambda_{L}$, which are known to be closely
related to the energy diffusion constant in strange metals [50, 51, 52, 66,
49, 44, 55, 43, 54, 53].

Figure 9: Temperature dependence of the energy diffusion constants (circles)
and the chaos bounds $\beta v_{B}^{2}/2\pi$ (solid lines) and
$v_{B}^{2}/\lambda_{L}$ (dashed lines) for different values of $q$. All three
quantities agree in the conformal limit.

In the remainder of this section, we will compute the value of $v_{B}$ and
show that although $\beta v_{B}^{2}/2\pi$ closely follows the energy diffusion
coefficient $D$ at low temperatures, the larger quantity
$v_{B}^{2}/\lambda_{L}$ provides a true upper bound on energy transport. Note
that we only have to find the momenta $p_{s,1}$. We follow the procedure in
Ref. [93]. For the SYK chain, the retarded kernel factorizes in momentum space
$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 [43] and $K^{R}$ is the kernel for a single cluster with
effective coupling $J$ defined in Eq. (93). Hence the eigenvalues of the
kernel simply get rescaled $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$. The location of the saddle can 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$. Note that for $q=2$, the model is non-
chaotic and the Lyapunov exponent vanishes [35, 36].

We numerically diagonalize the kernel and compute $\beta v_{B}^{2}/2\pi$ and
$v_{B}^{2}/\lambda_{L}$ as a function of temperature. Our results are shown in
Fig. 9. At infinite temperature, both $v_{B}$ and $\lambda_{L}$ approach a
constant. Therefore, their ratio $v_{B}^{2}/\lambda_{L}$ also approaches a
constant, while $\beta v_{B}^{2}/2\pi$ decays to zero. At high temperatures,
we have $D>\beta v_{B}^{2}/2\pi$, which was previously observed in a large $q$
expansion of this model [54]. At low temperatures, on the other hand, we find
$\lambda_{L}=2\pi/\beta$ and all three quantities converge to the same
diffusion constant $D=v_{B}^{2}/\lambda_{L}=\beta v_{B}^{2}/2\pi$, as was
previously shown analytically [43]. This indicates that the SYK chain is
maximally chaotic in the conformal limit. This remarkable result is a
consequence of the fact that the same reparameterization degrees of freedom
are responsible both for energy diffusion and the OTOC chaos dynamics [17,
43].

Our results also show that $D\leq v_{B}^{2}/\lambda_{L}$ at all temperatures,
suggesting that chaos upper bounds energy diffusion. Moreover, we see that the
weaker bound $D\leq\beta v_{B}^{2}/2\pi$ is violated at high temperatures
[54]. In was previously conjectured that there is a fundamental lower bound on
transport in incoherent metals $D\geq v^{2}\tau$, where $v$ and $\tau$ are the
characteristic velocity and relaxation time respectively [49]. It was later
suggested that many-body chaos provides a natural choice for these quantities:
$v=v_{B}$ and $\tau=1/\lambda_{L}$ [50, 51, 52]. However, we see that the
inequality is precisely reversed in our case, as was also found in other
systems [55, 65, 44, 54]. This indicates that $v_{B}$ and $\lambda_{L}$ are
not always the appropriate scales, and there is no simple bound relating
transport and chaos in all incoherent metals [53].

##  Appendix C Exact calculations of the diffusion constant

In Sec. IV.3 we introduced a simple ansatz for the non-equilibrium
contribution to the Green’s function
$F(t)=\mathop{}\\!\mathrm{d}G^{(0)}(t)/\mathop{}\\!\mathrm{d}\beta$. Both the
current and energy gradient have a relatively simple form in terms of $F(t)$.
In what follows, we will consider the special cases where it is possible to
compute the equilibrium Green’s function $G^{(0)}(t)$, and hence $F(t)$,
analytically. This will lead to closed-form expressions for the diffusion
constant $D$ in various limits. In particular, we will consider the case of
$q=2$, as well as the large $q$ approximation, in the limit of zero and
infinite temperatures. For the zero-temperature limit, we will recover
precisely the conformal answer in Eq. (55). Our results are summarized in Tab.
1.

Table 1:  Diffusion constant computed exactly in various limits. Limit | $q=2$ | large $q$ | $q=\infty$
---|---|---|---
$\beta\to 0$ | $\frac{32}{15\pi}\frac{J_{1}^{2}}{J}$ | $\frac{q}{q-2}\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}$ | $\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}$
$\beta\to\infty$ | $\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}$ | $\frac{q}{q-2}\frac{\pi}{6}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}$ | $\frac{\pi}{6}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}$

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

The SYK Hamiltonian for $q=2$ is equivalent to a random hopping model of free
Majorana fermions, which can be solved exactly [17, 37]. The spectral function
is given by

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

In equilibrium, the greater Green’s function can be obtained from the
fluctuation-dissipation theorem

$G^{(0)}(\omega)=\frac{A(\omega)}{1+e^{-\beta\omega}},$ (102)

followed by an inverse Fourier transform

$\begin{split}G^{(0)}(t)&=\int_{-\infty}^{\infty}\frac{\mathop{}\\!\mathrm{d}\omega}{2\pi}e^{-i\omega
t}G^{(0)}(\omega)\\\ &=-\frac{i}{2Jt}B(2Jt)-\frac{1}{\pi
J}\int_{-2J}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{\sin(\omega
t)}{1+e^{-\beta\omega}}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}},\end{split}$
(103)

where $B$ is the Bessel function of the first kind. Finally, we arrive at

$F(t)=-\frac{2}{\pi
J}\int_{0}^{2J}\mathop{}\\!\mathrm{d}\omega\frac{\omega\sin(\omega
t)e^{-\beta\omega}}{(1+e^{-\beta\omega})^{2}}\sqrt{1-\left(\frac{\omega}{2J}\right)^{2}}.$
(104)

For $q=2$, the energy gradient in Eq. (59) takes the form

$\nabla E=J\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\frac{F(t)}{t}B(2Jt),$
(105)

and the current in Eq. (61) becomes

$j=-J_{1}^{2}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\frac{F(t^{\prime})}{t(t^{\prime}-t)}B(2Jt)B(2J(t^{\prime}-t)).$
(106)

These integrals can be performed analytically in the $\beta\to\infty$ limit.
To leading order in $1/\beta$, we find $\nabla E=-\pi/(6J\beta^{3})$ and
$j=\pi J_{1}^{2}/(6J^{2}\beta^{3})$, concluding that

$D=\frac{J_{1}^{2}}{J}\quad(q=2,\beta\to\infty).$ (107)

This matches exactly the conformal answer in Eq. (55), since $\alpha_{K}=\pi$
for the $q=2$ theory [17].

On the other hand, in the $\beta\to 0$ limit, we find $\nabla E=-J^{2}/8$ and
$j=4J_{1}^{2}J/15\pi$. This leads to

$D=\frac{32}{15\pi}\frac{J_{1}^{2}}{J}\quad(q=2,\beta\to 0),$ (108)

in agreement with the results presented in Fig. 4(a).

###  C.2 $q\to\infty$ limit

We now return to the large $q$ analysis of Sec. IV.4 and take the limit of
infinite $q$, while keeping $\beta$ arbitrary. We already have an expression
for $\nabla E$ in Eq. (66), so we only have to compute the current. When we
take $q\to\infty$, while keeping $\mathcal{J}_{0}$ and $\mathcal{J}_{1}$
constant, Eq. (67) simplifies to

$j=\frac{\mathcal{J}_{1}^{2}\mathcal{J}^{2}}{2q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(\,\mathrm{Re}\left[e^{g^{(0)}(t)}f(t)\right]\,\mathrm{Im}\left[e^{g^{(0)}(t^{\prime})}\right]-\,\mathrm{Re}\left[e^{g^{(0)}(t)}\right]\,\mathrm{Im}\left[e^{g^{(0)}(t^{\prime})}f(t^{\prime})\right]\right).$
(109)

We will use the ODE for $f(t)$ (Eq. (64)) and the explicit form of
$g^{(0)}(t)$ (Eq. (42)) repeatedly to simplify the equation above. First,
using the differential equation, we find

$\displaystyle\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\,\mathrm{Im}\left[e^{g^{(0)}(t^{\prime})}\right]$
$\displaystyle=-\frac{1}{2\mathcal{J}^{2}}\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\,\mathrm{Im}\left[g^{(0)}(t^{\prime})^{\prime\prime}\right]=\frac{\,\mathrm{Im}\left[g^{(0)}(t)^{\prime}\right]}{2\mathcal{J}^{2}},$
(110)
$\displaystyle\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\,\mathrm{Im}\left[e^{g^{(0)}(t^{\prime})}f(t^{\prime})\right]$
$\displaystyle=-\frac{1}{2\mathcal{J}^{2}}\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\,\mathrm{Im}\left[f^{\prime\prime}(t^{\prime})\right]=\frac{\,\mathrm{Im}\left[f^{\prime}(t)\right]}{2\mathcal{J}^{2}}.$
(111)

Next, we perform an integration by parts using the fact that
$\,\mathrm{Re}\left[g^{(0)}(0)^{\prime}\right]=0$ and
$\,\mathrm{Im}\left[f^{\prime}(0)\right]=0$

$\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Re}\left[e^{g^{(0)}(t)}\right]\,\mathrm{Im}\left[f^{\prime}(t)\right]=-\frac{1}{2\mathcal{J}^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Re}\left[g^{(0)}(t)^{\prime\prime}\right]\,\mathrm{Im}\left[f^{\prime}(t)\right]=-\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Re}\left[g^{(0)}(t)^{\prime}\right]\,\mathrm{Im}\left[e^{g^{(0)}(t)}f(t)\right].$
(112)

Plugging this back into Eq. (109) we get

$\begin{split}j&=\frac{\mathcal{J}_{1}^{2}}{4q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\left(\,\mathrm{Re}\left[e^{g^{(0)}(t)}f(t)\right]\,\mathrm{Im}\left[g^{(0)}(t)^{\prime}\right]+\,\mathrm{Re}\left[g^{(0)}(t)^{\prime}\right]\,\mathrm{Im}\left[e^{g^{(0)}(t)}f(t)\right]\right)\\\
&=\frac{\mathcal{J}_{1}^{2}}{4q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[e^{g^{(0)}(t)}g^{(0)}(t)^{\prime}f(t)\right]=-\frac{\mathcal{J}_{1}^{2}}{4q^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,\mathrm{Im}\left[e^{g^{(0)}(t)}f^{\prime}(t)\right],\end{split}$
(113)

where in the last step we used integration by parts and
$e^{g^{(0)}(\infty)}=f(0)=0$. The last integral can be evaluated by
substituting the explicit formula for $f(t)$ from Eq. (65)

$j=\frac{\mathcal{J}_{1}^{2}}{12q^{2}}\left(\frac{\pi v}{2}\sin\left(\pi
v\right)+2\cos^{2}\left(\frac{\pi v}{2}\right)\right).$ (114)

Finally, combining this with Eq. (66) yields

$D=\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\left(\frac{\pi v}{2}\sin(\frac{\pi
v}{2})+\cos(\frac{\pi v}{2})\right)\quad(q\to\infty).$ (115)

Furthermore, we can investigate the different temperature limits. If $\beta\to
0$, then $v\to 0$ and

$D=\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\quad(q\to\infty,\beta\to 0).$
(116)

On the other hand, if $\beta\to\infty$, then $v\to 1$ and

$D=\frac{\pi}{6}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}\quad(q\to\infty,\beta\to\infty).$
(117)

This again agrees with the conformal limit in Eq. (55), since $\alpha_{K}=3$
at infinite $q$ [17].

###  C.3 Finite $q$ corrections

In our derivation above, we took the infinite $q$ limit first, followed by a
temperature limit. However, we can also reverse the order to find finite $q$
corrections to the zero- and infinite-temperature diffusion constants. For
instance, consider taking the limit $\beta\to 0$ first. In this case, $v\to 0$
as well, but their ratio approaches a constant $v/\beta\to\mathcal{J}/\pi$. To
leading order, we approximate

$\displaystyle e^{g^{(0)}(t)}$
$\displaystyle=\frac{1}{\cosh^{2}(\mathcal{J}t)},$ (118) $\displaystyle f(t)$
$\displaystyle=-i\tanh(\mathcal{J}t).$ (119)

We now have an elementary solution for $f(t)$, which allows us to directly
evaluate the integral in Eq. (67)

$j=\frac{\mathcal{J}_{1}^{2}}{2q(q-2)}\left(1-\frac{\sqrt{\pi}}{2}\frac{\Gamma(2-2/q)}{\Gamma(5/2-2/q)}\right),$
(120)

where $\Gamma$ denotes the gamma function. Combining this with $\nabla
E=\mathcal{J}/2q^{2}$, we arrive at

$D=\frac{q}{q-2}\frac{\mathcal{J}_{1}^{2}}{3\mathcal{J}}\quad(\beta\to 0),$
(121)

to leading order in $1/q$. A similar calculation for the $\beta\to\infty$
limit results in the same prefactor

$D=\frac{q}{q-2}\frac{\pi}{6}\frac{\mathcal{J}_{1}^{2}}{\mathcal{J}}\quad(\beta\to\infty).$
(122)

These results match the numerical values in Fig. 4(b).

##  Appendix D Solution of the KB equations with non-equilibrium AC ansatz

In this section, we present the KB equations for the non-equilibrium
contributions to the Green’s functions in the case of AC transport and discuss
how to solve them numerically. We start by substituting the ansatz in Eq. (72)
into Eq. (33), and expanding to first order in $F_{1,2}(t)$. All the functions
can be written in terms of the new time variables $t$ and $T$. Moreover, after
a change of variables, all the dependence on $T$ factors out into $\cos(\omega
T)$ and $\sin(\omega T)$ terms. Since the KB equations hold for arbitrary $T$,
the prefactors in front of $\cos(\omega T)$ and $\sin(\omega T)$ must vanish
independently. Therefore, each KB equation leads to a set of two new equations
as follows:

$\displaystyle\partial_{t}F_{1}(t)+\frac{\omega}{2}F_{2}(t)$
$\displaystyle=2i^{q}J^{2}I_{1}(F_{1},F_{2}),$ (123)
$\displaystyle\partial_{t}F_{2}(t)-\frac{\omega}{2}F_{1}(t)$
$\displaystyle=2i^{q}J^{2}I_{1}(F_{2},-F_{1}),$ (124)
$\displaystyle-\partial_{t}F_{1}(t)+\frac{\omega}{2}F_{2}(t)$
$\displaystyle=2i^{q}J^{2}I_{2}(F_{1},F_{2}),$ (125)
$\displaystyle-\partial_{t}F_{2}(t)-\frac{\omega}{2}F_{1}(t)$
$\displaystyle=2i^{q}J^{2}I_{2}(F_{2},-F_{1}),$ (126)

where we define the integrals

$\displaystyle I_{1}(F_{1},F_{2})$
$\displaystyle=-\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\cos(\omega
t^{\prime}/2)+b\sin(\omega
t^{\prime}/2)\right)G^{(0)}(t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t-t^{\prime})^{q-2}F_{1}(t-t^{\prime})\right]$
$\displaystyle+\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\cos(\omega
t^{\prime}/2)-b\sin(\omega
t^{\prime}/2)\right)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})G^{(0)}(t+t^{\prime})^{q-2}F_{1}(t+t^{\prime})\right]$
$\displaystyle-\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\cos(\omega
t^{\prime}/2)F_{1}(t-t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}\right]+\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\cos(\omega
t^{\prime}/2)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}F_{1}(t^{\prime}-t)\right]$
$\displaystyle-\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\sin(\omega
t^{\prime}/2)-b\cos(\omega
t^{\prime}/2)\right)G^{(0)}(t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t-t^{\prime})^{q-2}F_{2}(t-t^{\prime})\right]$
$\displaystyle-\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\sin(\omega
t^{\prime}/2)+b\cos(\omega
t^{\prime}/2)\right)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})G^{(0)}(t+t^{\prime})^{q-2}F_{2}(t+t^{\prime})\right]$
$\displaystyle+\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\sin(\omega
t^{\prime}/2)F_{2}(t-t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}\right]-\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\sin(\omega
t^{\prime}/2)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}F_{2}(t^{\prime}-t)\right],$
(127)

$\displaystyle I_{2}(F_{1},F_{2})$
$\displaystyle=\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\cos(\omega
t^{\prime}/2)-b\sin(\omega
t^{\prime}/2)\right)G^{(0)}(t-t^{\prime})^{q-2}F_{1}(t-t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t^{\prime})\right]$
$\displaystyle-\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\cos(\omega
t^{\prime}/2)-b\sin(\omega
t^{\prime}/2)\right)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})G^{(0)}(t^{\prime}-t)^{q-2}F_{1}(t^{\prime}-t)\right]$
$\displaystyle+\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\cos(\omega
t^{\prime}/2)G^{(0)}(t^{\prime})^{q-1}\,\mathrm{Im}\left[F_{1}(t-t^{\prime})\right]-\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\cos(\omega
t^{\prime}/2)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}F_{1}(t+t^{\prime})\right]$
$\displaystyle-\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\sin(\omega
t^{\prime}/2)+b\cos(\omega
t^{\prime}/2)\right)G^{(0)}(t-t^{\prime})^{q-2}F_{2}(t-t^{\prime})\,\mathrm{Im}\left[G^{(0)}(t^{\prime})\right]$
$\displaystyle+\int_{t}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\left(a\sin(\omega
t^{\prime}/2)+b\cos(\omega
t^{\prime}/2)\right)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})G^{(0)}(t^{\prime}-t)^{q-2}F_{2}(t^{\prime}-t)\right]$
$\displaystyle+\int_{0}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\sin(\omega
t^{\prime}/2)G^{(0)}(t^{\prime})^{q-1}\,\mathrm{Im}\left[F_{2}(t-t^{\prime})\right]+\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t^{\prime}\sin(\omega
t^{\prime}/2)\,\mathrm{Im}\left[G^{(0)}(t^{\prime})^{q-1}F_{2}(t+t^{\prime})\right].$
(128)

Note that the second equation in each set can be obtained from the first by
changing $F_{1}\to F_{2}$ and $F_{2}\to-F_{1}$. To shorten the notation, we
also introduced the constants

$\displaystyle a$
$\displaystyle=(q-1)+q\frac{J_{1}^{2}}{J^{2}}\left(\frac{1}{2}\left(c+\frac{1}{c}\right)\cos\delta-1\right),$
(129) $\displaystyle b$
$\displaystyle=\frac{q}{2}\frac{J_{1}^{2}}{J^{2}}\left(c-\frac{1}{c}\right)\sin\delta.$
(130)

Figure 10: Frequency dependence of the (a) amplitude decay and (b) phase shift
for the AC NESS solution. Filled circles represent numerical values extracted
directly from the boundary-driven NESS. Empty circles are obtained by
iteratively solving the system of equations in Appendix D. The insets show
high-frequency power-law fits to $c\sim\omega^{-2}$ and
$\delta\sim\omega^{-1}$ respectively.

We solve this system of equations iteratively. Recall that the parameters
$(c,\delta)$ are not known a priori and have to be calculated during each
iteration. In order to achieve convergence, we group the equations in a
specific way. Our algorithm can be summarized as follows:

  * •

Solve the linear system of equations for $a$ and $b$ (or equivalently $c$ and
$\delta$)

$\displaystyle\omega F_{1}(t)$
$\displaystyle=-2i^{q}J^{2}\left(I_{1}(F_{2},-F_{1})+I_{2}(F_{2},-F_{1})\right),$
(131) $\displaystyle\omega F_{2}(t)$
$\displaystyle=2i^{q}J^{2}\left((I_{1}(F_{1},F_{2})+I_{2}(F_{1},F_{2})\right).$
(132)

  * •

Solve the integral-differential equations for $F_{1,2}(t)$

$\displaystyle\partial_{t}F_{1}(t)$
$\displaystyle=i^{q}J^{2}\left(I_{1}(F_{1},F_{2})-I_{2}(F_{1},F_{2})\right),$
(133) $\displaystyle\partial_{t}F_{2}(t)$
$\displaystyle=i^{q}J^{2}\left(I_{1}(F_{2},-F_{1})-I_{2}(F_{2},-F_{1})\right).$
(134)

We start with an initial guess for $F_{1,2}(t)$ and $(c,\delta)$, and repeat
the procedure above with a weighted update at each iteration until convergence
(usually within $300$ iterations). It is worth mentioning that we can exactly
recover the DC solution $F(t)=F_{1}(t)=F_{2}(t)$ by setting $\omega=\delta=0$
and $c=1$ in the previous equations.

The iterative procedure is a lot faster to compute than the full time
evolution of a boundary-driven chain and the results are more accurate. A
sample solution for $\omega=5$ is shown in Fig. 7(d). The values of $c$ and
$\delta$ at different frequencies are plotted in Fig. 10. We see a good
agreement with the numerical results extracted from the boundary-driven NESS.
Surprisingly, we find these parameters to be independent of $q$. At high
frequencies, we find a quadratic dependence for the amplitude decay
$c\sim\omega^{-2}$. This can be easily deduced from the scaling properties of
Eq. (123); the left-hand side scales as $\omega$, while the right-hand side
scales as $1/\omega c$ after integration. For the two sides to match, we must
have $c\sim\omega^{-2}$. On the other hand, the phase scales as
$\delta\sim\omega^{-1}$. At infinite frequency, we expect the consecutive
sites to oscillate exactly out of phase with $\delta=\pi$.