SYK Superconductivity: Quantum Kuramoto and Generalized Richardson Models

Sachdev-Ye-Kitaev (SYK) model has emerged as a new paradigm of the non-Fermi-liquid behavior. Here we investigate a possibility of having a superconducting off-diagonal long-range order (ODLRO) and a pseudogap phase within the SYK framework. We found that ODLRO may be established in spin-1/2 version of the model with the time-reversal invariance and an extra attractive interaction. If the latter is taken as the on-site negative-$U$ Hubbard term, it leads to the pseudogap phase at $U<U_c$ dominated by quantum fluctuations of local phases. These fluctuations are described by a quantum version of the Kuramoto model, traditionally employed to illustrate synchronization of classical non-linear oscillators. In the opposite limit of large $U$, the SYK+Hubbard model is approaching a certain generalization of the integrable Richardson model. We present exact diagonalization studies, along with analytic solutions of the aforementioned limiting cases. We also discuss possible holographic interpretations of the model, ODLRO and the pseudogap.

Although the initial SYK model is zero-dimensional (0D) with all-to-all random interactions, it was soon generalized to include D-dimensional arrays of connected SYK grains [3,7,22,[29][30][31][32][33][34][35][36].Such models were shown to exhibit T -linear resistivity, making them attractive candidates for description of strongly correlated materials [37].An account of quantum fluctuations in such arrays reveals [36] a quantum phase transition (QPT) between a gapless (thermal) insulator and the Fermi liquid at certain critical inter-grain coupling.In these picture the T -linear metallic phase appears as the quantum critical region [38] of the aforementioned QPT.
Success of the SYK model in describing the non Fermi liquid state raises the question if superconductivity may be included in the same framework.A number of models were suggested with this goal in mind both in 0D [39,40] and in the array [41,42] context.All of them found that the original SYK model must be upgraded to complex spin-full fermions with an extra mechanism of attraction, such as phonons [39,40], pair hopping [41], or special correlations between matrix elements [42].Such upgraded SYK-like models indeed exhibit superconducting correlations, which may be treated within the large N mean-field approach.Similarly to the Fermi liquid BCS mechanism, an infinitesimal attraction is sufficient to develop the superconductivity.
In this paper we consider a different 0D model, where the attraction is provided by a negative U Hubbard term.Contrary to the mechanisms mentioned above, the mean-field treatment completely fails to describe the SYK+Hubbard model even in the N → ∞ limit.This is because the Hubbard term does not inhibit on-site phase fluctuations, which invalidate the mean-field approach.Such quantum phase fluctuations result in an insulating pseudogap phase at U < U c , where U c is a critical attraction strength.For U > U c there is a superconducting "dome" on the U vs. temperature plane.The superconducting phase under the dome is also strongly affected by the quantum fluctuations and does not conform to the mean-filed description.
In view of the failure of the mean-field, one needs to develop alternative means, capable of treating strongly fluctuating phases.Fortunately, within the 0D framework this can be achieved.In the limit of large U we found that the model may be mapped onto a certain generalization of the exactly solvable Richardson model [43][44][45].It's solution reveals a superconducting low temperature state with the first order transition to the normal non Fermi liquid state at T c ∝ U −1 .The first order transition between a superconductor and a non Fermi liquid has been already noticed in Refs.[41,46].It's possible that SYK+Hubbard and the associated Richardson models provide the simplest cartoon for this phenomenon.
In the opposite limit of a weak attraction, the phase fluctuations may be described by an effective model, which we call the quantum Kuramoto model.The classical Kuramoto model is a paradigm for synchronization of non-linear stochastic oscillators [47][48][49][50][51][52][53][54][55][56][57].It's quantum counterpart provides a description of a continuous QPT between the pseudogap state with unsynchronized phases and the phase-coherent superconductor.We found it remarkable that the SYK framework is capable to exhibit the pseudogap physics.
To verify validity of this theory we resort to an exact diagonalization of spin-1/2 SYK+Hubbard model.To detect superconductivity numerically in a finite size sys-tem, we employ the notion of the off-diagonal long-range order (ODLRO) [58,59].It allows for a sharp definition of the condensate fraction and its dependence on temperature and the attraction strength for a large, but finite N (number of sites).Numerical results are in a qualitative (and in cases where numerical coefficients may be evaluated, a quantitative) agreement with the theory.
The paper is organized in the following way.In section II we discuss the models and the notations.Section III is devoted to the notion of ODLRO.It is followed by section IV, where we outline mean-field treatment and expectations for the models at hand, the latter are then compared with the results of the exact diagonalization in section V.In section VI we explain numerical observations by mapping onto Richardson and quantum Kuramoto models in the regimes of strong and weak attraction correspondingly.In section VII we discuss a possible holographic interpretation of the fluctuation-dominated SYK superconductivity in terms of the "bulk" description.Finally, section VIII briefly summarizes our findings and lists some open problems.

II. NOTATIONS AND MODELS
We consider 0D models, consisting of N 1 orbitals (or sites), labeled as i, j, . . .= 1, 2, . . ., N .Each orbital may be occupied by a complex spin-1/2 fermion annihilated with the operator c iσ , where σ =↓, ↑ is the spin index.In the spirit of the SYK model, we assume that all orbitals are exactly degenerate with the on-site energy taken to be zero.The orbitals interact through the four-fermion interaction with real spin-independent matrix elements.These interactions are summarized by the SYK part of the Hamiltonian: where J ij;kl is a real tensor with the following symmetry properties: We also demand that non-zero elements must have all four indexes i, j, k, l distinct.Up to these symmetries, the matrix elements J ij;kl are assumed to be real independent random variables, drawn from the Gaussian distribution with the zero mean, J ij;kl = 0, and the variance We'll show below (both numerically and analytically) that the pure SYK Hamiltonian (1) does not lead to ODLRO [60].For ODLRO to develop, one needs to supplement SYK Hamiltonian with an attractive term, facilitating fermion pairing.One possibility is a site-local negative U Hubbard term: Another option is all-to-all pair hopping [41]: which annihilates a pair at an orbital j and creates at, in general, different orbital i.Both Hamiltonians contain a chemical potential to adjust the occupation fraction.The three Hamiltonians, written above, conserve particle number and are symmetric under the time-reversal transformations.States of these models are governed by temperature, T , fermion occupation number, N f , and the dimensionless parameter, U/J, characterizing the attraction strength.
In the absence of the SYK term the ground state of the pure Hubbard model, Eq. ( 4), consists of localized pairs and does not exhibit ODLRO.Its energy is obviously −U per fermion pair and its degeneracy is given by the number of combinatorial possibilities of distributing a given number of pairs among N orbitals.Excited states are formed by breaking some of the pairs and creating single occupied orbitals with zero energy.As we show below, ODLRO may be established, mediated by the SYK interactions.
The pure pair-hopping Hamiltonian, Eq. ( 5), is somewhat different.It constitutes a limiting case of the Richardson model [43][44][45] (see section VI A and Appendix C for details).The latter predicts a nondegenerate ground state with ODLRO separated by the gap, ∝ U , from the first excited state, which is (N − 1)fold degenerate.We'll show that SYK interactions do not destroy ODLRO, but weaken it substantially if J > U .
Numerically we first block diagonalize the 2 2N × 2 2N matrix Hamiltonian in the many-body space, using particle number conservation and other symmetries (e.g.particle-hole symmetry for the half-filled case).We then exactly diagonalize the relevant blocks to extract their spectrum and eigenfunctions.

III. THE OFF-DIAGONAL LONG-RANGE ORDER
The standard definition of the superconductivity implies a finite anomalous expectation value, ∆i ∝ c † i↑ c † i↓ .It is clear however, that for a finite size system with a particle conserving Hamiltonian such expectation value is bound to vanish.One thus needs another measure of the superconducting order.The corresponding concept of ODLRO is well known from, e.g., the theory of cold atom Bose condensates in optical or magnetic traps [59].
Let us define the bosonic pair creation operator as Since there can't be more than one such boson per orbital, we are dealing with the hard-core bosonic particles.One then defines the reduced single-particle bosonic density matrix as where . . .implies the exact many-body ground state (or thermal) expectation value.Defined this way, ρ ij , is an N × N positive-definite matrix.Its trace is a total number of local pairs, which is less or equal than N f /2 (we typically consider half-filled systems with N f = N ).
One is interested in the spectrum of eigenvalues of ρ ij : λ α , where The absence of the pair condensate corresponds to all N eigenvalues λ α being of order one, O(1).On the other hand, the pair condensate corresponds to the largest eigenvalue λ 0 being O(N ), while the remaining N − 1 eigenvalues being O(1). Figure 1 shows T = 0 spectrum of ρ ij for SYK +Hubbard model with U/J = 2 for N = N f = 8, 10, 12.One can clearly see the largest eigenvalue splits from the rest and approaches N/4 + 1/2.The remaining eigenvalues coalesce towards ≈ 1/4.This behavior may be understood with the help of the generalized Richardson model, as explained in Section VI A. The presence of the single eigenvalue with the O(N ) scaling is the hallmark of ODLRO [59].Indeed, admitting a nonzero anomalous average ∆i ∝ c † i↑ c † i↓ , one finds ρ ij ∝ ∆i ∆ j .This is the rank-1 matrix with the single non-zero eigenvalue, given by its trace (∝ N ).
Figure 2 shows temperature dependence of the condensate density, (λ 0 − 1/2)/N , (subtraction of 1/2 is motivated by the expectation that, in the absence of ODLRO, all λ's approach 1/2).One notices the approximate crossing point at T c ≈ 0.1W , where W is the energy scale of the Richardson model, W = 3J 2 /32U , (see Eq. (12)  Section VI A).Such crossing point indicates a phase transition in the N → ∞ limit between phases with a finite and zero condensate density.

IV. MEAN-FIELD TREATMENT
To develop a large N mean-field treatment, one follows the standard root [2,61] of averaging over the random SYK matrix elements J ij;kl and deriving the so-called GΣ action.There is a peculiarity, though, associated with the matrix elements being real.It is coming from the fact that there are two distinct terms in the square bracket on the right hand side of Eq. (1), see Appendix A. Upon averaging over the Gaussian distribution of J ij;kl , one obtains two types of terms which are expressed through the normal and anomalous two-point fields: The normal component is spin-diagonal and independent of the spin-projection.Here we have suppressed replica indexes for brevity.The normal and anomalous components may be combined in the Nambu matrix field Ĝτ,τ .The definitions (8) are enforced by conjugate non-local fields, which may be also combined into the Nambu space matrix Στ,τ , playing the role of the self-energy.
The Hubbard term, Eq. ( 4), may be decoupled in the Cooper channel with the help of the local fields ∆ i (τ ), leading to the effective action of the form: where ∆i = ∆ i σ + + ∆i σ − is the off-diagonal Nambu matrix.For the pair-hopping model, Eq. ( 5), one needs a single field ∆(τ ) to decouple it.One thus arrives at the same action (9) with the constraint ∆ i = ∆.In the latter case there is a large factor N in front of the entire action, justifying the mean-field saddle point approximation.
The mean-field equations, obtained upon variation of the action over the matrix fields Ĝ, Σ as well as over ∆ are specified in Appendix A. Their numerical analysis [60] shows that in the absence of attraction (U = 0 and thus ∆ i = 0) the lowest free energy solution is purely normal, i.e.F τ,τ = 0, while Ĝτ,τ ∝ |τ − τ | −1/2 , same as in conventional complex-J SYK model.
One can investigate now stability of such nonsuperconducting SYK solution against a small attractive U perturbation.The corresponding self-consistency equation for ∆ takes the form U −1 = C(∆), where the Cooper channel polarization C = dτ G 2 τ .In the normal phase of SYK, G τ ∝ (Jτ ) −1/2 , and therefore C is given by the logarithmic integral.In the IR limit the latter is cut by either temperature or |∆|, leading to U −1 ∝ J −1 ln(J/|∆|) and thus |∆| ∼ Je −const•J/U for U J. Thus the mean-field treatment predicts that, similarly to BCS case, an arbitrarily weak attraction results in a finite superconducting order parameter, albeit an exponentially small one.
A detailed calculation, presented in Appendix A, leads to the following mean-field solution for the absolute value of the order parameter It is worth mentioning that the energy gap in the manybody spectrum scales as |∆| 2 /J for U J and as |∆| for U J, Appendix A. As mentioned above, one expects the mean-field treatment to be accurate for the SYK+pair hopping model in N → ∞ limit.It is not clear a priory if SYK+Hubbard is also accurately described by this theory.Indeed, in the latter case the order parameters, ∆ i , on individual orbitals fluctuate independently (first line in Eq. ( 9)) and such fluctuations are not necessarily decreasing as N → ∞.To check this we perform finite-size exact diagonalization study, summarized below.

V. EXACT DIAGONALIZATION
Figure 3 shows the exact diagonalization results for the SYK+pair hopping Hamiltonian, Eqs.(1), (5), for the half-filled N = 12 case -the largest size accessible in our simulations.The top panel shows ODLRO, defined as the difference between the largest and the second largest eigenvalues of ρ ij , Eq. ( 7), as a function of U/J.The bottom panel shows the gap in the many-body spectrum, defined as the difference between the energies of the first excited and the ground-state, also as a function of U/J.At U J the ODLRO saturates to N/4, while the many-body gap approaches U -in agreement with the mean-field.Due to finite size effects, it is hard to draw definitive conclusions about small U behavior.Qualitatively it is also consistent with the mean-field expectations, Eq. ( 10).This behavior should be contrasted with the results of the exact diagonalization of the SYK+Hubbard, Eqs. ( 1), (4), presented in Fig. 4. One notices a critical value U c ≈ 0.24J, below which there is no any evidence of neither ODLRO, nor the many-body gap (beyond a finite-size effect of the SYK model).As indicated in the inset, U c does not decrease with increasing N and thus it's unlikely to be a finite-size artifact.Another marked difference is the behavior of the many-body gap at large U .Unlike the pair-hopping model, where the many-body gap increases with U , the Hubbard model exhibits a non-monotonic dependence of the gap with U , with the maximum gap reached at U ≈ 0.4J.The finite-temperature behavior of the SYK+Hubbard model is illustrated in Fig. 5, where we present the color plot of the logarithm of ODLRO on the temperature vs. U/J plane.Once again, one notices the non-monotonic behavior of the critical temperature, where ODLRO is suppressed.
The presence of the critical interaction strength, U c , and the non-monotonic behavior of the gap and T c are contrary to the mean-field predictions, Eq. ( 10).We attribute both phenomena to the strong quantum fluctuations in the SYK+Hubbard model.To account for such large N , non-mean-field phenomenology, we investigate below the SYK+Hubbard model in the two limiting cases of strong and weak attraction.In both cases we are able to account for the quantum fluctuations and show that they indeed explain the observed behavior.
In the case of the strong attraction this is achieved by mapping onto an exactly solvable generalized Richardson model.It provides an asymptotically exact description of the low-energy part of the SYK+Hubbard model in the limit U √ N J/13.In the opposite limit of the weak attraction we reduce the problem to the quantum version of the Kuramoto model.It's classical counterpart [47][48][49][50][51][52][53][54][55][56][57] provides a paradigm for synchronization of nonlinear oscillators.We show that the quantum Kuramoto model provides description of the pseudogap phase for U < U c and the continuous superconducting QPT at U = U c .

VI. QUANTUM FLUCTUATIONS IN SYK +HUBBARD MODEL A. Generalized Richardson Model
The many-body spectrum of the SYK+Hubbard model with U = 2J and N = 8 is shown in Fig. 6 as a function of the fermion number, N f .One notices strong alternation of the ground state energies between even and odd fermion number.The low-energy part of the spectrum, which is not resolved in the main plot is shown in the inset for even N f .These low-energy bands are separated by the gap ∼ U from the rest of the spectrum.Number of many-body states in these low-energy bands is exactly N f /2 N , i.e. the number of ways to place N f /2 indistinguishable pairs over N orbitals.Therefore the low-energy bands are described by models of hard-core bosons, Eq. ( 6).In the absence of the SYK term, bosons are localized and all N f /2 N bosonic states are degenerate with the energy −U per boson.The SYK term induces an effective bosonic hopping and thus leads to a formation of the low-energy bands.
To gain an insight in the physics of the corresponding bosonic model, consider a state with N f /2 hard-core bosons occupying a subset of N orbitals.Acting with a given term of the SYK Hamiltonian, (1), say J ij;kl , on such a state produces a non-zero result only if orbitals k and l are occupied, while i and j are empty (or vise versa).It leads to a state with N f /2−2 bosons and 2 broken pairs (i.e. 4 unpaired fermions on orbitals i, j, k, l).Such a state costs energy 2U and resides outside of the low-energy bosonic sector.From the point of view of an effective bosonic model, it is a virtual state, which ought to be integrated out.To bring the system back to the bosonic sector one has to act on it with the same SYK term, J ij;kl .This either brings the system back to the initial state (generating an uninteresting on-site energy shift), or results into hopping of two bosons from the orbitals k, l to i, j.The latter option gives rise to the effective bosonic Hamiltonian: where the factor of 6 = 2 + 4 is coming from the opposite and same spin terms in the SYK Hamiltonian, correspondingly.There is also a one boson hopping term of the form jk M jk b † j b k , where M jk ∝ − N il J ij;kl J lj;ki /U .Since the two matrix elements here are uncorrelated, the corresponding sum includes N 2 sign alternating terms, implying for a typical matrix element ).This makes one boson hopping insignificant at large N .
Hamiltonian (11) represents a version of the bosonic SYK model [62,63].Specifics of our model is that we work with real matrix elements J ij;kl and thus there is a non-random sign-definite part of the Hamiltonian (11), which we call a generalized Richardson model: where W = 3J 2 /32U and all indexes i, j, k, l must be distinct.We introduced operator B 0 = These operators form the su(2) algebra upon identification L+ = B † 0 , L− = B 0 , Lz = Nb − N/2.One thus finds that: B † 0 B 0 = L2 − L2 z + Lz .This observation allows one to solve the Richardson Model [43][44][45] with degenerate on-site energies, H R = − W N B † 0 B 0 .Let us focus for simplicity on the half-filled model, with N b = N/2 and thus L z = 0.The spectrum of the half-filled Richardson model is thus given by E R (L) = −W L(L+1)/N , where the total angular momentum runs L = 0, 1, . . ., N/2.The unique ground state corresponds to L = N/2.The degeneracies of the excited states are given by the multiplicity of the corresponding representations: with the total number of states: N , which is the Hilbert space dimensionality for the halffilled hard-core particles.
In the same way one finds the spectrum of the halffilled generalized Richardson model, Eq. ( 12), to be: with the same set of degeneracies, Eq. ( 14).The manybody gap between the ground state, L = N/2, and the first excited band with L = N/2 − 1 and degeneracy D(N/2 − 1) = N − 1 is approaching W/2 at large N .The ground state is |GS ∝ (B † 0 ) N/2 |0 .The corresponding single-particle density matrix ρ ij , Eq. ( 7), has diagonal elements ρ ii = 1/2 and off-diagonal ones Thus its largest eigenvalue is λ 0 = N/4 + 1/2 (dashed line in Fig. 1).The fact that it scales as N signals the presence of ODLRO in the ground state of the generalized Richardson model.The remaining N − 1 eigenvalues are degenerate at λ α = 1 4 N N −1 .These features are qualitatively consistent with the exact diagonalization results of SYK+Hubbard shown in Fig. 1 for To access distraction of ODLRO at elevated temperature one considers the partition function: where we introduced l = L/N , substituting summation with the integration, and the free energy density, Fig. 7, is defined as where γ = 4 for the generalized Richardson model.
In the large N limit, the integral in Eq. ( 16) is dominated by the minima of f (l).The latter changes from being l = 1/2 at T = 0 to l = 0 at T c /W = 1/16 ln 2 ≈ 0.09, where the model undergoes the first order transition to a state with no ODLRO.This behavior is illustrated in Fig. 8, which shows results of the exact diagonalization for the generalized bosonic Richardson model, Eq. ( 12).The crossing point at T /W ≈ 0.1 marks the first order transition, were ODLRO jumps from 1/4 to zero in the N → ∞ limit.This should be compared with the exact diagonalization of the SYK+Hubbard model shown in Fig. 2.
It is instructive to compare this behavior with that of the traditional Richardson model, H R = − W N B † 0 B 0 , whose partition function is again given by Eqs. ( 16), (17) with γ = 2.The latter model may be seen to undergo a continuous phase transition at T c = W/2.This model with W = U is exactly the pure pair hopping model, Eq. ( 5).
One may worry if the generalized Richardson model, Eq. ( 12), is a reasonable approximation for the lowenergy bosonic model (11).To answer this question one needs to examine the role of the random part of J 2 ij;kl in Eq. ( 11).This random part removes degeneracies, Eq. ( 14), between excited states with L < N/2, transforming them into the bands.Let's focus on the lowest such band with L = N/2 − 1, consisting of N − 1 states.One can write an effective model for this band as (N − 1) × (N − 1) matrix Hamiltonian with the random elements h rs .Their variance can be estimated from the fact that a matrix element h rs is given by a sum of N 4 random sign terms each of the order W/N 3 .As a result, The density of states of such random matrix is given by a semicircle with the bandwidth √ N W/N = W/ √ N .Since the gap between the band and the ground-state scales as W , the latter remains well separated as long as N 1 even for the random model, Eq. ( 11), see Fig. 9.
We thus conclude that the generalized Richardson model, Eqs. ( 12)-( 17), provides an accurate description of the low-energy sector of the SYK+Hubbard model for U J. It predicts ODLRO at low temperature.The many-body gap and critical temperature both scale as J 2 /U with the large ratio between the two, 8 ln 2 ≈ 5.55 (cf.with the BCS gap to T c ratio of 3.53).An enhancement of this ratio is also known in the context of quantum critical models [64], holographic superconductors [65] and other SYK-like models [39,41].These features are qualitatively consistent with the exact diagonalization results for the moderate N SYK+Hubbard model.The singleparticle fermionic excitations are separated by a larger gap ∼ U .It is important to notice that the full bandwidth of the bosonic states is N W/16 = 3N J 2 /512U .The requirement for the Richardson model to be quantitatively accurate is U > 3N J 2 /512U , i.e.U √ N J/13.This condition is satisfied for Figs. 1, 2 and 6.

B. Pseudogap and the Quantum Kuramoto Model
We turn now to the opposite limit of U J, where there is no separation between bosonic and fermionic sectors.To describe this limit, we notice that the action (9) exhibits a non-trivial saddle point with 2U ) , Eq. (10).However, the phases, φ i , of the local order parameters, ∆ i = |∆|e iφi , are not fixed by the saddle point equations.They constitute thus the soft degrees of freedom, which are (almost) free to fluctuate.Such fluctuations are capable of destroying ODLRO, despite presence of the non-zero |∆|, even in the N → ∞ limit.
The action which governs the low-energy dynamics of the local phases is given by: The two constants here, m and g, are both related to the thermodynamic susceptibilities of the model.In principal, they are site specific, m i and g ij , however in the large N limit they may be substituted by the corresponding ensemble averages: m = m i and g = g ij .The local compressibility m i = −∂ 2 E GS /∂µ 2 i is the susceptibility of the ground-state energy, E GS , to a local chemical potential, entering the Hamiltonian as −µ i c † iσ c iσ .In the |∆| = 0 case it was evaluated in Ref. [66] and found to be m ≈ 1.04/J.We do no expect it to be significantly affected by the presence of small |∆|.The off-diagonal Cooper susceptibility g ij /N = |∆| 2 ∂ 2 E GS /∂ ∆i ∂∆ j is a response to an extra term in the Hamiltonian of the form ∆i c i↑ c i↓ + h.c..It is evaluated in Appendix B and shown to be g ∼ |∆| 2 /J, with the mean-field pairing field |∆| given by Eq. (10).
The action (18) describes a quantum version of the celebrated classical Kuramoto model [47][48][49][50][51][52][53][54][55][56][57].The latter was proposed [47] to describe synchronization of coupled non-linear oscillators.It's quantum version, Eq. ( 18), may be interpreted as N -body quantum mechanics of particles with mass m and coordinates φ i , residing on the unit circle and interacting via all-to-all cospotential.The synchronized phase of the classical Kuramoto model is analogous to a φ-localized ground state wavefunction of this quantum mechanics.Within the SYK+Hubbard model such synchronized phase means globally phase-coherent superconductivity with ODLRO.Below we show that the synchronized phase of the quantum Kuramoto model, Eq. ( 18), emerges above some critical coupling g > g c (i.e. at U > U c ) as a continuous QPT.
Since the ground state is expected to be symmetric with respect to particle permutations, it may be thought off as a Bose condensate.Due to all-to-all nature of the interactions, the Bose condensation in the large N limit is accurately described by the Gross-Pitaevskii equation.
In the present context it takes the non-local form: where the condensate wave-function is normalized as 2π 0 dφ |Ψ(φ)| 2 = N and obeys the periodic boundary conditions, Ψ(2π) = Ψ(0).Employing separability of the exponential potential, e ±i(φ −φ) , one may reduce the non-linear equation (19) to the linear Matheiu equation: supplemented with the self-consistency condition where ρ 1 the first Fourier harmonics of the normalized condensate density, |Ψ(φ )| 2 /N .The strategy is to find a ground state wave-function of the Matheiu equation (20) for a given amplitude of the cos-potential, gρ 1 , and substitute it into the self-consistency condition (21) to find ρ 1 .A trivial solution, ρ 1 = 0, with the uniform condensate, Ψ = N/2π, and µ = 0 exists for any g.A non-trivial solution with µ < 0 requires g > g c .To find the non-trivial solution, one notices that the right hand side of Eq. ( 21) is an odd function of gρ 1 .Its behavior at small gρ 1 may be found from the first order perturbation theory for the Matheiu equation (20), yielding the linear slope 2mgρ 1 .On the other hand, at large mgρ 1 1 the ground state wave function of Eq. ( 20) is a narrow Gaussian, centered at φ = 0.This implies that the right hand side of Eq. ( 21) saturates to one for mgρ 1 1.As a result, Eq. ( 21) is the standard mean-field equation for a second order transition with the order parameter ρ 1 .It yields a finite order parameter ρ 1 ∝ √ g − g c for g g c with g c = 1/2m ≈ 0.48J, Fig. 10.An alternative way to determine g c is to investigate a spectrum of linearized fluctuations on top of the uniform solution, Ψ(φ, t) = N/2π + l ψ l e ilφ−iω l t , where l = ±1, ±2, . . .labels angular momentum components.Substituting this into the time-dependent Gross-Pitaevskii equation, Eq. ( 19) with i∂ t Ψ on the right hand side, and linearizing it with respect to ψ l , one finds the spectrum: Therefore for g > g c = 1/2m the frequency of the l = ±1 components becomes imaginary, indicating instability towards a non-uniform condensate.This expression shows that the continuous QPT is indeed associated with the time-scale ω −1 ±1 ∝ |g c − g| −zν , which is divergent at the transition with the Gaussian exponent zν = 1/2.
We thus conclude that the quantum Kuramoto model exhibits the synchronized phase for mg > 1/2, where the local phases, φ i , are coallesing.In the large N limit this spells spontaneous breaking of the U (1) symmetry.In terms of the SYK+Hubbard model these observations translate into formation of ODLRO for U > U c , where, employing Eqs. ( 10 For U < U c the on-sites phases φ i fluctuate freely and prevent formation of the global ODLRO.This phenomenon renders the mean-field treatment of Sec.IV grossly inadequate for U < U c and leads to creation of the pseudogap phase.The latter is characterized by the evenodd alternation in the ground state energies, cf.Fig. 6, thus exhibiting a single-particle energy gap (i.e. a finite energy to add or subtract a single fermion).However, there is neither ODLRO nor the many-body gap within a sector with a fixed N f .Therefore from the transport perspective, the pseudogap state is characterized as an insulator.Correspondingly the Kuramoto QPT should be termed an insulator-superconductor one.
The line 2πT ≈ ω ±1 , Eq. ( 22), spells the boundary of the quantum critical regime.If 2πT < ω ±1 , the quantum Kuramoto phase fluctuations, governed by e iφi(τ ) e −iφi(0) = e −ω±1|τ | , are averaged out to zero.This leads thus to the familiar SYK non Fermi liquid fermionic correlations.However, for ω ±1 < 2πT |∆| the imaginary time circle is too short to completely wash out the superconducting correlations.This creates an interesting quantum critical scenario, where superconducting correlations show up as a finite temperature effect.

VII. TOWARDS A HOLOGRAPHIC INTERPRETATION
In this section we briefly comment on a possible holographic interpretation of our findings.Recall that at T = 0 we have seen formation of local Cooper pairs at arbitrary small attraction between fermions.Their phases are incoherent at intermediate U , separated by the continuous QPT from the superfluid phase with ODLRO at large U .The complex SYK dot we work with is now used as a toy model for "near AdS/almost CFT" correspondence in quantum mechanics.From a higherdimensional perspective the Reissner-Nordstrom (RN) black hole (BH) is considered as the bulk whose geometry involves a long AdS 2 throat near the horizon.The large N SYK quantum mechanics lives at the boundary of the throat and conjecturally is dual to the AdS 2 near horizon, flavored with some matter.
The effective low energy boundary action for the unperturbed complex SYK involves two Goldstone modes -the reparameterization of time, governed by the Schwarzian action, and the U (1) phase field, φ(τ ), governed by the kinetic term φ2 (see, for example [66] and references therein).The bulk theory in addition to Jackiw-Teitelboim (JT) gravity involves the gauge fields and some matter.If we focus at T = 0 the extremal RN BH is unstable under small perturbations.The mode of instability can be interpreted as the Schwinger pair creation (see [67] for a review).Usually the bulk instability is treated as formation of the homogeneous condensate described in terms of the boundary behavior of the bulk complex scalar or the bulk fermion.
The individual local Cooper pairs play an important role in our analysis hence their holographic meaning needs to be clarified.Let us start with the bulk identification of the Goldstone phase field.To this aim consider for a moment the U (1) bulk 2d field (A τ , A r ) with the boundary behavior involving chemical potential and density In the boundary theory the density, ρ, and the phase, φ, Hence the phase has to be canonically conjugated to F τ r in the bulk.To get the correct conjugated variable recall the canonical pair in 2d gauge theory which allows to identify the phase field, φ(τ ), as the gauge holonomy along the radial direction, r.
Note that if we choose A r = 0 gauge, the holonomy factor appears in the boundary conditions.A similar identification of the Goldstone phase modes has been developed in holographic QCD [68] and in the holographic hydrodynamics [69].In QCD the bulk flavor gauge group SU (N f ) L ×SU (N f ) R is broken by the Higgs mechanism down to the diagonal SU (N f ) and the pions π a , which are non-abelian Goldstone phases, of the chiral (excitonic) condensate are identified as exp(iπ a t a ) = P exp drA r (r, x).In the holographic hydrodynamics a similar identification of the Goldstone phase is emerging upon breaking of U (1) × U (1) symmetry to the diagonal U (1) [69].
We turn now to the interpretation of the Hubbard U .Fortunately, the Hubbard model has been treated in the holographic approach for Bose [70] and Fermi systems [71], where it was realized that the Hubbard U is to be identified with the radial position of the hard wall r U = U .Therefore the control parameter, U/J, tells how close to the horizon the hard wall is placed.Small U corresponds to IR near horizon region, while large U corresponds to the hard-wall at UV near the boundary of AdS 2 .As follows from our analysis, a perturbation induced by the IR wall at small U amounts to the instability of the extremal RN geometry and formation of the Cooper pairs.At large U the gap of an individual Cooper pair, ∆ ∼ U , fits the length of two strings extended up to the U scale, representing two fermions at the boundary.
The most subtle question concerns the identification of the bulk counterpart of Cooper pairs.To formulate the conjecture let us remind universal aspects of the instanton solutions in different dimensions.Consider, for instance, instanton in the SU (N ) gauge theory on the R 3 × S 1 geometry.If there is non-vanishing holonomy around S 1 the instanton with unit topological charge is split into N constituents with the topological charges 1/N and non-vanishing monopole charges [72] (the total monopole charge is zero, of course).This is the "caloron" solution, known for the theory with finite temperature, or with one compact space coordinate.The positions of the constituents (x 1 , . . ., x N ) are fixed by the eigenvalues of the SU (N ) holonomy, V , around S 1 , V = diag(e ix1 . . ., e ix N ).
The holographic Skyrmion solution in D = 4 is the instanton in the gauge theory with the flavor gauge group [73], which involves three space coordinates and the radial coordinate A inst (x 1 , x 2 , x 3 , r).In the string theory framework equivalently it is a baryon vertex identified with the particular D-brane wrapped around the internal sphere [74].Due to the anomaly, N strings are attached to the baryon vertex located at some point in the bulk which amounts to N fermions at the boundary.In holographic QCD the baryonic vertex is placed dynamically nearby the effective IR wall [68], which at small U corresponds exactly to our r U = U scale.
Let's assume for a moment that our large N SYK+Hubbard dot is a kind of Skyrmion-instanton state that is a baryon vertex placed at r U in the throat region like in holographic QCD.N strings are attached to the vertex hence this picture at the first glance has nothing to do with the ensemble of Cooper pairs.However, we have identified the Goldstone phases as the holonomies in the radial coordinate and, as was established here, these holonomies do not vanish.Since the Skyrmioninstanton solution involves the radial coordinate, the nontrivial SU (N ) radial holonomy splits the Skyrmioninstanton into Skyrmion-caloron with at most N/2 fractional Skyrmion constituents with a fractional "topological charge".Such fractional Skyrmions host now two strings instead of N strings and therefore amount to the pair of fermions at the boundary.Hence the fractional Skyrmion is a candidate for a bulk counterpart of an individual Cooper pair.Similarly to the standard caloron, the holonomies of the Cooper pairs correspond to the positions of constituents on the "dual circle" providing the Kuramoto-like picture for the individual phases.
Establishing a potential for phases or the fractional Skyrmions is a dynamical issue.The potential for the phases of constituents of the caloron solution involves perturbative and non-perturbative contributions and typically reads as where δ is a gap in the model.I.e.typically potential for the phases involves only nearest neighbors interaction.In the case of Skyrmion-caloron the all-to-all SYK Hamiltonian apparently induces an all-to-all interaction between phases of individual components.
Summarizing, we conjecture that there is the baryonic vertex placed nearby r U radial scale.The non-trivial phases, corresponding to the radial holonomies exp(iφ i ) of the U (1) i bulk gauge field from the total U (1) N flavor gauge group, result in the splitting of the Skyrmion supporting N strings into the fractional Skyrmions supporting two strings.The disorder SYK Hamiltonian induces the all-to-all Kuramoto potential for the phases of the Cooper pairs and the phases become synchronized at some position of the hard wall specified by U c .Of course, many aspects of this conjecture deserves clarification.
Note some analogy with QCD at non-vanishing density.It is well-known that at large baryonic density QCD is in the color-flavor locking phase with the Cooper condensate of quarks.However it was argued in [75] that at smaller chemical potential there is a transition from Skyrmions into half-Skyrmions.It is assumed that at the transition the common gap and exciton(chiral) condensate disappears.Still there are "islands" of gapped phase with disordered chiral phases.This resembles the behavior of our model near the QPT.
Two additional remarks are in order.The insulatorsuperfluid QPT in 2 + 1 has been discussed in the holographic framework in [76] and has clear parallels with our 0 + 1 case.The insulator phase was related with the AdS soliton background while the superfluid phase with the AdS BH background.The AdS soliton solution has the effective IR cut-off at a tip of the cigar, which is an analogue of our small U regime, since U provides the IR cut-off as well.When U is large it no longer serves as an IR parameter, yielding the UV scale instead.The BH physics starts to dominate in the superfluid phase in IR similar to our case.
Another point concerns the origin of the Hubbard perturbation of the SYK model.We have chosen it by hand, but it has appeared in the holographic setup in an interesting manner in large N N =4 SYM theory [77].Namely, consider a string moving in S 3 or equivalently study the anomalous dimensions of the particular scalar operators with large conformal dimension in the boundary theory.In this scenario, the dilatation operator at three loop exactly coincides with the Hubbard Hamiltonian.The latter plays the role of a conventional Hamiltonian of the discretized string, propagating in the nontrivial background.Parameter U in the Hubbard model is identified with an inverse coupling in the boundary theory.It is unclear at the moment if these ideas provide an additional intuition for our model.The discussion in this section is clearly only qualitative and tentative.We postpone a more detailed analysis of the holographic picture for a separate study.

VIII. CONCLUSION AND OUTLOOK
Following the earlier studies [39][40][41][42], we found that the spin-full version of the SYK model with an extra attractive interactions may exhibit ODLRO and superconduc-tivity.Furthermore we found that details of this extra attraction are crucially important in dictating the global phase diagram of the model.The previous studies focused on an effective all-to-all attraction, which conform to the large N mean-field treatment.The latter calls for superconducting instability of the non Fermi liquid groundstate at an arbitrarily weak attraction.This is indeed the case for the SYK+pair hopping model briefly considered here.
Our main finding is that a local attraction, such as onsite negative U Hubbard term, leads to a qualitatively different scenario of the superconducting transition.In this case the physics is dictated by quantum fluctuations of local phases.They destroy ODLRO in a sizable part of the phase diagram, confining the superconductivity to a dome-like region, Fig. 5.In particular, they lead to the pseudogap phase at small U and the continuous QPT to the superconducting phase at U = U c .These features are described by the quantum version of the celebrated Kuramoto model.At strong attraction, the local nature of the attractive interactions is also of crucial importance, resulting in T c ∼ U −1 scaling of the critical temperature.This limit is mapped on the Richardson-like model with two-boson hopping.It's exact solution predicts the first order transition at T = T c from ODLRO into a bosonic insulator state.The latter consists of fermions, paired with the binding energy U T c , forming a gas incoherent bosons.Fermion transport in this state is suppressed as e −U/T .
We list now some of the open questions raised by our study: (i) What are fermionic correlation functions in the pseudogap phase at U < U c ?The naive answer is that they are the same as in the non Fermi liquid SYK model.Yet, contrary to SYK, fermions interact with the dynamical phases as |∆|e iφi(τ ) c i↓ c i↑ + h.c., where the phases, φ i (τ ), are governed by the Kuramoto quantum mechanics, Eq. (18).Close to the QPT this dynamics becomes increasingly slow, Eq. ( 22), and may significantly alter the fermionic correlation functions.
(ii) What are the implications of our 0D treatment for the array geometry?In particular, is the dome-like phase diagram, Fig. 5, applicable to arrays and how it depends on the coupling (hopping) strength between the dots in the array?(iii) Is there an interaction and an interplay between the phases, governed by the Kuramoto, and the reparametrization modes [2,61], governed by the Schwarzian action?The latter modes are described by the Liouville quantum mechanics [61], which predicts metal-insulator crossover at the energy scale J/N .For a finite N this energy scale may compete with the many-body gap |∆| 2 /J, possibly affecting the insulatorsuperconductor QPT [36].
(iv) An interesting generalization is a model with a weak time reversal symmetry breaking parameter.In the Richardson model such generalization leads to the Russian Doll (RD) model, Appendix C, which is known to be integrable.One may expect that deformed in this manner the large U generalized Richardson is also integrable.SYK corresponds to the completely degenerate local Richardson parameters, i = 0, which means that holographically all flavor branes are sitting on the top of each other in the IR and the SU (N ) symmetry is classically unbroken.Generic values of i correspond to displacements of flavor branes in the radial coordinate in the holographic treatment of Richardson or RD models.It would be interesting to elucidate the role of non-vanishing local parameters, i , in the generalized Richardson model.
(v) The quantum Kuramoto mechanism of the condensate formation could fit within a more general framework.In particular, an intermediate pseudogap phase is believed to exist in the thermal QCD below the deconfinement phase transition, where the local phases of the chiral condensate are disordered.The synchronization of the chiral phases leading to formation of the homogeneous chiral condensate may occur in a Kuramoto-like way.Indeed as shown above, at the 1/N order the nearhorizon gravity (RG) dynamics induces the Kuramoto potential for phases of the local Cooper pairs.Formation of the chiral condensate in the holographic QCD, being also a near-horizon effect, may thus lead to a nonabelian generalization of the Kuramoto potential for the exciton pairs.Furthermore, we perform Hubbard-Stratonovich transformation to decouple the Hubbard term in the Cooper channel, introducing site-local complex fields ∆ i , exp After the decoupling procedure, the action reads: (A9)

Saddle point ansatz
We assume the fields ∆ i to be time-and siteindependent at the saddle point.Then we integrate out fermion fields, and obtain the action in the form Variation of the action Eq.(A10) results in the following set of saddle point equations Hereafter we restrict ourselves to the case of the halffilling, where, due to the particle-hole symmetry, the normal components are odd, while anomalous are even functions of time, eg., Ξ τ τ = Ξ τ τ , Σ τ τ = −Σ τ τ .

Approximate solution of the mean-field equations
The anomalous fields Ξ and F , entering the saddle point equations, admit non-zero solutions only in the presence of ∆.Similarly to the BCS case, we'll find that F ∝ ∆.According to Eq. (A13), Ξ ∝ F 3 ∝ ∆ 3 .Therefore in the limit of (exponentially) small ∆ one may consider dropping Ξ from the set of the mean-field equations and restricting them down to: where we fixed the phase of ∆ to make the latter real.We'll see below that neglecting Ξ is not, strictly speaking, justified, even for the small ∆.Nevertheless Eqs.(A17)-(A19) will be shown to be a qualitatively (if not quantitatively) accurate representation of the full set.Since all propagators are site-diagonal, correlations between distinct sites appear in the order 1/N in expansion of the action.The physical mechanism of correlations between the superconducting fluctuations at different sites consists of correlated hopping of Cooper pairs facilitated by SYK-interactions J ik;lj c † iσ c † kσ c l,σ c j,σ .Because all four sites here are distinct, no direct hopping of a Cooper pair is possible.Rather, a transfer of a single Cooper pair from a site i to another site j involves at least two correlated acts of interaction that cause transfer of two Cooper pairs from the sites i, k to the sites , j. Thereby the second Cooper pair plays the role of an assisting agent for the hopping of the first one.The amplitude of an elementary jump of a Cooper pair from a site i to a different site j, assisted by a hopping of another Cooper pair from a site k to a site is represented by the diagram in Fig. 12a).The hopping of the assisting Cooper pair is depicted by the insertion of an anomalous loop between the normal Green's functions.Thus, insertion of anomalous loops is a necessary ingredient of diagrammatic representation of interaction between superconducting fluctuations at different sites.
Taking into account summation over the spin indexes and over the intermediate sites k, , one obtains contribution to the average susceptibility, κ ij , from the lowest order diagram in Fig. 12a): Substituting the variational solutions, Eqs.(A20), (A21), for normal and anomalous propagators and introducing dimensionless time-variables t where C (1) 2 is given by a convergent integral, which does not depend on any parameters, This numerical constant should not be taken too seriously.Indeed, our variational ansatz for the propagators, Eqs.(A20), (A21), is not exact but only interpolates between correct short and long time asymptotics.The reason of presenting this calculation is to point out the absence of logarithmic factors.The latter may be naively expected, due to the presence of two runs of the Cooper ladder in the diagram of Fig. 12a).If the anomalous loop in the middle would be confined in time to some scale τ 0 τ ∆ , the diagram would be ∝ ln 2 (τ ∆ /τ 0 ) 1.This is because the two integrals over τ and τ 1 ≈ τ 2 would be logarithmic.In this case summation of the entire Cooper ladder, Fig. 12b), would be of a crucial importance.However, our case happens to be different.The reason is that the anomalous loop has the same characteristic time scale, τ ∆ , as the normal Green functions, which form runs of the Cooper ladder.As a result, logarithms are not present and all the terms of the ladder have the same order of magnitude as the first diagram, Fig. 12a).Therefore the ladder summation only changes the numerical coefficient, C 2 , rather than the large logarithmic factor.Let us note in passing that Hubbard U , being time-local, induces the conventional logarithmic Cooper ladder and thus Eq. (A27).This ladder, however, is strictly diagonal in the site index (and is already incorporated in diagonal G and F , Eqs. (A20), (A21)).The off-diagonal ladder and thus the off-diagonal susceptibility, κ ij , needed in the quantum Kuramoto action, requires long-range anomalous loops inserted in each run of the ladder as in Fig. 12b).
Another consequence of the long-range nature of the anomalous loop is that the susceptibility, κ ij , Eq. (B2), is not proportional to ∆, despite each of the two anomalous propagators, F , being proportional to ∆, Eq. (A21).The reason is that the integrations range, given by τ ∆ , is inversely proportional to ∆ 2 , Eq. (A24).In the absence of other long time cutoffs, e.g. a finite size gap, this leads to ∆-independent susceptibility, Eq. (B2).
Finally, the interaction term in the quantum Kuramoto action, Eq. ( 18), is given by κ ij ∆ i ∆j + h.c.∼ κ ij |∆| 2 cos(φ i − φ j ).As a result, the interaction constant, g, in the quantum Kuramoto action, Eq. ( 18), is where d is a mean value of ( i+1 − i ) sequence.The RD model reduces to the Richardson model in the limit η → 0.
The RD model turns out to be integrable as well.Now instead of the Gaudin model a proper spin chain counterpart is the generic quantum twisted inhomogeneous XXX spin chain [82].The equation defining a spectrum of the RD model reads as and coincides with the BA equations for the spin chain.It reduces to the BA equation of the Richardson model (C5) in the limit η → 0.
The RD model enjoys the gap equation, which reads as follows: , ∆j = ∆ j e iφi , (C15) where V ij is a scattering potential, which depend on the parameters G, α.In the thermodynamical limit it becomes an integral equation with multiple solutions for the gaps.Different solutions to the gap equation yield different superconducting states.
Solutions of the gap equation in the large N limit are parametrized as follows: ∆ n = ω sinh t n , t n = t 0 + πn θ n = 0, 1 . . ., (C16) where t 0 is a solution to the following equation: tan(θt 0 ) = θ g 0 < t 0 < π θ (C17) and ω = dN for equal spacing ( i+1 − i ) = const.This behavior can be derived in the mean-field approximation [81].In the limit θ → 0 the gaps ∆ n>0 → 0 and This way the standard BCS expression for the gap is recovered.At a weak coupling the gaps behave as For Cooper pair degeneracies on orbitals, d i , the RD model is modified a bit and is related to the higher spin XXX spin chain.The local spins s i are determined by the corresponding pair degeneracy, d i , of the i-th orbital and the corresponding BA equations read as exp(−2iα) N E j − k − iη/2 + iηs i The RD model involves an interesting RG behavior of couplings with respect to RG time s = log N , [81].The coupling constant exhibit the cyclic RG flow (a recent review on the cyclic RG can be found in [83]), while the TRS parameter does not renormalize g(s + λ) = g(s), g(e −λ N ) = g(N ).(C23) The RG period reads as and the total number of the independent gaps in the model is The cyclic RG behavior reflects the breaking of the scale invariance down to the discrete subgroup and the spectrum of gaps manifests in the Efimov scaling The sizes of the Cooper pairs in the n-th condensates also have the Efimov-like scaling.

