Unconventional superconductivity as a quantum Kuramoto synchronization problem in random elasto-nuclear oscillator networks

We formulate the problem of unconventional d − wave superconductivity, with phase fluctuations, pseudogap phenomenon, and local Cooper pairs, in terms of a synchronization problem in random, quantum dissipative, elasto-nuclear oscillator networks. The nodes of the network correspond to localized, collective quadrupolar vibrations of nuclei-like, elastic inhomogeneities embedded in a dissipative medium. Electrons interacting with such vibrations form local Cooper pairs, with a superfluid d − wave pseudogap Δ PG , due to an effective, short range attractive interaction of dx2−y2 character. Phase coherent, bulk superconductivity, with a d − wave gap Δ, is stabilized when the oscillator network is asymptotically entangled in a nearly decoherence-free environment. Phase coherence will in turn be destroyed, at T c , when the thermal noise becomes comparable to the coupling between oscillators, the superfluid density K. The 2Δ/k B T c ratio is a function of Kuramoto’s order parameter, r=1−Kc/K , for the loss of synchronization at K c , and is much larger than the nonuniversal 2Δ PG /k B T * ratio, where T * is the temperature at which Δ PG is completely destroyed by thermal fluctuations. We discuss our findings in connection to the available data for various unconventionally high-temperature superconductors.


