SYK2-FL and impurity-like regimes, and where the gap
$\Delta_{0}\approx\Delta(i\omega_{1})$ is calculated at each $g$ from the
self-consistent imaginary-axis solution at the first Matsubara frequency; see
Fig. 10. (b) Example of scaling collapse of
$\omega_{0}\rho_{S}(T)/(\Theta_{L}\omega_{0}zt_{0}^{2})\propto g^{-4}$ as a
function of $g$ in the strong-coupling limit: the curves for
$t_{0}=\left\\{0,1,2\right\\}$ all collapse to the single-dot result at large
$g$. All other parameters are the same as in panel (a).

In SYK2-FL regime $\tilde{\Omega}\mapsto\sqrt{2z}t_{0}$, and we have

$\rho_{S}(T)=\Theta_{L}\Delta\left[1-2f_{FD}(\Delta)\right].$ (68)

Therefore, the stiffness is proportional to the superconducting gap in the
regime $\Delta/(k_{B}T)\gg 1$. Indeed, at $T=0$ using Eq. (43) we obtain

$\rho_{S}(0)=\Theta_{L}\Delta_{0}=\Theta_{L}2\omega_{0}e^{-\frac{1}{\bar{\lambda}}},$
(69)

with coupling constant (36), which shows that $\rho_{S}(0)$ is exponentially
suppressed as $g$ is decreased in SYK2-FL regime. Such exponential suppression
is shared by the quasiparticle weight (56). Let us remark that the result (69)
is consistent with the dirty limit of BCS theory: there, the stiffness is
affected by nonmagnetic disorder, but the pairing temperature and the zero-
temperature gap are not. This difference is reflected by our findings for
$T_{c}$ and $\Delta_{0}$ in SYK2-FL regime – see Eqs. (35) and (43) – which
are in agreement with the respective BCS formulae.

The result (69) corresponds to the dashed light-green curves in Fig. 16(a),
and it is there compared with the same numerical data points reported in Fig.
15(b). Comparing Eq. (69) with Eqs. (43) and (35), we deduce that the ratio
between the energy scale given by the zero-temperature stiffness, and the
mean-field transition temperature $k_{B}T_{c}$, is a universal constant in
SYK2-FL regime:

$\frac{\rho_{S}(T)}{\Theta_{L}k_{B}T_{c}}=\frac{\Delta_{0}}{k_{B}T_{c}}=\frac{\pi}{e^{\gamma}}=0.5669329586555488,$
(70)

where $\gamma$ is again the Euler-Mascheroni constant. A constant ratio of the
order of Eq. (70) between $k_{B}T_{c}$ and stiffness is indeed retrieved in
the SYK2-FL regime, from the full numerical solution of the Eliashberg
equations. This is shown in Fig. 3 of the companion paper Ref. 159, and there
further discussed.

In the impurity-like regime $\tilde{\Omega}\mapsto\Omega_{0}$, and Eq. (66)
gives

$\displaystyle\rho_{S}(T)$
$\displaystyle=\frac{\Theta_{L}zt_{0}^{2}}{2\left[8g^{2}\omega_{0}/(3\pi)\right]^{2}}\Delta\left[1-2f_{FD}(\Delta)\right]$
$\displaystyle=\frac{\Theta_{L}}{2}zt_{0}^{2}\left(\frac{3\pi}{8}\right)^{2}\frac{1}{\omega_{0}^{2}g^{4}}\Delta\left[1-2f_{FD}(\Delta)\right].$
(71)

In particular, at $T=0$ we have

$\rho_{S}(0)=\frac{\Theta_{L}}{2}zt_{0}^{2}\left(\frac{3\pi}{8}\right)^{2}\frac{1}{\omega_{0}^{2}g^{4}}\Delta_{0}.$
(72)

In the impurity-like regime, the gap asymptotically reaches a constant in the
infinite-interaction limit, in accordance with Eq. (45). Therefore, Eqs. (72)
and (45) tell us that the stiffness goes like $\rho_{S}(0)\propto\bar{g}^{-4}$
at large interaction $g\gg 1$. We can also make sense of this conclusion by
realizing that $Z(i\omega_{n})\propto\bar{g}^{2}$ in the impurity-like regime
[105], and the stiffness contains
$\left[Z(i\omega_{n})\right]^{-2}\propto\bar{g}^{-4}$ as seen from Eq. (65a).
Therefore, Eq. (72) predicts that the normalized quantity
$\rho_{S}(0)/(\Theta_{L}zt_{0}^{2})$ is independent from hopping in impurity-
like regime: this is exactly what we see in Fig. 16(b), where the numerical
data points for $t_{0}=\left\\{1,2\right\\}\omega_{0}$ (red squares and light-
blue circles) collapse onto the $t_{0}=0$ single-dot results (golden
triangles) at large $g$.

In summary, we find that the low-temperature phase stiffness is controlled by
the dependence of $Z(i\omega_{n})$ and $\Delta_{0}$ on coupling: in the Fermi-
liquid SYK2-FL regime, $\rho_{S}(0)\propto\Delta_{0}$ is exponentially
diminished as we decrease $g$, because $\Delta_{0}$ follows the BCS formula;
in the impurity-like regime,
$\rho_{S}(0)\propto\sum_{i\omega_{n}}\left[Z(i\omega_{n})\right]^{-2}\propto
g^{-4}$, because the gap $\Delta_{0}$ is approximately constant and the
imaginary-axis weight $Z(i\omega_{n})\propto g^{2}$. In the intermediate
coupling regime, the stiffness reaches a maximum. In fact, this maximum occurs
precisely at the SYK2-FL/SYK-NFL crossover energy $\omega_{c}$, previously
found in Sec. V and marked by the red and blue arrows in Fig. 15(b) for
$t_{0}=\left\\{1,2\right\\}\omega_{0}$. This is the second fundamental result
of this section: the phase stiffness is maximum at the crossover between the
SYK-NFL phase (non-Fermi liquid) and the SYK2-FL regime (disordered Fermi
liquid). Such non-monotonic evolution is reminiscent of the doping evolution
of the superfluid density extracted from muon relaxation experiments in
cuprate high-temperature superconductors [14, 22]. In the same experimental
works [14, 22], a correlation between the superfluid density and the
condensation energy, extracted from calorimetric measurements, was discovered.
In the same spirit, in the next section we conclude our investigation of the
superconducting properties of our model (III) with a discussion of the grand
potential and the condensation energy in the paired phase.

##  XII Condensation energy

Figure 17:  Grand potential and condensation energy as a function of coupling
$g$ for different hoppings $t_{0}$, numerically computed from Eqs. (73) and
(VI). (a) Normalized grand potential $\Omega/(\mathscr{N}\omega_{0})$ in the
normal state (circles) and in the superconducting state (squares) at
temperature $k_{B}T=0.01\omega_{0}$, for hopping
$t_{0}/\omega_{0}=\left\\{0,1,2\right\\}$. (b) Condensation energy
$\Delta\Omega/(\mathscr{N}\omega_{0})$, evaluated from the difference between
the normal-state and superconducting-state grand potentials in panel (a), with
values reported on the left $y$ axis, as indicated by arrows. The dashed lines
are the SYK2-FL/SYK-NFL crossover energies $\omega_{c}/\omega_{0}$ at a given
hopping, according to Sec. V, with values shown on the right $y$ axis, as
indicated by arrows.

The superconducting transition in the Yukawa-SYK model entails a lowering of
the total energy of the system, expressed by the grand potential $\Omega$,
with respect to the normal-state grand potential (VI). In a strict sense is
the condensation energy an ill-defined concept, as the rigorous solution of a
statistical mechanics problem for a superconductor knows nothing about the
behavior of the unstable normal state. The exception may be a normal state
stabilized by a magnetic field, which may however dramatically alter the
system on its own right. On the other hand, the condensation energy can be
well-defined and is physically insightful within mean-field theory. Thus,
keeping the formal caveat in mind we will, in what follows, analyze the energy
gain of superconducting solutions if compared to the corresponding normal-
state case. In the single-dot limit, such condensation energy was analyzed in
Ref. 105. To calculate the same quantity on the lattice, we have to derive an
expression for the grand potential $\Omega^{{\rm sc}}$ for $T<T_{c}$. We
report the technical derivation in Appendix F.2, while the final result is
quoted here below:

$\displaystyle\frac{\Omega^{{\rm sc}}}{\mathscr{N}}$ $\displaystyle=$
$\displaystyle
k_{B}T\sum_{i\omega_{n}}\left\\{2\log\left[\frac{\sqrt{G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})}}{G_{0}(i\omega_{n})}\right]\right.$
(73) $\displaystyle-$
$\displaystyle\left.\frac{zt_{0}^{2}}{2}\left[G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})\right]\right\\}$
$\displaystyle-$ $\displaystyle
k_{B}T\sum_{i\Omega_{n}}\left\\{\frac{1}{2}\log\left[\frac{D(i\Omega_{n})}{D_{0}(i\Omega_{n})}\right]+D(i\Omega_{n})\Pi(i\Omega_{n})\right\\}$
$\displaystyle-$ $\displaystyle
2k_{B}T\log\left[1+e^{\mu/(k_{B}T)}\right]-\frac{\omega_{0}}{4}$
$\displaystyle-$
$\displaystyle\frac{k_{B}T}{2}\log\left[1-e^{-\omega_{0}/(k_{B}T)}\right].$

As in the normal-state case, the grand potentials of free fermions and bosons
have been simultaneously added and subtracted in Eq. (73), to improve the
numerical convergence of the Matsubara sums. Fig. 17(a) shows the grand
potential in the superconducting state (colored squares) and in the normal
state (colored circles), as a function of coupling and at
$k_{B}T=0.01\omega_{0}$, for different hoppings. The data points are
numerically calculated from Eqs. (73) and (VI) respectively, and using the
saddle-point propagators self-consistently evaluated from the saddle-point
equations (4). We see that the difference between the normal-state and the
condensed-state results is very small on the scale of the grand potentials
themselves.

Subtracting the normal-state result (VI) from the superconducting grand
potential (73) gives the condensation energy:

$\frac{\Delta\Omega}{\mathscr{N}}=\frac{\Omega^{{\rm
sc}}}{\mathscr{N}}-\frac{\Omega^{{\rm ns}}}{\mathscr{N}}.$ (74)

Eq. (74) is employed together with the grand potentials in Fig. 17(a), to
obtain the numerical results shown in Fig. 17(b). We observe that the almost
linear evolution $\Delta\Omega\propto g$ at small $g$, in the single-dot case
(golden circles) and in the SYK-NFL phase, is replaced by an exponential
suppression of the condensation energy in the SYK2-FL regime (red and light-
blue circles), for finite hopping. Indeed, $\Delta\Omega$ starts to
significantly increase in magnitude only for couplings above the SYK2-FL/SYK-
NFL crossover, signaled by $\omega_{c}$; see Sec. V. The crossover energy
$\omega_{c}$ is marked by the dashed red and blue curves for hopping
$t_{0}=\left\\{1,2\right\\}\omega_{0}$, respectively. Therefore, the non-Fermi
liquid phase is characterized by a significantly higher condensation energy
than the disordered Fermi-liquid regime; such phenomenon might bear observable
consequences for the specific-heat jump at the superconducting transition in
non-Fermi liquid superconductors [199, 200, 201].

On the other hand, the data points in the impurity-like regime at finite
hopping are superimposed to the single-dot results, since $zt_{0}^{2}\ll
g^{2}\omega_{0}^{2}$ and hopping is an irrelevant perturbation. Then, we find
the same slow decrease of $\Delta\Omega$ with increasing $g$ as found in the
single-dot limit [105].

Notice that the exponential-like suppression of $\Delta\Omega$ in the SYK2-FL
phase, as well as the decrease with $g$ in the impurity-like state, remarkably
correlate with the coupling evolution of the quasiparticle weight and of the
stiffness, as summarized in Fig. 2. This is the fundamental result of the
present section: in our Yukawa-SYK superconductor on a lattice, the superfluid
stiffness, the condensation energy, and the quasiparticle weight are all
correlated with each other and share similar behaviors as functions of
fermion-boson coupling. Given the analogy with the experimental findings on
cuprate superconductors [14, 22], the qualitative trends proposed in the
present model may transcend the specific aspects of SYK physics, and be
generic for non-Fermi liquid superconductors. We discuss these speculations,
together with future useful extensions of our lattice model, in the conclusive
remarks.

##  XIII Conclusions and perspectives

In summary, we formulated and solved a model for spinful fermions with
multiple flavors, interacting at local sites through all-to-all random
couplings to many Einstein phonon modes and through a random coherent single-
particle hopping between nearest-neighbor sites. This forms a lattice of
Yukawa-SYK dots. Both couplings and hoppings follow a distribution function
with zero mean and finite variance. Our mean-field equations on the saddle
point of the disorder-averaged action (2) are exact in the limit of a large
number of fermion and boson flavors; they assume a diagonal structure in
replica space for the disordered interactions, an equal number of fermions and
boson, and particle-hole symmetry (zero chemical potential). They are
equivalent to self-consistent Eliashberg equations, where fermion-boson
interaction is at once responsible for the destruction of fermionic
quasiparticles (NFL state), the softening of the renormalized boson frequency
to criticality, and a superconducting state at low temperatures if every
disorder realization of the couplings preserves time-reversal symmetry.

In the normal state, the competition between coupling and hopping leads to
crossovers between FL and NFL phases. In addition to the high-temperature
classical-gas phase, the low-temperature SYK-NFL phase, and the intermediate-
temperature impurity-like phase, all identified and studied in the single-dot
limit in Refs. 105, 106, 107, we find crossovers to a disordered Fermi liquid,
the SYK2-FL state, when hoppings prevail over on-site coupling. The
characteristic crossover energies, estimated in Sec. V, are decreasing
functions of coupling at fixed hopping, and they correspond to observable
changes in the fermionic and bosonic spectral functions: in the SYK-NFL phase,
the fermion and boson propagators follow a power-law form governed by the
exponent $\Delta$, and by increasing hopping this evolves into the
semicircular Wigner spectral function of SYK2-FL fermions with bandwidth
controlled by hopping, with almost free bosons and weak fermion-boson
coupling. The crossover from the impurity-like regime is more subtle, as here
the fermions are also broadened into a Wigner spectral function controlled by
coupling, but bosons behave like static impurities, with a sharp but soft
renormalized frequency. This frequency stiffens towards the SYK2-FL regime.
The crossovers leave qualitative differences in the normal-state entropy,
which is finite at $T=0$ in the single-dot limit of negligible hopping (SYK-
NFL and imp.-like phases), and vanishes for any finite hopping.
Superconductivity emerges at low temperature in all regimes. The critical
temperature saturates to a fraction of the bare boson frequency $\omega_{0}$,
$k_{B}T_{c}\sim 0.11\omega_{0}$ [202, 203], in the strong-coupling imp.-like
regime where fermions are fully incoherent [105, 106]. Conversely, in the
weak-coupling SYK2-FL phase we retrieve a disordered version of BCS theory,
where $T_{c}$ and the gap at zero temperature and energy follow the BCS
formulae with a coupling constant $\bar{\lambda}\propto g^{2}\omega_{0}/t_{0}$
that depends on both coupling and hopping. The tell-tale signature of the
crossovers in the superconducting spectral functions is the alteration and
eventual suppression of peak-dip-hump features, in the SYK-NFL and imp.-like
states, with increasing hopping towards the SYK2-FL regime. Such resonances
are self-trapping bound states of the interacting Cooper-pair fluid [105].

Our most essential result is that the superconducting state of our model has a
finite phase stiffness. The latter is calculated through the electromagnetic
linear response function to an external vector potential, which couples to
phase fluctuations of the order parameter away from the saddle-point through
electric charge. At vanishing temperature and weak coupling (SYK2-FL regime),
the stiffness follows the exponential evolution of the gap with coupling
$\bar{\lambda}$, so that the ratio between $T_{c}$ and the stiffness
corresponds to the “universal” BCS ratio $\pi/e^{\gamma}\approx 0.567$. On the
contrary, at strong coupling (imp.-like phase) the gap saturates, but the
stiffness decreases as the inverse fourth power of coupling due to the
decrease of the quasiparticle weight. Remarkably, we find a correlation in the
evolution with coupling of the stiffness, the quasiparticle weight, and the
condensation energy. Their magnitudes all peak at the crossover between NFL
and FL behavior. Such correlation is reminiscent of the experimental
measurements of the same quantities in cuprate superconductors [14, 22]. More
generally, our results demonstrate that the low-temperature condensed phase of
the Yukawa-SYK model is indeed superconducting, as it exhibits perfect
diamagnetism according to the Meißner effect. Therefore, our calculations
offer a suitable platform to compare the behavior of the stiffness within
similar SYK-derived models, as well as in different classes of models for
quantum-critical superconductors [66, 204, 67, 205, 185, 206, 76, 207, 208,
209, 210, 203, 211]. For the Yukawa-SYK model in particular, the single-dot
limit corresponds to a superconductor with negligible interdot hopping and
vanishingly small stiffness per flavor, which represents a special case of our
lattice embedding.

The difference between a saturating mean-field $T_{c}$ and a decreasing
stiffness at strong coupling, as well as a decreasing condensation energy in
the same regime, suggest strong superconducting phase fluctuations in the
system [19, 105]. If we assume that phase fluctuations are the driving cause
of the disappearance of superconducting coherence, the phase stiffness
obtained in our theory would suggest that the fluctuation-corrected transition
temperature $T_{\Theta}$ vanishes as $T_{\Theta}\sim g^{-4}$ for
$g\rightarrow+\infty$. However, performing a quantitative analysis of
superconducting phase fluctuations requires to include $1/\mathscr{N}$
corrections with respect to our saddle-point results, which hold for
$\mathscr{N}\rightarrow+\infty$. The effect of such corrections has recently
been studied in an SYK-type system which superconducts due to an attractive
local Hubbard interaction: the $1/\mathscr{N}$ terms corresponding to phase
fluctuations lead to a pseudogap behavior [132]. A similar phenomenon, that
has been pointed out in various microscopic models for quantum critical
superconductors [66, 99], could also be revealed at order $1/\mathscr{N}$ in
our model. In this respect, the mean-field $T_{c}$ calculated in this work
would not be the real superconducting temperature where zero resistivity
appears, but rather a pseudogap temperature. A future analysis of such
pseudogap phase is in order, in the light of the analogous phenomenon found in
strange metals [10].

Since the stiffness is connected to the electrodynamic response, it is
interesting to investigate whether the evolution of the stiffness with
coupling and hopping leaves distinctive traces in the long-wavelength optical
conductivity $\sigma(\omega)$ of our lattice model in the superconducting
state as a function of frequency $\omega$. This is the subject of a companion
paper [159], in which we demonstrate that the removed spectral weight for
energies $\omega<2\Delta_{S}$, where $\Delta_{S}$ is the spectroscopic gap, is
a direct signature of the different behavior of the stiffness in the SYK2-FL,
SYK-NFL and imp.-like regimes. Moreover, the Cooper-pair bound states,
appearing in the spectral functions of the SYK-NFL and imp.-like phases, are
reflected in resonances of $\sigma(\omega)$, which disappear in the crossover
towards the SYK2-FL regime. Then, it is crucial to compare our results with
optical spectroscopy measurements in strongly correlated superconductors such
as cuprates, which show a Planckian-like power-law behavior of the
conductivity in energy and temperature [212, 213, 214, 215, 216].

Numerous future developments and generalizations of our theory could provide
further connections with experiments and with complementary theoretical
perspectives.

