Pairing Instabilities of the Yukawa-SYK Models with Controlled Fermion Incoherence

The interplay of non-Fermi liquid and superconductivity born out of strong dynamical interactions is at the heart of the physics of unconventional superconductivity. As a solvable platform of the strongly correlated superconductors, we study the pairing instabilities of the Yukawa-Sachdev-Ye-Kitaev (Yukawa-SYK) model, which describes spin-1/2 fermions coupled to bosons by the random, all-to-all, spin-independent and dependent Yukawa interactions. In contrast to the previously studied models, the random Yukawa couplings are sampled from a collection of Gaussian ensembles whose variances follow a continuous distribution rather than being fixed to a constant. By tuning the analytic behaviour of the distribution, we could control the fermion incoherence to systematically examine various normal states ranging from the Fermi liquid to non-Fermi liquids that are different from the conformal solution of the SYK model with a constant variance. Using the linearised Eliashberg theory, we show that the onset of the unconventional spin-triplet pairing is preferred with the spin-dependent interactions while all pairing channels show instabilities with the spin-independent interactions. Although the interactions shorten the lifetime of the fermions in the non-Fermi liquid, the same interactions also dress the bosons to strengthen the tendency to pair the incoherent fermions. As a consequence, the onset temperature $T_c$ of the pairing is enhanced in the non-Fermi liquid compared to the case of the Fermi liquid.