Introduction
Quadrupolar vibrations are ubiquitous in nature and often constitute the fundamental normal modes of vibration in a plethora of different physical systems. They arise as the lowest frequency modes of vibration in wineglasses [1], and also as the most relevant tectonic field for density perturbations induced by earthquake ruptures [2]. Accelerating masses moving through spacetime produce ripples that propagate as transverse, quadrupolar gravitational waves [3], and quadrupolar vibrations of the inner crust in a neutron star controls its transient cooling, when coupled to a dissipative, outter thermal bath [4]. At smaller scales, quadrupolar surface vibrations in finite nuclei are known to contribute to the giant quadrupolar resonance [5]. Most importantly, these very same quadrupolar surface vibrations contribute also to a remarkable emergent phenomenon in nuclear matter: nuclear superfluidity [6]. For finite nuclei, this is the mechanism behind the opening of superfluid gaps, Δ, in nuclear spectra, when nucleons minimize their energy, in the presence of a short range attractive nuclear potential, by moving in Cooper paired, time-reversed orbits [7].
The spontaneous emergence of collective behaviour in large oscillator networks is also a phenomenon that has applications in many branches of science. These include the description of the synchronous flash of fireflies [8] and the generation of alpha rhythms in the brain [9]. Huygen's pendulum clocks, weakly coupled through a wooden beam [10], or metronomes sharing a common base [11] are examples of anti-phase and in-phase classical synchronization, respectively. The dynamics of fast spins coupled to slow exchange interactions in XY spin glasses [12] and the frequency locking in superconducting Josephson junction arrays [13] are, on the other hand, examples of emergent collective behaviour in quantum oscillator networks. In all those cases, the model  [19], Copyright (1998) by the American Physical Society.; (f) topographic image of a Bi 2 Sr 2 CaCu 2 O 8+δ crystal, with inset showing the resonance associated to interstitial oxygen, O i , defects, corresponding to a −0.9V resonance in the local density of states; (g) distribution of −0.9 V resonances throughout the sample (map of O i impurities) [18].
According to this scenario the pseudogap phase would then be explained by the so called phase fluctuation superconductivity with pre-formed Cooper pairs, and this is the scenario behind the theory developed in this work.
We merge the ideas of quadrupolar vibrations, synchronization, and phase fluctuations, discussed in the previous three paragraphs, and we formulate a novel approach to unconventional d − wave superconductivity in terms of a Kuramoto synchronization problem in random, quantum dissipative, elasto-nuclear quadrupolar oscillator networks. Our basic requirement is the correspondence between self-organized, elastic inhomogeneities and enhanced local superconductivity, as evidenced in figures 2(a)-(c), obtained both from STM and micro X-Ray spectroscopy (μXRS) spectra for different high-T c cuprates [18,22]. Figure 2(a) shows the inhomogeneous pseudogap intensity heatmap extracted from the STM spectra in underdoped Bi 2 Sr 2 CaCu 2 O 8+δ [18] overlaid by the position map for the distribution of the interstitial oxygen dopants, O i . The large (yellow) and small (white) circles suggest a correlation between O i agglomerates and the amplitude of the pseudogap, since regions with larger groups of dopants (yellow circle) are observed to correspond to regions of larger gap amplitudes (darker regions). This correspondence has in fact been quantified in figure 2(b) [18], in terms of a mathematical correlation between Δ and O i . Parallely, and within a different experimental context, O i dopants have been observed to self-organize into nanosize regions, or puddles, via μXRS in HgBa 2 CuO 4+δ [22], as seen in figure 2(c), as well as in other compounds YBa 2 Cu 3 O 6+δ [23] and La 2 CuO 4+δ [24,25]. Most importantly, it has also been observed that spatial variations in the self-organization of such nanosize, O i − rich puddles have a direct effect on superconductivity, through variations in T c [22][23][24][25]. The O i − rich nanopuddles have different elastic properties than their surroundings, and can therefore be considered as elastic insertions in an otherwise homogeneous medium.
With that in mind, we henceforth shall use elasticity theory to calculate the collective, normal modes of vibration at these elastic insertions (nanopuddles) and show that particles couple to their low-lying quadrupolar mode of vibration. As a result of such particle-vibration coupling, an effective, short range, --d x y 2 2 wave attractive, twoparticle interaction emerges, leading to the formation of time-reversed, local Cooper pairs at each nanopuddle, in close analogy to nuclear superfluidity. Here the pseudogap amplitudes, Δ PG,i , are determined by the elastic properties at the i-th nanopuddle, through their eigenfrequencies, Ω i , while their phases, θ i , are inherited from the phases of the corresponding quadrupolar oscillators. We then calculate the overlap K ij between the quantum wavefunctions for localized, collective quadrupolar vibrations forming a connected, nuclear oscillator network and show that, for K ij stronger than the decoherence noise of the environment, either quenched or thermal, δ, k B T, phase locking occurs, θ i → Θ, and, as a consequence, bulk superconductivity emerges. Finally, we calculate the homogeneous bulk gap, Δ, the transition temperature, T c , and compare their ratio, 2Δ/k B T c , given as a function of Kuramoto's order parameter, r, to the available data for various compounds.

