\section*{Abstract}
{\bf
We study an array of strongly correlated quantum dots of complex SYK type and account for the effects of
quadratic terms added to the SYK  Hamiltonian; both local terms and inter-dot tunneling are considered
in the non-Fermi-liquid temperature range $T \gg T_{FL}$.  We take into account soft-mode fluctuations
and demonstrate their relevance for physical observables.
Electric $\sigma(\omega,p)$ and thermal $\kappa(\omega,p)$  
conductivities are calculated as functions of frequency and momentum  for arbitrary values of 
the particle-hole asymmetry parameter $\mathcal{E}$.
At low-frequencies $\omega \ll T$ we find
the Lorenz ratio $L = \kappa(0,0)/T\sigma(0,0)$ to be non-universal and temperature-dependent.
At  $\omega \gg T$ the  conductivity $\sigma(\omega,p)$ contains a pole with nearly linear dispersion
$\omega \approx sp\sqrt{\ln\frac{\omega}{T}}$ reminiscent of the "zero-sound", known for Fermi-liquids.
We demonstrate also that the developed approach makes it possible to understand the origin of heavy Fermi liquids
with anomalously large Kadowaki-Woods ratio.



}


\vspace{10pt}
\noindent\rule{\textwidth}{1pt}
\tableofcontents\thispagestyle{fancy}
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}


\section{Introduction}
\label{sec:intro}
The  Sachdev-Ye-Kitaev (S) model~\cite{Kitaev2015} describes a zero-dimensional system of $N$  interacting fermions. 
The energy levels of electrons are degenerate, and the interaction matrix elements are random  with zero mean and variance with typical energy scale $J$.
On the one hand, the SYK  model presents a fully analytical description of a strongly interacted fermionic system \cite{maldacena2016remarks,bagrets2016sachdev,kitaev2018soft,kitaev2019statistical}. On the other hand, the model has "asymptotic symmetry" (at the mean-field level), which leads to the existence of Goldstone mode with strong infra-red fluctuations.
 

From the condensed matter theory viewpoint, two  questions seem to be natural:
1) how could one extend  the SYK model to describe  fermionic systems in dimensions $d \geq 1$  ?
2) how the properties of this model are modified under reasonable perturbations? 
The simplest way to extend SYK model to higher dimensions is to consider an array of SYK quantum dots coupled by
weak tunnelling~\cite{song2017strongly}.
Using the saddle-point approximation valid at $N \to \infty$, 
the authors of Ref.~\cite{song2017strongly} predicted non-Fermi-liquid temperature
behavior of \textit{dc} conductivity $\sigma(T) \propto 1/T$ \cite{varma1989phenomenology,gegenwart2008quantum,cooper2009anomalous,bruin2013similarity,martin1990normal} in the temperature range $T \gg T_{FL}$, while at $T$ 
below the crossover temperature $T_{FL}$ the usual Fermi-liquid behavior with $\sigma(T\to 0)  \to \mathrm{const}$ was found.
A number of other papers have also considered similar issues~\cite{altland2019sachdev,gnezdilov2018low,khveshchenko2020connecting,davison2017thermoelectric}.



The answer to the second question becomes non-trivial once the analysis goes beyond the mean-field limit $N\to \infty$.
At finite values of $N$ quadratic perturbations of the strength $\Gamma \ll J$ added to the Hamiltonian of a
single SYK quantum dot with basic energy scale $J$ leave intact~\cite{lunkin2018} the non-Fermi-liquid 
ground-state for sufficiently small $\Gamma < \Gamma_c \approx J/N$.
The behavior of the mixed SYK$_4$+SYK$_2$ quantum dot at stronger perturbation magnitudes  $\Gamma_c < \Gamma \leq J/\sqrt{N}$
was studied in Ref.~\cite{lunkin2020perturbed}  where unexpected "polaron" bound-state solution was found.
Physical significance of the corresponding bound-state energy was clarified in Ref.~\cite{lunkin2022non},
where non-equilibrium generalization of the approach~\cite{lunkin2020perturbed} was developed.
We demonstrated that high-frequency parametric modulation of the SYK$_2$ part of the Hamiltonian leads to
$T_{FL} = \Gamma^2/J < T \ll \Gamma$ to the energy dissipation rate that is peaked at the frequency 
corresponding to the polaron energy $\sim \Gamma$.  As a byproduct, we demonstrated that physical 
susceptibilities at  $T_{FL} \ll T \leq \Gamma$ acquire important fluctuational contributions while single-fermion 
Green's function is described by pure SYK$_4$ saddle-point solution.




In the present paper we study non-equilibrium properties of an array of SYK$_4$+SYK$_2$ quantum dots and
extend substantially the results of Ref.~\cite{song2017strongly}. On technical level, our extension amounts to
the account for the effects of fluctuation of the soft reparametrization mode (absent in Ref.~\cite{song2017strongly} 
due to the limit $N \to \infty$ employed there). Clearly, any real quantum dot may contain finite number of electrons only,
and this number does not need to be very large  generically, thus the account for the effects arsing due to finite $N$
is essential. Moreover, our analysis demonstrate the existence of  important temperature region $T_{FL} < T < \Gamma $
where the said fluctuation effects contribute considerably to thermal conductivity even in the $N \to \infty$ limit,
while static electric conductivity is not affected by those effects and is given correctly by the results of
Ref.~\cite{song2017strongly}.
Another direction of our study where it differs from Ref.~\cite{song2017strongly} and similar papers is that we consider
frequency- and momentum-dependent electric and thermal conductivities; moreover, we consider general situation 
with an arbitrary particle-hole asymmetry. Last but not least, we consider more general model where quadractic terms
in the Hamiltonian are present both as tunneling terms (with strength $w$) and
 \textit{inside} quantum dots, with  strength $\Gamma$. It will allow us to describe the situation of "heavy Fermi liquid"
with anomalously large Kadowaki-Woods ratio known for $^3$He under pressure and some HTSC materials, 
see~\cite{Behnia2022dynamic}.


In our model SYK dots are coupled by tunnelling matrix elements $t_{ij}$ with the typical 
energy scale $w^2 = N \overline{t_{ij}^2}$.  Intra-dot quadratic terms are characterized,
 like in Ref.~\cite{lunkin2022non}, by random matrix elements $\Gamma_{ij}$ with
effective strength $\Gamma $. In the absence of intra-dot interactions, 
the system's bandwidth would be of the order of $w$, like in usual granular metal~\cite{GranularMetalReview}.
Our goal is to study the opposite situation when intra-dot interactions are  very strong, with typical energy scale
$J \gg w,\Gamma$. At the same time, we do not include long-range (inter-dot Coulomb) interactions into our model.
For this reason, any response to electric field we calculate below should be understood as
response to the local field, but not to the external one.


The role of quadratic terms of the Hamiltonian depends crucially on the temperature scale considered.
At relatively high temperatures $ J \gg T \geq \tilde\Gamma$ charge and energy transport are described by
the large-$N$ saddle-point approximation developed in Ref.~\cite{song2017strongly}.
Here $\tilde\Gamma \sim \max(\Gamma,  w)$ is  the effective strength of quadratic perturbations, and we will assume it
to be large in comparison to the critical magnitude $\Gamma_c \approx J/N$ determined in Ref.~\cite{lunkin2018}.
We will be most interested in the moderate-$T$ region  $\tilde\Gamma^2/J = T_{FL} \leq T  \leq \tilde\Gamma$ which is
analogous to the region $ \Gamma^2/J \leq  T \leq \Gamma$  studied in Ref.~\cite{lunkin2022non}.
It will be shown that in this temperature region  transport properties are strongly affected by fluctuations
of the soft reparametrization mode known to determine the behavior of SYK model~\cite{kitaev2018soft,kitaev2019statistical}.
The condition $\tilde\Gamma \gg \Gamma_c$ mentioned above will allow us to consider these fluctuations Gaussian.
While single-Fermion Green's function is well-described by the $N \to \infty$ approximation, fluctuations crucially affect 
kinetic properties at $T \leq \tilde\Gamma$.


The method we employ is based on the use of continuity equations for charge and energy currents in presence of
external sources: electric field and gradient of the Luttinger potential~\cite{luttinger1964theory}.
We show that reparametrization mode is intrinsically related to the energy current.
We  consider frequency- and
momentum-dependent conductivities, paying major attention to the asymptotic regions of  low frequencies $\omega \ll T$
and high frequnecies $\omega \gg T$. In the static limit $\omega \to 0$ we confirm the results of Ref.~\cite{song2017strongly} 
for electric conductivity $\sigma$,
while thermal conductivity $\kappa$ is found to contain additional terms that become important at $T \leq \tilde\Gamma$.
As a result, the Lorenz ratio $L = \kappa/T\sigma$ appears to be
non-universal and strongly temperature-dependent.  At low frequencies and low momenta, the dispersion of 
electric and thermal conductivities is characterized by diffusion laws. Charge diffusion constant  $D_e$ scales 
$\propto 1/T$ like conductivity, while energy diffusion constant $D_T$ demonstrates rather nontrivial $T$-dependence
with a maximum at $ T \sim \tilde\Gamma$.