A thermodynamic quantity which is equally found to correlate with the
stiffness in cuprate superconductors is the heat-capacity jump at $T_{c}$ [14,
22]. This quantity can be computed within our formalism, from the second
derivative of the condensation energy at $T=T_{c}$, and if the correlation
with the stiffness persists in our model, the latter would see its position
strengthen as a suitable effective toy model for quantum critical
superconductors and strange metals. Nevertheless, other ingredients related to
strange-metal behavior, such as Mott-insulator physics [217], are not included
in our model. In this respect, it would be interesting to quantitatively
compare predictions for the superfluid density evolution stemming from Mott-
based pictures with our SYK-based approach.

Concerning the still unexplored parameter space, it would be interesting to
extend the lattice calculations, in particular for $T_{c}$, the gap, and the
phase stiffness, at finite chemical potential $\mu\neq 0$. A thorough study of
this kind in the single-dot limit became recently available [107], and could
serve as a springboard reference to assess whether our found correlation among
stiffness, quasiparticle weight and condensation energy persists as a function
of $\mu$. Indeed, valence transitions were found by varying $\mu$ in the
Yukawa-SYK dot [153, 104] and in purely fermionic versions of the SYK model
[218, 219, 220, 221]. These phenomena could serve as an ideal playground for
comparison of the thermodynamic and optical observables of the Yukawa-SYK
lattice model with experiments on one hand, and with numerical methods such as
Functional Renormalization Group (FRG) [221] and DMFT [141, 142] on the other
hand. Even more precisely, a self-consistent calculation of the chemical
potential [222, 223, 224, 225], $T_{c}$, and the stiffness at fixed density
$\mathscr{Q}\in\left[0,1\right]$ per fermion flavor [153] would allow one to
explore the phase diagram as a function of $\mathscr{Q}$, which is analogous
to chemical doping; this exploration would permit a direct comparison of our
theoretical results with the observed experimental correlations among
thermodynamic and spectroscopic data in cuprates [14, 22]. Another insightful
additional variable to tune is the pair-breaking parameter
$\alpha\in\left(0,1\right)$, which permits a continuous interpolation between
complex and real fermion-boson coupling constants, that violate or preserve
time-reversal symmetry, respectively [106, 119, 226]. Specifically as a
function of hopping, how the crossover between FL and NFL physics affects the
QCP, where $T_{c}(\alpha)\rightarrow 0$ with BKT-like evolution, remains an
open question. The numerical results for thermodynamic and dynamical
observables as a function of $\alpha$ also offer a direct comparison with the
AdS/CFT correspondence: recently, a holographic dual of the Yukawa-SYK model
in the single-dot limit was explicitly constructed, and shown to be equivalent
to the Eliashberg formulation of the same model in the regime
$T_{c}(\alpha)\rightarrow 0$ [119]. The comparison of quantities like the
pairing susceptibility, computed from the numerical saddle-point solution and
from holography, could extend the newly found proof of the correspondence
between the dual theories [119].

Another crucial generalization of our formalism involves the extension to
anisotropic pairing symmetries, for instance through coupling with a spin-1
boson $\phi_{\bm{k}}$ [105]. Anisotropic pairing is also relevant for
dispersive fermions, for instance Dirac fermions which could model
quasiparticles in flat-band systems like graphene [227, 228, 229, 230] at the
charge neutrality point, or nodal quasiparticles in cuprate superconductors
[231].

The foundational assumptions on which our results rest also necessitate
further clarification, in particular the replica-diagonal ansatz for disorder-
averages assumed in the large-$\mathscr{N}$ limit [232]. For instance, in a
bosonic variant of the SYK model it was shown that anharmonic boson-boson
interactions are able to break the replica symmetry and give rise to a glassy
phase [233, 234].

Finally, since our results are valid on the saddle point and at equilibrium,
nonequilibrium numerical studies of the Yukawa-SYK model are crucial to assess
the boundaries in parameter space where our results apply, and viceversa where
new physics arises, e.g., for the electrodynamic response and the pairing
susceptibility [235, 236, 237, 238, 134].

All in all, our lattice theory offers an exactly solvable and analytically
controlled scenario to study the interplay between NFL and FL phases, the
thermodynamics, and the electromagnetic response in quantum-critical
superconducting systems. It represents a platform based on the Eliashberg
formalism, to incorporate fluctuations beyond the mean-field level and to
investigate the properties of unconventional superconductors.

##  XIV Acknowledgments

D. V. acknowledges insightful discussions with C. Berthod, J. Zaanen, L.
Benfatto, B. Goutéraux, R. Riva, N. Stegani, F. Lombardi, S. Caprara, R.
Arpaia, G. Ghiringhelli, D. van der Marel, A. Chubukov, C. Bernhard, P.
Werner, M. Müller, R. Morin, V. Stangier, D. Hauck, G. Mazza, D. Kennes, L.
Grunwald, R. Ojajärvi, R. Willa. D. V. and J. S. also acknowledge fruitful
discussions with S. Sachdev, I. Esterlis, C. Li, A. Patel, and H. Guo. D. V.
acknowledges partial support by the Swiss National Science Foundation (SNSF)
through the SNSF Early Postdoc.Mobility Grant P2GEP2_181450, and by the
Mafalda cluster of the University of Geneva, on which some of the calculations
were performed. G. A. I. and J. S. were supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - TRR 288 - 422213477
Elasto-Q-Mat (project A07).

##  Appendix A Disorder-averaged effective action for the SYK model

In this appendix we derive the disorder-averaged effective action of the SYK
model, Eq. (III). The partition function of the $m$-replicated system is given
by

$\displaystyle\overline{\mathscr{Z}^{m}}=\int\mathcal{D}\Psi
e^{-\sum_{a}\mathscr{S}_{0}[\Psi]}\int
dg_{ij,k}\varrho(g_{ij,k})e^{-\sum_{a}\mathscr{S}_{g}[\Psi]}\int
dt_{ij,\bm{x}\bm{x}^{\prime}}\varrho(t_{ij,\bm{x}\bm{x}^{\prime}})e^{-\sum_{a}\mathscr{S}_{t}[\Psi]},$
(75)

where $\Psi=\\{\hat{c}_{i\sigma a\bm{x}},\hat{c}_{i\sigma
a\bm{x}}^{\dagger},\phi_{ka\bm{x}}\\}$ denotes the fields of the model and
$a=\left\\{1\cdots m\right\\}$ is the replica index. The coupling constants
$g$ and $t$ are drawn from the Gaussian orthogonal ensemble, with probability
distribution functions $\varrho(g)$ and $\varrho(t)$ defined in terms of the
second moments
$\overline{g_{ijk}g_{i^{\prime}j^{\prime}k^{\prime}}}=\frac{\bar{g}^{2}}{2\mathscr{N}^{2}}\delta_{kk^{\prime}}\left(\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}}\right)$
and
$\overline{t_{ij\bm{x}\bm{x}^{\prime}}t_{i^{\prime}j^{\prime}yy^{\prime}}}=\frac{t_{0}^{2}}{2\mathscr{N}}\left(\delta_{ii^{\prime}}\delta_{jj^{\prime}}\delta_{\bm{x}\bm{y}}\delta_{\bm{x}^{\prime}\bm{y}^{\prime}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}}\delta_{\bm{x}\bm{y}^{\prime}}\delta_{\bm{x}^{\prime}\bm{y}}\right)$.
The imaginary-time action is given by

$\displaystyle\mathscr{S}_{0}$ $\displaystyle=$
$\displaystyle\sum_{i\sigma\bm{x}}\int d\tau\,\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)(\partial_{\tau}-\mu)\hat{c}_{i\sigma
a\bm{x}}(\tau)+\frac{1}{2}\sum_{ia\bm{x}}\int
d\tau\phi_{ia\bm{x}}(\tau)(-\partial_{\tau}^{2}+m_{0})\phi_{ia\bm{x}}(\tau),$
(76a) $\displaystyle\mathscr{S}_{g}$ $\displaystyle=$
$\displaystyle\sum_{\sigma\bm{x};ijk}\int
d\tau\,\left[g_{ij,k}+g_{ji,k}^{\star}\right]\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\hat{c}_{j\sigma a\bm{x}}(\tau)\phi_{ka\bm{x}}(\tau),$
(76b) $\displaystyle\mathscr{S}_{t}$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}\langle\bm{x},\bm{x}^{\prime}\rangle\\\
ij\sigma\end{subarray}}\int
d\tau\,t_{ij,\bm{x}\bm{x}^{\prime}}\hat{c}_{j\sigma\bm{x}^{\prime}a}^{\dagger}(\tau)\hat{c}_{i\sigma\bm{x}a}(\tau),$
(76c)

where angled brackets denote the sum over nearest-neighbor sites. Performing
the Gaussian integration over the coupling constants $g$ and $t$ – see Sec.
6.1 and Apps. C1-C3 of Ref. 226 – yields the disorder-averaged effective
action

$\displaystyle\mathscr{S}^{{\rm eff}}$ $\displaystyle=$
$\displaystyle\sum_{i\sigma a\bm{x}}\int d\tau\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)(\partial_{\tau}-\mu)\hat{c}_{i\sigma
a\bm{x}}(\tau)+\frac{1}{2}\sum_{ia\bm{x}}\int
d\tau\phi_{ia\bm{x}}(\tau)(-\partial_{\tau}^{2}+m_{0})\phi_{ia\bm{x}}(\tau)$
(77) $\displaystyle-$
$\displaystyle\frac{\bar{g}^{2}}{2\mathscr{N}}\sum_{ijk}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\left[\hat{c}_{j\sigma
a\bm{x}}^{\dagger}(\tau)\hat{c}_{j\sigma^{\prime}b\bm{x}}^{\dagger}(\tau^{\prime})\,\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\hat{c}_{i\sigma
a\bm{x}}(\tau)-\hat{c}_{j\sigma
a\bm{x}}^{\dagger}(\tau)\hat{c}_{j\sigma^{\prime}b\bm{x}}(\tau^{\prime})\,\hat{c}_{i\sigma^{\prime}b\bm{x}}^{\dagger}(\tau^{\prime})\hat{c}_{i\sigma
a\bm{x}}(\tau)\right]\phi_{ka\bm{x}}(\tau)\phi_{kb\bm{x}}(\tau^{\prime})$
$\displaystyle-$
$\displaystyle\frac{t_{0}^{2}}{4\mathscr{N}}\sum_{ij}\sum_{\begin{subarray}{c}\langle\bm{x},\bm{x}^{\prime}\rangle\\\
ab\sigma\sigma^{\prime}\end{subarray}}\int_{\tau\tau^{\prime}}\,\,\left[\hat{c}_{j\sigma\bm{x}^{\prime}a}^{\dagger}(\tau)\hat{c}_{i\sigma\bm{x}a}(\tau)\hat{c}_{j\sigma^{\prime}\bm{x}^{\prime}b}^{\dagger}(\tau^{\prime})\hat{c}_{i\sigma^{\prime}\bm{x}b}(\tau^{\prime})+\hat{c}_{j\sigma\bm{x}^{\prime}a}^{\dagger}(\tau)\hat{c}_{i\sigma^{\prime}\bm{x}b}(\tau^{\prime})\hat{c}_{i\sigma\bm{x}a}^{\dagger}(\tau)\hat{c}_{j\sigma^{\prime}\bm{x}^{\prime}b}(\tau^{\prime})\right],$

where $\int_{\tau\tau^{\prime}}\equiv\int d\tau\int d\tau^{\prime}$. We now
introduce the following collective bi-local fields

$\displaystyle
G_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=\frac{1}{\mathscr{N}}\sum_{i}\hat{c}_{i\sigma^{\prime}b\bm{x}}^{\dagger}(\tau^{\prime})\hat{c}_{i\sigma
a\bm{x}}(\tau),\quad
F_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=\frac{1}{\mathscr{N}}\sum_{i}\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\hat{c}_{i\sigma
a\bm{x}}(\tau),$ (78a) $\displaystyle
D_{\bm{x}}^{ab}(\tau,\tau^{\prime})=\frac{1}{\mathscr{N}}\sum_{k}\phi_{kb\bm{x}}(\tau^{\prime})\phi_{ka\bm{x}}(\tau).$
(78b)

In order to insert such variables into the action we make use of the following
identity

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\mathcal{D}G\prod_{\begin{subarray}{c}ab\tau\tau^{\prime}\\\
\sigma\sigma^{\prime}\bm{x}\end{subarray}}\delta\left[\mathscr{N}G_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime},\tau)-\sum_{i}\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\right]=\int\mathcal{D}G\mathcal{D}\Sigma
e^{\left[\mathscr{N}G_{\sigma\sigma^{\prime},\bm{x}}^{ba}(\tau^{\prime},\tau)-\sum_{i}\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\right]\,\Sigma_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})},$

where in the second equality we have introduced the Lagrange-multiplier
$\Sigma$ of the field $G$ and the sum over repetated indices is kept implicit.
Alternatively, one can introduce such fields through a Hubbard-Stratonovich
transformation – see Refs. 169, 226. Similarly, for the other fields we use:

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\mathcal{D}F^{\dagger}\mathcal{D}\Phi
e^{\frac{1}{2}\left[\mathscr{N}F_{\sigma^{\prime}\sigma,\bm{x}}^{\dagger
ba}(\tau^{\prime},\tau)-\sum_{i}\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\,c^{\dagger}{}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\right]\,\Phi_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})},$
(80a) $\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\mathcal{D}F\mathcal{D}\Phi^{\dagger}e^{\frac{1}{2}\left[\mathscr{N}F_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime},\tau)-\sum_{i}\hat{c}_{i\sigma
a\bm{x}}(\tau)\,\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})\right]\,\Phi_{\sigma\sigma^{\prime},\bm{x}}^{\dagger
ab}(\tau,\tau^{\prime})},$ (80b) $\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\mathcal{D}D\mathcal{D}\Pi
e^{\frac{1}{2}\left[\mathscr{N}D_{\bm{x}}^{ba}(\tau^{\prime},\tau)-\sum_{i}\phi_{ia\bm{x}}(\tau)\,\phi_{ib\bm{x}}(\tau^{\prime})\right]\,\Pi_{\bm{x}}^{ab}(\tau,\tau^{\prime})}.$
(80c)

The partition function Eq.(75) then becomes
$\overline{\mathscr{Z}^{m}}=\int\mathcal{D}G\mathcal{D}\Sigma\mathcal{D}F\mathcal{D}\Phi\mathcal{D}F^{\dagger}\mathcal{D}\Phi^{\dagger}\mathcal{D}c\mathcal{D}c^{\dagger}\mathcal{D}\phi\,e^{-\mathscr{S}^{{\rm
eff}}}$ with action

$\displaystyle\mathscr{S}^{{\rm eff}}$ $\displaystyle=$
$\displaystyle\sum_{iab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\,\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\,[(\partial_{\tau}-\mu)\delta_{ab}\delta_{\sigma\sigma^{\prime}}\delta(\tau-\tau^{\prime})+\Sigma_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})]\,\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})$
(81) $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{iab\bm{x}}\int_{\tau\tau^{\prime}}\,\phi_{ia\bm{x}}(\tau)[(-\partial_{\tau}^{2}+m_{0})\delta_{ab}\delta(\tau-\tau^{\prime})-\Pi^{ab}(\tau,\tau^{\prime})]\phi_{ib\bm{x}}(\tau^{\prime})$
$\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{iab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}[\hat{c}_{i\sigma
a\bm{x}}^{\dagger}(\tau)\,\Phi_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})\hat{c}_{i\sigma^{\prime}b\bm{x}}^{\dagger}(\tau^{\prime})+\hat{c}_{i\sigma
a\bm{x}}(\tau)\,\Phi_{\sigma\sigma^{\prime},\bm{x}}^{\dagger
ab}(\tau,\tau^{\prime})\,\hat{c}_{i\sigma^{\prime}b\bm{x}}(\tau^{\prime})]$
$\displaystyle-$
$\displaystyle\mathscr{N}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}G_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime},\tau)\Sigma_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})+\frac{\mathscr{N}}{2}\int_{\tau\tau^{\prime}}D_{\bm{x}}^{ba}(\tau^{\prime},\tau)\Pi_{\bm{x}}^{ab}(\tau,\tau^{\prime})$
$\displaystyle-$
$\displaystyle\frac{\mathscr{N}}{2}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\left[F_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime}\tau)\Phi_{\sigma\sigma^{\prime},\bm{x}}^{\dagger
ab}(\tau,\tau^{\prime})+F_{\sigma^{\prime}\sigma,\bm{x}}^{\dagger
ba}(\tau^{\prime}\tau)\Phi_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})\right]$
$\displaystyle+$
$\displaystyle\frac{\bar{g}^{2}}{2}\mathscr{N}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\left[G_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau,\tau^{\prime})G_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)-F_{\sigma^{\prime}\sigma,\bm{x}}^{\dagger
ba}(\tau,\tau^{\prime})F_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)\right]D_{\bm{x}}^{ab}(\tau,\tau^{\prime})$
$\displaystyle+$
$\displaystyle\frac{t_{0}^{2}}{4}\mathscr{N}\sum_{\begin{subarray}{c}\langle\bm{x},\bm{x}^{\prime}\rangle\\\
ab\sigma\sigma^{\prime}\bm{x}\end{subarray}}\int_{\tau\tau^{\prime}}\left[G_{\sigma^{\prime}\sigma,\bm{x}^{\prime}}^{ba}(\tau,\tau^{\prime})G_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)-F_{\sigma^{\prime}\sigma,\bm{x}^{\prime}}^{\dagger
ba}(\tau,\tau^{\prime})F_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)\right].$

The fermionic part can be reorganized in the Nambu representation as:

$\mathscr{S}_{{\rm
ferm}}=-\frac{1}{2}\sum_{iab}\int_{\tau\tau^{\prime}}\,\,\bar{\psi}_{ia}^{\dagger}(\tau)\left[\bar{G}_{0,ab}^{-1}(\tau,\tau^{\prime})-\bar{\Sigma}_{ab}(\tau,\tau^{\prime})\right]\bar{\psi}_{ib}(\tau^{\prime}),$
(82)

with Nambu-spinor $\bar{\psi}_{ia}(\tau)=(\hat{c}_{i\uparrow
a}(\tau),\hat{c}_{i\downarrow a}(\tau),\hat{c}_{i\uparrow
a}^{\dagger}(\tau),\hat{c}_{i\downarrow a}^{\dagger}(\tau))^{T}.$ Bars denote
$4\times 4$ matrices in Nambu space:

$\displaystyle\bar{G}_{0}^{-1}(\tau,\tau^{\prime})=\begin{pmatrix}\hat{G}_{0}^{-1}(\tau,\tau^{\prime})&0\\\
0&-\hat{\tilde{G}}_{0}^{-1}(\tau^{\prime},\tau)\end{pmatrix},\quad\bar{\Sigma}(\tau,\tau^{\prime})=\begin{pmatrix}{\hat{\Sigma}}(\tau,\tau^{\prime})&\hat{\Phi}(\tau,\tau^{\prime})\\\
\hat{\Phi}^{\dagger}(\tau,\tau^{\prime})&-\hat{\Sigma}(\tau^{\prime},\tau)\end{pmatrix}.$
(83)

Hats denote $2\times 2$ matrices in the spin subspace. The bare propagator is
given by
$\hat{G}_{0}^{-1}(\tau,\tau^{\prime})=-(\partial_{\tau}-\mu)\delta(\tau-\tau^{\prime})\,\hat{\sigma}_{0},$
$\hat{\tilde{G}}_{0}^{-1}(\tau,\tau^{\prime})=(\partial_{\tau}+\mu)\delta(\tau-\tau^{\prime})\,\hat{\sigma}_{0}$,
with $\hat{\sigma}_{0}$ the identity matrix in the 2 dimensional spin
subspace. We now perform the Grassmann integral over the Nambu-fermions
$\bar{\psi}$ and the Gaussian one over the bosons $\phi$, and obtain:

$\displaystyle\frac{\mathscr{S}^{{\rm eff}}}{\mathscr{N}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\rm
Tr}\log\left(\bar{G}_{0}^{-1}-\bar{\Sigma}\right)+\frac{1}{2}{\rm
Tr}\log\left(D_{0}^{-1}-\Pi\right)$ (84) $\displaystyle-$
$\displaystyle\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}G_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime},\tau)\Sigma_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})+\frac{1}{2}\int_{\tau\tau^{\prime}}D_{\bm{x}}^{ba}(\tau^{\prime},\tau)\Pi_{\bm{x}}^{ab}(\tau,\tau^{\prime})$
$\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\left[F_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau^{\prime}\tau)\Phi_{\sigma\sigma^{\prime},\bm{x}}^{\dagger
ab}(\tau,\tau^{\prime})+F_{\sigma^{\prime}\sigma,\bm{x}}^{\dagger
ba}(\tau^{\prime}\tau)\Phi_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})\right]$
$\displaystyle+$
$\displaystyle\frac{\bar{g}^{2}}{2}\sum_{ab\sigma\sigma^{\prime}\bm{x}}\int_{\tau\tau^{\prime}}\left[G_{\sigma^{\prime}\sigma,\bm{x}}^{ba}(\tau,\tau^{\prime})G_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)-F_{\sigma^{\prime}\sigma,\bm{x}}^{\dagger
ba}(\tau,\tau^{\prime})F_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)\right]D_{\bm{x}}^{ab}(\tau,\tau^{\prime})$
$\displaystyle+$
$\displaystyle\frac{t_{0}^{2}}{4}\sum_{\begin{subarray}{c}\langle\bm{x},\bm{x}^{\prime}\rangle\\\
ab\sigma\sigma^{\prime}\bm{x}\end{subarray}}\int_{\tau\tau^{\prime}}\left[G_{\sigma^{\prime}\sigma,\bm{x}^{\prime}}^{ba}(\tau,\tau^{\prime})G_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)-F_{\sigma^{\prime}\sigma,\bm{x}^{\prime}}^{\dagger
ba}(\tau,\tau^{\prime})F_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau^{\prime},\tau)\right].$

In addition, we work within the spin-singlet and replica-diagonal ansatz

$\displaystyle{G}_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=G_{\bm{x}}(\tau,\tau^{\prime})\hat{\sigma}_{0}\delta_{ab},\,\,\,\,$
$\displaystyle\quad{\Sigma}_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=\Sigma_{\bm{x}}(\tau,\tau^{\prime})\hat{\sigma}_{0}\delta_{ab},$
$\displaystyle{F}_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=F_{\bm{x}}(\tau,\tau^{\prime})\,i\hat{\sigma}_{2}\delta_{ab},\,\,\,$
$\displaystyle\quad{\Phi}_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=\Phi_{\bm{x}}(\tau,\tau^{\prime})\,i\hat{\sigma}_{2}\delta_{ab},$
$\displaystyle{F}_{\sigma\sigma^{\prime},\bm{x}}^{\dagger
ab}(\tau,\tau^{\prime})=-F_{\bm{x}}^{\dagger}(\tau,\tau^{\prime})\,i\hat{\sigma}_{2}\delta_{ab},$
$\displaystyle\quad{\Phi^{\dagger}}_{\sigma\sigma^{\prime},\bm{x}}^{ab}(\tau,\tau^{\prime})=-\Phi_{\bm{x}}^{\dagger}(\tau,\tau^{\prime})\,i\hat{\sigma}_{2}\delta_{ab},$
(85)

with Pauli matrices $\hat{\sigma}_{i}$. Noting that ${\rm Tr}{\rm
log}(\bar{G}_{0}^{-1}-\bar{\Sigma})=2{\rm Tr}{\rm
log}(\hat{G}_{0}^{-1}-\hat{\Sigma})$, we finally get Eq. (2) of the main text.

##  Appendix B Numerical methods

This Appendix sketches the numerical algorithms employed throughout this work
to exactly solve the large-$\mathscr{N}$ saddle-point Eliashberg equations (4)
on the imaginary axis, in the normal and superconducting states. The
propagators for fermions and bosons $G(i\omega_{n})$, $F(i\omega_{n})$, and
$D(i\Omega_{n})$, and their associtated self-energies $\Sigma(i\omega_{n})$,
$\Phi(i\omega_{n})$, and $\Pi(i\Omega_{n})$, are written in terms of fermionic
and bosonic Matsubara frequencies, $\omega_{n}$ and $\Omega_{n}$ respectively.
These quantities can be transformed to imaginary time
$\tau\in\left[0,\beta\right]$ with $\beta=(k_{B}T)^{-1}$, which is especially
useful for convolutions like the ones contained in Eqs. (4a) or (4c), since in
the $\tau$-representation such convolutions become simple products between the
convolved functions.

###  B.1 Matsubara Fourier transforms on the imaginary axis

Numerically, we discretize the imaginary-time interval in equal steps
according to $\tau_{l}=l/(2N_{f}k_{B}T)$ with $l\in\left[0,2N_{f}-1\right]$,
such that the discretized Matsubara frequencies are $\omega_{n}=(2n+1)\pi
k_{B}T$ and $\Omega_{n}=2n\pi k_{B}T$, with $n\in\left[-N_{f},N_{f}-1\right]$.
$N_{f}$ is the high-frequency cutoff, that must be chosen to lie well into the
ultraviolet regime of the theory where the propagators have already decayed
from their low-frequency behavior. Then, the propagators and self-energies
become finite discrete lists of values. After algebraic manipulations
(circshifts) of these lists, the FFT protocol lends itself to the built-in
implementation provided by the optimized `Fourier[]` and `InverseFourier[]`
functions of Mathematica. Such implementation is similar to many other self-
consistent loops used to solve Shwinger-Dyson saddle-point equations of SYK-
like models [218, 102, 219, 239, 132, 220, 88, 81, 94, 221, 134].

###  B.2 Self-consistent loops for the saddle-point equations on the imaginary
axis

At the first iteration $j=0$, our self-consistent loop to solve Eqs. (4)
starts with a guess on $\Sigma(i\omega_{n})$, which for $g^{2}<k_{B}T<1/g^{2}$
is assumed to be given by the NFL-SYK low-energy solution (11), while for
$k_{B}T>1/g^{2}$ we adopt the impurity-like solution (IV.3). We find these
guesses to be sufficiently accurate that the loop converges both in the normal
and superconducting states, although one could refine the algorithm to assume
more accurate low-energy solutions in the superconducting state, such as Eq.
(144) for the fermionic propagator in the SYK2-FL and impurity-like regimes.
At the end of the current iteration $j>0$, the list of the normal self-energy
$\Sigma_{j}(i\omega_{n})$ is updated with a weighted sum of the solution
$\bar{\Sigma}_{j}(i\omega_{n})$ of Eq. (4a) and the solution
$\Sigma_{j-1}(i\omega_{n})$ at the previous iteration $j-1$, according to

$\Sigma_{j}(i\omega_{n})=\alpha_{s}\bar{\Sigma}_{j}(i\omega_{n})+(1-\alpha_{s})\Sigma_{j-1}(i\omega_{n}).$
(86)

The mixing factor $\alpha_{s}\in\left(0,1\right)$ helps to stabilize the
convergence of the self-consistent loop [134]. We heuristically choose
$\alpha_{s}\in\left(0.004,0.1\right)$ depending on the values of parameters
like $g^{2}/(zt_{0}^{2}\omega_{0})$. In practice, we find that smaller values
of $\alpha_{s}$ are required in the superconducting state, and especially at
strong coupling (impurity-like regime). The error between the current
iteration $j>0$ and the previous one $j-1$ is monitored by the sum over the
normal self-energy
$\epsilon_{\Sigma}=\sum_{i\omega_{n}}\left|\Sigma_{j}(i\omega_{n})-\Sigma_{j-1}(i\omega_{n})\right|$.
Convergence is reached when $\epsilon_{\Sigma}$ falls below a user-imposed
threshold. In the superconducting state, the gap function
$\Delta_{j}(i\omega_{n})$ is also weighted and summed to its value at the
previous iteration:

$\Delta_{j}(i\omega_{n})=\alpha_{d}\Delta_{j}(i\omega_{n})+(1-\alpha_{d})\Delta_{j-1}(i\omega_{n}),$
(87)

where $\alpha_{d}\in\left(0.004,0.5\right)$ is another heuristically chosen
weighting factor, which may be different from $\alpha_{s}$. We empirically
find that the lowest values of $\alpha_{d}$ are required in the strong
coupling, impurity-like regime.

###  B.3 Self-consistent calculation of the critical temperature

In order to find the superconducting critical temperature $T_{c}$, we solve
the linearized gap equation (34). Here we define $T_{c}$ at the mean-field,
_gap-closing_ temperature at which $\Delta(i\omega_{n})\rightarrow
0\,\forall\omega_{n}$. As described in Sec. XI, the superconducting transition
temperature is further decreased by phase fluctuations at finite
$\mathscr{N}$. Finding $T_{c}$ requires us to know the saddle-point normal-
state converged lists for $G(i\omega_{n})$, $\Sigma(i\omega_{n})$,
$D(i\Omega_{n})$, and $\Pi(i\Omega_{n})$, that we obtain through the self-
consistent loop previously sketched in Sec. B.2. Indeed, after linearization
of the Eliashberg equations, $G(i\omega_{n})$ does not depend on
$\Phi(i\omega_{n})$, so it coincides with the normal-state solution. Then, the
gap equation (34) has the structure of an infinite-dimensional eigenvalue
problem for the anomalous self-energy $\tilde{\Phi}(i\omega_{n})$, which can
be written schematically as
$\tilde{\Phi}(i\omega_{n})=\sum_{m}A_{nm}\tilde{\Phi}(i\omega_{m})$, where
$\underline{\underline{A}}=\left\\{A_{nm}\right\\}$ is the kernel matrix which
depends on $D(i\Omega_{n})$, $Z(i\omega_{n})$, and temperature $T$. The
dimensionality of $\underline{\underline{A}}$ is truncated to $2N_{f}\times
2N_{f}$ due to the high-frequency cutoff of the Matsubara propagators. At
$T<T_{c}$, the largest eigenvalue of $\underline{\underline{A}}$ is greater
than one, and consequently the self-consistent loop tends to increase the
value of $\tilde{\Phi}(i\omega_{m})$. Conversely, at $T>T_{c}$ the self-
consistency drives $\tilde{\Phi}(i\omega_{m})$ to zero and all eigenvalues of
$\underline{\underline{A}}$ are therefore smaller than one. At $T=T_{c}$, the
largest eigenvalue of $\underline{\underline{A}}$ is unitary [161, 223, 134].
Monitoring the $T$-dependent eigenvalues of $\underline{\underline{A}}$ is
implemented by a root-finding loop on temperature, which allows us to find
$T_{c}$. The accuracy on the final $T_{c}$ value is set to $dT=10^{-7}$.

Figure 18:  (a) Numerical solution for the fermionic propagator on the
imaginary axis, showing the crossover from the SYK-NFL phase to the SYK2-FL
regime, at temperature $T=0.002\omega_{0}/k_{B}$, coupling $g=1$, coordination
number $z=1$, and for different hoppings $t_{0}$. The dashed and dotted gray
lines are given by the low-energy expansions of Eqs. (11) and (IV.2),
respectively. (b) Numerical solution for the bosonic propagator, for the same
parameters as in panel (a). The dashed gray line stems from Eq. (12).

###  B.4 Estimation of the quasiparticle residue

To extract the quasiparticle residue $Z_{\rm qp}$, we use its definition (50)
on the real axis and its relation to the imaginary-axis quantity
$Z(i\omega_{n})$. From the Kramers-Kronig relation we have

$\mathrm{Re}\Sigma^{R}(\omega)=-\frac{1}{\pi}\mathscr{P}\int_{-\infty}^{+\infty}d\omega^{\prime}\frac{\mathrm{Im}\Sigma^{R}(\omega^{\prime})}{\omega-\omega^{\prime}},$
(88)

which in the $\omega\rightarrow 0$ limit yields ($\mathscr{P}$ stands for the
principal part of the integral)

$\displaystyle\left.\frac{\partial\mathrm{Re}\Sigma^{R}(\omega)}{\partial\omega}\right|_{\omega\rightarrow
0}=\mathscr{P}\int_{-\infty}^{+\infty}\frac{d\omega^{\prime}}{\pi}\frac{\mathrm{Im}\Sigma^{R}(\omega^{\prime})}{(\omega^{\prime})^{2}}.$
(89)

On the other hand, from the spectral representation

$\Sigma(i\omega_{n})=\int_{-\infty}^{+\infty}\frac{d\epsilon}{\pi}\frac{\mathrm{Im}\Sigma^{R}(\epsilon)}{\epsilon-i\omega_{n}}.$
(90)

In the $T\rightarrow 0$ limit, and taking the imaginary part of Eq. (90) we
have

$\mathrm{Im}\Sigma(i\omega)=\int_{-\infty}^{+\infty}\frac{d\epsilon}{\pi}\frac{\mathrm{Im}\Sigma^{R}(\epsilon)\omega}{\epsilon^{2}+\omega^{2}}.$
(91)

Combining equations (91) and (89), we have

$\frac{1}{Z_{\rm
qp}}=1-\left.\frac{\mathrm{Im}\Sigma(i\omega)}{\omega}\right|_{\omega\rightarrow
0}.$ (92)

Using the relation (6) valid for our model on the imaginary axis, we finally
achieve Eq. (51). Therefore, the real-axis zero-energy quasiparticle weight
can be estimated through Eq. (51) by the inverse of the imaginary axis
quantity $Z(i\omega)$, in the zero-temperature limit. In practice, we use the
imaginary-axis code described in Sec. B.2 at temperatures $T\ll T_{c}$ to
approximate Eq. (51). As a consistency check, we also evaluate $Z_{\rm qp}$
for selected couplings and hoppings through the definition (50) using the
real-axis code described in Sec. B.5, which allows us to go to the zero-energy
limit but still requires finite temperatures.

Figure 19:  (a) Numerical solution for the fermionic propagator on the
imaginary axis, showing the crossover from the impurity-like phase to the
SYK2-FL regime, at temperature $T=0.05\omega_{0}/k_{B}$, coupling $g=4$,
coordination number $z=1$, and for different hoppings $t_{0}$. The dashed
curve and dotted gray line are given by Eq. (IV.3) and by the low-energy
expansion of Eq. (IV.2), respectively. (b) Numerical solution for the bosonic
propagator, for the same parameters as in panel (a). The dashed gray line is
generated by Eqs. (21) and (22).

###  B.5 Saddle-point equations on the real axis

Our numerical solution of the saddle-point equations (4) on the real axis is
based on the spectral (Lehmann) representation of the fermionic and bosonic
propagators [161], which can be implemented through Laplace transforms [240].
In essence, if we have to perform a convolution on the imaginary axis of the
form

$\alpha(i\Omega_{n})=-k_{B}T\sum_{m=-\infty}^{+\infty}g_{1}(i\omega_{m}+i\Omega_{n})g_{2}(i\omega_{m}),$
(93)

on the real axis it follows that

$\displaystyle\alpha(\omega)$
$\displaystyle=\int_{-\infty}^{+\infty}\frac{d\epsilon}{\pi}f_{FD}(\epsilon)\left[g_{1}(\omega+\epsilon)\mathrm{Im}\left\\{g_{2}(\epsilon)\right\\}\right.$
$\displaystyle\left.+\mathrm{Im}\left\\{g_{1}(\epsilon)\right\\}g_{2}^{\ast}(\epsilon-\omega)\right].$
(94)

Defining the Laplace transform $\alpha(t)$ as

$\alpha(\omega)=\int_{0}^{+\infty}dt\alpha(t)e^{i\omega t},$ (95)

Eq. (B.5) translates as

$\alpha(t)=i(2\pi)^{2}\left[\rho_{1}(t)a_{2}^{\ast}(t)-\rho_{2}^{\ast}(t)a_{1}(t)\right],$
(96)

where

$\rho_{i}(t)=-\int_{-\infty}^{+\infty}\frac{d\epsilon}{2\pi}\frac{\mathrm{Im}\left\\{g_{i}(\epsilon)\right\\}}{\pi}e^{-i\epsilon
t}$ (97)

and

$a_{i}(t)=-\int_{-\infty}^{+\infty}\frac{d\epsilon}{2\pi}f_{FD}(\epsilon)\frac{\mathrm{Im}\left\\{g_{i}(\epsilon)\right\\}}{\pi}e^{-i\epsilon
t}.$ (98)

Employing Eq. (96), one can transform convolutions of the form (B.5) on the
real axis into simpler products between the functions (97) and (98). Then, an
inverse Laplace transform of the resulting $\alpha(t)$ yields $\alpha(\omega)$
as a function of real frequency $\omega$. We apply the above approach to the
fermionic self-energy (4a), the anomalous self-energy (4d), and the boson
self-energy (4c), written on the real axis by the means of the spectral
representation (B.5). The respective Laplace transforms $\Sigma^{R}(t)$,
$\Phi^{R}(t)$, and $\Pi^{R}(t)$, where the superscript $R$ denotes the
retarded quantities, are complemented by the Dyson equations

$G^{R}(\omega)=\frac{\omega+i\gamma-\Sigma^{R}(\omega)}{\left[\omega+i\gamma-\Sigma^{R}(\omega)\right]^{2}-\left[\Phi^{R}(\omega)\right]^{2}},$
(99)
$F^{R}(\omega)=\frac{\Phi^{R}(\omega)}{\left[\omega+i\gamma-\Sigma^{R}(\omega)\right]^{2}-\left[\Phi^{R}(\omega)\right]^{2}},$
(100)

and

$D^{R}(\omega)=\frac{1}{\omega_{0}^{2}-(\omega+i\gamma)^{2}-\Pi^{R}(\omega)},$
(101)

where $\gamma=i0^{+}$, and iterated self-consistently starting from either a
BCS-like solution of the form [241]

$G^{R}(\omega)=\frac{\omega+i\gamma}{(\omega+i\gamma)^{2}-\Delta^{2}},$ (102a)
$F^{R}(\omega)=\frac{\Delta}{(\omega+i\gamma)^{2}-\Delta^{2}},$ (102b)

or from the results of a previously converged self-consistent loop. In
practice, we employ a finite damping $\gamma=10^{-5}\omega_{0}$. The fermionic
self-energy $\Sigma(\omega)$ is weighted at each self-consistent iteration
with the result of the previous iteration, similarly to Eq. (86), with a
mixing factor $\alpha_{s}\in\left(0,1\right)$. The error between the current
iteration $j>0$ and the previous one $j-1$ is monitored by the sum over the
normal self-energy
$\epsilon_{\Sigma}=\sum_{\omega}\left|\Sigma^{R}_{j}(\omega)-\Sigma_{j-1}^{R}(\omega)\right|$:
when $\epsilon_{\Sigma}$ decreases below a user-imposed value, convergence is
reached.

##  Appendix C Derivation of the normal-state results

###  C.1 Fermion and boson propagators on the imaginary axis

In this section, for completeness we report the graphs of the normal-state
fermion $G(i\omega_{n})$ and boson $D(i\Omega_{n})$ propagators on the
imaginary axis. These propagators stem from the self-consistent solution of
the saddle-point equations (8), obtained with the methods described in
Appendix B.2. Notice that, on the imaginary axis, $G(i\omega_{n})$ and
$D(i\Omega_{n})$ are purely imaginary and purely real, respectively, as
deduced from Eqs. (8) and (6).