Quadrupolar vibrations
We begin by showing that, already at the classical level, the lowest energetic mode of vibration for an elastic sphere embedded in a infinite homogeneous medium, with different elasto-mechanical properties such as the Lamè parameters, λ and μ, and the density, ρ, is the quadrupole mode. For that we need to solve Navier's  [18]. Reprinted with permission from AAAS. ) -(a) gap inhomogeneities and oxygen interstitial impurities in Bi 2 Sr 2 CaCu 2 O 8+δ as seen by STM, and (b) the study of the correlation between gap intensity and the location of the O i impuritiesthe top curve (light blue) shows a clear and strong correlation between Δ and O i dopants, that is evident from (a) -the yellow (larger) circle contains several impurities and the gap is large, while the white (smaller) circle contains no impurities and the gap is small; (c) (Reprinted (abstract/excerpt/figure) with permission from [22], Copyright (2020) by the American Physical Society.) -self-organized O i − rich nanopuddles as seen by μ x-ray spectroscopy-these are the result of the clustering of O i interstitial dopants forming elastic inhomogeneities such as the ones we shall consider in this paper.
In equation (2) F represents all volume and surface equilibrium external forces that are provided by the crystalline host enclosing the spherical elastic insertion. The stationary, homogeneous, F = 0, solutions to (2) are of the form In what follows we shall be interested in considering radial spheroidal solutions to equation (2) of the kind [27] å å å q q f , with c s and c p being the velocities of the transverse and longitudinal elastic waves, respectively [27]. The coefficients A ℓ,m and B ℓ,m are fixed by the boundary conditions of continuity for the displacement, u, and by the radial, σ rr , and shearing, τ rθ and τ rf , stresses, at the radius of the inhomogeneous insertion, r = R 0 , relevant when F ≠ 0. As usual, are the spherical harmonics, and we have defined the function b ℓ (z) = h ℓ (z), for vibrations of the crystalling host, outside the sphere, and the function b ℓ (z) = j ℓ (z) for vibrations of the embedded elastic inhomogeneity, with j ℓ (z) and h ℓ (z) being the spherical Bessel and Hankel functions, respectively. The use of spherical Bessel and Hankel functions ensures that the displacement, u, vanishes both at the origin and at infinity, which are natural boundary conditions for any finite elastic deformation. Finally, we have also defined d 1 (z) = ℓb ℓ (z) − zb ℓ+1 (z), and d 2 (z) = (ℓ + 1)b ℓ (z) − zb ℓ+1 (z), for compactness.
The normal frequencies, Ω, are found, as usual, from the zeroes of the determinant for the system of coupled equations defined by equation (2), and their real or complex characters are determined by the ratio between the shear moduli of the enclosing medium and the inclusion, p = β 2 η [27]. Here η = ρ e /ρ i is the ratio between the external (e) and internal (i) densities, and β = c se /c si is the ratio between the transverse elastic velocities outside and inside the sphere. For the case of freely vibrating elastic objects (or for a very soft enclosing medium) then p → 0 and all normal frequencies are real, W Î : once excited the insertion is set into permanent, stationary vibration. For radial-spheroidal elastic vibrations embedded in an equally ellastic, homogeneous crystalline host, however, then p ≈ 1 and we end up with the transcendental equation [27] h z c h a z c - where χ = ΩR 0 /c se is a dimensionless frequency written in terms of the transverse elastic velocity in the exterior, c se , while α = c pe /c pi and ζ = c se /c pe are ratios between the transverse and longitudinal elastic velocities outside and inside the sphere. The functions g i (z) = zj ℓ+1 (z)/j ℓ (z), and g e (z) = zh ℓ+1 (z)/h ℓ (z), are also defined inside and outside the sphere, respectively [27]. In this case, equation (3) produces, instead, complex solutions, W Î , with a real part, , that sets the natural frequency of vibration and an imaginary part, , that provides damping, which can be naturally understood as the decay of the localized, normal vibrations into divergent spherical elastic waves of the enclosing host [27]. Finally, we find that the lowest W ¹ e 0 [ ] mode of vibration corresponds to ℓ = 2, localized, collective quadrupolar vibrations, just like for the vibrations at the border of wineglasses [1], at the core of neutron stars [4], or at the surface of nuclei [7]. The monopole, ℓ = 0, is a breathing mode associated to volume changes at the insertion, see bottom of figure 3, while dipole, ℓ = 1, corresponds to translations of the center of mass of the insertion, see bottom of figure 3, and both ℓ = 0, 1 solutions are thus too costly energetically.