In the high-frequency region $\omega \gg T$ we found a zero-sound - like
mode which provides a pole in the frequency- and momentum-dependent conductivity $\sigma(\omega,p)$ 
when temperature is low compared to the bandwidth, $T \ll w$.
Decay rate of this mode $1/\tau(\omega) \approx \omega/\ln\frac{\omega}{T}$ is somewhat smaller than its frequency
at $T \ll \omega \ll w$. Dispersion of this mode is $\omega(p) = sp \sqrt{\ln(sp/T)}$ where velocity $s \propto w$ 
is temperature-independent. Remarkably, neither the largest energy scale of the problem $J$, nor the number $N$ of
fermions in each dot  enter these results, which  form  a new "universality class" of strongly interacting Fermi-system.
High-frequency thermal conductivity
$\kappa(\omega,p)$ contains, in general, two contributions.  One of them is proportional to electric conductivity
$\sigma(\omega,p)$ and to the square of particle-hole asymmetry parameter  $\mathcal{E}$; the second contribution
is present also at $\mathcal{E}=0$, we will call it "intrinsic" and denote as $\tilde\kappa(\omega,p)$ below.
This contribution to thermal conductivity contains a weakly-dispersive resonance at the frequency $\omega_0 \sim \tilde{\Gamma}$,
similar to the result of Ref.~\cite{lunkin2022non} for a single SYK dot.




The rest of  the paper is organzied as follows: in  Sec.\ref{sec:The model} we introduce our model and  discuss the 
source terms needed  to calculate susceptibilities. In  Sec.~\ref{sec:Around SYK solution} we describe mean-field 
solution of the model.  The effective action for small fluctuations in the presence of the sources  is derived  in 
Sec.~\ref{sec:Quadratic Action for soft modes and sources}. Section~\ref{sec: Kinetic} contains our main results
for electric and thermal conductivities. Our conclusions are provided in Sec.~\ref{sec:conclusions}. Some technical details
are present in the Appendix.






\section{The model}
\label{sec:The model}
 We study an array  of SYK quantum dots with $N$ fermions each. The Hamiltonian  is as follows:
 \begin{eqnarray}
 & H =\sum\limits_{{\textbf r}} \left\{H_{{\textbf r}}+\sum\limits_{\delta{\textbf r}}  \sum\limits_{i,j}\left(t_{{\textbf r},i;{\textbf r}+\delta{\textbf r},j}  \psi^\dagger_{{\textbf r},i} \psi_{{\textbf r}+\delta{\textbf r},j}+h.c.\right)\right\} \nonumber \\
 &  H_{\textbf r}=\frac{1}{2!} \sum\limits_{i,j,k,l} J_{i,j;k,l;{\textbf r}} ~\mathcal{A}\left\{\psi_{{\textbf r},i}^\dagger \psi_{{\textbf r},j}^\dagger  \psi_{{\textbf r},k} \psi_{{\textbf r},l}\right\}+\sum\limits_{0<i,j\le N} \Gamma_{i,j,{\textbf r}}  \psi_{{\textbf r},i}^\dagger  \psi_{{\textbf r},j}-\mu \hat{N}
 \label{eq:Hamiltonian}
 \end{eqnarray}
 Here $\psi_{{\textbf r},i}$ is an annihilation operator  in the dot with coordinate ${\textbf r}$ on the $i$th site inside the dot. Here $\hat{N}$ is the operator of the total number of particles and $\mu$ is a chemical potential.
The first term in the first line represents the Hamiltonian of an individual dot. The second term describes tunneling between neighboring dots. The dots are labeled by the index ${\textbf r}$ and notation $\delta{\textbf r}$ is used to label the
 neighbors of a dot ${\textbf r}$. 
There are $N\gg 1$ positions for fermions in each dot (labeled by index $i,j,k\ldots$). The  Hamiltonian $ H_{\textbf r}$ describes fermions in each single dot ${\textbf r}$.  The first term in the expression for $ H_{\textbf r}$ is the Hamiltonian of the SYK models for the complex fermions \cite{gu2020notes} whereas the second term in $ H_{\textbf r}$ is a perturbation which lifts the degeneracy of their spectrum. The sign $\mathcal{A}\left\{\ldots\right\}$ denotes the antisymmetrized product of operators \cite{gu2020notes}. Tensors $J\ldots$ and $\Gamma\ldots$ satisfy following symmetry properties: $J_{i,j;k,l;{\textbf r}}^*=J_{l,k;j,i;{\textbf r}}^*$\,, $\Gamma_{i,j,{\textbf r}}^*=\Gamma_{j,i,{\textbf r}}$ and $J_{j,i;k,l;{\textbf r}}=J_{i,j;l,k;{\textbf r}}=-J_{i,j;k,l;{\textbf r}}$ \,.  Components of tensors $t_{\ldots}$, $J_{\ldots}$, $\Gamma_{\ldots}$ are independent random Gaussian variables  with zero mean value and the following variances:
\begin{eqnarray}
\overline{|J_{i,j;k,l;{\textbf r}}|^2} = \frac{2 J^2}{N^3}\qquad \overline{\Gamma_{i,j,{\textbf r}}^2} = \frac{ \Gamma^2 }{N} \qquad  \overline{t_{{\textbf r},i;{\textbf r}+\delta{\textbf r},j}^2} = \frac{ w^2 }{N}
\label{eq:variance}
\end{eqnarray}
We assume that the SYK term is dominant: $J\gg\Gamma,w$. 

The operators of charge density  and energy density are defined as follows:
\begin{eqnarray}
&\mathcal{Q}_{\textbf r} = \sum\limits_{i} \psi_{{\textbf r},i}^\dagger  \psi_{{\textbf r},i} -\frac{N}{2} \nonumber \\ & \mathcal{H}_{\textbf r} =  H_{\textbf r}+\frac{1}{2}\sum\limits_{\delta{\textbf r}}  \sum\limits_{i,j}\left[\left(t_{{\textbf r},i;{\textbf r}+\delta{\textbf r},j}  \psi^\dagger_{{\textbf r},i} \psi_{{\textbf r}+\delta{\textbf r},j}+h.c.\right)+{\textbf r}\mapsto {\textbf r}-\delta {\textbf r} \right]
\label{eq: charges}
\end{eqnarray}
Note that the latter one contains symmetrized terms responsible for inter-dot tunneling.
In thermal equilibrium, the average charge at a dot is equal to $\langle\mathcal{Q}_{\textbf r}\rangle =\mathcal{Q}$. 

To calculate charge and heat conductivity, we  introduce the  source terms that describe electric voltage and  Luttinger potential.
The Hamiltonian with the sources reads:
\begin{eqnarray}
H_{S}(t) = H + \sum_{{\textbf r}} \left(\mathcal{Q}_{\textbf r} U({\textbf r},t) + L({\textbf r},t) \mathcal{H}_{\textbf r}  \right)
\end{eqnarray} 
Our goal is to derive expressions for the charge and energy currents originating from the presence of these sources.

\section{Around SYK solution}
\label{sec:Around SYK solution}
\subsection{Action in terms of the Green's function $G$ and self-energy $\Sigma$}
We will use path integral representation of the problem in the Keldysh form\cite{kamenev2011field}, 
since we are interested in non-equilibrium properties. To simplify discussion of the  derivation to the soft-mode action  we will omit terms with sources.  We will restore them during the discussion of the continuity equations. The full action with sources could be found in the appendix. 

The action is defined as an integral over Keldysh contour and has the form:
\begin{eqnarray}
S=\int\limits_{\mathcal{C}} \left(\sum_{{\textbf r},i} \psi_{{\textbf r},k}^\dagger  i\partial_t \psi_{{\textbf r},k}-H_{S}(t)\right) dt
\label{eq: action}
\end{eqnarray}
where $\mu$ is the chemical potential of electrons.

First of all we need to take average over disorder. Since the Keldysh functional integral is normalized to unity (in the absence of source fields), the averaging of the functional with the action
(\ref{eq: action})  can be performed in the  straightforward  way, making use of Eqs.~(\ref{eq:variance}).


We would like to  rewrite the action in terms of the matrix field $G(t_1,t_0|{\bf r})$ defined as:
\begin{eqnarray}
G(t_1,t_0|{\bf r}) = -\frac{i}{N}\sum_{l} \psi_{{\textbf r},l}(t_1) \psi^\dagger_{{\textbf r},l} (t_0)
\label{eq: G defenition}
\end{eqnarray}
and axiliary  matrix field $\Sigma(t_1,t_0|{\bf r})$ which plays the role of Lagrangian multiplier conjugated to
$G(t_1,t_0|{\bf r})$. Then the action becomes quadratic in terms of fermions, and they can be integrated out.
The action acquires then the following form:
\begin{eqnarray}
&S=-iN\sum\limits_{{\textbf r}} \int\limits_{\mathcal{C}} dt_1 dt_0 (\Sigma(t_1,t_0|{\bf r})+\Sigma_{free}(t_1,t_0) )G(t_0,t_1|{\bf r})-iN~tr\log\{\hat{\Sigma}\}\nonumber\\
&iN  \int\limits_{\mathcal{C}} dt_1 dt_0 \sum\limits_{{\textbf r}}  \left\{\frac{J^2}{4}\left[G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r})\right]^2+\frac{\Gamma^2}{2}G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r})+w^2 \sum\limits_{\delta {\textbf r}}  G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r}+\delta {\bf r})\right\}\nonumber \\
\label{eq:action and sources}
\end{eqnarray}
We introduced  new notation $\Sigma_{free}(t_1,t_0)\equiv  \left\{ i\partial_{t_1}+\mu\right\} \delta(t_1-t_0)$. 
since  $N\gg1$  we can start from the mean-field (saddle-point) solution and then add fluctuation effects.