Figure 20:  Normal-state squared renormalized boson frequency
$\omega_{r}^{2}/\omega_{0}^{2}$ as a function of normalized temperature
$k_{B}T/(\omega_{0}g^{2})$ and fermion hopping parameter $t_{0}$, calculated
from the full numerical solution of the Yukawa-SYK Eliashberg equations in
accordance with Eq. (5). The dashed gray line is the power-law solution
$\omega_{r}^{2}/\omega_{0}^{2}\sim\left[(k_{B}T)/(\omega_{0}g^{2})\right]^{4\Delta-1}$
found in the SYK-NFL regime of the single Yukawa-SYK dot [105]. Panels (a),
(b), and (c) show the numerical results for $g=\left\\{0.25,0.5,1\right\\}$,
respectively, for different hopping parameters
$t_{0}/\omega_{0}=\left\\{0,0.5,1\right\\}$. Arrows indicate the analytical
estimation (18) of the zero-temperature value for
$\omega_{r}^{2}/\omega_{0}^{2}$ in the SYK2-FL regime, for the same parameters
as the corresponding numerical data points of the same color.

Fig. 18(a) shows the results for
$-g^{2}\mathrm{Im}\left\\{G(i\omega_{n})\right\\}$ across the SYK-NFL/SYK2-FL
crossover, at temperature $T=0.002\omega_{0}/k_{B}$, coupling $g=1$,
coordination number $z=1$, and for different hoppings $t_{0}$. At the lowest
values of $t_{0}$, the fermionic propagator displays the power-law dependence
(11) at small energies, which is characteristic of the SYK-NFL regime, and is
indicated by the dashed gray line. The crossover to the disordered FL regime
is realized as higher hoppings, with visible differences in the propagator
appearing for $t_{0}\approx\omega_{0}$, consistently with the crossover energy
$\omega_{c}$ estimated in Sec. V and with the spectral functions on the real
axis in Fig. 4(c). In the SYK2-FL regime, the low-energy part of the
propagator tends to the constant $\sqrt{2}/(\sqrt{z}t_{0})$ (the additional
factor of $g^{2}$ is due to the normalization of the $y$ axis in the figure),
as deduced from Eq. (IV.2) at small $\omega_{n}$; this behavior is shown by
the dotted gray line. Fig. 18(b) reports the results for $D(i\Omega_{n})$, for
the same parameters as in panel (a). In the SYK-NFL regime, i.e., at low
hopping, we recognize the power-law scaling of Eq. (12), which is highlighted
by the dashed gray line. Such scaling evolves into almost free bosons in the
high-hopping SYK2-FL regime, since the boson self-energy becomes negligible
when hopping dominates over fermion-boson coupling; this stiffening of the
renormalized boson frequency towards $\omega_{r}\approx\omega_{0}$ is
consistent with the analogous trend of the imaginary part of the real-axis
boson propagator, shown in Fig. 4(d).

Figure 21:  Normal-state squared renormalized boson frequency
$\omega_{r}^{2}/\omega_{0}^{2}$ as a function of normalized temperature
$k_{B}T/(\omega_{0}g^{2})$, for coupling $g=4$ and for different hoppings
$t_{0}$. Data points are calculated from the full numerical solution of the
Yukawa-SYK Eliashberg equations Eq. (5). The dashed gray line is the power-law
solution (22) found in the impurity-like regime of the single Yukawa-SYK dot
[105].

The results for $-g^{2}\mathrm{Im}\left\\{G(i\omega_{n})\right\\}$ across the
impurity-like/SYK2-FL crossover are shown in Fig. 19(a), at temperature
$T=0.05\omega_{0}/k_{B}$, coupling $g=4$, coordination number $z=1$, and for
different hoppings. At low hopping, the curves follow the analytical result
(IV.3), proper of the impurity-like regime, and indicated by the dashed gray
curve. Increasing hopping, we cross over to the SYK2-FL regime, and the low-
energy part of the propagator approaches $\sqrt{2}/(\sqrt{z}t_{0})$, as marked
by the dotted gray line. The crossover of fermions to the disordered FL regime
occurs at $t_{0}\approx 19\omega_{0}$ at $T=0.05k_{B}/\omega_{0}$,
consistently with the criterion (29) in Sec. V, and with the real-axis
spectral functions shown in Fig. 5(c). Fig. 19(b) shows the boson propagator,
which is consistent with the impurity-like solution given by Eqs. (21) and
(22) at low hoppings. At the highest value $t_{0}=30\omega_{0}$, there are
modifications which appear at low frequency, due to the stiffening of the
renormalized boson frequency. These features are most easily seen by looking
at the analogous real-axis calculations in Fig. 5(d).

###  C.2 Normal-state renormalized boson frequency

The renormalized boson frequency is given by Eq. (5), while the boson self-
energy in the normal state follows Eq. (8b). We first present the full
numerical results stemming from Eqs. (8), and then report the derivations of
analytical approximations. Fig. 20 shows $\omega_{r}^{2}/\omega_{0}^{2}$ as a
function of normalized temperature $k_{B}T/(\omega_{0}g^{2})$, for different
hopping parameters and different couplings. In the single-dot limit $t_{0}=0$,
all curves follow the power-law scaling (13) found in SYK-NFL regime.
Increasing hopping at fixed coupling, we cross over to the SYK2-FL regime and
$\omega_{r}^{2}$ stiffens because the static boson self-energy is negligible
for $g^{2}\omega_{0}\ll zt_{0}^{2}$, as seen from Eq. (18). The latter
analytical expression, valid at $T=0$ and derived in App. C.2.1, yields the
values indicated by arrows in Fig. 20. Fig. 21 displays
$\omega_{r}^{2}/\omega_{0}^{2}$ as a function of normalized temperature
$k_{B}T/(\omega_{0}g^{2})$, for $g=4$ and different hoppings. The single-dot
curve $t_{0}=0$ is consistent with Eq. (22), characteristic of the impurity-
like regime: the bosons are critical at zero temperature. As in Fig. 20, the
crossover to the SYK2-FL regime is realized by increasing hopping, which makes
$\omega_{r}^{2}$ stiffen: the bosons are no longer critical in the disordered
FL phase.

We now derive analytical expressions for $\omega_{r}$ at $T=0$ in the SYK2-FL
and impurity-like regimes.

####  C.2.1 Fermi-liquid SYK2-FL regime

In the SYK2-FL regime, for $\bar{g}^{2}/\omega_{0}^{2}\ll zt_{0}^{2}$, the
dynamical quasiparticle weight is approximately given by Eq. (26), i.e.,

$Z(i\omega_{n})=\frac{1}{2}\left(1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega_{n}^{2}}}\right)$
(103)

in the normal state, a result which follows from Eq. (IV.2) and the normal-
state Dyson equation (9). Let us first calculate the normal-state polarization
bubble at zero boson frequency and at $T=0$. Using Eqs. (8b) and (103), we
have

$\Pi(0)=2\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\omega_{m}^{2}}\frac{1}{\left[Z(i\omega_{m})\right]^{2}}\\\
\approx
2\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\omega_{m}^{2}}\frac{4}{\left[1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega_{m}^{2}}}\right]^{2}}\\\
=8\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\left[\omega_{m}+\sqrt{\omega_{m}^{2}+2zt_{0}^{2}}\right]^{2}}.$
(104)

At $T=0$, the sums over Matsubara frequencies become integrals, so that

$\displaystyle\Pi(0)$
$\displaystyle=8\bar{g}^{2}\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{1}{\left(\omega+\sqrt{\omega^{2}+2zt_{0}^{2}}\right)^{2}}$
$\displaystyle=\frac{8\sqrt{2}\bar{g}^{2}}{3\pi\sqrt{z}t_{0}}.$ (105)

The renormalized boson frequency thus has an analytical expression at $T=0$,
which stems from Eqs. (5) and (C.2.1): it is Eq. (18). Of course Eq. (18) is
valid only as long as $\omega_{r}(0)>0$, which means at sufficiently low
interaction $\bar{g}^{2}/\omega_{0}^{2}\ll 3\pi\sqrt{z}t_{0}/(8\sqrt{2})$. The
estimation (18) favorably compares with the full numerical solution of the
Eliashberg equations in the SYK2-FL regime, as shown by the arrows in Fig. 20.

####  C.2.2 Strong-coupling impurity-like regime

For the dynamical weight $Z(i\omega_{n})$ in the impurity-like regime we can
employ the parametrization (140) that interpolates between the SYK2-FL and
impurity-like regions of the phase diagram:

$Z(i\omega_{n})\approx\frac{1}{2}\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{(\omega_{n})^{2}}}\right],$
(106)

where $\tilde{\Omega}$ follows Eq. (54): it is a characteristic frequency
which interpolates between the SYK2-FL and impurity-like regimes.

Using Eqs. (8b) and (106), we achieve

$\Pi(0)=2\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\omega_{m}^{2}}\frac{1}{\left[Z(i\omega_{m})\right]^{2}}\\\
\approx
2\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\omega_{m}^{2}}\frac{4}{\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega_{m}^{2}}}\right]^{2}}\\\
=8\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{\omega_{m}^{2}}\frac{1}{\left[1+\sqrt{1+\tilde{\Omega}^{2}/\omega_{m}^{2}}\right]^{2}}.$
(107)

At $T=0$, the sum over Matsubara frequencies becomes an integral, therefore

$\displaystyle\Pi(0)$
$\displaystyle=8\bar{g}^{2}\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{1}{\omega^{2}}\frac{1}{\left(1+\sqrt{1+\tilde{\Omega}^{2}/\omega^{2}}\right)^{2}}$
$\displaystyle=\frac{16\bar{g}^{2}}{3\pi\tilde{\Omega}}.$ (108)

The renormalized boson frequency stems from Eqs. (5) and (C.2.2):

$\displaystyle\left[\omega_{r}(0)\right]^{2}$
$\displaystyle\approx\omega_{0}^{2}-\frac{16\bar{g}^{2}}{3\pi\tilde{\Omega}}$
$\displaystyle=\omega_{0}^{2}-\frac{16\bar{g}^{2}}{3\pi}\left[2zt_{0}^{2}+4\left(\frac{8g^{2}\omega_{0}}{3\pi}\right)^{2}\right]^{-\frac{1}{2}}.$
(109)

In the SYK2-FL limit $\bar{g}^{2}/\omega_{0}^{2}\ll\sqrt{z}t_{0}$, we retrieve
Eq. (18) from the just obtained result (C.2.2). On the other hand, in the
impurity-like regime $\bar{g}^{2}/\omega_{0}^{2}\gg\sqrt{z}t_{0}$, neglecting
hopping altogether yields

$\left[\omega_{r}(0)\right]^{2}\approx\omega_{0}^{2}-\frac{16\bar{g}^{2}}{3\pi}\frac{3\pi\omega_{0}^{2}}{16\bar{g}^{2}}=0,$
(110)

which is consistent with Eq. (22): the bosons are fully massless at zero
temperature in the limit of zero hopping.

##  Appendix D Derivation of the superconducting-state results

In this Appendix, we show that analytical results are reachable for various
quantities in the superconducting state, in both the Fermi-liquid and
impurity-like regimes. The general strategy in the SYK2-FL regime involves
approximating $D(i\Omega_{n})$ with the propagator for free bosons,

$D(i\Omega_{n})\approx\frac{1}{\omega_{0}^{2}},$ (111)

since the boson self-energy is negligible in such regime, where
$g^{2}\omega_{0}\ll zt_{0}^{2}$. Furthermore, in both the SYK2-FL and
impurity-like regimes, a useful approximation for the gap function
$\Delta(i\omega_{n})$ is the piecewise-constant

$\Delta(i\omega_{n})=\begin{cases}\Delta\,:\,\left|\omega_{n}\right|\leq\omega_{0},\\\
0\,:\,\left|\omega_{n}\right|>\omega_{0}.\end{cases}$ (112)

In the impurity-like regime, the rather crude approximation (112) neglects all
strong-coupling retardation effects on the gap $\Delta(i\omega_{n})$, which
manifest themselves in the oscillations observed in the superconducting
spectral functions of the interacting Cooper-pair fluid; see Fig. 12.
Nevertheless, Eq. (112) is sufficient to achieve qualitatively valid
expressions for thermodynamic quantities even in the impurity-like regime,
knowing the zero-energy gap $\Delta_{0}$ as a function of coupling and
hopping. One example of such quantities is the phase stiffness, analyzed in
Sec. XI. Given the above assumptions, we proceed with the analysis of the
superconducting state in the SYK2-FL and impurity-like regimes.

###  D.1 Weak-coupling disordered Fermi liquid: the SYK2-FL regime

####  D.1.1 Dynamical quasiparticle weight

In the SYK2-FL regime, we can assume the bosons to obey Eq. (111), i.e., for a
negligible boson self-energy $\Pi(i\Omega_{n})\approx 0$ and in the static
limit, as verified with a full numerical solution of the Eliashberg equations
for $\bar{g}^{2}/\omega_{0}^{2}\ll zt_{0}^{2}$. Using Eqs. (4a) and (6) in
this regime, where we can neglect the term in $\Sigma(i\omega_{n})$ that
depends on fermion-boson coupling, we find

$Z(i\omega_{n})\approx
1+\frac{zt_{0}^{2}}{2}\frac{1}{Z(i\omega_{n})\left\\{\omega_{m}^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}.$
(113)

We now make the assumption (112) of a constant gap up to $\omega_{0}$, which
is numerically confirmed by the full self-consistent solution of the lattice-
SYK Eliashberg equations. Then, at low energies
$\left|i\omega_{n}\right|\lessapprox\omega_{0}$, we can write

$Z(i\omega_{n})\approx
1+\frac{zt_{0}^{2}}{2}\frac{1}{Z(i\omega_{n})\left(\omega_{m}^{2}+\Delta^{2}\right)}.$
(114)

The solution to the quadratic equation (114) is

$Z(i\omega_{n})=\frac{1}{2}\left(1\pm\sqrt{1+\frac{2zt_{0}^{2}}{\Delta^{2}+\omega_{n}^{2}}}\right),$
(115)

where we have to take the positive root in order for $Z(\omega)$ to have
physically sound asymptotics, as in the normal state. We thus conclude that in
the superconducting state, for $\bar{g}^{2}/\omega_{0}^{2}\ll zt_{0}^{2}$,

$Z(i\omega_{n})=\frac{1}{2}\left(1+\sqrt{1+\frac{2zt_{0}^{2}}{\Delta^{2}+\omega_{n}^{2}}}\right).$
(116)

In the normal state, $\Delta\rightarrow 0$ and we retrieve the quasiparticle
weight of the SYK2-FL model, stemming from Eq. (IV.2): it is given by Eq.
(103).

####  D.1.2 Fermion propagators and self-energies in the superconducting state

Using Eq. (6), we can calculate the fermion self-energy from the previously
obtained dynamical quasiparticle weight (116) in the superconducting state:

$\displaystyle\Sigma(i\omega)$
$\displaystyle=\frac{1}{2}\left[i\omega_{n}-i\mathrm{sign}(\omega_{n})\right.$
$\displaystyle\times\left.\sqrt{\omega_{n}^{2}+\frac{2zt_{0}^{2}\omega_{n}^{2}}{\omega_{n}^{2}+\Delta^{2}}}\right].$
(117)

In the normal state, i.e., for $\Delta\rightarrow 0$, Eq. (D.1.2) correctly
yields the SYK2-FL self-energy stemming from Eqs. (IV.2) and (16). The fermion
propagator results from the normal-state Dyson equation (16), which is also
approximately valid in the SYK2-FL superconducting regime, since fermion-boson
coupling is negligible with respect to hopping in such regime. We thus achieve

$G(i\omega_{n})=\frac{1}{zt_{0}^{2}}\left[i\omega_{n}-i\mathrm{sign}(\omega_{n})\sqrt{\omega_{n}^{2}+\frac{2zt_{0}^{2}\omega_{n}^{2}}{\omega_{n}^{2}+\Delta^{2}}}\right].$
(118)

In the normal state we consequently obtain the SYK2-FL propagator (IV.2).
Performing the analytic continuation $i\omega_{n}\rightarrow\omega+i0^{+}$ on
Eq. (118), we arrive at the (retarded) fermion propagator on the real axis:
Eq. (47). In the superconducting phase, the anomalous propagator
$F(i\omega_{n})$ appears as well, in accordance with Eq. (4f). Combining the
latter with Eq. (7), we obtain

$F(i\omega_{n})=-\frac{\Delta(i\omega_{n})}{Z(i\omega_{n})\left\\{\omega_{n}^{2}+\left[\Delta(i\omega_{n})\right]^{2}\right\\}}.$
(119)

Making the assumption (112) of an energy-independent gap in the SYK2-FL state,
and using the previously determined dynamical weight (116), we arrive at

$F(i\omega_{n})=-\frac{2\Delta}{\left[\omega_{n}^{2}+\Delta^{2}\right]\left[1+\sqrt{1+\frac{2z{t_{0}^{2}}}{\omega_{n}^{2}+\Delta}^{2}}\right]}.$
(120)

The real-axis retarded anomalous propagator stems from setting
$i\omega_{n}\rightarrow\omega+i0^{+}$ in Eq. (120):

$\displaystyle F^{R}(\omega)$
$\displaystyle=-\frac{2\Delta}{\Delta^{2}-(\omega+i0^{+})^{2}}$
$\displaystyle\times\frac{1}{1+\sqrt{1+\frac{2zt_{0}^{2}}{\Delta^{2}-(\omega+i0^{+})^{2}}}}.$
(121)

####  D.1.3 Superconducting-state renormalized boson frequency

In the superconducting state, the boson frequency (5) is altered with respect
to the normal-state result (18), due to the opening of the gap
$\Delta(i\omega_{n})$ in the spectral function. The boson self-energy follows
Eq. (4c), which can be recast in terms of the dynamical quasiparticle weight
$Z(i\omega_{n})$ by the means of Eqs. (6) and (9):

$\Pi(i\Omega_{n})=-2\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\sum_{m=-\infty}^{+\infty}\left[\frac{1}{i\omega_{m}Z(i\omega_{m})}\right.\\\
\left.\times\frac{1}{(i\omega_{m}+i\Omega_{n})Z(i\omega_{m}+i\Omega_{n})}\right.\\\
\left.-F(i\omega_{m})F(i\omega_{m}+i\Omega_{n})\right].$ (122)

We now approximate the dynamical weight $Z(i\omega_{n})$ by Eq. (116) in the
superconducting state, with the piecewise-constant gap (112). One can verify
both numerically and analytically that the $Z(i\omega_{n})$-dependent term in
Eq. (122) dominates over the $F(i\omega_{n})$-dependent term for
$\Delta\lessapprox 1$. Therefore, in the following we neglect the
$F(i\omega_{n})$-dependent term in Eq. (122). Using Eqs. (122) and (116), and
going to the $T=0$ limit, we are left with

$\Pi(0)=2\bar{g}^{2}\left\\{\int_{-\infty}^{-\omega_{0}}\frac{d\omega}{2\pi}\frac{1}{\omega^{2}}\frac{4}{\left[1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega^{2}}}\right]^{2}}\right.\\\
\left.+\int_{\omega_{0}}^{+\infty}\frac{d\omega}{2\pi}\frac{1}{\omega^{2}}\frac{4}{\left[1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega^{2}}}\right]^{2}}\right.\\\
\left.+\int_{-\omega_{0}}^{\omega_{0}}\frac{d\omega}{2\pi}\frac{1}{\omega^{2}}\frac{4}{\left[1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega^{2}+\Delta^{2}}}\right]^{2}}\right\\}.$
(123)

The three integrations in Eq. (123) can be performed analytically, with the
latter one yielding a rather lengthy expression in terms of imcomplete
elliptic integrals of the second kind. Expanding the explicit result of Eq.
(123) at leading order for $zt_{0}^{2}\rightarrow+\infty$, we achieve