Introduction
Understanding unconventional superconductivity of strongly correlated electrons is a long-standing goal of modern condensed matter physics [1][2][3][4][5][6][7][8][9][10]. It is generally believed that dynamical interactions mediated by collective charge or spin fluctuations are responsible for the Cooper pairing in the correlated superconductors [11,12]. Major challenges of the problem come from the emergence of non-Fermi liquid normal states due to the same dynamical interactions [13][14][15][16][17][18]. Since both superconductivity and non-Fermi liquid are stabilized by the same physical origin, systematic investigations of two competing effects are necessary [19][20][21][22]. However, the strongly coupled nature of the problem makes it difficult to draw a concrete theoretical conclusion as no small parameter exists to control the theory perturbatively.
In this work, we circumvent such difficulty by examining a variant of the Sachdev-Ye-Kitaev (SYK) model [23][24][25][26][27], so-called the Yukawa-SYK model [28][29][30][31], which is exactly solvable and supports non-Fermi liquid ground states. While the previously studied models use the fixed variance of the random coupling, we introduce a continuous distribution of variances. The model consists of N number of fermions (c i=1,...,N ) strongly coupled to M number of bosons (φ k=1,...,M ) via the random all-to-all Yukawa couplings (g i j,k ): where σ a is the Pauli matrix acting on the spin space α, β =↑, ↓, and the summation is assumed for the repeated Greek indices. Physically, the scalar bosons φ k represent the collective charge (or spin) fluctuations of the fermion bilinear i, j c † iα c jα or i, j c † iα σ 3 αβ c jβ . The recurring theme of the SYK model and its variants is that the disorder averaging over the random coupling constants (g i j,k ) suppresses almost all quantum fluctuations except one or a few families of the Feynmann diagrams in the large M and N limits [23,24,32]. With those handful number of diagrams, we can solve the model self-consistently and identify the leading pairing instabilities without any perturbative approximations.
It is important to note that the Yukawa-SYK model is defined by not only the Hamiltonian but also the statistical properties of the random couplings (g i j,k ). Most of previous studies on various families of the SYK model focused on the random couplings with zero mean (g i j,k = 0) and constant variance ((g i j,k ) 2 = λ) [26-31, 31, 33-38]. However, we can also consider the random couplings whose variances obey a well-defined distribution, i.e., (g i j,k ) 2 = λ k has the k dependence such that the set of the variances {λ k } forms a continuous distribution ρ(λ) in the large M limit. Pioneering work on the low-rank SYK models [39], which are equivalent to the Yukawa-SYK models with the extensive (M /N ∼ O(1)) number of nondynamical massive bosons, first notices the significance of the distribution ρ(λ); depending on the singular behaviour of the distribution ρ(λ) near the maximum variance λ max , the low-rank SYK models show a rich variety of the low energy states ranging from the Fermi liquid to non-Fermi liquids [39,40]. By tuning the distribution ρ(λ), we can systematically control the fermion incoherence and push the system toward the non-Fermi liquid. Therefore the current variant of the Yukawa-SYK model is an excellent solvable platform to examine the interplay of non-Fermi liquid and superconductivity with the distribution ρ(λ) as a theoretical handle to control the incoherence of fermions.
While the flourishing papers discussed the SYK superconductivity, they focused on the fast scrambling conformal solution of the SYK model (and its variants) with a fixed constant variance [28,29,[34][35][36][37][38]. The pairing instabilities of the Fermi liquid and the nonconformal non-Fermi liquid states of the low-rank SYK models are not examined yet [39]. Since the variance distribution ρ(λ) opens up a new direction of the controllability for the "non-Fermi-liquidness", we would like to understand whether the strong interaction, which makes the fermions more incoherent but the bosons to glue the fermions stronger, is an ally or a foe of the Cooper pair formation. The enhanced transition temperatures T c of the non-Fermi liquid state ( Figure 3) demonstrate that the highly incoherent fermions can prefer the pairing more than the well-defined quasiparticles of the Fermi liquid due to the significant enhancement of the bosonic glue in the Yukawa-SYK model. Furthermore, to understand the distinct contributions of the collective charge and spin fluctuations to the pairing, we examine both the spin-singlet and triplet pairing instabilities with the linearized Schwinger-Dyson equations. The unconventional dynamical pairing between the equal-spin fermions at distinct times, i.e., 〈c † is found to occur. The remaining part of the paper is organized as follows. In Sec. 2, we introduce the Yukawa-SYK model and its effective action in terms of the Green functions and self-energies. Sec. 3 discusses the Schwinger-Dyson equations, which are the saddle point equations of the effective action. We first consider the high-temperature normal state solutions in Sec. 4, which demonstrate how the distribution of variances can result in both the Fermi liquid and non-Fermi liquids. Then, in Sec. 5, we discuss the leading pairing instabilities of the Fermi liquid and the non-Fermi liquid normal states by solving the linearized Schwinger-Dyson equations. At last, we summarize and conclude our work in Sec. 6.

Model
We consider spin-1/2 fermions (c) coupled to real scalar fields (φ) by all-to-all random Yukawa couplings (g), S = S c + S φ + S g : We use the natural unit ħ h = k B = 1 so that β = 1/T is the inverse temperature. Since S c and S φ are invariant under spin rotation, it is sufficient to investigate two cases: a = 0 and a = 3.
The real symmetric Yukawa couplings g i j,k = g ji,k ∈ are sampled from the Gaussian orthogonal ensemble (GOE) for each k, i.e., g i j,k follows the Gaussian distribution with zero mean g i j,k = 0 and variance g i j,k g i j ,k = λ k δ k,k (δ ii δ j j + δ i j δ ji ) for λ k > 0. Assuming that the model is self-averaging, we can derive the effective action from the disorder average of the partition function Z: where Note that the pairing correlations among fermions are generated because the random Yukawa couplings are averaged over GOE [36].
With the bilocal fields we can rewrite the interacting part of the effective action S λ defined in Eq. (5): where γ = M /N ∼ O(1) is the ratio between the number of bosons and fermions. Note that D(τ, τ ) is the bilocal field that becomes the sum of the bosonic propagators weighted by the variances λ k at the saddle point of the action. By introducing the Lagrange multipliers Σ, Φ + , Φ, and Π, we can relate the dynamics of the fermions and bosons with the bilocal fields G, F , F + , and D, respectively (see Appendix A). Physically, the bilocal fields become the fermionic (G, F, F + ) and bosonic (D) Green functions, and the Lagrange multipliers become the corresponding fermion (Σ, Φ + , Φ) and boson (Π) self-energies, at the saddle point of the action. In this model, the bosonic part of the action S φ = S φ + S Π (see Appendix A for the definition of the bosonic self-energy action S Π ) needs special attention because the bosons may condense at low temperatures. After the Fourier transformations, we split S φ into the normal [Eq. (11)] and condensed parts [Eq. (12)]: where ν n = 2πn/β are the bosonic Matsubara frequencies. The bosons are condensed when the quadratic potential for some bosonic modes is no longer convex. As the zero frequency modes φk(0) with λk = λ max = max[{λ k }] first become unstable when m 2 − λ max Π(0) = 0, they are condensed at T < T BEC [39]. Then can be treated as a classical degree of freedom. By integrating out the fermions and remaining uncondensed bosons, we can obtain the large N effective action S eff in terms of the bilocal fields and the Lagrange multiplier fields (see Appendix A).

Schwinger-Dyson Equations
In the large M and N limits, the saddle point of S eff precisely describes the low-energy dynamics of the Yukawa-SYK model. Hence, we derive the Schwinger-Dyson equations from the condition δS eff = 0. The normal (G) and anomalous (F ) Green functions for the fermions are where the spin indices are suppressed for notational convenience. We assume that the set of variances {λ k } forms a well-defined distribution ρ(λ) in the large M limit: which is regular at λ = λ max for η > 0 (class I) but diverges algebraically as λ → λ max for −1 < η < 0 (class II) ( Figure 1) [39]. Then the bosonic propagator is The first part of Eq.
where the spin indices for G, F , Σ, and Φ are suppressed for simpler notation, and "tr" is the trace over the spin indices. After we plug in Eqs. (14) and (15) to Eqs. (18 -20), we can get a set of nonlinear equations for the normal and anomalous fermionic self-energies Σ(iω n ) and Φ(iω n ).
Since the fermions would not be paired at high temperatures, our analysis starts from the normal state with F (iω n ) = Φ(iω n ) = 0.

Normal State Analysis
Without the pairing among fermions, Eq. (14) gives the fermion Green function where Σ α (iω n ) ≡ Σ αα (iω n ). For both a = 0, 3 in S g , the fermion Green function is spin-diagonal (G ↑↓ = G ↓↑ = 0) and independent of spin polarization (G ↑↑ = G ↓↓ ). Therefore we write and the effective condensateφ = ϕ + γD N (0)/βλ max . Then the fermion Green function solves the Schwinger-Dyson equation with J(iω n ) = ω n + iΣ N (iω n ) [39]. With our model distribution ρ η (λ) in Eq. (16), the propagator for the uncondensed bosons where The function f η is positive and monotonic, and Figure 2 (a)]. Note that the distribution ρ η (λ) shows larger value near λ = λ max when η is smaller, i.e., , it is enhanced when there is higher chance to sample the Yukawa couplings g i j,k with large variance λ k . The asymptotic expansion of the hypergeometric function 2 F 1 (a, b; c; x) gives The class II bosonic propagator (η < 0) is larger than the class I propagator (η > 0) for all frequency range. Since the distribution ρ η<0 (λ) is mostly concentrated around λ ∼ λ max , there is high chance to sample strong Yukawa coupling g i j,k . Hence, the bosonic propagator is more strongly enhanced by the interactions between fermions and bosons.
A qualitatively important distinction between class I (η > 0) and class II (η < 0) is the boundedness of the Green function for uncondensed bosons, D N (iν n ). While f η (x) ≤ 1/η is bounded from above for class I (η > 0), f η (x) diverges algebraically as x → 1 − for class II (−1 < η < 0). Thus, the boson Green function D(iν n ) for the class I model is bounded from above if there were no Bose-Einstein condensation. When the boson Green function is bounded from above, the bosons can yield limited quantum corrections to the fermion self-energy. Hence, the fermion self-energy also becomes bounded from above. Without large fermion self-energy, the strong Yukawa interactions result in large boson self-energy which can make the bosons unstable, i.e., the renormalized squared mass of the bosons m 2 ren = m 2 − λΠ(0) can be negative due to the large boson self-energy at zero frequency Π(0). Hence, the zero frequency bosons need to be condensed when x = λ max Π(0)/m 2 = 1. If the bosons are condensed, the total boson Green function D(iν n ) = (β/γ)λ max ϕδ n,0 + D N (iν) is no longer bounded from above because the Bose condensate ϕ is not bounded from above. Therefore, the Bose-Einstein condensation is essential for the class I model (η > 0). On the other hand, the class II model (η < 0) or the model with a fixed variance coupling ρ(λ) ∼ δ(λ − λ 0 ) do not show the Bose-Einstein condensation. More detailed mathematical discussion about the Bose-Einstein condensation in the Yukawa-SYK model can be found in Appendix B.
In the absence of the pairing F = Φ = 0, the same Schwinger-Dyson equations are solved in the context of the low-rank SYK models, which can be obtained from the Yukawa-SYK models by integrating out the massive bosons. Since the asymptotic expansion of our bosonic propagator, Eq. (29), coincides with the bosonic propagator in Ref. [39], thermodynamics of the Yukawa-SYK models are equal to that of the low-rank SYK models. For the self-contained discussion, we will briefly review the frequency scaling and thermodynamics of the normal state Green functions and self-energies [39].
The distribution of variances ρ η (λ) plays important role in the frequency scaling of the fermion Green function and the self-energy. For small frequencies at low temperatures, The self-consistency of the above frequency scaling properties can be checked as follows. Since the overall argument for class I (η > 0) and class II (−1 < η < 0) are similar, let us focus on class II because the frequency scaling of the self-energy will play important role in non-Fermi liquid behaviours. If Eq. (31) is true, then the fermion self-energy Σ N (iω n ) → 0 as ω n → 0 for η > −1. Hence, from Eq. (24), G 0 (iω n ) ∼ isgn(ω n )/ λ maxφ ∼ isgn(ω n ) for small frequencies at low temperatures. This implies G 0 (τ) ∼ 1/τ after the Fourier transformation. Then the boson self-energy Π(τ) ∼ G 0 (τ)G 0 (−τ) ∼ 1/|τ| 2 . Hence, From the asymptotic expansion in Eq. (29), the Green function for the uncondensed bosons for −1 < η < 0 in the low frequency limit ω m. Also, we used the numerical observation that Π(0) ≈ m 2 /λ max at low temperatures. Since 1+η . Therefore, the frequency scalings in Eqs. (30) and (31) are self-consistent. More detailed numerical and analytical investigations about the scaling properties can be found in Ref. [39].
At low temperatures, the low-frequency dynamics of bosons are important to determine the thermodynamics of the normal state. By approximating the free energy as a saddle point action, i.e., β F = S eff such that δS eff = 0, one can show that the leading contribution to the energy density comes from the effective condensateφ, i.e., E/N ≈ −φ/2 [39]. The effective condensateφ can be approximately calculated from the condition In the third line, the Euler-McLaurin formula was used to approximate the series. We can solve the integral equation with the successive approximation. Letφ =φ 0 + δφ such that i.e.,φ 0 is the effective condensate at zero temperature. At low temperatures, δφ must be small. So we expand Eq. (35) with respect to δφ: where g(φ, ω) = g(φ 0 , ω) + g 1 (φ 0 , ω)δφ + O (δφ) 2 . Then, one can easily read out the temperature dependence of δφ ∼ J (π/β)/β 2 . Although the fermion self-energy |iΣ N (iω n )| ∼ |ω n | 1+η ω n is negligible for class I (η > 0) at low frequencies, |iΣ N (iω n )| ∼ |ω n | 1+η ω n is significant in class II (−1 < η < 0). Hence, in the low frequency limit (|ω n | → 0), which implies Thus, the heat capacity demonstrates the Fermi-liquid-like behaviour for class I and a non-Fermi liquid property of the class II Yukawa-SYK model [39]. While the class I (η > 0) shows conventional linear temperature dependence, the class II (−1 < η < 0) exhibits anomalously large heat capacity at low temperatures because of algebraically diverging ρ(λ) → ∞ as λ → λ max . Note that these two new classes of the normal states are qualitatively different from the quantum critical SYK non-Fermi liquid (SYK-NFL) normal state of the Yukawa-SYK model with a fixed variance coupling, ρ(λ) ∼ δ(λ − λ 0 ) [28,36,41]. The SYK-NFL state is the fast scrambling conformal solution of the fixed variance model, and its non-Fermi liquid nature originates from the strong boson-fermion interactions which dynamically tune the mass of the bosons to zero [28]. The scaling property of the fermion Green function G(iω n ) for the SYK-NFL state is dominated by the fermion self-energy Σ(iω n ) ∼ isgn(ω n )|ω n | 1−2∆ at low frequencies, and the scaling dimension, which lies between 1 4 < ∆ < 1 2 , depends on the ratio between the number of bosons (M ) and fermions (N ), γ = M /N .
On the contrary, the normal states of our model with the distribution of the variances ρ(λ) show different scaling behaviour. Both "Fermi liquid" (η > 0) and "non-Fermi liquid" (−1 < η < 0) states show the "impurity-like" behaviour, which has been observed from the fixed-variance model at intermediate temperature window [36,41]. The strong interactions do not tune the mass of the boson to zero. Instead, the bosons act as impurity centres that scatter fermions. Especially, as we can see from Eq. (24), the impurity-like scattering due to the Bose condensate ϕ and the static uncondensed bosons D N (iν 0 = 0) lead to the leading scaling dimension of the fermions ∆ = 1 2 , i.e., G(iω n ) ∼ isgn(ω n ), for all γ and η > −1. While the impurity-like non-Fermi liquid fixed point is not stable in the fixed-variance models, our model supports the impurity-like Fermi liquid and non-Fermi liquid states as stable infrared fixed points at low temperatures.