The
Once we established that the lowest mode of vibration is the ℓ = 2 quadrupole mode, let us now proceed and narrow the true vibrational ground state down to thed x y 2 2mode, which will be always valid for layered, anisotropic systems such as the high temperature cuprates. For that, we recall that quadrupolar deformations from a spherical equilibrium can be parametrized in terms of collective coordinates α 2,m [28] å a The equilibrium configuration requires that α 2,m = α 2,−m and, since the radius of a sphere is always a real quantity, Î  R , one must also impose that a a = - , eliminating also deformations associated to combinations of + Y 2 2 and -Y 2 2 containing the imaginary unit, i. Altogether, these constraints exclude completely any deformations associated to the three t 2g vibration orbitals constraints leave us then with only two symmetry allowed, real, e g vibration orbitals where s, c stand for sine and cosine, respectively. For isotropic or weakly anisotropic systems the two e g deformations in equation (5) are degenerate (or nearly) and the lowest elastic quadrupolar vibrations should contain admixtures between the two Y 2,0 and Y 2,2c vibration orbitals. This might, perhaps, be relevant to the description of the two-gap structure observed by ARPES in some cuprates [29,30]. For strongly anisotropic, weakly-coupled, layered structured systems, however, such as the majority of the high-temperature cuprates, the degeneracy between the two e g vibration orbitals is lifted because of the large energy cost for planar stretches required by the Y 2,0 oblate and prolate deformations (see figure 4, left). In this case, the true, lowest order, radial-spheroidal vibration mode is the quadrupolar --d x y 2 2 mode, associated to the Y 2,2c vibration orbital (see figure 4, right). It is worth point out that such elasticd x y 2 2mode breaks the C 4 rotational symmetry of the lattice and may be related to the nematicity recently observed in several layered cuprates, in the pseudogap phase [31][32][33].
The classical analysis performed so far, derived from elasticity theory, is valid for describing the quadrupolar vibrations in classical, macroscopic systems such as wineglasses [1] or the normal modes of vibrations of Earth's crust after an earthquake rupture [2]. Before we are able to write down a particle vibration coupling (PVC) Hamiltonian for particles at nanopuddles, however, we must first extended the classical analysis of quadrupolar vibrations to include quantum mechanical systems of different sizes. The phenomenon of nuclear superfluidity [6], which we shall claim to take place also at each of the individual nanopuddles depicted in figure 2(c), occurs in systems as tiny as the interior of atomic nuclei up to systems as gigantic as the core neutron stars [6]. For this purpose, a quantum mechanical Hamiltonian for the two, lowest order quadrupolar deformations Y 2,0 and Y 2,2c can be written in terms of the Bohr shape variables, a b g = cos where B is some inertial parameter for the quadrupolar vibration of normal frequency W = C B. This is a five dimensional problem in which the separable eigenfunctions of ( ), are written not only in terms of the two Bohr shape variables, β and γ, but also in terms of three Euler angles, θ i=1... 3 . In what follows we fix the crystal-to-laboratory reference frames, through a convenient choice for the Euler angles, q, and we find for the radial part with spectrum E n,τ = ÿΩ(N + 5/2), and N = 2n + τ [34]. Here F n,τ is a normalization constant, n is the index of the radial solution, =  s BC 2 1 4 ( ) is the oscillator stiffness, t+ L n 3 2 is a Laguerre polynomial, and τ is the so called seniority [34].
Particles of effective mass, m * , moving in the vicinity of an elastic inhomogeneity are subject to a local potential, U(r), entering the single-particle Schrödinger's equation from which one obtains the single-particle wave functions, j k (r), and energies, ò(k). The linear deformation potential for small displacements, α 2 = α, reads [6]  (5). For layered, anisotropic systems, such as the cuprates, the configuration which minimizes the elastic energy is Y 2,2c (right), and the true ground state for elasto-nuclear quadrupolar vibrations corresponds to the and the matrix element for a process where an electron, at initial state k, scatters off an elastic insertion, setting it into vibration of the Y 2,2c type, and goes into a final state ¢ k can then be written straightforwardly as [6] ò j ) are the single-particle wave functions solutions to equation (8). Here is the reduced matrix element for the quadrupolar deformations, a c 2,2 , which, in second quantized form and in the Heisenberg picture, are given as where b c 2,2 † and b 2,2c are bosonic, phonon creation and anihilation operators for the Y 2,2c quadrupolar vibrational modes with frequencies W = C B, and θ(t) = Ωt + f. We recognize equation (10) as the second quantized, quantum mechanical version of equation (1).
Within the adiabatic, Born-Oppenheimer approximation the particle motion and quadrupolar vibrations occur at very different time scales and to all purposes the phases e ± i θ( t) appearing in equation (10) can be omitted. We are then ready to write down the local PVC Hamiltonian between conduction electrons and the , † are the usual fermionic, creation and anihilation operators for electrons with wavevector k and spin σ, associated to the single-particle wavefunctions j k (r), and with dispersion given by the single-particle energies ò(k).