Possible generalizations
Here we consider generalizations of the Richardson model, involving four-boson interactions.The Hamiltonian (12), appropriate for large U , is which follows from the hidden algebraic structure of the Gaudin model.Therefore R 2 i yield two-boson interaction term only.
To obtain the four-boson interaction term we can consider the quadratic form with arbitrary matrix, A ij .The integrable Hamiltonians, H 4 , involve the desired four-boson interactions.In general, if i = 0, the resulting interaction coupling constants are site-and i -dependent.In our case all i = 0 and hence the Hamiltonian (C27) can be considered as the peculiar limit of the generic quadratic form, Eq. (C29).Moreover, all Bethe states creation operators, B i , at i = 0 reduce to the single operator B 0 = i b † i .

12 FIG. 1 .
FIG. 1. Spectrum of ρij, i.e. λα vs. number of orbitals N for the ground state of SYK+Hubbard model with U/J = 2 and N f = N .Dashed line is a result from generalized Richardson model (Section VI A).

FIG. 3 .
FIG.3.ODLRO (top) and many-body energy gap in units of J (bottom) vs. U/J for the SYK+pair-hopping model with N f = N = 12.

FIG. 6 .
FIG.6.Many-body energy spectrum (in units of J) of SYK+Hubbard model vs. number of fermions N f for N = 8 and U = 2J.Chemical potential, µ, is set µ = −U/2 to preserve particle-hole symmetry between N f and 2N − N f sectors.The inset shows the lowest bands of the spectrum for even N f .Band-width for the lowest energy sector at halffilling is consistent with N W/16 0.023, predicted by the generalized Richardson solution, Eq. (15).
N i b i and the boson number operator Nb = N i b † i b i .Employing the (anti)commutation relations for the hard-core bosons: b

FIG. 9 .
FIG. 9. Spectra of n|b †i bj|n for each many-body state |n vs. it's energy, E − EGS, for (a) SYK+Hubbard with N = 12 and (b) effective low-energy bosonic theory with the Hamiltonian(11).Black dashed lines are energies of the generalized Richardson model, Eq. (15).
) and (B3), U c ≈ J √ π/(4 √ 2 log C 2 ), see Appendix B. The quantum Kuramoto model synchronization transition is indeed seen in the exact exact diagonalization of the SYK+Hubbard model, Fig. 4, as the continuous QPT at U = U c .

FIG. 11 .
FIG. 11.Schematic picture of the conjectured bulk representation of Cooper pairs.The SYK and the hard wall are placed in the throat region of the BH geometry.A baryonic vertex is near the radial position of hard wall and gets fractionalized into constituents supporting two strings each.

FIG. 12 .
FIG. 12. Diagrammatic representation of the off-diagonal Cooper susceptibility: a) the lowest order diagram; b) the ladder series.
[84]Hence one may question if Hamiltonians with four-boson interactions can be derived from the commuting set, R i .Such representation would prove the integrability of the model.It is known that the Hamiltonians, R i , obey a nontrivial algebraic relation[84]