Pairing Instabilities of Fermi and Non-Fermi Liquids
We are interested in pairing instabilities of fermions in the presence of the singlet (a = 0) and the triplet (a = 3) Yuakawa interactions [Eq. (4)]. Hence, we consider not only singlet pairing but also triplet pairings. Let us expand the anomalous part of the Green function and the self-energy in the singlet (µ = 0) and the triplet channels (µ = 1, 2, 3): Then Eq. (19) becomes where ζ = 1 if σ a and σ 2 σ µ commutes and ζ = −1 if σ a and σ 2 σ µ anticommutes. Hence, ζ = 1 for all pairing channels (µ = 0, 1, 2, 3) in case of the singlet Yukawa coupling (a = 0). However, ζ = 1 for µ = 1, 2 and ζ = −1 for µ = 0, 3 in case of the triplet Yukawa coupling (a = 3). At the critical temperatures T c , we consider a continuous phase transition to a paired state. Near T c , the anomalous part of the self-energy Φ(iω n ) and the Green function F (iω n ) must be very small. Hence, we linearize the Schwinger-Dyson equations to estimate T c and identify the leading pairing instability. Then we can approximate the anomalous Green function F (iω n ) with the normal state Green function G 0 (iω n ) near T c :

Non-Fermi Liquid Fermi Liquid
In the second line, we used the odd parity of G 0 (iω n ) = −G 0 (−iω n ). Then we get the linearized Schwinger-Dyson equations for the paring channels: Since the bosonic propagator is in the numerator while the fermionic self-energy is in the denominator of Eq. (46), strong Yukawa couplings lead to two competing effects: enhancement of the bosonic propagator D, which is the pairing glue of fermions, and decoherence of fermions due to large fermionic self-energy Σ 0 . Using the bosonic propagator and fermionic self-energy of the normal state, we can calculate the transition temperature T c from the condition that the linearized equation, Eq. (46), has the nontrivial solution. Figure 3 shows the phase diagram of the Yukawa-SYK model for the various distribution parameter η. The phase boundary implies the leading pairing instabilities of the model. While the non-Fermi liquid states with η < 0 (class II) are known to have large fermionic self-energy Σ N (iω n ) ∼ |ω n | 1+η (compared to the free Green function G free (iω n ) −1 ∼ ω n ) due to the uncondensed bosons [39], their transition temperatures are greater than those of the Fermi liquid states with η > 0. Our result implies that the enhancement of the pairing glue D(iν n ) ( Figure 2) plays a more important role in the pairing than the decoherence of fermions in the Yukawa-SYK model.
While the singlet coupling (a = 0) yields the same linearized equations for both singlet (µ = 0) and triplet pairing channels (µ = 1, 2, 3), the triplet Yukawa coupling (a = 3) turns out to have the attractive pairing channels only for the spin-preserving triplet pairing (µ = 1, 2). Note that the spin-preserving triplet pairings are Due to the Pauli exclusion principle, these pairings must be vanishing in the static limit τ → 0.
Only the dynamical pairing among fermions at distinct times can be finite. Therefore, the leading pairing instabilities with the triplet Yukawa coupling (a = 3) correspond to the dynamical pairing of fermions. Such a feature is distinguished from the conventional pairing in the BCS theory. Apart from the nature of the paired states, the transition temperature T c for both the singlet and triplet Yukawa-SYK models are the same for a given value of η. Hence, the phase diagrams for the singlet and triplet couplings are the same although the paired states' nature is different.