$\Pi(0)\approx\frac{8\bar{g}^{2}}{6\pi\omega_{0}zt_{0}^{2}}\left(-3\Delta^{2}+2\sqrt{2}\omega_{0}\sqrt{z}t_{0}\right).$
(124)

Notice that the normal-state limit of Eq. (124), i.e., $\Delta\rightarrow 0$,
consistently yields Eq. (C.2.1). Using the definition (5) and Eq. (124), the
renormalized boson frequency becomes Eq. (48).

####  D.1.4 Gap equation

The gap equation stems from Eq. (4d) for the anomalous self-energy
$\Phi(i\omega_{n})$. Here we assume to be in the SYK2-FL regime, where
$\bar{g}^{2}/\omega_{0}^{2}\ll z{t_{0}^{2}}$. For the boson propagator, we
neglect the boson self-energy $\Pi(i\Omega_{n})\approx 0\,\forall i\omega_{n}$
altogether, which nevertheless yields a satisfactory agreement of the ensuing
analytical expressions, for the transition temperature and the zero-
temperature gap at vanishing energy, with the full numerical solutions of the
model. Hence, we assume Eq. (111) in the present regime. Using Eqs. (6) and
(7), together with Eq. (111), the gap equation (4d) becomes

$\displaystyle Z(i\omega_{n})\Delta(i\omega_{n})$
$\displaystyle\left\\{1-\frac{zt_{0}^{2}}{2}\right.$
$\displaystyle\left.\times\frac{1}{\left[Z(i\omega_{n})\right]^{2}\left\\{\omega_{n}^{2}+\left[\Delta(i\omega_{n})\right]^{2}\right\\}}\right\\}$
$\displaystyle=\frac{\bar{g}^{2}k_{B}T}{\omega_{0}}\sum_{m=-\infty}^{+\infty}Z(i\omega_{m})\Delta(i\omega_{m})$
$\displaystyle\times\frac{1}{\left[Z(i\omega_{m})\right]^{2}\left\\{\omega_{m}^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}.$
(125)

Further approximating the gap function with Eq. (112), and using the SYK2-FL
result (116) for the dynamical weight in the superconducting state, the left-
hand side of Eq. (D.1.4) simplifies to $\Delta$, and we are left with

$\displaystyle\Delta$
$\displaystyle=\frac{\bar{g}^{2}k_{B}T}{\omega_{0}}\sum_{n=-\infty}^{+\infty}2\Delta$
$\displaystyle\times\frac{1}{\left(1+\sqrt{1+\frac{2zt_{0}^{2}}{\omega_{n}^{2}+\Delta^{2}}}\right)\left(\omega_{n}^{2}+\Delta^{2}\right)}.$
(126)

In the large-hopping regime, we can approximate
$1+\sqrt{1+\left(2zt_{0}^{2}\right)/\left(\omega_{n}^{2}+\Delta^{2}\right)}\approx\sqrt{2z}t_{0}/\sqrt{\omega_{n}^{2}+\Delta^{2}}$,
so that

$1=\frac{\sqrt{2}\bar{g}^{2}k_{B}T}{\omega_{0}\sqrt{z}t_{0}}\sum_{n=-\infty}^{+\infty}\frac{1}{\sqrt{\omega_{n}^{2}+\Delta^{2}}}.$
(127)

Eq. (127) has now assumed the form of a BCS-like gap equation [172, 173, 174],
which we can employ to find the critical temperature and the zero-temperature
zero-energy gap in Fermi-liquid regime.

####  D.1.5 Zero-temperature gap: BCS formula

In the $T=0$ limit, where the Matsubara sums become continuous integrals, the
gap equation (127) in SYK2-FL regime translates as

$1=\frac{\sqrt{2}\bar{g}^{2}}{\pi\omega_{0}\sqrt{z}t_{0}}\int_{-\infty}^{+\infty}\frac{d\omega}{2}\frac{1}{\sqrt{\omega^{2}+\Delta^{2}}}.$
(128)

We use a ultraviolet cutoff $\omega_{0}$ for the integral in Eq. (128), since
pairing occurs at much lower energies than the bare boson energy $\omega_{0}$;
this feature is completely analogous to BCS theory [172, 173, 174]. Then we
have

$\frac{1}{\bar{\lambda}}=\int_{-\omega_{0}}^{+\omega_{0}}\frac{d\omega}{2}\frac{1}{\sqrt{\omega^{2}+\Delta^{2}}}=\mathrm{arcsinh}\left(\frac{\omega_{0}}{\Delta}\right),$
(129)

where we retrieve the superconducting coupling constant (36). Eq. (129) is the
BCS formula (43).

####  D.1.6 Superfluid phase stiffness

Approximating the gap to a constant with Eq. (112), and using the
superconducting dynamical weight (116) in SYK2-FL regime, together with Eq.
(7), Eq. (65a) becomes

$\displaystyle\rho_{S}(T)$
$\displaystyle=\Theta_{L}zt_{0}^{2}k_{B}T\sum_{i\omega_{n}}\frac{4\Delta^{2}}{\left[Z(i\omega_{n})\right]^{2}\left(\omega_{n}^{2}+\Delta^{2}\right)^{2}}$
$\displaystyle=\Theta_{L}zt_{0}^{2}k_{B}T\sum_{i\omega_{n}}\frac{1}{\left(\omega_{n}^{2}+\Delta^{2}\right)^{2}}$
$\displaystyle\times\frac{4\Delta^{2}}{\left(1+\sqrt{1+\frac{2zt_{0}^{2}}{\Delta^{2}+\omega_{n}^{2}}}\right)^{2}}.$
(130)

For $zt_{0}^{2}\gg\left|\omega_{n}^{2}+\Delta^{2}\right|$, we can expand as

$\rho_{S}(T)\approx
2\Theta_{L}\Delta^{2}k_{B}Tk_{B}T\sum_{i\omega_{n}}\frac{1}{\omega_{n}^{2}+\Delta^{2}}.$
(131)

Using the Matsubara sum

$k_{B}T\sum_{i\omega_{n}}\frac{1}{\omega_{n}^{2}+\mathscr{E}^{2}}=\frac{1}{2\mathscr{E}}\left[1-2f_{FD}(\mathscr{E})\right]$
(132)

with $\mathscr{E}\equiv\Delta$ in Eq. (131), we finally achieve Eq. (68).

###  D.2 Strong-coupling non-Fermi liquid: the impurity-like regime

####  D.2.1 Dynamical quasiparticle weight

To analyze the dynamical quasiparticle weight in the impurity-like regime
$\bar{g}^{2}/\omega_{0}\gg zt_{0}^{2}$, we start again from the normal self-
energy (4a) written in terms of $Z(i\omega_{n})$, using Eq. (6). We now employ
the Eliashberg equation (4e) for $G(i\omega_{n})$, and we approximate the
boson propagator with its impurity-like expression in the superconducting
state [105]

$D(i\omega_{n}-i\omega_{m})\approx\frac{1}{(\omega_{n}-\omega_{m})^{2}+(\omega_{r}^{SC})^{2}},$
(133)

with the finite renormalized boson frequency (49). The latter appears as an
in-line equation before Eq. (40) in Ref. 105. We also use Eq. (7) for
$\Delta(i\omega_{n})$, with the result

$\displaystyle Z(i\omega_{n})$
$\displaystyle=1+\frac{\bar{g}^{2}k_{B}T}{\omega_{0}^{2}}\sum_{m=-\infty}^{+\infty}\frac{1}{(\omega_{n}-\omega_{m})^{2}+(\omega_{r}^{SC})^{2}}\frac{\omega_{m}}{\omega_{n}}$
$\displaystyle\times\frac{1}{Z(i\omega_{n})\left\\{\omega_{m}^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}$
$\displaystyle+\frac{zt_{0}^{2}}{2}\frac{1}{Z(i\omega_{n})\left\\{\omega_{m}^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}.$
(134)

In the impurity-like regime, the $t_{0}$-dependent term in Eq. (D.2.1) is
negligible with respect to the fermion-boson term that depends on
$\bar{g}^{2}$. However, as we will see in the following, we can obtain an
analytical result for $Z(i\omega_{n})$ even at finite hopping, i.e., retaining
the $t_{0}$-dependent part of Eq. (D.2.1). While formally valid for
$\bar{g}^{2}/\Omega_{0}\gg z{t_{0}^{2}}$, the resulting expression – see Eqs.
(D.2.1) and (53), at zero and finite temperature respectively – allows us to
smoothly interpolate between the SYK2-FL and impurity-like regimes, which
occurs because the two regimes essentially share the same physics: fermions
randomly interacting with static scattering centers. At $T=0$, we convert the
Matsubara sum in Eq. (D.2.1) into an integral and we achieve

$\displaystyle Z(i\omega)$
$\displaystyle=1+\frac{\bar{g}^{2}}{\omega_{0}^{2}}\int_{-\infty}^{+\infty}\frac{d\omega^{\prime}}{2\pi}\frac{1}{(\omega-\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}\frac{\omega^{\prime}}{\omega}$
$\displaystyle\frac{1}{Z(i\omega^{\prime})\left\\{(\omega^{\prime})^{2}+\left[\Delta(i\omega^{\prime})\right]^{2}\right\\}}$
$\displaystyle+\frac{zt_{0}^{2}}{2}\frac{1}{Z(i\omega)\left\\{\omega^{2}+\left[\Delta(i\omega)\right]^{2}\right\\}}.$
(135)

Now, in the impurity-like regime, we have $\omega_{r}^{SC}\ll k_{B}T$. Then,
the term depending on $\omega_{r}^{SC}$ in Eq. (D.2.1) behaves like a Dirac
delta function [105],

$\lim_{\omega_{r}^{SC}\rightarrow
0}\frac{1}{(\omega-\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}\approx\frac{\pi}{\omega_{r}^{SC}}\delta(\omega-\omega^{\prime}).$
(136)

Inserting Eq. (136) into the quasiparticle weight (D.2.1) yields

$Z(i\omega)=1+\frac{\bar{g}^{2}}{\omega_{0}^{2}\omega_{r}^{SC}}\int_{-\infty}^{+\infty}\frac{d\omega^{\prime}}{2\pi}\pi\delta(\omega-\omega^{\prime})\\\
\times\frac{1}{Z(i\omega^{\prime})\left\\{(\omega^{\prime})^{2}+\left[\Delta(i\omega^{\prime})\right]^{2}\right\\}}\\\
+\frac{zt_{0}^{2}}{2}\frac{1}{Z(i\omega)\left\\{\omega^{2}+\left[\Delta(i\omega)\right]^{2}\right\\}}\\\
=1+\left[\left(\frac{8g^{2}\omega_{0}}{3\pi}\right)^{2}+\frac{z{t_{0}^{2}}}{2}\right]\\\
\times\frac{1}{Z(i\omega)\left\\{\omega^{2}+\left[\Delta(i\omega)\right]^{2}\right\\}},$
(137)

where we used $g^{2}=\bar{g}^{2}/\omega_{0}^{3}$. Therefore, Eq. (137) gives
the solutions

$Z(i\omega)=\frac{1}{2}\left[1\pm\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega^{2}+\left[\Delta(i\omega)\right]^{2}}}\right].$
(138)

To comply with the normal-state result, we choose the positive root of Eq.
(138). Actually, one can numerically confirm that the gap $\Delta(i\omega)$
can be approximated by a constant at low energies, as in Eq. (112), also in
the impurity-like regime. Therefore, we can write

$\displaystyle Z(i\omega)$
$\displaystyle=\frac{1}{2}\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega^{2}+\Delta^{2}}}\right]$
$\displaystyle\approx\frac{1}{2}\left(1+\sqrt{1+\frac{\left[16g^{2}\omega_{0}/(3\pi)\right]^{2}}{\omega^{2}+\Delta^{2}}}\right),$
(139)

where the last step is valid for $\bar{g}^{2}/\omega_{0}\gg zt_{0}^{2}$.

At low but nonzero temperatures, we can substitute $i\omega\mapsto
i\omega_{n}$ in Eq. (D.2.1), which gives Eq. (53).

The normal-state limit of Eq. (53) is obtained by taking $\Delta\rightarrow
0$, with the result

$Z(i\omega_{n})\approx\frac{1}{2}\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega_{n}^{2}}}\right].$
(140)

At low frequencies $\left|\omega_{n}\right|\rightarrow 0^{+}$, we have

$Z(i\omega_{n})\approx\frac{\tilde{\Omega}}{2\left|\omega_{n}\right|}\approx\frac{8g^{2}\omega_{0}/(3\pi)}{\left|\omega_{n}\right|},$
(141)

The last step of Eq. (141), again valid for $\bar{g}^{2}/\omega_{0}\gg
zt_{0}^{2}$, is consistent with Eq. (26) of Ref. 105, through Eq. (6).

####  D.2.2 Fermion propagators and self-energies in the superconducting state

In the superconducting state, in general from Eqs. (D.2.1) and (6) we have

$\Sigma(i\omega)=\frac{i\omega}{2}\left(1-\sqrt{1+\frac{\tilde{\Omega}}{\omega^{2}+\Delta^{2}}}\right)$
(142)

at $T=0$. Instead, at finite temperature we can employ

$\Sigma(i\omega_{n})=\frac{i\omega_{n}}{2}\left(1-\sqrt{1+\frac{\tilde{\Omega}}{\omega_{n}^{2}+\Delta^{2}}}\right).$
(143)

The fermion propagator which corresponds to the self-energy (143) is

$G(i\omega_{n})=\frac{2}{i\omega_{n}+i\mathrm{sign}(\omega_{n})\sqrt{\omega_{n}^{2}+\frac{\omega_{n}^{2}}{\omega_{n}^{2}+\Delta^{2}}\tilde{\Omega}^{2}}},$
(144)

which reduces to Eq. (IV.3) in the normal state, where $\Delta=0$. Performing
the analytic continuation $i\omega_{n}\rightarrow\omega+i0^{+}$ on Eq. (144),
we obtain the (retarded) fermion propagator on the real axis:

$G^{R}(\omega)=\frac{2}{(\omega+i0^{+})\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\Delta^{2}-(\omega+i0^{+})^{2}}}\right]}.$
(145)

In addition, in the superconducting state the anomalous propagator (4f)
emerges. Performing the same steps leading to Eq. (120), but now substituting
the impurity-like dynamical weight (53) in Eq. (4f), we obtain

$F(i\omega_{n})=-\frac{2\Delta}{\left(\omega_{n}^{2}+\Delta^{2}\right)\left(1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega_{n}^{2}+\Delta^{2}}}\right)},$
(146)

where $\tilde{\Omega}$ satisfies Eq. (54). The retarded anomalous propagator
is obtained from the analytic continuation of Eq. (146), and it reads

$F^{R}(\omega)=-\frac{2\Delta}{\left[\Delta^{2}-(\omega+i0^{+})^{2}\right]\left[1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\Delta^{2}-(\omega+i0^{+})^{2}}}\right]}.$
(147)

Notice that, in the limit $\bar{g}^{2}/\omega_{0}\ll zt_{0}^{2}$, we directly
retrieve the SYK2-FL result (D.1.2) from Eq. (147).

####  D.2.3 Superconducting-state renormalized boson frequency

In the impurity-like superconducting state, the renormalized boson frequency
does not vanish at $T=0$ as would be predicted by the normal-state result
(110), due to the gapped nature of the bosons at $T<T_{c}$ [105]. Instead, at
$k_{B}T\ll\omega_{0}g^{2}$ we have Eq. (49). The latter can be obtained from
the high-energy behaviour of the spectral function in the superconducting
state [105]. More precisely, we exploit the approximate Dirac-delta form (136)
of the bosonic propagator (133) in the impurity-like regime, and we insert it
in Eq. (4a) for the fermionic self-energy in the $T=0$ limit. Due to the Dirac
delta sampling only the $\Omega=0$ component of the bosonic propagator
$D(i\Omega)$, the integral over bosonic energy $\Omega$ yields

$\Sigma(i\omega)=\frac{\bar{g}^{2}}{2\pi}\frac{\pi}{\omega_{r}^{SC}(0)}G(i\omega)+\frac{zt_{0}^{2}}{2}G(i\omega).$
(148)

Assuming $\Sigma(\omega)\approx-1/G(i\omega)$ in Eq. (148) leaves a relation
between $\omega_{r}^{SC}$ and $G(i\omega)$. Using $G(i\omega)=1/\left[i\omega
Z(i\omega)\right]$ and the high-energy leading-order expansion of the
dynamical weight (D.2.1), which gives

$\left|\omega\right|Z(i\omega)\approx\frac{\tilde{\Omega}}{2}$ (149)

even in the superconducting state (the gap does not contribute at leading
order for $\omega\rightarrow+\infty$), we finally achieve

$\displaystyle\omega_{r}^{SC}(0)$
$\displaystyle\approx-\frac{\bar{g}^{2}\left[G(i\omega)\right]^{2}}{2+zt_{0}^{2}\left[G(i\omega)\right]^{2}}$
$\displaystyle=\frac{\omega_{0}}{2}\left(\frac{3\pi}{8g}\right)^{2}.$ (150)

Eq. (D.2.3) is consistent with Eq. (49). Notice that Eq. (D.2.3) has been
derived in the impurity-like regime, but in the presence of small hopping
$zt_{0}^{2}\ll\bar{g}^{2}/\omega_{0}$. This means that small hopping
perturbations with respect to fermion-boson coupling are not expected to
significantly alter the renormalized boson frequency. This is indeed observed
in Fig. 11(b).

####  D.2.4 Gap equation

We consider again Eq. (4d) for the anomalous self-energy, with the anomalous
propagator (4f), in the impurity-like regime where
$\bar{g}^{2}/\omega_{0}^{2}\gg zt_{0}^{2}$. The boson propagator can be
approximated by Eq. (133). Then, using Eqs. (6) and (7), and approximating the
gap function with Eq. (112), we translate Eq. (4d) into

$\displaystyle
Z(i\omega_{n})\left\\{1-\frac{zt_{0}^{2}}{2}\frac{1}{\left[Z(i\omega_{n})\right]^{2}\left(\omega_{n}^{2}+\Delta^{2}\right)}\right\\}$
$\displaystyle=\bar{g}^{2}k_{B}T\sum_{m=-\infty}^{+\infty}\frac{1}{(\omega_{n}-\omega_{m})^{2}+(\omega_{r}^{SC})^{2}}$
$\displaystyle\times\frac{1}{Z(i\omega_{m})\left(\omega_{m}^{2}+\Delta^{2}\right)}.$
(151)

####  D.2.5 Zero-temperature gap: asymptotic strong-coupling limit

Let us consider the gap equation (D.2.4) at $T=0$. Then, the Matsubara sum
becomes an integral and