Local quadrupolar superconductivity
Consider now the problem of a single band of conduction electrons interacting with an isolated, localized, collective quadrupolar vibration, as described by the total Hamiltonian and ξ k = ò(k) − μ, with dispersion relation ò(k) relative to the chemical potential μ. The particle-vibration-coupling in (11) produces, in second order perturbation theory, an effective two-particle interaction The explicit form of such two-particle interaction is obtained through a canonical transformation to eliminate the phonons producing a BCS Hamiltonian [17] å å Since by construction ¢ V k k , is separable we may write [35][36][37] h h 2anisotropy form factor inherited from the Y 2,2c deformation and is a momentum dependent form factor [37,38], obtained from the Fourier transform of a Gogny-type short range interaction [35], in terms of a material dependent parameter b∼R 0 , where R 0 is the typical radius of the elastic insertion. At zero temperature, T = 0, the pseudogap equation reads , where 〈 L 〉 FS stands for angular average over the Fermi surface, and since ÿΩ = ò F we can set w 2 (k) → w 2 (k F ) to arrive at the local (pseudo) gap at T = 0 where λ = N(ò F )V 0 is proportional to the density of states at the Fermi level. At the temperature T * the local pseudogap, Δ PG , is completely destroyed by thermal fluctuations. In order to calculate T * one needs to solve the gap equation (12) for whose solution is, after performing the same averaging procedures as done for Δ PG , given by We are now ready to calculate the pseudo-gap-to-T * ratio which is a nonuniversal ratio, due to the factor w(k F ), larger than BCS's 3.53, and can, in fact, be as large as 8, as observed experimentally in [20].