Conclusion
In summary, we present a solvable strongly coupled theory of spin-half fermions c iσ interacting with scalar bosons φ k by the all-to-all random Yukawa couplings g i j,k . For each boson φ k , the Yukawa coupling constant g i j,k is sampled from the Gaussian orthogonal ensemble of zero mean, g i j,k = 0, and finite variance, (g i j,k ) 2 = λ k . With a large number of fermions and bosons, we assume that the theory is self-averaging and the set of the variances {λ k } forms a continuous distribution ρ(λ) (Figure 1). An important aspect of the theory is the systematic controllability of the fermionic incoherence with the distribution ρ(λ) responsible for the statistical nature of the Yukawa interaction g i j,k . The model can realize both the Fermi liquid normal state when ρ(λ) is regular at the maximum variance λ max and the non-Fermi liquid normal state when ρ(λ) diverges algebraically at λ max . These Fermi and non-Fermi liquid normal states correspond to the low-energy states of class I and class II low-rank SYK models in Ref. [39].
Starting from these normal states, we examined the leading pairing instabilities in both spinsinglet and triplet channels by solving the linearized Schwinger-Dyson equations. The spin independent Yukawa interactions g i j,k (c † i↑ c j↑ φ k + c † i↓ c j↓ φ k ), which model the charge fluctuations of correlated metals, show the pairing instabilities from both spin singlet and spin triplet channels. However, the spin dependent Yukawa interactions g i j,k (c † i↑ c j↑ φ k − c † i↓ c j↓ φ k ), which represent the spin fluctuations, yield the leading pairing instabilities from the spin triplet channels Although both the spin-independent and dependent Yukawa interactions result in the same normal states, the resulting pairing instabilities are not the same. Furthermore, it is interesting to note that the critical temperature for the pairing state arising from the non-Fermi liquid is higher than that of the Fermi liquid ( Figure 3). Although conventional wisdom may expect that the pairing would be eventually suppressed due to incoherence of the fermions, our theory demonstrates an example that the enhancement of the boson propagator, which glues the fermion pair, dominates the effect of the large fermion self-energy, which shortens each dressed fermion's lifetime. In this theory, there is no ad hoc parameter to control the relative contributions of the boson propagator and fermion self-energy to the pairing instabilities. The control knob of our theory ρ(λ) influences both the enhancement of the pairing glue and the incoherence of the fermions, revealing a concrete physical meaning of the distribution ρ(λ).
Since the Yukawa-SYK model is zero-dimensional, the natural follow-up question is the extension of our work to higher dimensions. If a quantum dot that consists of a large number of bosons and fermions realizes the paired state of the Yukawa-SYK model, we can consider an array of the coupled quantum dots as a higher dimensional generalization of our theory. Then, the leading spin-triplet pairing instabilities from the spin-dependent Yukawa interactions raise an interesting question: can the array of the coupled Yukawa-SYK quantum dots realize any unconventional (topological) superconductor? Furthermore, our analysis is based on the linearized Schwinger-Dyson equations. To examine the thermodynamic properties of the strongly interacting paired states below T c , it would be interesting to explore the solutions of the full nonlinear Schwinger-Dyson equations.