$\displaystyle
Z(i\omega)\left\\{1-\frac{zt_{0}^{2}}{2}\frac{1}{\left[Z(i\omega)\right]^{2}\left(\omega^{2}+\Delta^{2}\right)}\right\\}$
$\displaystyle=\bar{g}^{2}\int_{-\infty}^{+\infty}\frac{d\omega^{\prime}}{2\pi}\frac{1}{(\omega-\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}$
$\displaystyle\times\frac{1}{Z(i\omega)\left(\omega^{2}+\Delta^{2}\right)}.$
(152)

Let us analyze Eq. (D.2.5) in the single-dot limit $t_{0}\rightarrow 0$. In
such limit, Eq. (D.2.5) for the dispersionless Yukawa-SYK model can be
formally brought into the same form as the Eliashberg equation for a wide-band
system interacting with an Einstein phonon of frequency $\Omega_{E}$ in the
asymptotic strong-coupling limit [187] – see also the discussion after Eq.
(32) in Ref. 105. We can then borrow the asymptotic results for the gap and
the quasiparticle weight from strong-coupling Eliashberg theory, as follows.
Labeling $\Delta_{0}$ as in Eq. (41), we employ Eq. (53) and for the
quasiparticle weight, in the very large-coupling limit $g\rightarrow+\infty$,
we obtain explicitly

$Z(\omega)\approx\frac{8g^{2}\omega_{0}}{(3\pi)\sqrt{\omega^{2}+\Delta^{2}}}.$
(153)

On the other hand, the dynamical quasiparticle weight satisfies Eq. (D.2.1),
where we can neglect the $t_{0}$-dependent term deep in the impurity-like
regime. Hence, inserting the result (153) into Eq. (D.2.1), we have

$Z(i\omega)=1+\frac{g^{2}\omega_{0}}{2\pi\omega}\int_{-\infty}^{+\infty}d\omega^{\prime}\frac{1}{(\omega-\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}\\\
\times\frac{\omega^{\prime}}{Z(i\omega)\left(\omega^{2}+\Delta^{2}\right)}\\\
=1+\frac{3\pi\omega_{0}^{2}}{8}\frac{1}{2\pi\omega}\int_{-\infty}^{+\infty}d\omega^{\prime}\frac{1}{(\omega-\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}\\\
\times\frac{\omega^{\prime}}{\sqrt{\omega^{2}+\Delta^{2}}}.$ (154)

Eq. (154) yields an expression for the weight
$\left.Z(\omega)\right|_{\omega\approx 0}$ [187]:

$\displaystyle\left.Z(\omega)\right|_{\omega\approx 0}$
$\displaystyle=1+2\frac{3\omega_{0}^{2}}{16\omega_{r}^{SC}}\frac{d}{d\omega_{r}^{SC}}\left[(\omega_{r}^{SC})^{2}\mathscr{J}(\omega_{r}^{SC})\right]$
$\displaystyle+\frac{3\omega_{0}^{2}}{16}\frac{1}{2\Delta_{0}^{2}},$ (155)

where

$\mathscr{J}(\omega_{r}^{SC})=\int_{0}^{2\Delta_{0}}\frac{d\omega^{\prime}}{(\omega^{\prime})^{2}+(\omega_{r}^{SC})^{2}}\frac{1}{\sqrt{(\omega^{\prime})^{2}+\Delta_{0}^{2}}}.$
(156)

In the same way, using Eq. (156) the gap equation (D.2.5) for $t_{0}=0$
reduces to

$\left.Z(\omega)\right|_{\omega\approx
0}=2\frac{3\omega_{0}^{2}}{16}\mathscr{J}(\omega_{r}^{SC}).$ (157)

Therefore, we now have Eqs. (D.2.5) and (157) for the quasistatic weight,
which constitute a closed problem for the gap $\Delta_{0}$. To obtain the
latter analytically, we use the series expansion of Eq. (156),

$\mathscr{J}(\omega_{r}^{SC})\approx\frac{1}{\omega_{r}^{SC}\Delta_{0}}\left[\frac{\pi}{2}-\frac{\sqrt{5}}{2}\frac{\omega_{r}^{SC}}{\omega_{0}}-\frac{\sqrt{5}}{24}\left(\frac{\omega_{r}^{SC}}{\omega_{0}}\right)^{3}\right],$
(158)

into Eqs. (D.2.5) and (157). Solving the latter system for $\Delta_{0}$, we
finally achieve

$\Delta_{0}\approx\sqrt{\frac{-1+2\sqrt{5}}{2}}\sqrt{\frac{3}{16}}\omega_{0}\approx
0.5705372431149527\omega_{0}.$ (159)

An analogous calculation for $T_{c}$ in the asymptotic limit
$g\rightarrow+\infty$ would lead to [187]

$k_{B}T_{c}=\frac{1}{2\pi}\sqrt{\frac{3}{8}}\omega_{0}\approx
0.09746210015420952\omega_{0},$ (160)

so that the gap to $T_{c}$ ratio would be

$\frac{2\Delta_{0}}{k_{B}T_{c}}\approx 11.708,$ (161)

in good agreement with Eq. (46) obtained from the numerical solution of the
Eliashberg equations (4), and with Eq. (38) in Ref. 105. Actually, the
approximated analytical expressions (159) and (160) underestimate the gap and
$T_{c}$ by a factor $C^{-1}\approx 1/1.1481028752388993$, so that more
accurate results are

$\Delta_{0}\approx\sqrt{\frac{-1+2\sqrt{5}}{2}}\sqrt{\frac{3}{16}}C\omega_{0}\approx
0.6550354492511521\omega_{0}$ (162)

and

$k_{B}T_{\rm c}=\frac{1}{2\pi}\sqrt{\frac{3}{8}}C\omega_{0}\approx
0.111897\omega_{0}.$ (163)

Eq. (162) and (163) are in good agreement with the gap value quoted above Eq.
(38) and with Eq. (30) in Ref. 105. Notice that, even with the correction $C$,
the gap to $T_{\rm c}$ ratio still has the value (161).

####  D.2.6 Superfluid phase stiffness

Recognizing the formal equivalence between the results for the dynamical
weight $Z(i\omega_{n})$ in the SYK2-FL and impurity-like regimes, in
accordance with Eqs. (116) and (53), we can readily calculate the phase
stiffness in the impurity-like regime as well, employing the results of Sec.
D.1.6. Starting from Eq. (65a), the derivation is exactly the same upon
mapping $2z{t_{0}^{2}}\rightarrow\tilde{\Omega}^{2}$ in $Z(i\omega_{n})$, in
accordance with Eq. (54). We then have

$\rho_{S}=\Theta_{L}zt_{0}^{2}k_{B}T\sum_{i\omega_{n}}\frac{4\Delta^{2}}{\left(1+\sqrt{1+\frac{\tilde{\Omega}^{2}}{\omega_{n}^{2}+\Delta^{2}}}\right)^{2}\left(\omega_{n}^{2}+\Delta^{2}\right)^{2}}.$
(164)

In the regime where $\tilde{\Omega}^{2}\gg\omega_{n}^{2}+\Delta^{2}$, we can
write

$\rho_{S}=\frac{4\Theta_{L}\Delta^{2}zt_{0}^{2}}{\tilde{\Omega}^{2}}k_{B}T\sum_{i\omega_{n}}\frac{1}{\omega_{n}^{2}+\Delta^{2}}.$
(165)

Using Eq. (132) for the Matsubara sum, we are left with Eq. (66).

##  Appendix E Analysis of the critical temperature

###  E.1 Weak-coupling single-dot limit: the NFL-SYK regime

To analyze the weak-coupling limit in the single dot, we start from the
linearized gap equation (32c). For the SYK-NFL fixed point, the fermionic and
bosonic low-energy propagators are given by Eqs. (11) and (12), respectively,
with the renormalized boson frequency (13).

With Eqs. (11)-(12), neglecting the term $\Omega_{n}^{2}$ in Eq. (12) at weak
coupling [105], and using the relation (7), the linearized gap equation (32c)
becomes

$\displaystyle\Delta(i\omega_{n})$
$\displaystyle=a_{0}\sum_{m=-\infty}^{+\infty}(T_{f}/T)^{2\Delta}\mathrm{sign}(\omega_{m})\times$
$\displaystyle\frac{1}{(T_{f}/T)^{-2\Delta}\left|m+1/2\right|+\left|m+1/2\right|^{1-2\Delta}}$
$\displaystyle\qquad\times\left[\frac{\Delta(i\omega_{m})}{i\omega_{m}}-\frac{\Delta(i\omega_{n})}{i\omega_{n}}\right]\frac{1}{m_{0}+\left|n-m\right|^{4\Delta-1}},$
(166)

where $m_{0}=c_{2}/\left[c_{3}(2\pi)^{4\Delta-1}\right]$,
$a_{0}=1/\left[2\pi(c_{1})^{2}c_{2}\right]$, and $2\pi
k_{B}T_{f}=c_{1}^{1/(2\Delta)}\bar{g}^{2}/\omega_{0}^{2}\approx
0.1888\bar{g}^{2}/\omega_{0}^{2}$. Since the only temperature scale in Eq.
(166) is $T_{f}/T$, the transition temperature $T_{c}$ is of the order of
$T_{f}$ and is numerically found to be $k_{B}T_{c}\approx
0.16\bar{g}^{2}/\omega_{0}^{2}$ [105].

Figure 22:  Comparison between the numerical solution of the linearized
Eliashberg equations (32) for the Yukawa-SYK model on a lattice (point
markers) and the BCS-like formula (35) with the coupling constant (36) (dashed
lines) as a function of dimensionless interaction $g$, for different hopping
parameters $t_{0}$.

###  E.2 Weak-coupling BCS formula from a disordered Fermi liquid: the SYK2-FL
regime

Let us consider the SYK2-FL regime at $T=T_{c}$. Then the gap equation (127)
for $\Delta\rightarrow 0$ becomes

$1=\frac{\sqrt{2}\bar{g}^{2}k_{B}T}{\omega_{0}\sqrt{z}t_{0}}\sum_{m=-\infty}^{+\infty}\frac{1}{\left|\omega_{n}\right|}.$
(167)

The result (167) is also consistent with the series expansion of the summation
term in Eq. (39) at leading order in small $\omega_{m}$. At this point, we can
exploit the analogy with standard BCS theory, and write Eq. (167) in analogy
with its BCS counterpart [172, 173, 174], using the result

$\sum_{m=-\infty}^{+\infty}\frac{1}{\left|\omega_{n}\right|}\approx\ln\left(\frac{\bar{\Omega}}{k_{B}T}\frac{2e^{\gamma}}{\pi}\right),$
(168)

where $\bar{\Omega}$ is an appropriate UV cutoff energy and $\gamma\approx
0.5772156649$ is the Euler-Mascheroni constant. The approximation (168) is
valid for $\bar{\Omega}\gg k_{B}T$. The proper cutoff energy of the Yukawa-SYK
problem is $\bar{\Omega}=\omega_{0}$, so in the SYK2-FL regime we use Eqs.
(167) and (168) to achieve

$\frac{1}{\bar{\lambda}}=\ln\left(\frac{\omega_{0}}{k_{B}T}\frac{2e^{\gamma}}{\pi}\right),$
(169)

which is an effective BCS equation for the critical temperature in the Fermi-
liquid regime, with the coupling constant (36). Inverting Eq. (169), we find
the critical temperature (35).

We find good agreement between Eq. (35) and the numerics for $T_{c}$,
employing the full-fledged linearized Eliashberg equations (32), as shown Fig.
22.

###  E.3 Strong-coupling single-dot limit: the impurity-like regime

For completeness, we start by reviewing the strong-coupling analysis in the
single-dot limit, i.e., $t_{0}=0$, reported in Ref. 105. We will exploit this
derivation to deduce the same infinite-coupling limit on the lattice. Using
the functions (6) and (7), the gap equation (4d) can be equivalently rewritten
as

$\displaystyle Z(i\omega_{n})\Delta(i\omega_{n})$
$\displaystyle=\bar{g}^{2}k_{B}T\sum_{m}\frac{D(\omega_{n}-\omega_{m})\Delta(i\omega_{m})}{Z(i\omega_{m})\left\\{(\omega_{m})^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}.$
(170)

On the other hand, the self-energy equation (4a) for $t_{0}=0$ translates as

$\displaystyle Z(i\omega_{n})$
$\displaystyle=1+\bar{g}^{2}k_{B}T\sum_{m}\frac{D(\omega_{n}-\omega_{m})}{Z(i\omega_{m})\left\\{(\omega_{m})^{2}+\left[\Delta(i\omega_{m})\right]^{2}\right\\}}$
$\displaystyle\times\frac{\omega_{m}}{\omega_{n}}.$ (171)

Combining Eq. (E.3) and (E.3) gives [105, 106]

$\displaystyle\Delta(i\omega_{n})$
$\displaystyle=\bar{g}^{2}k_{B}T\sum_{m}\frac{D(\omega_{n}-\omega_{m})}{Z(i\omega_{m})\sqrt{(\omega_{m})^{2}+\left[\Delta(i\omega_{m})\right]^{2}}}$
$\displaystyle\times\frac{1}{\sqrt{(\omega_{m})^{2}+\left[\Delta(i\omega_{m})\right]^{2}}}\left[\Delta(i\omega_{m})-\frac{\omega_{m}}{\omega_{n}}\Delta(i\omega_{n})\right].$
(172)

The linearized version of Eq. (172) is

$\displaystyle\Delta(i\omega_{n})$
$\displaystyle=\bar{g}^{2}k_{B}T\sum_{m}\frac{D(\omega_{n}-\omega_{m})}{Z(i\omega_{m})(\omega_{m})^{2}}$
$\displaystyle\times\left[\Delta(i\omega_{m})-\frac{\omega_{m}}{\omega_{n}}\Delta(i\omega_{n})\right].$
(173)

Using the normal-state boson propagator (21) and the low-frequency strong-
coupling quasiparticle residue $\left|\omega_{n}\right|Z(i\omega_{n})\approx
8\bar{g}^{2}/\left(3\pi\omega_{0}^{2}\right)$, which results from Eqs. (IV.3)
and (6), we obtain

$\Delta(i\omega_{n})=\frac{3\pi\omega_{0}^{2}k_{B}T}{8}\sum_{m}\frac{1}{(\omega_{n}-\omega_{m})^{2}+(\omega_{r})^{2}}\left[\frac{\Delta(i\omega_{m})}{\omega_{m}}-\frac{\Delta(i\omega_{n})}{\omega_{n}}\right]\mathrm{sign}(\omega_{m}).$
(174)

Since $k_{B}T\gg\omega_{r}$ in the strong-coupling limit, and since the zeroth
Matsubara frequency does not participate to pairing [105], we set
$\omega_{r}=0$ in Eq. (174). Hence, as the fermionic Matsubara frequencies are
$\omega_{n}=(2n+1)\pi k_{B}T$, Eq. (174) becomes

$\Delta(i\omega_{n})=\frac{3\omega_{0}^{2}}{8(\pi k_{B}T)^{2}}\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\left[\frac{\Delta(i\omega_{m})}{2m+1}-\frac{\Delta(i\omega_{n})}{2n+1}\right]\mathrm{sign}\left(m+\frac{1}{2}\right).$
(175)

Now, we can rewrite

$\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\left[\frac{\Delta(i\omega_{m})}{\Delta(i\omega_{n})}\frac{1}{2m+1}-\frac{1}{2n+1}\right]\mathrm{sign}\left(m+\frac{1}{2}\right)\\\
=\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\left[\frac{\Delta(i\omega_{m})}{\Delta(i\omega_{n})}\frac{1}{\left|2m+1\right|}-\frac{\mathrm{sign}\left(m+\frac{1}{2}\right)}{2n+1}\right]\\\
=\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\left[\frac{r_{mn}}{\left|2m+1\right|}-\frac{r_{nn}}{2n+1}\mathrm{sign}\left(m+\frac{1}{2}\right)\right],$
(176)

where $r_{mn}=\Delta(i\omega_{m})/\Delta(i\omega_{n})$. Hence, Eq. (175) is an
eigenvalue problem for the ratios $r_{mn}$. Additionally, in Eq. (175) we have

$-\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\frac{\mathrm{sign}\left(m+\frac{1}{2}\right)}{2n+1}=\\\
\frac{1}{2n+1}\left\\{\sum_{m=-\infty}^{-1}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}-\sum_{m=0}^{n-1}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}-\sum_{m=n+1}^{+\infty}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\right\\}=\\\
\frac{1}{2n+1}\left[\frac{\psi(1+n)}{4}-\frac{\pi^{2}}{24}+\frac{\psi(1+n)}{4}-\frac{\pi^{2}}{24}\right]=\left[-\frac{\pi^{2}}{12}+\frac{\psi(1+n)}{2}\right]\frac{1}{2n+1},$
(177)

where $\psi(z)=\lim_{M\rightarrow+\infty}\left[\ln
M-\sum_{n=0}^{M}1/(n+z)\right]$ is the digamma function. Therefore

$1=\alpha\left\\{\underbrace{\sum_{m\neq
n}\frac{1}{\left[2n+1-(2m+1)\right]^{2}}\frac{r_{mn}}{\left|2m+1\right|}+\left[-\frac{\pi^{2}}{12}+\frac{\psi(1+n)}{2}\right]\frac{r_{nn}}{2n+1}}_{\sum_{n,m}\Lambda_{nm}r_{mn}}\right\\},$
(178)

with $\alpha=3\omega_{0}^{2}/[8(\pi k_{B}T)^{2}]$ and $r_{nn}=1\,\forall n$.
The strategy to solve Eq. (178) is to diagonalize the matrix
$\underline{\underline{\Lambda}}$ formed by the coefficients of the ratios
$r_{mn}$: at $T=T_{c}$, the largest eigenvalue of
$\alpha\underline{\underline{\Lambda}}$ has to be unitary. Of course, this is
the same criterion we adopt to numerically find $T_{c}$ for the full
linearized problem (34), as described in Appendix B.3. Without loss of
generality, we take $0\leq n<+\infty$ and $0\leq m<+\infty$, and the diagonal
matrix elements can be written as

$\Lambda_{nn}=\frac{1}{2n+1}\frac{1}{2^{2}(2n+1)^{2}}+\left[-\frac{\pi^{2}}{12}+\frac{\psi(1+n)}{2}\right]\frac{1}{2n+1},$
(179)

where the first term comes from $m=-n-1$ in Eq. (178), when we rewrite the sum
$\sum_{m\neq n}$ with $-\infty<m<\infty$ as a sum $\sum_{m\neq n}$ with $m>0$
exploiting the even parity of $\left|2m+1\right|$ and of $r_{nm}$ under the
exchange $m\mapsto-n-1$. The off-diagonal matrix elements are

$\displaystyle\Lambda_{nm}$
$\displaystyle=\frac{1}{2n+1}\left\\{\frac{1}{[2n+1-(2m+1)]^{2}}\right.$
$\displaystyle\left.+\frac{1}{[2n+1+(2m+1)]^{2}}\right\\},$ (180)

where the first term comes from $m>0$ and the second term from $m<0$. The
largest eigenvalue of the matrix $\underline{\underline{\Lambda}}$ is

$\lambda_{\Lambda}=0.32953505303295694.$ (181)

Therefore, the $T_{c}$ equation amounts to

$1=\alpha_{c}\lambda_{\Lambda},\,\alpha_{c}=\frac{3(\omega_{0})^{2}}{8(\pi
k_{B}T)^{2}}.$ (182)

Eq. (182) is solved for $T_{c}>0$, to finally give with reasonable accuracy

$\displaystyle k_{B}T_{c}$
$\displaystyle=\frac{\sqrt{\alpha_{c}}}{\pi}\sqrt{\frac{3}{8}}=0.18272624777228844\sqrt{\frac{3}{8}}\omega_{0}$
$\displaystyle=0.11189651741386951\omega_{0},$ (183)

which is Eq. (40). Hence, in the strong-coupling regime $T_{c}$ does not
depend on $\bar{g}^{2}$ for an isolated dot coupled to an Einstein boson.
Notice that Eq. (40) is in satisfactory agreement with the analytical
estimation (163).

###  E.4 Impurity-like regime on the lattice

The estimation of the asymptotic limit for $T_{c}$ on the lattice closely
follows the developments in Appendix E.3, modulo the alterations to the
fermion self-energy, the gap equation, and the renormalized boson frequency
due to hopping. The gap equation becomes Eq. (32c), while the normal part of
the self-energy follows Eq. (4a). However, combining Eqs. (32c) and (4a),
together with Eqs. (6) and (7), we realize that the hopping-dependent part
disappears from the resulting gap equation. Then, using the low-frequency
quasiparticle residue
$\left|\omega_{n}\right|Z(i\omega_{n})\approx\sqrt{zt_{0}^{2}/2+\left[8\bar{g}^{2}/\left(3\pi\omega_{0}^{2}\right)\right]^{2}}$,
which results from Eqs. (144) and (6), neglecting the $t_{0}$-dependent part
of $Z(i\omega_{n})$ in the $g\rightarrow+\infty$ limit, we retrieve the
single-dot Eq. (174). One can confirm that even on the lattice, in the
impurity-like regime we can approximate $\omega_{r}\approx 0$, which yields
again Eq. (175). Therefore, the limit $\lim_{g\rightarrow+\infty}T_{c}$ on the
lattice is the same as the one in the single dot, namely, Eq. (40).

##  Appendix F Derivation of the lattice grand potential

###  F.1 Normal-state grand potential

The grand potential in the normal state stems from the disorder-averaged
Yukawa-SYK action (2), where the anomalous propagator $F(i\omega_{n})$ and
self-energy $\Phi(i\omega_{n})$ are both null. Assuming translation invariance
in time and space after the disorder average, the effective action per lattice
site $\bm{x}$ in imaginary time is given by

$\displaystyle\frac{\mathscr{S}}{\mathscr{N}}=$
$\displaystyle-\text{Tr}\log\left\\{\left(\hat{G}_{0}^{-1}-\hat{\Sigma}\right)\right\\}$
$\displaystyle+\frac{1}{2}\text{Tr}\log\text{$\left\\{\left(D_{0}^{-1}-\Pi\right)\right\\}$
}$ $\displaystyle-2\int_{\tau_{12}\tau_{21}}G(\tau_{21})\Sigma(\tau_{12})$
$\displaystyle+\frac{1}{2}\int_{\tau_{12}\tau_{21}}D(\tau_{21})\Pi(\tau_{12})$
$\displaystyle+\bar{g}^{2}\int_{\tau_{12}\tau_{21}}G(\tau_{12})G(\tau_{21})D(\tau_{12})$
$\displaystyle+\frac{zt_{0}^{2}}{2}\int_{\tau_{12}\tau_{21}}G(\tau_{12})G(\tau_{21}).$
(184)

Here $\hat{\Sigma}(i\epsilon_{n})={\rm
diag}\left\\{\Sigma(i\epsilon_{n}),-\Sigma(-i\epsilon_{n})\right\\}$, while
$\hat{G}_{0}^{-1}(\tau,\tau^{\prime})=-(\partial_{\tau}-\mu)\delta(\tau-\tau^{\prime})\,\hat{\sigma}_{0}$
– see also Appendix A.

On the normal-state saddle point the fields are homogeneous in time and space,
and are given by

$\displaystyle\hat{G}^{-1}(i\omega_{n})=\hat{G}_{0}^{-1}(i\epsilon_{n})-\hat{\Sigma}(i\epsilon_{n}),$
(185a) $\displaystyle
D^{-1}(i\Omega_{n})=D_{0}^{-1}(i\Omega_{n})-\Pi(i\Omega_{n}),$ (185b)
$\displaystyle\Sigma(\tau)=\bar{g}^{2}G(\tau)D(\tau)+\frac{zt_{0}^{2}}{2}G(\tau),$
(185c) $\displaystyle\Pi(\tau)=-2\bar{g}^{2}G(\tau)G(-\tau).$ (185d)

We now go on-shell and insert Eqs. (185) into the action (F.1). Focusing on
the integrands in the last two lines of the effective action, and using Eqs.
(185c) and (185d), we obtain

$\displaystyle-2G(\tau_{21})\Sigma(\tau_{12})+\frac{1}{2}D(\tau_{21})\Pi(\tau_{12})$
$\displaystyle\qquad+\bar{g}^{2}G(\tau_{12})G(\tau_{21})D(\tau_{12})+\frac{zt_{0}^{2}}{2}G(\tau_{12})G(\tau_{21})$
$\displaystyle\qquad\qquad\equiv
D(\tau_{21})\Pi(\tau_{12})-\frac{zt_{0}^{2}}{2}G(\tau_{12})G(\tau_{21}).$
(186)

The trace-log term can be written in Matsubara space as

$\displaystyle-\text{Tr}\log(\hat{G}^{-1})$ $\displaystyle=$
$\displaystyle\sum_{i\omega_{n}}\text{Tr}_{{\rm
N}}\log\left\\{\hat{G}(i\omega_{n})/\beta\right\\},$ (187)

where $\text{Tr}_{{\rm N}}$ denotes the trace over the Nambu subspace. Using
also the Dyson equations (185a) and (185b), we end up with

$\displaystyle\frac{\mathscr{S}^{{\rm sp}}}{\mathscr{N}}=$
$\displaystyle\sum_{n}\text{Tr}_{{\rm
N}}\log\left\\{\hat{G}(i\omega_{n})/\beta\right\\}-\frac{1}{2}\sum_{n}\text{$\log\left\\{D(\Omega_{n})/\beta^{2}\right\\}$
}$ $\displaystyle+$
$\displaystyle\int_{\tau_{12}\tau_{21}}D(\tau_{21})\Pi_{x}(\tau_{12})-\frac{zt_{0}^{2}}{2}\int_{\tau_{12}\tau_{21}}G(\tau_{12})G(\tau_{21}).$
(188)

The fermionic trace-log term can be simplified by using the explicit
expression of the Nambu propagator:

$\displaystyle\text{Tr}_{{\rm
N}}\log\left\\{\hat{G}(i\omega_{n})/\beta\right\\}$
$\displaystyle\qquad=\log{\rm
det_{N}}\left(\begin{array}[]{cc}G(i\omega_{n})/\beta&0\\\
0&-G(-i\omega_{n})/\beta\end{array}\right)$ (191)
$\displaystyle\qquad\qquad=2\log[G(i\omega_{n})/\beta],$ (192)

where we have used the odd parity of normal-state fermionic propagator with
respect to $\omega_{n}$. We Matsubara-transform the terms depending on
$\tau_{12}$ in Eq. (F.1), then we add and subtract the grand potential density
(per flavor $\mathscr{N}$) for free fermions and bosons in equilibrium:

$\displaystyle\frac{\beta\Omega_{c}^{(0)}}{\mathscr{N}}$
$\displaystyle=-\ln(1+e^{\beta\mu})$
$\displaystyle=\sum_{i\omega_{n}}\log\left[G_{0}(i\omega_{n})/\beta\right],$
(193a) $\displaystyle\frac{\beta\Omega_{b}^{(0)}}{\mathscr{N}}$
$\displaystyle=\frac{\beta\omega_{0}}{2}+\log(1-e^{-\beta\omega_{0}})$
$\displaystyle=\sum_{i\Omega_{n}}\log\left[D_{0}(i\Omega_{n})/\beta^{2}\right],$
(193b)

with bare propagators $G_{0}^{-1}(i\omega_{n})=i\omega_{n}+\mu$ and
$D_{0}^{-1}(i\Omega_{n})=\Omega_{n}^{2}+\omega_{0}^{2}$. Inserting Eqs. (193)
in the Matsubara-transformed version of the action (F.1) gives Eq. (VI).

###  F.2 Superconducting-state grand potential

To calculate the grand potential in the superconducting state, we have to
consider the full Yukawa-SYK action (2), including the anomalous Green’s
function and self-energy. Working within the spin-singlet ansatz, and assuming
translation invariance in time and space, the effective action per lattice
site $\bm{x}$ is given by

$\displaystyle\frac{\mathscr{S}}{\mathscr{N}}=$
$\displaystyle-\text{Tr}\log\left(\hat{G}_{0}^{-1}-\hat{\Sigma}\right)+\frac{1}{2}\text{Tr}\log\left({D}_{0}^{-1}-{\Pi}\right)$
$\displaystyle-2\int_{\tau_{12}\tau_{21}}G(\tau_{21})\Sigma(\tau_{12})+\frac{1}{2}\int_{\tau_{12}\tau_{21}}D(\tau_{21})\Pi(\tau_{12})$
$\displaystyle-\int_{\tau_{12}\tau_{21}}\left[F(\tau_{21})\Phi^{\dagger}(\tau_{12})+F^{\dagger}(\tau_{21})\Phi(\tau_{12})\right]$
$\displaystyle+\bar{g}^{2}\int_{\tau_{12}\tau_{21}}\left[G(\tau_{12})G(\tau_{21})-F^{\dagger}(\tau_{12})F(\tau_{21})\right]D(\tau_{12})$
$\displaystyle+\frac{zt_{0}^{2}}{2}\int_{\tau_{12}\tau_{21}}\left[G(\tau_{12})G(\tau_{21})-F^{\dagger}(\tau_{12})F(\tau_{21})\right].$
(194)

The saddle-point equations (4), entirely written in imaginary time, are

$\displaystyle\hat{G}(\tau)$ $\displaystyle=$
$\displaystyle\left[\hat{G}_{0}^{-1}(\tau)-\hat{\Sigma}(\tau)\right]^{-1},$
(195a) $\displaystyle D(\tau)$ $\displaystyle=$
$\displaystyle\left[{D}_{0}^{-1}(\tau)-{\Pi}(\tau)\right]^{-1},$ (195b)
$\displaystyle\Sigma(\tau)$ $\displaystyle=$
$\displaystyle\bar{g}^{2}G(\tau)D(\tau)+\frac{zt_{0}^{2}}{2}G(\tau),$ (195c)
$\displaystyle\Phi(\tau)$ $\displaystyle=$
$\displaystyle-\bar{g}^{2}F(\tau)D(\tau)-\frac{zt_{0}^{2}}{2}F(\tau),$ (195d)
$\displaystyle\Pi(\tau)$ $\displaystyle=$
$\displaystyle-2\bar{g}^{2}\left[G(\tau)G(-\tau)-F^{\dagger}(\tau)F(-\tau)\right].$
(195e)

Introducting the saddle-point equations (195) in Eq. (F.2), we find the on-
shell action

$\displaystyle\frac{\mathscr{S}^{{\rm sp}}}{\mathscr{N}}$ $\displaystyle=$
$\displaystyle-\text{Tr}\log(\hat{G}^{-1})+\frac{1}{2}\text{Tr}\log(D^{-1})$
(196) $\displaystyle-$ $\displaystyle
2\bar{g}^{2}\int_{\tau_{12}\tau_{21}}\left[G(\tau_{21})G(\tau_{12})-F^{\dagger}(\tau_{21})F(\tau_{12})\right]D(\tau_{12})$
$\displaystyle-$
$\displaystyle\frac{zt_{0}^{2}}{2}\int_{\tau_{12}\tau_{21}}\left[G(\tau_{12})G(\tau_{21})-F^{\dagger}(\tau_{12})F(\tau_{21})\right],$

where we have used the notation $\tau_{12}=\tau_{1}-\tau_{2}.$ Now, if we
further use Eq. (195) we can rewrite the second line in terms of the boson
self-energy:

$\displaystyle\frac{\mathscr{S}^{{\rm sp}}}{\mathscr{N}}$ $\displaystyle=$
$\displaystyle-\text{Tr}\log(\hat{G}^{-1})+\frac{1}{2}\text{Tr}\log(D^{-1})$
$\displaystyle+$
$\displaystyle\int_{\tau_{12}\tau_{21}}\Pi(\tau_{21})D(\tau_{12})$
$\displaystyle-$
$\displaystyle\frac{zt_{0}^{2}}{2}\int_{\tau_{12}\tau_{21}}\left[G(\tau_{12})G(\tau_{21})-F^{\dagger}(\tau_{12})F(\tau_{21})\right].$

Let us manipulate the fermionic trace-log term. In analogy with the normal-
state calculation, we go to Matsubara space and employ Eq. (187), where now
the trace is performed in Nambu subspace. We recall that, in our conventions,
$\hat{G}$ is a $2\times 2$ matrix. We have

$\displaystyle\log\left[\frac{1}{\beta^{2}}\text{det}_{{\rm
N}}\left\\{\hat{G}(i\omega_{n})\right\\}\right]$
$\displaystyle\qquad=\log\left[\frac{1}{\beta^{2}}\text{det}_{{\rm
N}}\left(\begin{array}[]{cc}G(i\omega_{n})&F(i\omega_{n})\\\
F^{\dagger}(i\omega_{n})&G(i\omega_{n})\end{array}\right)\right]$ (200)
$\displaystyle\qquad\qquad=\log\left\\{\frac{1}{\beta^{2}}\left[G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})\right]\right\\}.$
(201)