The quadrupolar oscillator network
Consider now a large number, N, of mutually interacting quadrupolar oscillators with distributed natural frequencies as shown in figure 5. We want to establish the conditions for the spontaneous emergence of synchronization, or phase locking, in our quadrupolar oscillator network, as well as its consequences to unconventional, bulk superconductivity. For that purpose, we write down the coupled Hamiltonian [39]. Here W = C B i i i are the quadrupolar natural frequencies for the i = 1...N elastic insertions, distributed according to a Lorentzian spectral distribution, that is symmetric with respect to some average frequency, W, with a quenched spread, δ, determining the amount of disorder in the material, while controls the asymptotic hybridization between oscillator wave functions (7) at positions r i and r j through the stiffness s. The model that has become the simplest paradigm for the synchronization phenomenon is the Kuramoto model [14,15]. In terms of the Hamiltonian (18), the spectral distribution (19), and the overlaps (20), the where θ i is the phase of the i − th oscillator, and ζ i are Gaussian noises, associated, via fluctuation-dissipation theorem, to the damping, , that controls the decay of the localized quadrupolar vibrations into thermal modes of the bath [40]. For N independent oscillators, K ij = 0, ∀i, j, the solution to Kuramoto's model in equation (21) is simply where f i is an arbitrary phase constant. In this case, the second-quantized displacements at the independent quadrupolar oscillators are written in Heisenberg's representation as Second order perturbation theory then gives rise to attractive, two-particle interactions at each independent oscillator, and to local superconductivity, as described earlier, being, however, incoherent and unable to produce bulk superconductivity. On the other hand, when K ij ≠ 0, Kuramoto's order parameter, defined as [15] quantifies the overall entanglement of the oscillator network, in terms of the relative strengths between fluctuations and hibridization. If the spread of g(w) is larger than K ij , which is the case of dilute, distant elastic insertions, s 2 |r i − r j | 2 ? 1, then all oscillators perform individual cycles, all phases θ i (t) are uniformly distributed over a circle [0, 2π], and r = 0: synchronization is not possible. In figure 5 (left) we show a dilute, random collection of quadrupolar oscillators whose overlap is not enough to overcome the spread of their distributed, individual natural frequencies. As a result, although at each elastic inhomogeneity attractive two-particle interactions emerge, their local Cooper pairs do not exhibit phase coherence and bulk superconductivity is not possible. For stronger overlaps, K ij , however, when elastic insertions are more abundant and closer together, s 2 |r i − r j | 2 = 1, larger number of oscillators start to group together into clusters of a definite phase, Θ, and then r ≠ 0: partial or full synchronization is achieved. In figure 5 (right) we show a larger, random collection of quadrupolar oscillators whose overlap is now enough to overcome the spread of the distributed, individual natural frequencies. In this case the attractive two-particle interactions at all elastic inhomogeneities give rise to coherent Cooper pairing and to bulk superconductivity. The surge of bulk superconductivity in random oscillator networks, as N → ∞ , can thus be described by Kuramoto's order parameter, r, which, by using a mean field approximation, K ij = K, ∀i, j, is found as [40]  The critical temperature, in turn, is obtained as solution to the gap equation (12) and the BCS interaction . Nevertheless, the ratio between Δ and T c can be calculated for 2γk B T c = r 2 (δ, and is exponentially larger than 2Δ PG /k B T * because of dissipation, γ ≠ 0, associated to the decay of the collective quadrupolar vibrations into thermal modes of the bath

Comparison to experiments
We are now ready to make a direct connection between the results obtained in the previous sections to the experimental situation discussed in the introduction. In particular, the transition between a severely inhomogeneous pseudogap phase towards a somewhat more uniform superconducting phase where no abrupt changes in the gap amplitude are observed but phases get synchronized as the superfluid density outgrows the quenched or thermal noises.
• Before phase coherence all phases θ i (t) evolve independently and, as we have seen from equation (15), the amplitude of the superfluid pseudogap at each nanopuddle, Δ PG,i , is proportional to its normal mode of vibration, Ω i , which is, in turn, determined by the elastic parameters η i , ζ i , α i of each individual nanopuddle through equation (3). Furthermore, the momentum dependence of the superfuid pseudogap in equation (13) has a Gogny factor w(k) shown in equation (14) that ensures its locality in the sample. Therefore, for a random distribution of nanopuddles as in figure 2(c) one expects a distribution of eigenfrequencies Ω i which produces a severely inhomogeneous landscape for the pseudogap as shown in figure 1(c). • Upon phase coherence, however, all phases become locked θ i (t) → Θ(t) and, as we have seen from equation (27), the amplitude of the bulk superconducting gap, Δ, is proportional to the collective, bonding normal mode of vibration, Ω σ , obtained after the globally connected quadrupolar oscillator network synchronizes and the entire system behaves as a single oscillator. Notice furthermore that now no locality form factors are present in equation (26) showing that the bulk gap is uniform throghout the sample, once full synchronization is achieved and r → 1, thus producing the nearly uniform landscape for the bulk gap as shown in figure 1(a).
• As for the gap-to-T c ratios before and after synchronization, in figure 7 we show the behaviour of (17) and (29) as functions of the superfluid density K. The Δ PG -to-T * ratio does not depend on K and, apart from a nonuniversal factor, w(k F ), it is given by h á ñ » = 3.53 3.53 2 4.99 FS 2 , which is pretty close to the expected mean field value for d − wave BCS superconductors. The Δ-to-T c ratio, in turn, depends exponentially on K. Even if in cuprates the Uemura plot suggests that the ratio appearing as r 4 (δ, 0) in the denominator of the exponential. For this reason the Δ-to-T c ratio is considerably larger at low K (within the underdoped regime) while it saturates at large K (within the overdoped regime). Rough estimates for the maximum values for T c , assuming full synchronization, r → 1, are given in table 1 for two layered, high-temperature superconducting cuprates. The fairly large values for T c obtained in our calculations can be traced back to the high normal frequencies of the collective quadrupolar vibrations, Ω σ , that correspond to the roots to equation (3), and lay far above Debye's.

Conclusion
Localized vibrations in cuprates have been observed with Raman spectroscopy [45], inelastic neutron scattering [46], and ultrafast optical coherent spectroscopy [47], but their role to the mechanism of superconductivity remains controversial. STM spectroscopy, on the other hand, show unambiguously that local Cooper pairs form at elastic inhomogeneities, as shown in figure 1 [18]. In this work, we introduce a novel formulation of unconventional, d − wave superconductivity in which collective quadrupolar vibrations, localized at random Figure 7. Comparison of the gap-to-T c ratios for different regimes. The 2Δ/k B T c ratio (top, left, blue) is exponentially large in the underdoped regime, small superfluid density, K. The reason is because while the amplitude of the bulk gap, Δ, remains robust (in the numerator) the transition temperature, T c , becomes vanishingly small at the synchronization transition r → 0 for K → K c . The pseudogap amplitude, Δ PG , (center, right, red) however, as well as its related temperature scale, T * , are independent on the superfluid density, K, and their ratio 2Δ PG w(k F )/k B T * is a nonuniversal constant. For the sake of comparison we also show BCS's universal ratio (bottom, dashed, black). The decay of vibrations into bath modes results in 2Δ/k B T c ? 2Δ PG w(k F )/k B T * (blue area). Reproduced with permission from [43]. © [1998] The Physical Society of Japan. Table 1. ζ values extracted from [44]. Rough estimates for T c from (28), assuming full synchronization, r → 1, nanometer-size insertions, and the weak-coupling parameter, l = 0.8. elastic inhomogeneities, are able to capture several key features of high-T c cuprates. Using concepts from elasticity theory, nuclear superfluidity, and the emergence of collective behaviour in oscillator networks, we propose that the phenomena of phase fluctuations, pseudogap formation, local --d x y 2 2 wave Cooper pairing, and phase coherent, bulk high-temperature superconductivity arise naturally from a Kuramoto synchronization problem in nuclear oscillator networks. In a related but different context, a quantum version of the Kuramoto model has been recently used to study the emergence of coexisting superconducting off diaginal long range order and a pseudogap phase within the Sachdev-Ye-Kitaev model with an attractive Hubbard interaction [48], in which the pseudogap phase is dominated by quantum fluctuations of local phases. As such, the work reported in [48] is closely related to the problem of frequency locking in Josephson arrays with the Kuramoto model [13]. Although at the first sight the ideas explored in the present work can be considered as similar to the problem of percolation between superconducting islands with the overlap K ij playing the role of Josephson couplings, the present approach incorporates also an extra ingredient: the close correspondence between disordered elastic inhomogeneities and fluctuating (both amplitude and phase) local superconductivity, thus leading to highly a nonlinear percolation problem as recently acknowledged in [49]. Finally and furthermore, we especulate that the quadrupolar nature of the node oscillators may be associated to the nematicity observed in the pseudogap phase [31][32][33], while its double degeneracy, with admixture betweend x y 2 2(Y 2,2c ) and d z 2 (Y 2,0 ) vibrational orbitals for weakly anisotropic systems, may be relevant to the understanding of the two-gap structure also observed in the pseudogap phase [29] that can be explored to further enhance the superconducting transition temperature [30]. We hope our formulation of unconventional superconductivityn may help to shed some light into the remarkable emergent phenomenon of unconventional high-temperature superconductivity in cuprates and related inhomogeneous systems.