A Derivation of the Effective Action
We derive the effective action by averaging over the random Yukawa couplings, g i j,k . Assuming that the model is self-averaging, we construct the large N effective action from the disorder average of the partition function Z instead of the free energy log Z. In the language of the replica field theory, we are assuming that the replica diagonal terms dominate the low-energy physics while the replica non-diagonal terms are suppressed by O (1/N ).
The summation is assumed for the repeated Greek indices. Therefore where "tr" is the trace over the spin indices. To impose the relationship between the bilocal fields and the fermions and bosons, we introduce the Lagrange multipliers: Let us define the Fourier transformations where ω n = (2n + 1)π/β and ν n = 2nπ/β are the fermionic and bosonic Matsubara frequencies, respectively. Since the model is time-translation invariant, the bilocal fields are functions of τ−τ . The consistent definition of the Fourier transformations for the bilocal fields is Then our modified action S = S c + S φ + S λ including the Lagrange multipliers in the Fourier space is where "Tr" is the trace over the indices for the four-component spinor By integrating out the fermions and bosons, we obtain the effective action S eff = S 0 + S λ in terms of the bilocal fields, where ϕ is the magnitude of the condensed bosons defined in Eq. (13). When the set of the variances {λ k } form a well-defined distribution in the large M limit, we can rewrite S eff as