\subsection{Pure SYK solution}
\label{subsec: Pure SYK solution }
Saddle-point equations for $\Sigma$ and $G$ can be obtained by minimization of the action (\ref{eq:action and sources})
over these matrix variables in the absence of sources. The terms describing quadratic perturbations to the Hamiltonian
are irrelevant at this stage as long as we consider temperatures $T \gg T_{FL}$, where 
$T_{FL}\equiv\max\{\frac{w^2}{J},\frac{\Gamma^2}{J}\}$.  It will be convenient to replace time variable $t$ by
its dimensionless analog $u \equiv 2\pi T t$. Then self-energy and Green's function will change according to
\begin{eqnarray}
 G(t_1,t_0|{\bf r}) \equiv \left(\frac{2\pi T}{J}\right)^{1/2} G(u_1,u_0|{\bf r})  \qquad  \Sigma(t_1,t_0|{\bf r}) \equiv J^2 \left(\frac{2\pi T}{J}\right)^{3/2} \Sigma(u_1,u_0|{\bf r})
\end{eqnarray}
 Using these new variables we   write the SYK action in the form:
\begin{eqnarray}
&S_{SYK}= -iN~tr\log\{\hat{\Sigma}\}-iN\sum\limits_{{\textbf r}} \sum\limits_{\sigma_1,\sigma_0}\sigma_1 \sigma_0\int\limits du_1 du_0 \Sigma_{free,\sigma_1\sigma_0}(u_1,u_0) G_{\sigma_1\sigma_0}(u_0,u_1|{\bf r})\nonumber \\&iN  \sum\limits_{\sigma_1,\sigma_0}\sigma_1\sigma_0\int du_1 du_0 \sum\limits_{{\textbf r}}  \left\{\frac{1}{4}\left[G_{\sigma_1\sigma_0}(u_1,u_0|{\bf r})G_{\sigma_0\sigma_1}(u_0,u_1|{\bf r})\right]^2-\Sigma_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) G_{\sigma_0\sigma_1}(u_0,u_1|{\bf r}) \right\} \nonumber \\ 
\label{eq:approximated action}
\end{eqnarray}
Here we split a Keldysh contour $\mathcal{C}$ into positive and negative branches~\cite{kamenev2011field}, 
so that summation  variables $\sigma_{0,1}$ take values $\pm 1$.
Minimization  of the action (\ref{eq:approximated action}) provides the saddle-point functional equations in the form:
\begin{eqnarray}
\Sigma_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) = \left[G_{\sigma_1\sigma_0}(u_1,u_0|{\bf r})G_{\sigma_0\sigma_1}(u_0,u_1|{\bf r})\right]^{\frac{q}{2}-1} G_{\sigma_0\sigma_1}(u_0,u_1|{\bf r}) \nonumber \\
\sum_{\sigma}\sigma \sigma_1\int \Sigma_{\sigma_1\sigma}(u_1,u|{\bf r}) G_{\sigma\sigma_0}(u,u_0|{\bf r}) du = -\delta_{\sigma_1,\sigma_0} \delta(u_1-u_0)
\label{eq:mean-field equation}
\end{eqnarray}
We omit the term  $\Sigma_{free}$ from the saddle-point equations (\ref{eq:mean-field equation}) 
 since it is parametrically  small at time scales $t \gg 1/J$ where we will consider these equations.
We introduce parameter $q$ which is equal to  $4$ in our case; this  form can be important in the context of generalized SYK model (with arbitrary number of fermions, \cite{kitaev2018soft}) we also used this parameter for the dimensional regularization  (see comment in \ref{sec: appendix_A})
Translationally invariant solution of  these equations is (see also \cite{sachdev2015bekenstein}):
\begin{eqnarray}
\Large G^{(s)}_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) = b^\Delta g_{\sigma_1\sigma_0}(u_1-u_0) \qquad b = \frac{(1-2\Delta)\sin(2\pi \Delta)}{2\pi (\cos(2\pi \Delta)+\cosh(2\pi\mathcal{E}))}
 \label{eq:Green function of fermions}
\end{eqnarray}
where we have used function $g$ defined as follows:
\begin{eqnarray}
\Large g_{\sigma_1\sigma_0}(u)\equiv i e^{-i\mathcal{E}u}\left(4\sinh^2\left(\frac{u}{2}\right)\right)^{-\Delta} 
\left[\theta(u)\left(\begin{smallmatrix}
-e^{-i\pi \left(\Delta+i\mathcal{E}\right)} & e^{i\pi \left(\Delta+i\mathcal{E}\right)} \\ -e^{-i\pi \left(\Delta+i\mathcal{E}\right)} & e^{i\pi \left(\Delta+i\mathcal{E}\right)}
\end{smallmatrix}\right)+\theta(-u)\left(\begin{smallmatrix}
e^{-i\pi \left(\Delta-i\mathcal{E}\right)} & e^{-i\pi \left(\Delta-i\mathcal{E}\right)} \\ -e^{i\pi \left(\Delta-i\mathcal{E}\right)} & -e^{i\pi \left(\Delta-i\mathcal{E}\right)} 
\end{smallmatrix}\right)\right]_{\sigma_1\sigma_0} \nonumber \\
\end{eqnarray}
Here $\Delta= \frac{1}{q}$ and $\mathcal{E}\in (-\infty,\infty)$ is the parameter of particle-hole asymmetry which can be 
related  (see \cite{gu2020notes,sachdev2015bekenstein}) to the average charge of the dot $\mathcal{Q}$ as follows:
\begin{eqnarray}
\mathcal{Q}= -\frac{\theta}{\pi}-\left(\frac{1}{2}-\Delta\right) \frac{\sin(2\theta)}{\sin(2\pi \Delta)}\qquad  e^{2\pi \mathcal{E} }\equiv \frac{\sin(\pi\Delta+\theta)}{\sin(\pi\Delta-\theta)}
\label{eq:charge and eps}
\end{eqnarray}
Actually the asymetry parameter $\mathcal{E}$ cannot be too large in its absolute value. The  stable solutions of the type provided by  Eq.(\ref{eq:Green function of fermions})  can be obtained for $|\mathcal{E}| \leq 0.14$ only, as it was found numerically
 in Refs.~\cite{azeyanagi2018phase,tikhanovskaya2021excitation}. Although these papers studied models  similar but not identical
to ours, they worked  with the same saddle-point equations (\ref{eq:mean-field equation}), so we expect that their results
for maximum value of $|\mathcal{E}|$
apply to our case as well.

\subsection{Selection of the soft-mode manifold}
\label{subsec: Selection of the soft-mode manifold}
The saddle-point solution (\ref{eq:Green function of fermions})  
breaks down large symmetry of the saddle-point equations (\ref{eq:mean-field equation}). Namely, 
symmetry transformations allowed for Eqs.(\ref{eq:mean-field equation}) are as follows:
 for any solution $G_{\sigma_1,\sigma_0}(u_1,u_0|{\bf r})$  we can produce a new solution $\tilde{G}$ defined as
\begin{eqnarray}
\tilde{G}_{\sigma_1,\sigma_0}(u_1,u_0|{\bf r})= \left[F_{\sigma_1}^{'}(u_1,{\bf r})F_{\sigma_0}^{'}(u_0,{\bf r})\right]^\Delta e^{i\left(\varphi_{\sigma_1}(u_1,{\bf r})-\varphi_{\sigma_0}(u_0,{\bf r})\right)} G_{\sigma_1,\sigma_0}(F_{\sigma_1}(u_1,{\bf r}),F_{\sigma_0}(u_0,{\bf r})|{\bf r}) \nonumber \\
\label{eq:symmetry transform}
\end{eqnarray}
Here $F(u)$ is a  monotonous function, and $\varphi(u)$ is an arbitrary function. 
On the other hand, the solution (\ref{eq:Green function of fermions}) allows for much smaller symmetry group $SL(2,R)$.
As a  result,  approximate action (\ref{eq:approximated action}) has a zero mode. Once previously neglected
term with $\Sigma_{free}$ is put back into the action, this zero mode transforms into soft Goldstone mode,
like it happens in the $\sigma$-model~\cite{efetov1999supersymmetry}.
The corresponding  manifold is parametrized by the functions $F(u)$ and $\varphi(u)$  by means of  transformations
 (\ref{eq:symmetry transform}) applied to the mean-field solution (\ref{eq:Green function of fermions}). 