Using Eq. (F.2) in the on-shell action (196), we obtain

$\displaystyle\frac{\mathscr{S}^{{\rm sp}}}{\mathscr{N}}$ $\displaystyle=$
$\displaystyle\sum_{i\omega_{n}}\log\left[\frac{1}{\beta^{2}}(G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n}))\right]$
(202) $\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{i\Omega_{n}}\text{$\log\left\\{D(i\Omega_{n})/\beta^{2}\right\\}$
}+\sum_{i\Omega_{n}}D(i\Omega_{n})\Pi(i\Omega_{n})$ $\displaystyle-$
$\displaystyle\frac{zt_{0}^{2}}{2}\sum_{i\omega_{n}}\left[G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})\right].$

Setting $F(i\omega_{n})=F^{\dagger}(i\omega_{n})=0$ in Eq. (202), we retrieve
Eq. (VI), consistently with the normal-state calculations. For numerical
reasons we normalize the ${\rm log}$ terms as follows:

$\displaystyle\frac{\mathscr{S}^{{\rm sp}}}{{\mathscr{N}}_{t}}$
$\displaystyle=$
$\displaystyle\sum_{i\omega_{n}}\left\\{\log\left[\frac{G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})}{G_{0}^{2}(i\omega_{n})}\right]\right.$
$\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{i\Omega_{n}}\log\left[\frac{D(i\Omega_{n})}{D_{0}(i\Omega_{n})}\right]+\sum_{i\Omega_{n}}D(i\Omega_{n})\Pi(i\Omega_{n})$
$\displaystyle-$
$\displaystyle\left.\frac{zt_{0}^{2}}{2}\sum_{i\omega_{n}}\left[G^{2}(i\omega_{n})-F^{\dagger}(i\omega_{n})F(i\omega_{n})\right]\right\\}$
$\displaystyle+$ $\displaystyle
2\sum_{i\omega_{n}}\log\left[G_{0}(i\omega_{n})/\beta\right]-\frac{1}{2}\sum_{i\Omega_{n}}\text{$\log\left[D_{0}(i\Omega_{n})/\beta^{2}\right]$},$

with the bare propagators $G_{0}^{-1}(i\omega_{n})$ and
$D_{0}^{-1}(i\Omega_{n})$ defined after Eq. (193). Using Eq. (193),
$\Omega^{{\rm sc}}/\mathscr{N}=\mathscr{S}^{\text{sp}}/(\beta{\mathscr{N}})$
and Eq. (F.2), we finally find Eq. (73).

##  Appendix G Action for charge fluctuations and superconducting kernel

In this Appendix we analyze phase fluctuations around the solution to the
stationary equations (4). At low energies, a generic solution can be easily
generated from the time and space translation-invariant one
$\left\\{\hat{G}_{{\rm s}}(\tau),\hat{\Sigma}_{{\rm s}}(\tau)\right\\}$
through a $U(1)$ transformation [226]:

$\displaystyle\hat{G}_{\bm{x}}\left(\tau_{1},\tau_{2}\right)$ $\displaystyle=$
$\displaystyle e^{i\varphi_{\bm{x}}(\tau_{1})\hat{\sigma}_{3}}\hat{G}_{{\rm
s}}\left(\tau_{1}-\tau_{2}\right)e^{-i\varphi_{\bm{x}}(\tau_{2})\hat{\sigma}_{3}},$
$\displaystyle\hat{\Sigma}_{\bm{x}}\left(\tau_{1},\tau_{2}\right)$
$\displaystyle=$ $\displaystyle
e^{i\varphi_{\bm{x}}(\tau_{1})\hat{\sigma}_{3}}\hat{\Sigma}_{{\rm
s}}\left(\tau_{1}-\tau_{2}\right)e^{-i\varphi_{\bm{x}}(\tau_{2})\hat{\sigma}_{3}}.$
(204)

Hats denote $2\times 2$ matrices in Nambu space, and $\varphi_{\bm{x}}(\tau)$
is a $U(1)$ phase field conjugated to charge fluctuations of the system [94,
133]. Away from the IR limit, the equations are no longer invariant under the
map Eq. (204) and the $U(1)$ symmetry is broken. The associated Goldstone mode
is precisely given by the field $\varphi_{\bm{x}}(\tau)$, and in the following
we derive the associated low-energy effective action.
Let us insert Eq. (204) into the disorder-averaged SYK action (2), and focus
on the phase degrees of freedom only. The terms which are not invariant under
the $U(1)$ transformation are the trace-log and hopping ones. This leads to
the following expression:

$\displaystyle\frac{\mathscr{S}[\varphi]}{\mathscr{N}}$
$\displaystyle=-\sum_{\bm{x}}{\rm Tr}{\rm
log}\left\\{[(\mu-i\partial_{\tau}\varphi)\hat{\sigma}_{3}-\hat{\sigma}_{0}\partial_{\tau}]\delta(\tau_{1}-\tau_{2})-\hat{\Sigma}_{\bm{x}}(\tau_{1},\tau_{2})\right\\}$
$\displaystyle+$
$\displaystyle\frac{t_{0}^{2}}{2}\sum_{\begin{subarray}{c}\left\langle\bm{x},\bm{x}^{\prime}\right\rangle\end{subarray}}\int_{\tau_{1}\tau_{2}}\left[G_{{\rm
s}}(\tau_{1}-\tau_{2})G_{{\rm
s}}(\tau_{2}-\tau_{1})\,e^{i\left[\partial_{\bm{x}}\varphi(\tau_{1})-\partial_{\bm{x}}\varphi(\tau_{2})\right]\cdot\bm{a}}-F_{{\rm
s}}^{\dagger}(\tau_{1}-\tau_{2})F_{{\rm
s}}(\tau_{2}-\tau_{1})\,e^{-i\left[\partial_{\bm{x}}\varphi(\tau_{1})+\partial_{\bm{x}}\varphi(\tau_{2})\right]\cdot\bm{a}}\right]$

with $\hat{\sigma}_{0}$ and $\hat{\sigma}_{3}$ Pauli matrices. Angled brackets
denote the sum over nearest-neighbor sites on the lattice. Space gradients are
defined via

$\displaystyle\bm{a}\cdot\partial_{\bm{x}}\varphi(\tau)$ $\displaystyle=$
$\displaystyle\varphi_{\bm{x}}(\tau)-\varphi_{\bm{x}^{\prime}}(\tau),$ (206)

where $\bm{a}=\bm{x}-\bm{x}^{\prime}$ gives the spacing between nearest-
neighbor sites, which we have assumed to be the same for any pairs of sites.
Time gradients enter the action by shifting the chemical potential as
$\mu\mapsto\mu-i\partial_{\tau}\varphi$. This is consequence of a gauge
transformation of the fermion fields $\hat{c}_{i\bm{x}\uparrow}(\tau)\mapsto
e^{-i\varphi_{\bm{x}}(\tau)}\hat{c}_{i\bm{x}\uparrow}(\tau),\,\,\hat{c}_{i\bm{x}\downarrow}^{\dagger}(\tau)\mapsto
e^{i\varphi_{\bm{x}}(\tau)}\hat{c}_{i\bm{x}\downarrow}^{\dagger}(\tau)$ that
eliminates phase factors from the self-energy term [226].
Expanding the action Eq. (G) to quadratic order in phase gradients, and
ignoring contributions independent of $\varphi$, yields the Gaussian action:

$\displaystyle\frac{\mathscr{S}[\varphi]}{{\mathscr{N}}}$
$\displaystyle=-\frac{k_{B}T}{a^{3}}\sum_{\bm{p},i\Omega_{n}}\varphi_{-\bm{p}}(i\Omega_{n})$
(207)
$\displaystyle\qquad\qquad\quad\times\left[\Omega^{2}_{n}\bar{\Pi}(i\Omega_{n})+\epsilon(p)\Lambda(i\Omega_{n})\right]\varphi_{\bm{p}}(-i\Omega_{n}),$

which is the lowest-order effective theory for phase fluctuations. Here we
have used the double Fourier transform of the phase field to the space of
momenta $\bm{p}$ and bosonic Matsubara frequencies $\Omega_{n}$:

$\varphi_{\bm{x}}(\tau)=\frac{k_{B}T}{a^{3}}\sum_{\bm{p},i\Omega_{n}}\varphi_{\bm{p}}(i\Omega_{n})e^{i(\bm{p}\cdot\bm{x}-\Omega_{n}\tau)},$
(208)

where $a^{3}$ is the volume of a unit cell, which enters into the definition
of the spatial Fourier transform $\bm{x}\mapsto\bm{p}$. The dispersion
relation of the phase fluctuations in the action (207) is given by
$\epsilon(p)=\frac{1}{2}\sum_{\left\langle
x,x^{\prime}\right\rangle}(\bm{p}\cdot\bm{a})^{2}$. In the case of isotropic
Bravais lattices with spacing $|\bm{a}|=a$, the dispersion is simply
$\epsilon(p)=z\,a^{2}p^{2}/2$, where $p=|\bm{p}|$ and $z$ is the coordination
number. The latter form of the dispersion $\epsilon(p)$ will be assumed from
now on. Here $\bar{\Pi}(i\Omega_{n})$ is the Matsubara transform of
$\bar{\Pi}(\tau)=\Pi_{G}(\tau)-\Pi_{F}(\tau)$, with standard and anomalous
bubbles $\Pi_{G}(\tau)=G(\tau)G(-\tau)$ and
$\Pi_{F}(\tau)=F^{\dagger}(\tau)F(-\tau)$. Moreover,
$2\Lambda(i\Omega_{n})=t_{0}^{2}\left[\Pi_{G}(0)-\Pi_{G}(i\Omega_{n})-\Pi_{F}(0)-\Pi_{F}(i\Omega_{n})\right]$.

In the following we analyze the effects of an external electromagnetic field
and compute the superconducting kernel (58). In order to introduce a vector
potential $\bm{A}$, we perform a Peierls substitution [242, 243] of the
fermions

$\hat{c}_{i\bm{x}\sigma}(\tau)\mapsto\hat{c}_{i\bm{x}\sigma}(\tau)e^{-ie\int_{\bm{x}_{0}}^{\bm{x}}d\bm{y}\cdot\bm{A}(\bm{y},\tau)},$
(209)

with $e$ the fermion electric charge and $\bm{x}_{0}$ a reference location.

Only the non-local-in-$\bm{x}$ hopping term of the SYK action Eq. (2) is
affected by the transformation (209). As a consequence, space gradients in Eq.
(G) are shifted as
$\bm{a}\cdot\partial_{\bm{x}}\varphi(\tau)\mapsto\bm{a}\cdot\partial_{\bm{x}}\varphi(\tau)-e\int_{\bm{x}^{\prime}}^{\bm{x}}d\bm{y}\cdot\bm{A}(\bm{y},\tau)$.
Furthermore, for slowly varying-in-space vector fields we can approximate
$\int_{\bm{x}^{\prime}}^{\bm{x}}d\bm{y}\cdot\bm{A}(\bm{y},\tau)\approx\bm{a}\cdot\bm{A}(\bm{x},\tau)$.
The effective theory in the presence of an external vector potential then
simply follows from Eq. (207) via the minimal substitution
$i\bm{p}\,\varphi_{\bm{p}}(i\Omega_{n})\mapsto
i\bm{p}\,\varphi_{\bm{p}}(i\Omega_{n})-a\,e\,\bm{A}_{\bm{p}}(i\Omega_{n})$.
The change in the action with respect to its unperturbed counterpart (G), due
to the presence of the vector potential, is given by

$\frac{\Delta\mathscr{S}[\varphi,\bm{A}]}{{\mathscr{N}}}=\frac{1}{2}\frac{k_{B}T}{a^{3}}\sum_{\bm{p},i\Omega_{n}}\sum_{\alpha,\beta}A_{\alpha}(\bm{p},-i\Omega_{n})m_{\alpha\beta}(\bm{p},i\Omega_{n})A_{\beta}(-\bm{p},i\Omega_{n})-\frac{k_{B}T}{a^{3}}\sum_{\bm{p},i\Omega_{n}}\sum_{\alpha}\,j_{\alpha}(\bm{p},-i\Omega_{n})A_{\alpha}(-\bm{p},i\Omega_{n}),$
(210)

where $\alpha$ and $\beta$ label the spatial directions. The mass and source
terms are given by

$m_{\alpha\beta}(\bm{p},i\Omega_{n})=-{\mathscr{N}}e^{2}za^{2}\Lambda(i\Omega_{n})\delta_{\alpha\beta}$
(211a)
$j_{\alpha}(\bm{p},-i\Omega_{n})=-i{\mathscr{N}}eza^{2}\Lambda(i\Omega_{n})\varphi_{\bm{p}}(-i\Omega_{n}),$
(211b)

where we have used $\Lambda(-i\omega_{n})=\Lambda(i\omega_{n})$, which enters
into the propagator for phase fluctuations that we defined in Eq. (207). The
partition function, including the effect of the vector potential $\bm{A}$, is
$\mathscr{Z}=\int D\varphi
e^{-(\mathscr{S}[\varphi]+\Delta\mathscr{S}[\varphi,\bm{A}])}$. The functional
derivative of the partition function with respect to the vector potential
gives the total current $\bm{J}_{\text{tot}}(\bm{q},\omega)$ of the system
[158, 226], which comprises the external current as well as the induced
internal current, in accordance with Eq. (58) upon the analytic continuation
$i\Omega_{n}\rightarrow\omega+i0^{+}$ [174, 194]. Formally,
$a^{3}/(k_{B}T)J_{\text{tot},\alpha}(\bm{p},-i\Omega_{n})=\delta\log\mathscr{Z}/\delta
A_{\alpha}(-\bm{p},i\Omega_{n})=\left.\left\langle
j_{\alpha}(\bm{p},-i\Omega_{n})\right\rangle\right|_{\bm{A}}-m_{\alpha\beta}(\bm{p},i\Omega_{n})A_{\alpha}(\bm{p},-i\Omega_{n})$.
Here the brackets
$\left\langle(\cdots)\right\rangle_{\bm{A}}=\mathscr{Z}^{-1}\int
D\varphi(\cdots)e^{-(\mathscr{S}[\varphi]+\Delta\mathscr{S}[\varphi,\bm{A}])}$.
Then, the electromagnetic kernel is defined as the ratio between the total
current $J_{\text{tot},\alpha}(\bm{p},-i\Omega_{n})$ and the imposed vector
potential $A_{\beta}(\bm{p},-i\Omega_{n})$, evaluated in the vanishing vector-
potential limit – see Eq. (58). This additional derivative yields two terms:
one involves the product between two one-point functions $\left\langle
j_{\alpha}(-\bm{p},i\Omega_{n})\right\rangle_{\bm{A}}\left\langle
j_{\alpha}(\bm{p},-i\Omega_{n})\right\rangle_{\bm{A}}$, which is assumed to be
null in the absence of the external source (in the unperturbed system, with
$\bm{A}=0$); the second term is the correlator $\left\langle
j_{\alpha}(\bm{p},-i\Omega_{n})j_{\beta}(-\bm{p},i\Omega_{n})\right\rangle_{\bm{A}}$.
Employing Eq.(211b), we thus obtain for the kernel

$\displaystyle\mathcal{K}_{\alpha\beta}(\bm{p},i\Omega_{n})$
$\displaystyle=\frac{1}{k_{B}T}\frac{\delta^{2}\log{\cal Z}}{\delta
A_{\beta}(\bm{p},-i\Omega_{n})\delta
A_{\alpha}(-\bm{p},i\Omega_{n})}\Big{\rvert}_{\bm{A}=0}$
$\displaystyle=\frac{\mathscr{N}e^{2}\,z\Lambda(i\Omega_{n})}{a}$
$\displaystyle\times\left[\delta_{\alpha\beta}+\frac{za^{2}\Lambda(i\Omega_{n})p_{\alpha}p_{\beta}}{\Omega_{n}^{2}\bar{\Pi}(i\Omega_{n})+za^{2}p^{2}\Lambda(i\Omega_{n})/2}\right].$
(212)

We immediately see that (minus) the static, zero-momentum limit of the kernel
(G) corresponds to the phase stiffness of the system, consistently with Eq.
(63): it is proportional to $-\Lambda(0)$, which can be equivalently obtained
from the phase action (207), written in real space of coordinates, by deriving
the logarithm of $\mathscr{Z}_{\varphi}=\int D\varphi
e^{-\mathscr{S}[\varphi]}$ twice with respect to spatial gradients
$\nabla\varphi(\bm{r},\tau)$ of the phase. Hence, the long-wavelength static
electromagnetic response function probes the rigidity of the system, i.e. the
second derivative of the grand potential, with respect to global phase twists
[191, 192, 193, 158]. As such, the stiffness has units
$\left[\mathscr{E}L^{2-d}\right]$, where $\mathscr{E}$ is an energy, $L$ is a
length, and $d$ is the system dimensionality. We rapidly comment on the
modifications to the various quantities in this section, when we employ SI
units instead: in this case, in the Peierls-substitution formula (209) we have
$ie\mapsto ie/\hbar$, with $\hbar$ reduced Planck’s constant; this constant
propagates to Eqs. (211), where $m_{\alpha\beta}(\bm{p},i\Omega_{n})\mapsto
m_{\alpha\beta}(\bm{p},i\Omega_{n})/\hbar^{2}$ and
$J_{\alpha}(\bm{p},-i\Omega_{n})\mapsto
J_{\alpha}(\bm{p},-i\Omega_{n})/\hbar^{2}$, and finally to the kernel
$\mathcal{K}_{\alpha\beta}(\bm{p},i\Omega_{n})\mapsto\mathcal{K}_{\alpha\beta}(\bm{p},i\Omega_{n})/\hbar^{2}$.
Going to the $i\Omega_{n}=0$ and $q\rightarrow 0$ limit of Eq. (G), and
employing Eq. (63), finally yields Eq. (64) for the superconducting phase
stiffness.

The electromagnetic kernel (G) allows us to calculate the conductivity tensor
of the Yukawa-SYK lattice through [174, 194, 158]

$\sigma_{\alpha\beta}(\bm{p},\omega)=-i\frac{\mathcal{K}_{\alpha\beta}(\bm{p},\omega)}{\omega+i0^{+}},$
(213)

where we have analytically continued the kernel to real frequencies as
$\mathcal{K}_{\alpha\beta}(\bm{p},\omega)=\lim_{i\Omega_{n}\rightarrow\omega+i0^{+}}\mathcal{K}_{\alpha\beta}(\bm{p},i\Omega_{n})$.
Taking the longitudinal part $\mathcal{K}_{L}(\bm{p},\omega)$ of the kernel –
see discussion after Eq. (59) – and going to the $p\rightarrow 0$ limit [226],
we obtain the zero-momentum conductivity
$\sigma(\omega)\equiv-i\lim_{p\rightarrow
0}\mathcal{K}_{L}(\bm{p},\omega)/(\omega+i0^{+})$. At $T>T_{c}$ (normal state)
the conductivity only depends on the retarded normal propagator
$G^{R}(\omega)$ through $\Lambda(\omega)$, while for $T<T_{c}$
(superconducting state) there is an additional contribution to the
conductivity stemming from $F^{R}(\omega)$; see also Appendix B.5. An analysis
of the conductivity in the normal and superconducting states is reported in
the companion paper [159].

## References

  * Zaanen [2004] J. Zaanen, Superconductivity: Why the temperature is high, Nature 430, 512 (2004).
  * Valla  _et al._ [1999] T. Valla, A. Fedorov, P. Johnson, B. Wells, S. Hulbert, Q. Li, G. Gu, and N. Koshizuka, Evidence for quantum critical behavior in the optimally doped cuprate Bi2Sr2CaCu2O8+δ, Science 285, 2110 (1999).
  * Bruin  _et al._ [2013] J. Bruin, H. Sakai, R. Perry, and A. P. Mackenzie, Similarity of scattering rates in metals showing T-linear resistivity, Science 339, 804 (2013).
  * Giraldo-Gallo  _et al._ [2018] P. Giraldo-Gallo, J. Galvis, Z. Stegen, K. A. Modic, F. Balakirev, J. Betts, X. Lian, C. Moir, S. Riggs, J. Wu, A. T. Bollinger, X. He, I. Božović, B. J. Ramshaw, R. D. McDonald, G. S. Boebinger, and A. Shekhter, Scale-invariant magnetoresistance in a cuprate superconductor, Science 361, 479 (2018).
  * Naqib and Islam [2019] S. H. Naqib and R. S. Islam, Possible quantum critical behavior revealed by the critical current density of hole doped high-Tc cuprates in comparison to heavy fermion superconductors, Sci. Rep. 9, 1 (2019).
  * Mandal  _et al._ [2019] P. R. Mandal, T. Sarkar, and R. L. Greene, Anomalous quantum criticality in the electron-doped cuprates, Proc. Natl. Acad. Sci. U.S.A. 116, 5991 (2019).
  * Michon  _et al._ [2019] B. Michon, C. Girod, S. Badoux, J. Kačmarčík, Q. Ma, M. Dragomir, H. Dabkowska, B. Gaulin, J.-S. Zhou, S. Pyon, T. Takayama, H. Takagi, S. Verret, N. Doiron-Leyraud, C. Marcenat, L. Taillefer, and T. Klein, Thermodynamic signatures of quantum criticality in cuprate superconductors, Nature 567, 218 (2019).
  * Legros  _et al._ [2019] A. Legros, S. Benhabib, W. Tabis, F. Laliberté, M. Dion, M. Lizaire, B. Vignolle, D. Vignolles, H. Raffy, Z. Li, P. Auban-Senzier, N. Doiron-Leyraud, P. Fournier, D. Colson, L. Taillefer, and C. Proust, Universal T-linear resistivity and Planckian dissipation in overdoped cuprates, Nat. Phys. 15, 142 (2019).
  * Grissonnanche  _et al._ [2021] G. Grissonnanche, Y. Fang, A. Legros, S. Verret, F. Laliberté, C. Collignon, J. Zhou, D. Graf, P. A. Goddard, L. Taillefer, and B. J. Ramshaw, Linear-in temperature resistivity from an isotropic Planckian scattering rate, Nature 595, 667 (2021).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Dessau  _et al._ [1991] D. S. Dessau, B. O. Wells, Z.-X. Shen, W. E. Spicer, A. J. Arko, R. S. List, D. B. Mitzi, and A. Kapitulnik, Anomalous spectral weight transfer at the superconducting transition of ${\mathrm{Bi}}_{2}$${\mathrm{Sr}}_{2}$${\mathrm{CaCu}}_{2}$${\mathrm{O}}_{8+\mathrm{\delta}}$, Phys. Rev. Lett. 66, 2160 (1991).
  * Campuzano  _et al._ [1996] J. C. Campuzano, H. Ding, M. R. Norman, M. Randeira, A. F. Bellman, T. Yokoya, T. Takahashi, H. Katayama-Yoshida, T. Mochiku, and K. Kadowaki, Direct observation of particle-hole mixing in the superconducting state by angle-resolved photoemission, Phys. Rev. B 53, R14737 (1996).
  * Loeser  _et al._ [1997] A. G. Loeser, Z.-X. Shen, M. C. Schabel, C. Kim, M. Zhang, A. Kapitulnik, and P. Fournier, Temperature and doping dependence of the Bi-Sr-Ca-Cu-O electronic structure and fluctuation effects, Phys. Rev. B 56, 14185 (1997).
  * Shen and Schrieffer [1997] Z.-X. Shen and J. R. Schrieffer, Momentum, temperature, and doping dependence of photoemission lineshape and implications for the nature of the pairing potential in high- ${T}_{c}$ superconducting materials, Phys. Rev. Lett. 78, 1771 (1997).
  * Fedorov  _et al._ [1999] A. V. Fedorov, T. Valla, P. D. Johnson, Q. Li, G. D. Gu, and N. Koshizuka, Temperature dependent photoemission studies of optimally doped ${\mathrm{Bi}}_{2}{\mathrm{Sr}}_{2}{\mathrm{CaCu}}_{2}{\mathrm{O}}_{8}$, Phys. Rev. Lett. 82, 2179 (1999).
  * Kaminski  _et al._ [2016] A. Kaminski, S. Rosenkranz, M. R. Norman, M. Randeria, Z. Z. Li, H. Raffy, and J. C. Campuzano, Destroying coherence in high-temperature superconductors with current flow, Phys. Rev. X 6, 031040 (2016).
  * Orenstein and Millis [2000] J. Orenstein and A. Millis, Advances in the physics of high-temperature superconductivity, Science 288, 468 (2000).
  * Hashimoto  _et al._ [2014] M. Hashimoto, I. M. Vishik, R.-H. He, T. P. Devereaux, and Z.-X. Shen, Energy gaps in high-transition-temperature cuprate superconductors, Nat. Phys. 10, 483 (2014).
  * Emery and Kivelson [1995a] V. J. Emery and S. A. Kivelson, Superconductivity in bad metals, Phys. Rev. Lett. 74, 3253 (1995a).
  * Hartnoll [2014] S. A. Hartnoll, Theory of universal incoherent metallic transport, Nat. Phys. 11, 54 (2014).
  * Zaanen [2019] J. Zaanen, Planckian dissipation, minimal viscosity and the transport in cuprate strange metals, SciPost Phys. 6, 61 (2019).
  * Feng  _et al._ [2000] D. Feng, D. Lu, K. Shen, C. Kim, H. Eisaki, A. Damascelli, R. Yoshizaki, J.-i. Shimoyama, K. Kishio, G. Gu, S. Oh, A. Andrus, O. J., E. J. N., and Z.-X. Shen, Signature of superfluid density in the single-particle excitation spectrum of Bi2Sr2CaCu2O8+δ, Science 289, 277 (2000).
  * Zaanen [2021] J. Zaanen, Lectures on quantum supreme matter, arXiv:2110.00961 (2021).
  * Sachdev [2000a] S. Sachdev, Quantum criticality: competing ground states in low dimensions, Science 288, 475 (2000a).
  * Sachdev [2000b] S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, Cambridge, 2000).
  * Coleman and Schofield [2005] P. Coleman and A. J. Schofield, Quantum criticality, Nature 433, 226 (2005).
  * Hartnoll  _et al._ [2018] S. A. Hartnoll, A. Lucas, and S. Sachdev,  _Holographic quantum matter_ (MIT press, Cambridge, Massachussets, 2018).
  * Hartnoll and Mackenzie [2022] S. A. Hartnoll and A. P. Mackenzie, Colloquium: Planckian dissipation in metals, Rev. Mod. Phys. 94, 041002 (2022).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, Rev. Mod. Phys. 75, 473 (2003).
  * Wang  _et al._ [2004] S.-C. Wang, H.-B. Yang, A. K. P. Sekharan, H. Ding, J. R. Engelbrecht, X. Dai, Z. Wang, A. Kaminski, T. Valla, T. Kidd, A. V. Fedorov, and P. D. Johnson, Quasiparticle line shape of ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$ and its relation to anisotropic transport, Phys. Rev. Lett. 92, 137002 (2004).