B Bose-Einstein Condensation for the η > 0 Model
In low-dimensional systems, the violent quantum fluctuations often prevent the presence of longrange order or the Bose-Einstein condensation (BEC). For example, in the Yukawa-SYK model with a fixed variance coupling, i.e., ρ(λ) ∼ δ(λ − λ 0 ), bosons are not condensed although the strong Yukawa interactions with fermions renormalize their mass to zero [28,36]. However, the low dimensionality does not always rule out the possibility of BEC. In this appendix, we demonstrate that our model with η > 0 (class I) can show the Bose-Einstein condensation despite strong quantum fluctuations due to the zero-dimensional nature of the all-to-all interactions. The Bose-Einstein condensation occurs when the number of excited states is bounded from above, i.e., N excited < N 0 , at some temperatures below T BEC . If there were N number of bosons, the remaining N −N 0 number of bosons are forced to be in the ground state since the Bose statistics limits the maximum number of the bosons in the excited states. Then, the macroscopic number (N − N 0 1) of bosons are said to be condensed in the ground state. For our model, it is difficult to calculate the maximum number of available excited states N 0 explicitly. Instead, we can demonstrate that the bosons inevitably become unstable (i.e., m 2 − λ max Π(0) < 0) without the Bose-Einstein condensation when η > 0. If there were no BEC (ϕ = 0), we will first show that the boson Green function D(iν n ) is bounded from above when η > 0. Next, we will prove that the fermion self-energy Σ(iω n ) is bounded from above because D(iν n ) is bounded from above. At last, when the fermion self-energy is bounded from above, we can show that the boson self-energy Π(0) at zero frequency is bounded from below. The lower bound of Π(0) is a function of temperature, and we will show that m 2 − λ max Π(0) inevitably becomes negative at some finite temperatures because of this lower bound. In short, we will prove that BEC must occur in the η > 0 model because "no BEC" implies the bounded boson Green function which results in the unstable bosons. Below, we mathematically demonstrate this idea.
Thus, m 2 −λ max Π(0) < 0 if T < λ max /2m 2 , i.e., the bosons become unstable at finite temperatures if there were no Bose-Einstein condensation. Suppose iΣ 0 (iω n ) for finite γ > 0 is continuously connected to the γ = 0 limit, i.e., if we increase γ from 0, then S also continuously increases from zero. By numerically solving the Schwinger-Dyson equations, we confirmed that the self-consistent Green functions and self-energies for finite γ > 0 is continuously connected to the self-consistent solution for γ = 0 (up to γ = 1. This is also previously confirmed both numerical calculations and analytical perturbation theory in [39]. Since ψ (1) (1/2 + S/2π) is an analytic function of S ≥ 0, the existence of the solution at γ = 0 implies the existence of the solution for some finite γ > 0. To sum up, the Bose-Einstein condensation is necessary for the η > 0 model to avoid the instability of bosons. If the bosons are condensed (ϕ = 0), the bosons can be either critical or stable because the total boson Green function D(iν n ) = (βλ max /γ)ϕδ n,0 + D N (iν n ) is no longer bounded from above. Then the quantum fluctuations of the bosons can yield arbitrarily large fermion self-energy. When the fermion self-energy is sufficiently large, it will not result in large boson self-energy which makes bosons unstable. Note that our conclusion is consistent with the previous work for the fixed-variance model. Since the boson Green function is not bounded from above for the fixed-variance model, the bosons can remain critical without the condensation [28]. We would also like to note that the physics of our model is different from that of free Bose gas at d = 0 because the bosons are strongly coupled to the massless fermions. The intertwined dynamics of the strongly correlated bosons and fermions can result in the Bose condensation even at d = 0 under certain conditions.