As a result, we obtain the following action for the soft modes:
\begin{eqnarray}
&S_{soft}= \int d{\textbf r} S_{SYK,{\textbf r}}^{(F,\phi)}+(S-S_{SYK})|_{G=G^{(F,\phi)}}\qquad G^{F,\phi}_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) = b^\Delta g^{(F,\phi)}_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) \nonumber \\
&g^{(F,\phi)}_{\sigma_1\sigma_0}(u_1,u_0|{\bf r}) =  \left[F_{\sigma_1}^{'}(u_1,{\bf r})F_{\sigma_0}^{'}(u_0,{\bf r})\right]^\Delta e^{i\left(\varphi_{\sigma_1}(u_1,{\bf r})-\varphi_{\sigma_0}(u_0,{\bf r})\right)}\times \nonumber \\& \times g_{\sigma_1,\sigma_0}(F_{\sigma_1}(u_1,{\bf r}),F_{\sigma_0}(u_0,{\bf r})) \nonumber \\
 \label{eq:soft mode action}
\end{eqnarray}
Here $ S_{SYK,{\textbf r}}^{(F,\phi)}$ is the original action for soft modes in the SYK model which has the form:
\begin{eqnarray}
 S_{SYK,{\textbf r}}^{(F,\phi)} =  \sum\limits_\sigma \sigma \int du \left\{-C_E Sch\left(e^{F_{\sigma}(u,{\bf r})},u\right)+C_Q\left(\varphi^{\prime}_{\sigma}(u,{\bf r})+\mathcal{E}F^\prime_{\sigma}(u,{\bf r})-\mathcal{E}\right)^2\right\}
  \label{eq:SYK soft mode action} 
\end{eqnarray}
Here 
\begin{eqnarray}
C_E = N \alpha_s \frac{2\pi T}{J V_0}\qquad C_Q = N K  \frac{2\pi T}{J  V_0 }
\label{CQ}
\end{eqnarray}
where $\alpha_s\approx 0.05$ and $K\approx  1.04$ see \cite{kitaev2018soft,gu2020notes} and $ V_0$ is the system
 volume per  single dot . The second term in $S_{soft}$
comes from quadratic perturbations in the Hamiltonian and contains terms $\propto \Gamma^2$ and  $w^2$.



Soft fluctuations play a crucial role in the theory of the SYK model: they change asymptotic behaviour of the Green function,
as shown in Ref.~\cite{bagrets2016sachdev}).  On the other hand, quadratic perturbations (terms with $\Gamma$ and $w$) partially suppress these fluctuations, as demonstrated in Ref.~\cite{lunkin2020perturbed}. 
These fluctuations are also responsible for the kinetic properties of the system \cite{song2017strongly,lunkin2022non}. Bellow we will consider a quadratic action for these fluctuations to find saddle-point solution in the presence of sources. We will discuss physical properties of these fluctuations in   section \ref{sec: Kinetic}. 



 
\section{Quadratic action for soft modes in presence of  source fields}
\label{sec:Quadratic Action for soft modes and sources}

 We  will consider soft-mode fluctuations  around the mean-field solution (\ref{eq:Green function of fermions}) at Gaussian level in this section; smallness of these fluctuations is controlled by
the following inequalities:
\begin{eqnarray}
F_{\sigma}(u,{\bf r}) = u+f_{\sigma}(u,{\bf r})\qquad  f^\prime_{\sigma}(u,{\bf r})\ll 1\qquad  \varphi^\prime_{\sigma}(u,{\bf r})\ll 1.
\end{eqnarray}
We  will expand  action  (\ref{eq:soft mode action}) up to the second order in small fluctuations 
$f$ and $\varphi$\, (similar procedure was used in Ref.~\cite{lunkin2022non}).

The result can be conveniently written in the Fourier representation:
\begin{eqnarray}
&S_{final}=\frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left(\hat{\varphi}_{{\bf p},\Omega}+\mathcal{E}\hat{f}_{{\bf p},\Omega}\right)^\dagger \left\{\left[\hat{\mathcal{G}}^{(\varphi)}_0(\Omega)\right]^{-1}- \hat{\Sigma}^{(\varphi)}_t(\Omega,{\bf p})\right\}\left(\hat{\varphi}_{{\bf p},\Omega}+\mathcal{E}\hat{f}_{{\bf p},\Omega}\right) \nonumber \\
& \frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \hat{f}_{{\bf p},\Omega}^\dagger\left\{\left[\hat{\mathcal{G}}^{(f)}_0(\Omega)\right]^{-1}-\hat{\Sigma}^{(f)}(\Omega)- \hat{\Sigma}^{(f)}_{w}(\Omega,{\bf p})\right\}\hat{f}_{{\bf p},\Omega}
\label{eq:soft mode action and sources}
\end{eqnarray}
Here $d$ is the space dimension of our problem.
We introduced vectors that describe fields and sources on positive ($\varphi^{+}_{{\bf p},\Omega}$) and negative  ($\varphi^{-}_{{\bf p},\Omega}$) parts of the Keldysh contour.  We also performed Keldysh rotation defined as follows:
\begin{eqnarray}
\varphi^{+}_{{\bf p},\Omega} = \varphi^{(cl)}_{{\bf p},\Omega}+\frac{\varphi^{(q)}_{{\bf p},\Omega}}{2}\qquad \varphi^{-}_{{\bf p},\Omega} = \varphi^{(cl)}_{{\bf p},\Omega}-\frac{\varphi^{(q)}_{{\bf p},\Omega}}{2}.
\end{eqnarray}
The same notations are used for other fields and sources. $2\times 2$ Keldysh space represented in terms of $(cl,q)$
variables is convenient for representations of matrices 
$\hat{\Sigma}$ and $\left[\hat{\mathcal{G}}^{(\varphi)}_0(\Omega)\right]^{-1} $.
 Since we consider fluctuations in
Gaussian approximation only, it will be sufficient to study the behavior of retarded (R) part of Green function and self-energy.

The corresponding explicit expressions are provided below:
\begin{eqnarray}
\label{GSR}
&\left[\hat{\mathcal{G}}^{(\varphi)}_0(\Omega)\right]^{-1}_{(R)}=\left[\mathcal{G}^{(\varphi)}_{0,(R)}(\Omega)\right]^{-1}=C_{Q} \Omega^2 \qquad \qquad \left[\mathcal{G}^{(f)}_{0,(R)}(\Omega)\right]^{-1} = C_{E} \Omega^2 (\Omega^2+1) \nonumber \\
 &\Sigma^{(f)}_{(R)}(\Omega) =  \frac{C_\Gamma}{2}  \Omega^2 \psi(\Omega)\qquad \qquad \Sigma^{(\varphi)}_{w,(R)}(\Omega,{\bf p}) = 8 C_w\sum\limits_{\delta{\bf r}} \left(\frac{{\bf p}\delta{\bf r}}{2} \right)^2 \left(\psi(\Omega)+2\right)\nonumber \\
&\Sigma^{(f)}_{w,(R)}(\Omega,{\bf p})=\frac{C_w}{2}\sum\limits_{\delta{\bf r}} \left\{\Omega^2 \psi(\Omega) +\frac{1}{2}\left(\frac{{\bf p}\delta{\bf r}}{2} \right)^2\left(2(1+\Omega^2)+(1+2\Omega^2)\psi(\Omega)\right) \right\} \nonumber \\
\end{eqnarray} 
Here we introduced the functions $\psi(\Omega)$, and constants $C_{\Gamma}$, $C_w$, which are defined as follows:
\begin{eqnarray}
\label{26}
\psi(\Omega)=\Psi\left(\frac{1}{2}-i\Omega\right)-\Psi\left(-\frac{1}{2}\right)\, ; \qquad \Psi(z)=\ln^\prime\left(\Gamma(z)\right)  \nonumber \\
C_{\Gamma} = \frac{N}{4\pi V_0}\frac{\Gamma^{2}}{J T} \sqrt{b} \qquad C_{w} = \frac{N}{2\pi V_0}\frac{w^{2}}{J T} \sqrt{b}.
\end{eqnarray}
The lattice will be considered as a generalized cubic one,  so we introduce the following notations:
\begin{eqnarray}
\tilde{C}_w \equiv C_w \sum\limits_{\delta{\bf r}}1=2d C_w\, ; \qquad \tilde{\sigma} \equiv C_w \frac1{p^2}\sum\limits_{\delta{\bf r}}\left(\frac{{\bf p}\delta{\bf r}}{2} \right)^2=C_w \frac{a^2}{2}  
\label{CtCG}
\end{eqnarray}
Here $d$ is a dimension of the lattice  and $a$ is a lattice constant; the above formulas are valid at $ap\ll 1$.


The expressions (\ref{eq:soft mode action and sources}) and (\ref{GSR}) for quadratic part of the soft-mode action in presence of sources is the main result of this Section.
These formulae  will be used below to determine kinetic properties of the system in both limits of  low ($\omega \ll T$) and high
($\omega \gg T$)  frequencies.


\section{Electric and thermal transport at arbitrary frequencies}
\label{sec: Kinetic}
\subsection{Noether's theorem}
\label{subsec: Noether's theorem} 

As we are interested in the kinetics of the system we would like to obtain expressions for currents of energy and charge in our system. We need to consider the following change of the initial fermionic fields to find these currents, according to Noether's theorem: 
\begin{eqnarray}
\psi_{{\textbf r},j}(t) \mapsto e^{i \delta\varphi(t,{\textbf{r}})} \psi_{{\textbf r},j}(t+\delta f(t,{\textbf{r}}))
\end{eqnarray}
On the other hand, we use fields $f$ and $\varphi$ instead of fermionic ones. To connect them let us write a transformation law for the field $G$ using its definition (\ref{eq: G defenition}):
\begin{eqnarray}
G_{\sigma_1,\sigma_2}(t_1,t_2|{\textbf{r}}) \mapsto  e^{i (\delta\varphi_{\sigma_1}(t_1,{\textbf{r}})-\delta\varphi_{\sigma_2}(t_2,{\textbf{r}}))}G_{\sigma_1,\sigma_2}(t_1+\delta f_{\sigma_1}(t_1,{\textbf{r}}),t_2+\delta f(t_2,{\textbf{r}})|{\textbf{r}})
\end{eqnarray}
This map is not equal to the symmetry map described in eq. (\ref{eq:symmetry transform}). As we want to limit ourselves by soft mode fluctuations only we need to project the result of the symmetry application to  the manifold of soft modes. Finally, we will obtain the following transformation law for fields $f$ and $\varphi$:
\begin{eqnarray}
\varphi^{\sigma}(u,{\textbf{r}}) \mapsto \varphi^{\sigma}(u,{\textbf{r}})+\delta\varphi_{\sigma}(u,{\textbf{r}}) \qquad \qquad f^{\sigma}(u,{\textbf{r}}) \mapsto f^{\sigma}(u,{\textbf{r}})+2\pi T \delta f_{\sigma}(u,{\textbf{r}})
\end{eqnarray}
This transformation rule reveals the physical origin of fields $\varphi$ and $f$. These fields are conjugated to the charge density fluctuations and energy density fluctuation  respectively. In other words they are related as:
\begin{eqnarray}
\delta\mathcal{Q}(t,{\textbf r }) = - \frac{\delta S}{\delta \partial_u \varphi^{(q)}(u,{\textbf{r}})} \qquad \qquad \delta E(t,{\textbf r }) = - \frac{1}{2\pi T} \frac{\delta S}{\delta \partial_u f^{(q)}(u,{\textbf{r}})}
\end{eqnarray} 
We can understand the meaning of coefficients in the action using this expression.  Let us slightly change the temperature and chemical potential of the system. As a result, the saddle-point solution will change. Green function of the new system could be obtained from the Green function of the original system using symmetry transformation (\ref{eq:symmetry transform}). For this transformation: $f^{(cl)}(u,{\textbf{r}})= \frac{\delta T}{T} u$ and $ \varphi^{(q)}(u,{\textbf{r}})=\frac{\delta \mu}{2\pi T} u$. As a result, the charge density and energy density will change as:
\begin{eqnarray}
\delta \mathcal{Q}= \frac{C_Q}{2\pi T}\left(\delta \mu +2\pi \mathcal{E} \delta T \right)\qquad \delta E=2\pi T \left(\left[C_E+C_\Gamma+\tilde{C}_w\right] \frac{\delta T}{T}+\mathcal{E}\delta Q\right) 
\label{eq:thermodanmyc relations}
\end{eqnarray}
Using these equations, we can  note that heat capacity of the system is:
\begin{eqnarray}
\left(\frac{\delta E}{\delta T}\right)_{Q}=  2\pi(C_E+C_\Gamma+\tilde{C}_w)
\label{eq: heat capacity}
\end{eqnarray} In the absence of quadratic terms, this result agrees with the heat capacity of the SYK model $C_{SYK}=2\pi C_E$ \cite{kitaev2018soft}.  Quadratic terms become important for heat capacity at low temperatures  when
\begin{eqnarray}
C_E \leq C_\Gamma+\tilde{C}_w \quad \Leftrightarrow \quad T \leq \tilde{\Gamma} =  \sqrt{\Gamma^2+4 d w^2}
\end{eqnarray}
As we see, quadratic terms play the crucial role here, even at $T \gg T_{FL}$.
 


 We also can find compressibility of the  system (see also \cite{gu2020notes}): 
\begin{eqnarray}
\left(\frac{\delta Q}{\delta \mu}\right)_T =  \frac{C_Q}{2\pi T}
\label{eq: compressibility}
\end{eqnarray}
Finally, we can verify Maxwell's relation (see \cite{davison2017thermoelectric,sachdev2015bekenstein,gu2020notes}):
\begin{eqnarray}
\left(\frac{\delta S}{\delta Q}\right)_{T}=-\left(\frac{\delta \mu}{\delta T}\right)_{\mathcal{Q}}= 2\pi \mathcal{E}
\end{eqnarray}
Here we have connected entropy and energy as $\delta E= T \delta S$.

In this subsection, we have considered the statistical mechanics of the system. These properties were found as a response of our system to static perturbation. In the following subsection we will find kinetic properties of the system as a response on dynamical sources.




\subsection{Continuity equations}
\label{subsec: Continuity equations} 
In the limit $N\gg1$ we can use the saddle point equation. In the absence of the quantum component of sources the quantum component of fields at saddle-point configuration will be zero. So we have only two saddle-point equation:
\begin{eqnarray}
\frac{\delta S}{\delta \varphi^{(q)}(u,{\textbf{r}})}=0 \qquad \qquad  \frac{\delta S}{\delta f^{(q)}(u,{\textbf{r}})}=0
\label{eq: soft saddle point equations}
\end{eqnarray} 

As fields $\varphi$ and $f$ are conjugated to charge and energy densities, we conclude that the above equations are continuity equations for charge and energy densities respectively.

We can also  obtain expressions for the densities of charge current  and energy  current  as:
\begin{eqnarray}
\label{currents-def}
j^{(Q),\alpha}_{{\bf p},\omega} = -\frac{\delta S_{final}}{\delta (i p^\alpha f_{{\bf p},\Omega}^{(q)})^{\dagger}} \,
\qquad \mathrm{and} \qquad
j^{(E),\alpha}_{{\bf p},\omega} = -\frac{\delta S_{final}}{\delta(ip^\alpha \varphi_{{\bf p},\Omega}^{(q)})^{\dagger}}
\end{eqnarray}

There are non-zero currents in the system at the presence of sources ($L$ and $U$). We can solve saddle-point equations (\ref{eq: soft saddle point equations}) to find fields $\varphi$ and $f$ in the presence of sources. Using this result we can connect currents and sources as: 

 
\begin{eqnarray}
\Large \left(\begin{smallmatrix}
j^{(Q),\alpha}_{{\bf p},\omega} \\  j^{(E),\alpha}_{{\bf p},\omega}
\end{smallmatrix}\right)= \left(\begin{smallmatrix}
\sigma({\bf p},\omega)e^2 && 2\pi\mathcal{E}e\sigma({\bf p},\omega) \\ 2\pi\mathcal{E}eT\sigma({\bf p},\omega) &&  \tilde\kappa({\bf p},\omega)+(2\pi\mathcal{E})^2 T\sigma({\bf p},\omega) 
\end{smallmatrix}\right) \left(\begin{smallmatrix}
-ip^{\alpha} U_{{\bf p},\omega} \\ -ip^{\alpha} T L^{(cl)}_{{\bf p},\omega}
\end{smallmatrix}\right)
\normalsize
\label{eq:currents}
\end{eqnarray}

Here  $e$ is the  charge of electron. We  have  also restored the dimensionalities of currents and sources. Left upper block of the matrix (\ref{eq:currents})  shows that electric conductivity at arbitrary frequency and momentum is given by  $e^2\sigma(\omega,p)$, where

\begin{eqnarray}
 \sigma({\bf p},\omega=2\pi T \Omega)= \frac{ 8 \tilde{\sigma}  \tilde{\psi}(\Omega)C_Q}{-iC_Q\Omega+8\tilde{\sigma}p^2\frac{\tilde{\psi}(\Omega)}{-i\Omega}} \qquad \text{where} \quad \tilde{\psi}(\Omega) = \psi(\Omega)+2 \approx -i\Omega \quad \text{at} \quad \Omega \to 0
\label{eq:conductivity}
\end{eqnarray}
Recall that notations $C_Q$ and $\tilde\sigma$ are  defined in Eqs.(\ref{CQ}) and (\ref{CtCG}) correspondingly.

The expression for heat conductivity (right-bottom matrix element in Eq. (\ref{eq:currents})) contains two terms. The first one describes thermal  current without transport of charge (we call it intrinsic thermal conductivity):
\begin{eqnarray}
\Large
\tilde\kappa({\bf p},\omega=2\pi T \Omega) =  \tilde{\sigma} T \pi^2 \times  \frac{ \left[C_{E} (\Omega^2+1) -  \tilde{C} \psi(\Omega)\right]\left((\Omega^2+1) \tilde{\psi}(\Omega)+5\Omega^2\psi(\Omega)\right)+ 2\Omega^2  \psi^2(\Omega)C_\Gamma }{- iC_{E} \Omega (\Omega^2+1)+ \frac{\tilde{C}}{2}  i\Omega \psi(\Omega)+ \frac{\tilde{\sigma}}{4}p^2  \left((\Omega^2+1) \frac{\tilde{\psi}(\Omega)}{-i\Omega}+i\Omega\psi(\Omega)\right)  } \nonumber \\
\normalsize
\label{eq:heat conductivity}
\end{eqnarray}
Here we introduced  a new constant $\tilde C$, see Eqs.(\ref{26}):
\begin{equation}
\tilde{C}\equiv C_\Gamma+\tilde{C}_w = 
\frac{N \sqrt{b}}{4\pi V_0 JT} \tilde{\Gamma}^2
\label{tildeC}
\end{equation}
This constant describes contribution of perturbations to the heat capacity (see also \ref{eq:thermodanmyc relations} ) 

The second term in the right-bottom block of (\ref{eq:currents}) is due to charge transport; 
the total thermal conductivity is given then by
\begin{equation}
\kappa({\bf p},\omega) = \tilde\kappa({\bf p},\omega) + (2\pi\mathcal{E})^2\,T\,\sigma({\bf p},\omega)
\label{kappa}
\end{equation}
Finally, off-diagonal elements of the conductivity matrix (\ref{eq:currents}) define a Seebeck (${\mathscr S}$) coefficient,
 which is equal to
\begin{eqnarray}
\mathscr{S}=\frac{2\pi  \mathcal{E} }{e}=\frac{\partial S}{\partial \mathcal{Q}}
\end{eqnarray}
where $S$ is the entropy of the system. This result was obtained in Refs.~\cite{sachdev2015bekenstein,gu2020notes}.

Below we consider simplified expressions for electric and thermal conductivities at low frequencies $\omega \ll T$ and
and high frequencies $\omega \gg T$.

\subsection{Hydrodynamic limit}
\label{subsec: Hydrodynamic limit}
In the hydrodynamic limit $\omega\ll T$, the system is in a local equilibrium. 
In this limit, electric conductivity and intrinsic heat conductivity acquire the following form:
\begin{eqnarray}
\Large \sigma({\bf p},\omega)=  (2\pi)^2 \tilde{\sigma} \frac{ -i \omega  \frac{C_Q}{2\pi T}}{ -i \omega  \frac{C_Q}{2\pi T}+(2\pi)^2\tilde{\sigma}p^2}
\end{eqnarray}
\begin{eqnarray}
\Large \tilde\kappa({\bf p},\Omega)=  \frac{\tilde{\sigma} T\pi^4}{2} \frac{-i 2\pi \left[C_{E}  +2 \tilde{C} \right]\omega }{- i 2\pi\left[C_{E}  + \tilde{C} \right]\omega +  \frac{\tilde{\sigma}T\pi^4}{2}p^2   } \normalsize
\label{eq:conductivity low frequency}
\end{eqnarray}
Uniform static electric conductivity is equal to
\begin{equation}
\sigma_0(T)=\sigma({\bf 0},\Omega \rightarrow 0;T) = \pi \sqrt{b} \frac{N}{ V_0}\frac{w^{2}}{J T} a^2.
\label{sigma0}
\end{equation}
Dispersions of  both conductivities as functions of frequency and momentum are characterized by diffusive dependencies.
The diffusion coefficients are as follows (see also  (\ref{eq: heat capacity}, \ref{eq: compressibility})):
\begin{eqnarray}
D_e=\sigma_0(T) \left[\frac{\delta Q}{\delta \mu}\right]^{-1} \propto \frac1{T} \qquad\text{charge diffusion} \\
D_T= \frac{\pi^2}{8} \sigma_0(T)\, T\, \left[\frac{\delta E}{\delta T}\right]^{-1} 
\propto \begin{cases}
& T \quad \textrm{at}\,\, T \ll \tilde{\Gamma}\\
& T^{-1} \quad \textrm{at}\,\, T \gg \tilde{\Gamma}`
\end{cases} ,
\qquad\text{energy diffusion}
\end{eqnarray}


Lorenz ratio  $L$ in the zero-$\omega$ and zero-$p$ limit is given (see Eqs.(\ref{eq:conductivity low frequency}), 
 (\ref{kappa}),(\ref{CQ}) and (\ref{CtCG}) ) by
\begin{eqnarray}
\label{L1}
L=\frac{\kappa({\bf 0},0)}{T\sigma({\bf 0},0)} +(2\pi\mathcal{E})^2 =  \frac{\pi^2}{8} \frac{C_{E}  +2  \tilde{C}  }{C_{E}+\tilde{C} }+(2\pi\mathcal{E})^2 = \frac{\pi^2}{8} \frac{  \alpha_s  +\sqrt{b}\frac{\tilde{\Gamma}^2}{   (2\pi T)^2}   }{  \alpha_s + \frac{ \sqrt{b}}{2}\frac{\tilde{\Gamma}^{2}}{ (2\pi T)^2} }+(2\pi\mathcal{E})^2
\end{eqnarray} 
At relatively high temperatures $T \geq \tilde{\Gamma}$ the heat capacity is dominated by pure SYK term, $C_E \gg C_\Gamma$,
and Lorenz ratio follows  the result of Ref.~\cite{davison2017thermoelectric},
$L = \pi^2/8+(2\pi \mathcal{E})^2$ . However, at lower temperatures $T \ll \tilde{\Gamma}$ the heat capacity is determined by the term coming from quadratic
perturbation in the Hamiltonian, $\tilde{C}\gg C_E$, and we find for Lorenz ratio another result:
\begin{eqnarray}
L = \frac{\pi^2}{4} +(2\pi\mathcal{E})^2 
\label{L2}
\end{eqnarray}
We emphasize that considerable variation of the Lorenz ratio with temperature occurs in domain $T \sim \tilde{\Gamma} \gg T_{FL}$
where modifications of the saddle-point solution are small.  The whole effect upon $L(T)$ which we found,
 comes due to soft-mode fluctuations. At lowest temperatures $T \leq T_{FL}$ the Lorenz  ratio returns to its Fermi-liquid
value $L_{FL} = \frac{\pi^2}{3}$ irrespectively of $\mathcal{E}$.






 
\subsection{High-frequency, low-$T$ limit}
\label{subsec: High-frequency, low-$T$ limit}
Relaxation rate in our system is bounded  as $1/\tau_{\mathrm{inel}} \leq T $ thus the high-frequency range $\omega \gg T$
corresponds to its non-equilibrium behaviour. Surprisingly, we find  electric current
response containing a pole in the complex plane of frequency with nearly linear dispersion and relatively small decay rate.
This resonance mode looks similar to the zero sound in Fermi-liquids~\cite{lifshitz2013statistical}
 in spite of completely different nature of our system.
Electric conductivity  in this regime can be written as follows (see Eq.(\ref{eq:conductivity}) ):
\begin{eqnarray}
\Large
\sigma({\bf p},\omega)=  8 \tilde{\sigma}C_Q \frac{  i\frac{\omega}{2\pi T} \psi(\frac{\omega}{2\pi T})}{C_Q\left[\frac{\omega}{2\pi T}\right]^2-8\tilde{\sigma}p^2 \psi(\frac{\omega}{2\pi T})} \qquad \psi(\Omega)\approx \ln(|\Omega|)-i\frac{\pi}{2}sign(\Omega) 
\label{sigmaHigh}
\end{eqnarray}
At $\Omega \gg 1$  the real part of the function $\psi(\Omega)$ is much larger than its imaginary part which makes
possible propagation of waves. The  dispersion of waves is determined by the position of the pole in in Eq.(\ref{sigmaHigh}):
\begin{eqnarray}
\omega^2(p)=s^2 p^2 \ln\left(\frac{ s p  }{T}\right)\qquad \, \quad  s=  2\sqrt{D_e T}
\label{0sound} 
\end{eqnarray}

The  characteristic velocity
$s$ is temperature-independent as $D_e\sim \frac{1}{T}$ and we can express  it as: 
\begin{equation}
 s = \frac{2b^{1/4}}{K^{1/2}} a w
\label{s}
\end{equation}
For the condition $sp \gg T$ to be satisfied, the temperature should be low enough, $T \ll w$.
Note that $s$ does not depend on the largest energy scale $J$, although the mere existence of this mode is related to
 strong electron-electron interaction. Similar situation was found in our previous study~\cite{lunkin2022non} of
the energy absorption in a single SYK$_4$+SYK$_2$ quantum dot.  Life-time $\tau(p)$ of the "zero-sound" modes (\ref{0sound})
 is relatively long at $\omega \gg T$, namely: $\omega(p)\tau(p) \approx \ln\frac{\omega}{T}$.


Thermal conductivity in the high-frequency regime contains two terms, one of them is proportional to $\sigma({\bf p},\omega)$
given by Eq.(\ref{sigmaHigh}), see general relation (\ref{kappa}).  The  second term, the intrinsic thermal conductivity,
is given by
\begin{eqnarray}
\label{kappa-high}
\tilde\kappa({\bf p},\omega) =  3\tilde{\sigma}  \pi \psi\left(\frac{\omega}{2\pi T}\right)   \frac{ i \omega\left[C_{E} \left[\frac{\omega}{2\pi T}\right]^2 -  (\tilde{C}-\frac{1}{3}C_\Gamma ) \psi(\frac{\omega}{2\pi T})\right]}{C_{E}  \left[\frac{\omega}{2\pi T}\right]^2- \frac{\tilde{C}}{2}  \psi(\frac{\omega}{2\pi T})- \frac{\tilde{\sigma}}{2}p^2 \psi(\frac{\omega}{2\pi T})   }
\end{eqnarray}

As a  function of frequency, $\tilde\kappa({\bf p},\omega)$ contains a pole at the  frequency $\omega_{th}(p)$ which is close
to a single resonance frequency $\omega_0$ defined below, since we always have the condition $\tilde\sigma p^2 \ll C_E$ fulfilled.
The equation for $\omega_0$ reads
\begin{eqnarray}
  \left[\frac{\omega_0}{2\pi T}\right]^2= \frac{\tilde{C}}{2 C_E} \ln\left(\frac{\omega_0}{2\pi T}\right) \quad \rightarrow \quad 
	\omega_0 = \frac{b^{1/4}}{2\sqrt{\alpha_S}} \left[\tilde{\Gamma}^2\ln\frac{\omega_0}{2\pi T}\right]^{1/2}
\label{omega01}
\end{eqnarray}
To characterize decay rate and dispersion of the intrinsic thermal mode $\omega_{th}(p)$ we set \\
$\omega_{th}(p) \equiv sign(\omega)\omega_0+ \delta\omega_{th}(p)$ and find:

\begin{eqnarray}
\delta \omega_{th}(p) =  - \frac{i\pi}{4}\frac{\omega_0}{\ln\left(\frac{\omega_0}{2\pi T}\right)} + \frac{sign(\omega)\sigma_0 \omega_0}{8\pi^2 C}p^2 
\end{eqnarray}
where $C=C_E+\tilde{C} $.  
The mode $\omega_{th}(p)$ is completely analogous to the  frequency of resonant absorption found in Ref.~\cite{lunkin2022non}.


\section{Conclusions}
\label{sec:conclusions}
We have studied a model of strongly correlated electron liquid described by SYK$_4$ + SYK$_2$ quantum dots coupled
by single-particle tunneling terms. Effective strength of quadratic perturbations 
$\tilde\Gamma = \sqrt{\Gamma^2 + 4d w^2}$  is small compared to interaction strength $J$, but it is still much larger
than the critical value $\Gamma_c \approx J/N$ found in Refs.\cite{lunkin2018,lunkin2020perturbed}.
Therefore non-Fermi-liquid behaviour of the system is limited from below by the temperature scale 
$T_{FL} \approx \tilde\Gamma^2/J$, and we consider the range of temperatures $T \geq T_{FL}$ only.
Our results do not depend on the number $N$ of electrons in each quantum dot, as long as we work at
 $\tilde\Gamma \gg \Gamma_c$. However, at smaller $\tilde\Gamma$ the parameter $N$ becomes important; this region
of parameters is still to be studied.

 The generic case of an electron - hole asymmetric model, $\mathcal{E} \neq 0$, is considered.
General expressions for frequency-dependent  and momentum-dependent electric and thermal conductivities 
$\sigma(\mathbf{p},\omega) $ and $\kappa(\mathbf{p},\omega)$ are derived.
For non-zero asymmetry parameter $\mathcal{E}$, thermal conductivity contains both the "intrinsic" 
(unrelated to charge transport) term, and the term proportional to electric conductivity.

Qualitatively new results are obtained for the intermediate temperature interval  $ T_{FL} \ll T \leq \tilde\Gamma $ where fluctuations of the reparametrization soft mode make strong effect upon kinetic properties. Physically this mode
can be understood as space-time fluctuations of the local energy density in strongly interacting system.
It is thus natural that these fluctuations affect  thermal conductivity, which is strongly modified
already in the static \textit{dc} limit, studied previously within saddle-point approximation~\cite{song2017strongly}.
Namely, the Lorenz ratio $L = \kappa/T\sigma$ is found to be temperature - dependent, 
it also depends on particle-hole asymmetry parameter $\mathcal{E}$, see Eqs.(\ref{L1},\ref{L2}).
Electric conductivity in the \textit{dc} limit is not sensitive to  soft-mode fluctuations and demonstrate
"strange-metal" scaling $\sigma \propto 1/T$ in agreement with Ref.~\cite{song2017strongly}, as well as kinetic 
relaxation time $\tau(T) \sim 1/T$.

Soft-mode fluctuations become even more important for the high-frequency $\omega \gg T$ transport that is strongly
non-equilibrium since $\omega\tau(T) \gg 1$. We found that
high-frequency electric conductivity contains a pole contribution (\ref{sigmaHigh}), with the dispersion relation of 
nearly-linear form, Eq.(\ref{0sound}), which seems to be a "strongly - correlated partner" of a zero-sound mode known for
Fermi-liquids. This mode exists in an intermediate temperature range $T_{FL} < T \ll w$. We emphasize that
"sound velocity" given by Eq.(\ref{s}) does not depend on  interaction parameter $J$, although existence
of such a mode is related to strong interactions. Lifetime of these nearly-sound modes $\tau_s(p)$ is relatively
long in the $\omega/T \gg 1$ limit: $\tau_s(p)\omega(p) \approx \ln\frac{sp}{T}$.

High-frequency thermal conductivity contains (at $\mathcal E \neq 0$) similar pole contribution; in addition,
the intrinsic contribution to high-frequency thermal conductivity $\tilde\kappa(\omega)$ contains a weakly dispersive
resonance at  frequency $\omega_0 \sim \Gamma$,  see Eqs.(\ref{kappa-high},\ref{omega01}). This resonance is 
similar to  the one found in Ref.~\cite{lunkin2022non} for parametric excitation of the SYK quantum dot;
its quality factor is $ Q \sim \ln\frac{\omega_0}{2\pi T}$.  In general, our results are fully compatible with
those obtained in Ref.~\cite{lunkin2022non}; the difference between Majorana and complex fermion models consists
 just in the presence of charge transport in the latter one.

\begin{figure}
	\centering
	\includegraphics[width=0.67\linewidth]{SYK_regions}
	\caption{Transport regimes}
	\label{fig:sykregions}
\end{figure}

Fig.~\ref{fig:sykregions} provides schematic view of the "phase diagram" of the model we studied, in the case $\Gamma \gg w$.
We used logarithmic axis  for the ratios $\omega/w$ and $T/w$
to make our phase diagram compact. Depending on relations between temperature $T$, 
frequency $\omega$ and transport bandwidth $w$, various transport regimes discussed above can found in the phase diagram;
they are indicated in the Figure.~\ref{fig:sykregions} in the self-explanatory way; blue region in the centre corresponds to crossover
between various specific regimes.


Recently we became aware of interesting preprint~\cite{Behnia2022dynamic} discussing phenomenology of several
 strongly interacting fermionic systems, like overdoped La$_{2-x}$Sr$_x$CuO$_4$ and $^3$He under strong pressure.
One of important observations made in Ref.\cite{Behnia2022dynamic} is that both these very different objects demonstrate
significant deviation from the  "universal"  Kadowaki-Woods ratio $A/\gamma^2$ known for many various electronic materials
including strongly correlated ones; here $A$ is the coeffcient in the $T^2$ law of resistivity known for Fermi-liquids,
$R(T) = R_0 + AT^2$, and $\gamma$ enters specific heat, $C(T)=\gamma T$. Namely, for La$_{1.67}$Sr$_{0.33}$CuO$_4$ the ratio
$A/\gamma^2$ is anomalously large; for the compressed $^3$He one cannot define conductivity, but the analysis of
transport relaxation rate leads to very similar conclusion as for La$_{1.67}$Sr$_{0.33}$CuO$_4$.

Kadowaki-Woods ratio was discussed in Ref.~\cite{song2017strongly} where it was found to be universal;
however in the model~\cite{song2017strongly} the only quadractic terms were the tunneling ones.
 Here we would like
to notice that "heavy Fermi-liquid"  with anomalously large Kadowaki-Woods ratio can be obtained within
SYK$_4$+SYK$_2$ model we consider here, if we assume $\Gamma \gg w$. Indeed, in this model conductivity
scales~\cite{song2017strongly} as $\sigma \sim \frac{w^2}{J T}$  at $T > T_{FL} = \Gamma^2/J$.  
Below $T_{FL}$  we expect Fermi-liquid behavior to set in, with $R(T) = AT^2$, thus by continuity arguments
 $A \approx (J/w\Gamma)^2$. Then, the specific heat coefficient $\gamma \approx J/\Gamma^2$, and the
Kadowaki-Woods ratio is $A/\gamma^2 \sim \Gamma^2/w^2 \gg 1$  and does not depend on $J$.
We conclude that  the model studied in this paper may be relevant for the description of extremely strongly
interacting Fermi systems, like those discussed in Ref.~\cite{Behnia2022dynamic}.




\section*{Acknowledgements}
We are grateful to Alexei Kitaev, Kamran Behnia and  Grigory Tarnopolsky for useful discussions. We are also grateful to Aleksandr Svetogorov for his comments to the text.


\paragraph{Funding information}


Research of A.V.L. was partially supported by the Basis Foundation,  by 
the Basic research program of the HSE and by the RFBR grant \# 20-32-90057.

\begin{appendix}
\label{sec: appendix_A}

\section{Action and continuity equations in the presence of the  sources}
In this appendix we will cover gaps of the main text. Particular, we have not written action with sources to simplify formulas. 


The action for fields $G$ and $\Sigma$ in the presence of sources has the following form:

\begin{eqnarray}
&S=-iN\sum\limits_{{\textbf r}} \int\limits_{\mathcal{C}} dt_1 dt_0 (\Sigma(t_1,t_0|{\bf r})+\Sigma_{free}(t_1,t_0) )G(t_0,t_1|{\bf r})-iN~tr\log\{\hat{\Sigma}\}\nonumber\\
&iN  \int\limits_{\mathcal{C}} dt_1 dt_0 \sum\limits_{{\textbf r}}  \left\{\frac{J^2}{4}\left[G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r})\right]^2+\frac{\Gamma^2}{2}G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r})\right\}(1+ L({\textbf r},t_1)+L({\textbf r},t_0))\nonumber \\
&+iN w^2  \int\limits_{\mathcal{C}} dt_1 dt_0 \sum\limits_{{\textbf r},\delta {\textbf r}}  G(t_1,t_0|{\bf r})G(t_0,t_1|{\bf r}+\delta {\bf r})(1+\frac{L({\textbf r},t_1)+L({\textbf r}+\delta{\bf r},t_1)}{2}+\frac{L({\textbf r},t_0)+L({\textbf r}+\delta{\bf r},t_0)}{2})\times \nonumber \\ & \exp\{i\varphi^{(U)}({\bf r}+\delta {\bf r},t_1)-i\varphi^{(U)}({\bf r},t_1)-i\varphi^{(U)}({\bf r}+\delta {\bf r},t_0)+i\varphi^{(U)}({\bf r},t_0)\}\nonumber \\
\label{eq:action G and sources}
\end{eqnarray}

Our aim will be to find the action for the soft modes in the presence of sources. Using this action we apply procedure \ref{eq:soft mode action} from the main text.  

After substitution of the field $G$  from the soft manifold we can find that terms proportional to $\Gamma$ or $w$ contains logariphimcaly divergent integral (for $q=4$ ). This devergency does not affect on the action of the soft mode but to find one we need to use the following regularization procedure. We should consider the problem with general $q$ and expand action up to the second order in soft mode fluctuation, after than we need to take the limit $q\rightarrow 4-0$. The similar procedure was also performed in \cite{lunkin2022non}. As a result we have the  quadratic action for soft mode, which can be written in the Fourier domain as:
\begin{eqnarray}
&S_{final}=\frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left(\hat{\varphi}_{{\bf p},\Omega}-\mathcal{E}(\hat{f}_{{\bf p},\Omega}-\frac{i}{\Omega}\hat{L}_{{\bf p},\Omega})\right)^\dagger \left[\hat{\mathcal{G}}^{(\varphi)}_0(\Omega)\right]^{-1}\left(\hat{\varphi}_{{\bf p},\Omega}-\mathcal{E}(\hat{f}_{{\bf p},\Omega}-\frac{i}{\Omega}\hat{L}_{{\bf p},\Omega})\right) \nonumber \\
& \frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left(\hat{f}_{{\bf p},\Omega}-\frac{i}{\Omega}\hat{L}_{{\bf p},\Omega}\right)^\dagger\left[\hat{\mathcal{G}}^{(f)}_0(\Omega)\right]^{-1}\left(\hat{f}_{{\bf p},\Omega}-\frac{i}{\Omega}\hat{L}_{{\bf p},\Omega}\right)^\dagger\nonumber\\
&-\frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left\{\left(\hat{f}_{{\bf p},\Omega}-\frac{2i}{\Omega}\hat{L}_{{\bf p},\Omega}\right)^\dagger  \left[\hat{\Sigma}^{(f)}(\Omega)+\hat{\Sigma}_w^{(f)}(\Omega)\right]\left(\hat{f}_{{\bf p},\Omega}-\frac{2i}{\Omega}\hat{L}_{{\bf p},\Omega}\right)\right. \nonumber \\
&\left. +\hat{f}_{{\bf p},\Omega}^\dagger   \left[\hat{\Sigma}^{(f)}_{w}(\Omega,{\bf p})-\hat{\Sigma}^{(f)}_{w}(\Omega,{\bf 0})\right]\hat{f}_{{\bf p},\Omega} \right\} \nonumber \\
&-\frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left\{ \hat{L}_{{\bf p},\Omega}^\dagger  \hat{\Sigma}^{(L,f)}_{t}(\Omega,{\bf p})\hat{f}_{{\bf p},\Omega} +\hat{f}_{{\bf p},\Omega}^\dagger  \hat{\Sigma}^{(f,L)}_{t}(\Omega,{\bf p})\hat{L}_{{\bf p},\Omega} \right\} \nonumber \\
&-\frac{1}{2}\int \frac{d^d {\bf p}d \Omega}{(2\pi)^{d+1}} \left(\hat{\varphi}_{{\bf p},\Omega}-\mathcal{E}\hat{f}_{{\bf p},\Omega}+\hat{\varphi}^{(U)}_{{\bf p},\Omega}\right)^\dagger  \hat{\Sigma}^{(\varphi)}_t(\Omega,{\bf p})\left(\hat{\varphi}_{{\bf p},\Omega}-\mathcal{E}\hat{f}_{{\bf p},\Omega}+\hat{\varphi}^{(U)}_{{\bf p},\Omega}\right)
\label{eq:soft mode action and sources full}
\end{eqnarray}
This action without sources was written in the main text as: (\ref{eq:soft mode action and sources}). Besides operators $\mathcal{G}$ and $\hat{\Sigma}$ mentioned in the main text we also introduced the following operators:
\begin{eqnarray}
\Sigma^{(f,L)}_{t,(R)}(\Omega,{\bf p}) = -\Sigma^{(L,f)}_{t,(R)}(\Omega,{\bf p}) = C_w\sum\limits_{\delta{\bf r}} i \Omega \left(\frac{{\bf p}\delta{\bf r}}{2} \right)^2 \psi(\Omega)
\end{eqnarray}

The variation of the above action leads us to continuity equations (see (\ref{subsec: Noether's theorem},\ref{subsec: Continuity equations})) which has the following form in the presence of sources:

\begin{eqnarray}
&0 =-\frac{\delta S_{final}}{\delta \varphi_{{\bf p},\Omega}^{(q)\dagger}}= -iC_{Q} \Omega\left(V_{{\bf p},\Omega}+\mathcal{E}L_{{\bf p},\Omega}\right) +   8 \tilde{\sigma}p^2  \frac{\tilde{\psi}(\Omega)}{-i\Omega}\left(V_{{\bf p},\Omega}+U_{{\bf p},\Omega}\right)\nonumber\\
&0=-\frac{\delta S_{final}}{\delta f_{{\bf p},\Omega}^{(q)\dagger}}= \mathcal{E}\frac{\delta S_{final}}{\delta \varphi_{{\bf p},\Omega}^{(q)\dagger}} - iC_{E} \Omega (\Omega^2+1) \left(\tau_{{\bf p},\Omega}-L_{{\bf p},\Omega}\right)  \nonumber \\
&+  \frac{C_\Gamma}{2}  i\Omega \psi(\Omega)\left(\tau_{{\bf p},\Omega}-2 L_{{\bf p},\Omega}\right) + \frac{\tilde{\sigma}}{4}p^2 \left( \left((\Omega^2+1) \frac{\tilde{\psi(\Omega)}}{-i\Omega}+i\Omega\psi(\Omega)\right) \tau_{{\bf p},\Omega} + i 4\Omega  \psi(\Omega)L_{{\bf p},\Omega}\right)\nonumber\\
\label{eq:continuity  equation}
\end{eqnarray}
These equations describes charge and energy conservations respectively. 
Here we also introduced notations:
\begin{eqnarray}
V_{{\bf p},\Omega} \equiv  -i \Omega\left(\varphi^{(cl)}_{{\bf p},\Omega}-\mathcal{E} f^{(cl)}_{{\bf p},\Omega}\right) \qquad \tau_{{\bf p},\Omega} \equiv  -i \Omega f^{(cl)}_{{\bf p},\Omega}
\end{eqnarray}

We need to solve continuity equations to find changes of fields $\varphi$ and $f$ caused by sources. Using this information we can find currents using (\ref{currents-def}). It leads us to the  expression (\ref	{eq:currents}) which is the main result of the text. 
 \end{appendix}











