Phase diagram of the spin-1/2 Yukawa-SYK model: Non-Fermi liquid, insulator, and superconductor
Wei Wang, Andrew Davis, Gaopei Pan, Yuxuan Wang, Zi Yang Meng
PPhase diagram of the Yukawa-SYK model: Non-Fermi liquid, insulator, and superconductor
Wei Wang,
1, 2, ∗ Andrew Davis, ∗ Gaopei Pan,
1, 2
Yuxuan Wang, † and Zi Yang Meng
4, 1, ‡ Beijing National Laboratory for Condensed Matter Physics and Institute of Physics,Chinese Academy of Sciences, Beijing 100190, China School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China Department of Physics, University of Florida, Gainesville, FL 32601 Department of Physics and HKU-UCAS Joint Institute of Theoretical and Computational Physics,The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China (Dated: February 23, 2021)We analyze the phase diagram of the Yukawa-Sachdev-Ye-Kitaev model, which describes complex fermionsrandomly interacting with real bosons via a Yukawa coupling, at finite temperatures and varying fermion density.In a recent work [1] it has been shown that, upon varying filling or chemical potential, there exists a first-orderquantum phase transition between a non-Fermi liquid (nFL) phase and an insulating phase. Here we showthat in such a model with time-reversal symmetry this quantum phase transition is preempted by a pairingphase that develops as a low-temperature instability. As a remnant of the would-be nFL-insulator transition, thesuperconducting critical temperature rapidly decreases beyond a certain chemical potential. On the other hand,depending on parameters the first-order quantum phase transition extend to finite-temperatures and terminateat a thermal critical point, beyond which the nFL and the insulator become the same phase, similar to thatof the liquid-gas and metal-insulator transition in real materials. We determine the pairing phase boundaryand the location of the thermal critical point via combined analytic and quantum Monte Carlo numeric efforts.Our results provide the model realization of the transition of nFL’s towards superconductivity and insulatingstates, therefore offer a controlled platform for future investigations of the generic phase diagram that hosts nFL,insulator and superconductor and their phase transitions.
Introduction.—
To understand the non-Fermi liquid (nFL)behavior of interacting electron systems is one of the centralissues in modern condensed matter physics. Widely believedto be relevant to the microscopic origin of the strange metalbehavior in unconventional superconudctors [2–9], its theo-retical description remains a challenging issue due to the lackof a small control parameter. Recently, the Sachdev-Ye-Kitaev(SYK) model have garnered widespread attention as it emergesas a new paradigm for the study of nFLs [10–13], which is dif-ferent from most previous research of nFLs, where the systemis usually realized in itinerant fermions coupled to soft bosonicmodes near a quantum critical point [14–21]. The nFL in SYKmodel is exactly solvable in the large- 𝑁 limit. Beyond the con-text of non-Fermi liquids, the SYK model also has been foundto have a hidden holographic connection to quantum blackholes and saturates the limiting rate of scrambling due to itsshort equilibration time [22, 23].Motivated by the aforementioned fermionic systems nearquantum-critical points, recently, a variant of the SYK model,dubbed the Yukawa-SYK model, has been proposed [1, 24–27]. Such a model describes strong random Yukawa couplingbetween 𝑀 𝑁 flavors of fermions and and 𝑁 flavors of bosons.Analytical investigation at large- 𝑁 has revealed a saddle pointsolution in which the Yukawa coupling “self-tunes" the mas-sive bosons to criticality and the fermions form a nFL, whichsaturates the bound on quantum chaos [28, 29]. This saddlepoint solution has been verified at finite 𝑁 and finite 𝑇 viaquantum Monte Carlo (QMC) simulations with an additionalantiunitary time-reversal symmetry [27] by making use of areparametrization symmetry of the large- 𝑁 solution.In this work we focus on the fate of the nFL upon vary-ing temperature and density. Similar to the complex SYK model [30, 31], it has been recently shown that there exist afirst-order quantum phase transition at finite chemical potentialseparating a nFL state and a gapped insulating state [1]. On theother hand, in several versions of the Yukawa-SYK model thenFL phase becomes unstable to a pairing phase [24–26, 32–34]. The natural question then is the corresponding phasediagram of the model spanned by the axes of temperature 𝑇 and chemical potential 𝜇 .Here, by means of large- 𝑁 calculation and unbiased largescale QMC simulation, we present the global phase diagram(see Fig. 1) of the Yukawa-SYK model. Up to a critical valuein the chemical potential 𝜇 , a finite temperature phase transi-tion from nFL to superconductivity is observed. We determinethe pairing transition from solving the linear Eliashberg equa-tion using large- 𝑁 result of the Green’s functions, as well asfinite-size scaling of the pairing susceptibility in QMC simula-tions. We obtain a good agreement between the two methods,indicating the pairing transition is mean-field like. In particu-lar, in the weak coupling limit, we analytically determine thethreshold value for 𝜇 for the superconductor-insulator transi-tion at zero temperature, which agrees well with numericalresults. On the other hand, by solving the Schwinger-Dysonequation, we found the first-order quantum phase transition ex-tends to low temperature and terminates at a (thermal) criticalpoint, which is a generic feature in many metal-insulator tran-sition in correlated materials [35, 36]. However, depending onthe strength of the first-order quantum transition (previouslyfound to be controlled by the ratio 𝑀 / 𝑁 [1]), we show thatthe thermal critical point may be masked by the superconduct-ing phase. The phase diagram obtained offers a controlledplatform for future investigations of phase transitions betweennFL, insulator and superconductor, at generic electron fillings. a r X i v : . [ c ond - m a t . s t r- e l ] F e b (a)(b) N=M Superconductor nFL
Insulator
FIG. 1. (a) Phase diagram of the Yukawa-SYK Model at 𝑁 = 𝑀, 𝜔 = , 𝑚 =
2. From the large- 𝑁 calculation, one sees the nFLbecome superconductor at low temperature in a wide range of chemi-cal potential until the system experiences a first order superconductor-to-insulator phase transition, as denoted by the red points and lines forthe first-order hysteresis region. The thermal critical point that ter-minates the first order transition locates at ( 𝜇 𝑐 = . , 𝑇 𝑐 = . ) .The blue triangles are the transition points from nFL to supercon-ductor obtained from QMC at finite 𝑁, 𝑀 (cf. Fig. 4), whichare consistent with the results obtained from large- 𝑁 calculations(black squares). (b) Phase diagram of the Yukawa-SYK Model at 𝑁 = 𝑀, 𝜔 = , 𝑚 = 𝑁 calculation. The QMC 𝑛 − 𝜇 curves in Fig. 2 are along blue dashed paths. The thermalcritical point at ( 𝜇 𝑐 = . , 𝑇 𝑐 = . ) is denoted by the blackstar. Model.—
The Yukawa-SYK Model we study is describedby the following Hamiltonian, 𝐻 = 𝑀 ∑︁ 𝑖, 𝑗 = 𝑁 ∑︁ 𝛼,𝛽 = ∑︁ 𝑚,𝑛 = ↑ , ↓ ( 𝑖 √ 𝑀 𝑁 𝑡 𝑖𝛼, 𝑗𝛽 𝜙 𝛼𝛽 𝑐 † 𝑖𝛼𝑚 𝜎 𝑧𝑚,𝑛 𝑐 𝑗𝛽𝑛 )+ 𝑁 ∑︁ 𝛼,𝛽 = ( 𝜋 𝛼𝛽 + 𝑚 𝜙 𝛼𝛽 ) − 𝜇 𝑀 ∑︁ 𝑖 = 𝑁 ∑︁ 𝛼 = ∑︁ 𝑚 = ↑ , ↓ 𝑐 † 𝑖𝛼𝑚 𝑐 𝑖𝛼𝑚 , (1)where 𝑐 𝑖𝛼𝑚 ( 𝑐 † 𝑖𝛼𝑚 ) is the annihilation (creation) operatorfor a fermion with flavor 𝛼 and spin 𝑚, 𝑛 ( ↑ or ↓ ).The random Yukawa coupling parameter between fermionand boson is realized as (cid:104) 𝑡 𝑖𝛼, 𝑗𝛽 (cid:105) = (cid:10) 𝑡 𝑖𝛼, 𝑗𝛽 𝑡 𝑘𝛾,𝑙 𝛿 (cid:11) = (cid:0) 𝛿 𝛼𝛾 𝛿 𝑖𝑘 𝛿 𝛽 𝛿 𝛿 𝑗𝑙 + 𝛿 𝛼𝛿 𝛿 𝑖𝑙 𝛿 𝛽𝛾 𝛿 𝑗𝑘 (cid:1) 𝜔 . We set 𝜔 = 𝜋 𝛼𝛽 is thecanonical momentum of 𝜙 𝛼𝛽 . Hermiticity of the model re-quires 𝜙 𝛼𝛽 = − 𝜙 𝛽𝛼 . ( 𝛼, 𝛽 ) are flavor indices which run from1 to 𝑁 and ( 𝑖, 𝑗 ) are quantum dot indices which run from 1 to 𝑀 . 𝜎 𝑧 represent the 𝑧 component of the fermion spin. Dueto time-reversal symmetry, this Hamiltonian is free from thefermion sign problem and can be simulated by QMC at finite 𝑀 and 𝑁 and at finite doping with 𝜇 ≠
0. We prove the absenceof the sign problem and discuss the QMC implementation inSupplemental Material (SM) [37].The main results of this work can be summarized by the tworepresentative phase diagrams. We found that in general thereexist a superconducting dome 𝑇 𝑠𝑐 ( 𝜇 ) in 𝑇 − 𝜇 plane shownin Fig. 1 (the phases for positive and negative 𝜇 are identi-cal by particle-hole symmetry). The vanishing of pairing forlarger 𝜇 is driven by the underlying nFL/insulator transition.In the absence of pairing, we obtain a hysteresis region (thewedge region marked by red lines in Fig. 1) in which bothnFL and insulator states are metastable, divided a first-orderphase transition inside the wedge, similar to that of the liquid-gas transition and metal-insulator transition in many corre-lated materials [35, 36]. The exact location of the first-ordertransition requires comparing the free energy for different so-lutions, which is beyond the scope of the current work. For 𝑁 = 𝑀, 𝜔 = , 𝑚 =
2, the superconducting dome com-pletely preempts the would-be nFL/insulator transition, whilefor 𝑁 = 𝑀, 𝜔 = , 𝑚 =
2, the first-order phase transition isstronger and the corresponding thermal critical point occursoutside of the superconducting phase. In Fig. 1(a) we have alsomarked the superconducting critical temperature obtained byfinite-size scaling from QMC data with blue triangles, whichare consistent with the results obtained from large- 𝑁 calcu-lations denoted by the black squares. The exact position ofthe hysteresis region inside the superconducting phase, andvise versa, requires solving the non-linear superconductinggap equation, and is qualitatively marked by dashed lines re-spectively. Phases in the normal state.—
Within large- 𝑁 , we map outthe 𝑇 - 𝜇 phase diagram by numerically solving the Schwinger-Dyson Eqs. (2). In terms of the propagators 𝐺 − 𝑓 ( 𝑖𝜔 ) = 𝑖𝜔 + 𝜇 + Σ ( 𝑖𝜔 ) for the fermions and 𝐺 − 𝑏 ( 𝑖 Ω ) = Ω + 𝑚 + Π ( 𝑖 Ω ) for the bosons, the Schwinger-Dyson equations are Σ ( 𝑖𝜔 ) = − 𝜔 ∫ 𝑑 Ω 𝜋 𝐺 𝑏 ( 𝑖 Ω ) 𝐺 𝑓 ( 𝑖𝜔 − 𝑖 Ω ) , Π ( 𝑖 Ω ) = 𝑀𝑁 𝜔 ∫ 𝑑𝜔 𝜋 𝐺 𝑓 ( 𝑖 Ω / + 𝑖𝜔 ) 𝐺 𝑓 (− 𝑖 Ω / + 𝑖𝜔 ) . (2)Since we work in the large- 𝑀, 𝑁 regime, and only the ratio 𝑀 / 𝑁 enters the equations [27].We solve Eqs. (2) iteratively, by starting with a simple ansatzfor Σ and Π on the right-hand-sides, obtaining updated val-ues Σ and Π on the right-hand-sides, and repeating until thesolutions { Σ 𝑛 } and { Π 𝑛 } saturate, where 𝑛 is the iteration stepnumber. Noting that Eqs. (2) is consistent with the assump-tions that Π ( 𝑖 Ω ) is even and that the real and imaginary parts FIG. 2. Filling 𝑛 versus 𝜇 for selected 𝑇 in the vicinity of the thermalcritical point in Fig. 1 (b) for 𝑁 = 𝑀, 𝜔 = , 𝑚 =
2. At higher 𝑇 = . 𝑛 ( 𝜇 ) curves both from large- 𝑁 and QMCare smooth. 𝑇 = . , 𝜇 = .
375 (marked by the star), is a thermalcritical point. At temperature 𝑇 = .
083 close to the thermal criticalpoint, there is a sharp turn in of 𝑛 ( 𝜇 ) signifying the divergence ofthe compressibility 𝑑𝑛 / 𝑑𝜇 . At lower 𝑇 , there is a range of 𝜇 —thehysteresis region—where the filling is double-valued. The lowerbranch represents the nFL behavior and the upper branch representsthe insulating behavior; a first order transition connects the two at achemical potential given by a Maxwell construction. The dashed linedelimits the region where solutions are unstable. of Σ ( 𝑖𝜔 ) are even and odd, respectively, we need only com-pute the self-energies at nonnegative frequencies. However,directly implementing this strategy leads to divergent behavior,especially at Σ 𝑛 (± 𝜋𝑇 ) and Π 𝑛 ( ) . This issue is related to thefact that Σ 𝑛 (± 𝜋𝑇 ) and Π 𝑛 ( ) are determined by the behaviorsof 𝐺 and 𝐷 at all frequency scales, rather than their “local"behavior at nearby low energies [1]. Indeed, in analytical solu-tions of Eq. (1) at 𝑇 = Σ ( ) and Π ( ) were used to determine the ultraviolet energy scale beyondwhich nFL behavior crosses over to that of a free system. Toavoid the instability at lowest frequency points in the iterativemethod, we artificially introduce the “stabilizers" for each stepof the iteration by uniformly shifting Σ 𝑛 and Π 𝑛 such that Σ (cid:48) 𝑛 ( 𝜋𝑇 ) = 𝑠, Π (cid:48) 𝑛 ( ) = 𝑝. (3)This prescription prevents Σ 𝑛 ( 𝜋𝑇 ) and Π 𝑛 ( ) from runningaway. Of course, in general after the iteration converges, thesolution we get is not the solution of the original SD equation,unless the updated value at the next step coincides with thestabilizers, i.e., Σ 𝑛 + ( 𝜋𝑇 ) → 𝑠 and Π 𝑛 + ( ) → 𝑝 . Using thiscriterion, we can find the correct values of the stablizers 𝑠 and 𝑝 . The necessary shifts are typically extremely smallcompared to the values of the self-energies over the frequencyrange where most of their support lies.At low 𝑇 and for some ranges of 𝜇 , we obtain two differ-ent choices of stabilizers { 𝑠 , 𝑝 } which cause the iteration toconverge. This signals the hysteresis behavior, and the result- ing two types of solutions physically correspond to nFL andinsulator behaviors that are local minimum of the free energy.This method does not reproduce the unstable ( 𝑑𝑛 / 𝑑𝜇 <
0) so-lutions, whose boundary are sketched qualitatively in Fig. 2.The filling 𝑛 is calculated from the imaginary-time fermionicGreen’s function 𝑛 = + ∑︁ 𝜔 𝑛 𝐺 𝑓 ( 𝑖𝜔 𝑛 ) (4)Alternatively, one can obtain the filling from 𝑛 = lim 𝛿 → 𝐺 ( 𝜏 = − 𝛿 ) . However, due to the finite number offrequency points we keep, 𝐺 ( 𝜏 ) obtained from a Fouriertransform exhibits strong oscillations at small 𝜏 (the Gibbsphenomenon).From the fermion Green’s functions, we extract the 𝑛 − 𝜇 curve, shown in Fig. 2. In particular, for certain values of ( 𝜇, 𝑇 ) the two solutions coexist, indicating the existence ofmetastable states. At low temperatures, there is in general arange of 𝑛 for which no solutions were found. Nevertheless, weexpect in the full solution the 𝑛 ( 𝜇 ) the curve to be smooth. Thismissing portion of solutions (see the dashed lines in Fig. 2)thus correspond to those that cannot be obtained from a sta-ble convergent iterative series. We thus identify this missingportion as thermodynamically unstable saddle points of thefree energy. Such a behavior is typical of first-order phasetransitions. Like water-vapor transition, the actual 𝑛 − 𝜇 curveconnecting the two branches is a straight line determined bythe Maxwell construction. Above a certain temperature 𝑇 𝑐 , thetwo types of solutions become smoothly connected at a chem-ical potential 𝜇 𝑐 . Here the compressibility 𝑑𝑛 / 𝑑𝜇 diverges,and thus ( 𝑇 𝑐 , 𝜇 𝑐 ) is a thermal critical point of the system.Qualitatively, the value of 𝑇 𝑐 up to which the first-orderphase transition survives is related to the strength of quantum first-order transition at 𝑇 =
0. In Ref. [1], it was obtainedanalytically that the first-order quantum phase transition isweaker for a larger ratio 𝑀 / 𝑁 and becomes continuous at 𝑀 / 𝑁 → ∞ . Indeed, we obtain that 𝑇 𝑐 for the case 𝑁 = 𝑀 ishigher than that with 𝑁 = 𝑀 . Pairing transition.—
The interaction mediated by the bosonexchange is attractive in the equal-index, spin-singlet Cooperchannel [27], and the system has a tendency toward a low-temperature pairing phase. The Eliashberg equation is givenby Φ ( 𝜔 𝑛 ) = 𝜔 𝑇 ∑︁ Ω 𝑚 𝐺 𝑏 ( 𝑖 Ω 𝑚 ) 𝐺 𝑓 ( 𝑖 ( 𝜔 𝑛 + Ω 𝑚 )) 𝐺 𝑓 (− 𝑖 ( 𝜔 𝑛 + Ω 𝑚 ))× Φ ( 𝜔 𝑛 + Ω 𝑚 ) . (5)At 𝜇 =
0, the pairing problem has been analyzed in Ref. [27].For 𝜇 ≠
0, due to the breaking of particle-hole symmetry,the mismatch between 𝐺 (± 𝑖𝜔 𝑛 ) leads to a reduced pairingtendancy, much like a Zeeman splitting in momentum spacereduces the spin-singlet pairing susceptibility. We can gleansome insight about the pairing transition by considering theweak-coupling limit 𝜔 (cid:28) 𝑚 and determine the value of 𝜇 𝑠𝑐 beyond which pairing vanishes. At 𝑇 =
0, the Schwinger-Dyson equations admit an insulating solution approximated [1]by Σ ( 𝑖𝜔 ) = − 𝜔 𝐹 /
2, where 𝜔 𝐹 ≡ 𝜔 / 𝑚 , and Π ( 𝑖 Ω ) = 𝜇 > 𝜔 𝐹 /
2. These self-energies become exact in thelimit. In this regime the pairing equation becomes Φ ( 𝜔 ) = 𝜔 ∫ 𝑑 Ω 𝜋 ( Ω − 𝜔 ) + 𝑚 Ω + ( 𝜇 − 𝜔 𝐹 / ) Φ ( Ω ) (6)Most of the support of the integral comes from frequencieson the order of 𝜔 𝐹 , so at very weak coupling the frequencydependence of the boson propagator can be ignored. With theansatz Φ ( 𝜔 ) = const, performing the integral reveals a pairingtransition at 𝜇 𝑠𝑐 = 𝜔 𝐹 (cid:16) ≡ 𝜔 / 𝑚 (cid:17) . (7)To verify this analytic result, we solved the pairing equationnumerically at very weak coupling: 𝑚 = 𝜔 , and we weobtained 𝜇 𝑠𝑐 ≈ . 𝜔 𝐹 for the maximum chemical potentialbeyond which pairing vanishes, and we also found that Φ ( 𝜔 ) is virtually constant, justifying our ansatz. Extrapolating toour case with 𝜔 / 𝑚 = .
5, we expect 𝜇 𝑠𝑐 = .
25. Thisindeed matches well with the numerical results from large- 𝑁 (Fig. 1(a,b)) and from QMC (Fig. 1(a)).Using the numerical solutions for Eqs. (2), the Eliashbergequation can be viewed as a matrix equation | Φ (cid:105) = ˆ 𝐾 ( 𝑇, 𝜇 )| Φ (cid:105) (after imposing a large enough cutoff in frequency) and thelargest eigenvalue of the kernel ˆ 𝐾 can be computed. TheEliashberg equation has a nontrivial solution when the largesteigenvalue reaches 1, indicating the onset of pairing, and wecan map out the boundary of the superconducting region inthe 𝑇 - 𝜇 phase diagram. The numerical results for 𝑁 = 𝑀 , 𝑚 = 𝜔 and 𝑁 = 𝑀 , 𝑚 = 𝜔 are shown in Fig. 1 (a) and(b), respectively. The thermal critical point may lie inside oroutside the superconducting region depending on the ratio of 𝑀 / 𝑁 , as discussed above. Results from QMC.—
To analyze the phase diagram inFig. 1 with QMC, we focus on the Green’s functions andpairing susceptibility obtained in simulations at finite
𝑀, 𝑁 .We first perform the QMC simulations at the parameters of 𝑁 = 𝑀, 𝜔 = , 𝑚 = 𝛽 ≡ / 𝑇 and 𝜇 .We show in Fig. 3 the QMC Green’s functionsfor large 𝛽 with different 𝜇 , with 𝐺 𝑓 ( 𝜏, ) = ( 𝑀 𝑁 ) (cid:205) 𝑀𝑖, 𝑗 = (cid:205) 𝑁𝛼,𝛽 = (cid:104) 𝑐 𝑖,𝛼,𝜎 ( 𝜏 ) 𝑐 † 𝑗,𝛽,𝜎 ( )(cid:105) and 𝐺 𝑏 ( 𝜏, ) = 𝑁 ( 𝑁 − ) (cid:205) 𝑁𝛼,𝛽 = ,𝛼 ≠ 𝛽 (cid:104) 𝜙 𝛼𝛽 ( 𝜏 ) 𝜙 𝛼𝛽 ( )(cid:105) . One can clearly see thatthey exhibit distinct behaviors for small and large 𝜇 , consistentwith the phase diagram in Fig. 1(a). At 𝜇 = . 𝐺 𝑓 and 𝐺 𝑏 decays slowly in imaginary time, similar to the resultsin Ref. [27] exhibiting power-law scaling. Note that, for 𝜇 ≠ 𝐺 𝑓 isnot symmetric with respect to 𝜏 = 𝛽 /
2, and we normalize thedata with respect to 𝐺 𝑓 ( 𝜏 = 𝛽 ) and 𝐺 𝑏 ( 𝜏 = 𝛽 ) . At largerdoping, with 𝜇 = .
35, both 𝐺 𝑓 and 𝐺 𝑏 decay exponentially,consistent with insulating behavior. Since in QMC simulationwith finite 𝑁, 𝑀 , the system does not develop superconduc-tivity, and it is sensible to compare with Green’s functions at (a)(b)
FIG. 3. The QMC fermonic Green’s functions (a) and the bosonicGreen’s functions (b) with different 𝜇 . 𝑀 = 𝑁 = 𝛽 = 𝜔 = 𝑚 =
2, plot with semi-log axes . For convenience boththe 𝐺 𝑓 and 𝐺 𝑏 have been normalized to 1 at 𝜏 = 𝛽 . The systembecome gapped with the increase of the chemical potential. The sharpdownturn of the large- 𝑁 result in panel (a) is an artifact of keepingfinite frequency points. large- 𝑁 for the normal state. As can be seen in Fig. 3, theagreement is excellent.For 𝑀 = 𝑁 case in Fig. 1 (b), the thermal critical pointis located around ( 𝜇 𝑐 = . , 𝑇 𝑐 = . ) , which is withinreach of our QMC simulations. We compute the 𝑛 ( 𝜇 ) curves(whose derivative is the charge compressibility) near and faraway from 𝑇 𝑐 . We can see in Fig. 2 that, in excellent agree-ment with the large- 𝑁 solution, the compressibility is con-stant when the temperature is much higher than the criti-cal point ( 𝑀 = 𝑁 = , 𝑇 = . 𝑛 ( 𝜇 ) when the temperature is close to the critical point( 𝑀 = 𝑁 = , , , 𝑇 = . 𝑁 .For 𝑁 = 𝑀 our QMC results further reveal that the nFLdevelop a superconductivity at low temperature in a widerange of chemical potential, reaching beyond the would-befirst order phase transition. To extract the superconduct-ing transition temperature, we measure the pairing correla-tion in our QMC simulation, and analyze its scaling behav-ior as system size. The pair susceptibility is expressed as (a)(b) (c)(d) (e)(f) (g)(h) (i)(j) FIG. 4. Pair susceptibility 𝑃 𝑠 measured at different chemical potential 𝜇 . The obtained 𝑇 𝑐 ( 𝛽 𝑐 ) are denoted by the blue triangles in Fig. 1(a). From the temperature dependence of the 𝑃 𝑠 with different system sizes ( 𝑀, 𝑁 ) we perform the data collapse using mean-field exponents 𝛾 = 𝜈 = 𝑇 𝑐 ( 𝛽 𝑐 ) are obtained accordingly. The parameters are: 𝑁 = 𝑀 , 𝜔 = 𝑚 = 𝜇 = 𝜇 = .
05 in (c) and (d); 𝜇 = .
125 in (e) and (f); 𝜇 = . 𝜇 = .
25 in (i); 𝜇 = .
35 in (j). The superconducting transitiontemperature reduces as 𝜇 increases. For 𝜇 = .
25 (i) and 𝜇 = .
35 (j) the pairing susceptibilities are not divergent at larger 𝑁 , and the systementer a gapped insulator phase. 𝑃 𝑠 = ∫ 𝛽 d 𝜏 (cid:104) Δ ( 𝜏 ) Δ † ( )(cid:105) , where Δ is the pairing field de-fined as Δ † = √ 𝑀 𝑁 (cid:205) 𝑀𝑖 = (cid:205) 𝑁𝛼 = 𝑐 † 𝑖,𝛼, ↑ 𝑐 † 𝑖,𝛼, ↓ . At finite 𝑁 , thepairing susceptibility 𝑃 𝑠 does not diverge, and can be writtenas 𝑃 ( 𝑁 ) 𝑠 ∼ 𝑁 𝑎 𝑓 (cid:2) 𝑁 / 𝜈 ( 𝑇 − 𝑇 𝑠𝑐 ) (cid:3) , in which 𝑁 (and 𝑀 for afixed ratio) plays the role of the system size [9, 38–40]. Forour large- 𝑁 system without the notion of space, the role ofcorrelation length is replaced with a correlation “cluster size" 𝑛 ∼ ( 𝑇 − 𝑇 𝑠𝑐 ) − 𝜈 , and hence the functional dependence of 𝑓 ( 𝑥 ) .In the large- 𝑁 limit, all fluctuation effects are suppressed by1 / 𝑁 , and such a phase transition is mean-field like [24]. Thismeans that for a fixed 𝑇 − 𝑇 𝑐 the exponent 𝜈 = 𝑓 ( 𝑥 ) ∼ / 𝑥 .Further requiring that in the large- 𝑁 (thermodynamic) limitthe susceptibility diverges independent of 𝑁 , we obtain that 𝑎 = /
2. Using these exponents, we indeed obtain decent fi-nite size scaling with a 𝛽 𝑠𝑐 by data collapse, see Fig. 4 (a) ∼ (h)with different fermion densities (different 𝜇 ). We can see thatwhen the 𝜇 increase the superconducting transition tempera-ture is moderately reduced until a sudden drop at larger 𝜇 . For 𝜇 > .
25 the pairing susceptibility no longer diverges withlarge 𝑁 , and the system does not form a pairing state (seeFig. 4(i) and (j)). The corresponding QMC 𝑇 𝑠𝑐 points areshown in the Fig. 1 (a). The values of 𝑇 𝑠𝑐 ( 𝜇 ) from QMC arelarger than their large- 𝑁 counterparts, but they are close. Inparticular, the values of 𝜇 𝑠𝑐 from QMC and large- 𝑁 are ingood agreement, consistent with analytical result Eq. (7). Discussion.—
With combined analytical and numerical ef-forts, we reveal the 𝑇 − 𝜇 phase diagram of the Yukawa-SYKmodel. We identified that an underlying first-order quantumphase transition between a non-Fermi liquid and an insulatorleads to a dome-like structure of the pairing phase, and de-pending on the parameter 𝑁 / 𝑀 , survives at finite- 𝑇 until asecond-order thermal tri-critical point between non-Fermi liq-uid and an insualtor. Our results provide the model realizationof the SYK-type nFL and its transition towards superconduc- tivity and insulating states, therefore offer a controlled platformfor future investigations of the generic phase diagram that hostsnFL, insulator and superconductor phases and their transitionsat generic fermion densities. Acknowledgments — We thank Andrey Chubukov, Ilya Es-terlis, Yingfei Gu, Grigory Tarnopolsky, Joerg Schmalian,Subir Sachdev, Steven Kivelson for insightful discussions. ADand YW are supported by startup funds at the University ofFlorida. WW, GPP and ZYM acknowledge support from theRGC of Hong Kong SAR of China (Grant Nos. 17303019and 17301420), MOST through the National Key Researchand Development Program (Grant No. 2016YFA0300502)and the Strategic Priority Research Program of the ChineseAcademy of Sciences (Grant No. XDB33000000). We thankthe Computational Initiative at the Faculty of Science and theInformation Technology Services at the University of HongKong and Tianhe-2 platform at the National SupercomputerCenters in Guangzhou for their technical support and gener-ous allocation of CPU time. ∗ These two authors contributed equally. † yuxuan.wang@ufl.edu ‡ [email protected][1] Y. Wang and A. V. Chubukov, Phys. Rev. Research , 033084(2020).[2] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, andJ. Zaanen, Nature , 179 (2015).[3] Z. Liu, Y. Gu, W. Zhang, D. Gong, W. Zhang, T. Xie, X. Lu,X. Ma, X. Zhang, R. Zhang, J. Zhu, C. Ren, L. Shan, X. Qiu,P. Dai, Y.-f. Yang, H. Luo, and S. Li, Phys. Rev. Lett. ,157002 (2016).[4] Y. Gu, Z. Liu, T. Xie, W. Zhang, D. Gong, D. Hu, X. Ma, C. Li,L. Zhao, L. Lin, Z. Xu, G. Tan, G. Chen, Z. Y. Meng, Y.-f. Yang,H. Luo, and S. Li, Phys. Rev. Lett. , 157001 (2017).[5] J. Custers, P. Gegenwart, H. Wilhelm, K. Neumaier, Y. Tokiwa, O. Trovarelli, C. Geibel, F. Steglich, C. Pépin, and P. Coleman,Nature , 524 (2003).[6] B. Shen, Y. Zhang, Y. Komijani, M. Nicklas, R. Borth, A. Wang,Y. Chen, Z. Nie, R. Li, X. Lu, H. Lee, M. Smidman, F. Steglich,P. Coleman, and H. Yuan, Nature , 51 (2020).[7] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda,K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero,Phys. Rev. Lett. , 076801 (2020).[8] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang,J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng,D. Shi, O. V. Yazyev, and G. Zhang, Nature Physics (2020),10.1038/s41567-020-0825-9.[9] C. Chen, T. Yuan, Y. Qi, and Z. Y. Meng, arXiv e-prints ,arXiv:2007.05543 (2020), arXiv:2007.05543 [cond-mat.str-el].[10] S. Sachdev and J. Ye, Phys. Rev. Lett. , 3339 (2015).[11] A. Kitaev, Entanglement in Strongly-Correlated Quantum Mat-ter .[12] S. Sachdev, Phys. Rev. X , 041025 (2015).[13] A. Kitaev and S. J. Suh, Journal of High Energy Physics ,183 (2018).[14] A. Abanov, A. Chubukov, and J. Schmalian, Advances inPhysics , 119 (2003).[15] M. A. Metlitski and S. Sachdev, Phys. Rev. B , 075127 (2010).[16] M. A. Metlitski and S. Sachdev, Phys. Rev. B , 075128 (2010).[17] Z. H. Liu, X. Y. Xu, Y. Qi, K. Sun, and Z. Y. Meng, Phys. Rev.B , 045116 (2018).[18] Z. H. Liu, G. Pan, X. Y. Xu, K. Sun, and Z. Y. Meng, Proceedingsof the National Academy of Sciences , 16760 (2019).[19] X. Y. Xu, Z. H. Liu, G. Pan, Y. Qi, K. Sun, and Z. Y. Meng,Journal of Physics: Condensed Matter , 463001 (2019).[20] X. Y. Xu, A. Klein, K. Sun, A. V. Chubukov, and Z. Y. Meng,npj Quantum Materials , 65 (2020).[21] J. A. Damia, M. Solís, and G. Torroba, Phys. Rev. B , 045147(2020).[22] H. Guo, Y. Gu, and S. Sachdev, Annals of Physics , 168202(2020).[23] Y. Gu, X.-L. Qi, and D. Stanford, Journal of High EnergyPhysics , 125 (2017).[24] Y. Wang, Phys. Rev. Lett. , 017002 (2020).[25] I. Esterlis and J. Schmalian, Phys. Rev. B , 115132 (2019).[26] D. Hauck, M. J. Klug, I. Esterlis, and J. Schmalian, Annals ofPhysics , 168120 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.[27] G. Pan, W. Wang, A. Davis, Y. Wang, and Z. Y. Meng, arXive-prints , arXiv:2001.06586 (2020), arXiv:2001.06586 [cond-mat.str-el].[28] J. Kim, X. Cao, and E. Altman, Phys. Rev. B , 125112(2020), arXiv:1910.10173 [cond-mat.str-el].[29] J. Kim, E. Altman, and X. Cao, arXiv e-prints ,arXiv:2010.10545 (2020), arXiv:2010.10545 [cond-mat.str-el].[30] T. Azeyanagi, F. Ferrari, and F. I. S. Massolo, Phys. Rev. Lett. , 061602 (2018).[31] R. Smit, D. Valentinis, J. Schmalian, and P. Kopietz, arXivpreprint arXiv:2010.01142 (2020).[32] H. Wang, A. L. Chudnovskiy, A. Gorsky, and A. Kamenev,Phys. Rev. Research , 033025 (2020).[33] C. Setty, Phys. Rev. B , 184506 (2020).[34] Y. Cheipesh, A. I. Pavlov, V. Scopelliti, J. Tworzydło, and N. V.Gnezdilov, Phys. Rev. B , 220506 (2019).[35] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. ,1039 (1998).[36] P. Limelette, A. Georges, D. Jérome, P. Wzietek, P. Metcalf,and J. M. Honig, Science , 89 (2003).[37] QMC implementation of the model and the proof of sign- problem free at arbitrary chemical potential are presented inthis Supplemental Material.[38] S. V. Isakov and R. Moessner, Phys. Rev. B , 104409 (2003).[39] T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer,Phys. Rev. B , 184501 (2004).[40] N. C. Costa, T. Blommel, W.-T. Chiu, G. Batrouni, and R. T.Scalettar, Phys. Rev. Lett. , 187003 (2018).[41] For a mean field theory in finite spatial dimensions, Josephson’sidentity states that 𝜈𝑑 =
2, and in our model the the dimensional-ity 𝑑 does not enter the theory. We have instead 𝜈 =
2. We thankIlya Esterlis and Joerg Schmalian for sharing their unpublishedresults with us on this.
Supplemental Material for "Phase diagram of the Yukawa-SYK model: Non-Fermi liquid, insulator, andsuperconductor"
I. MODEL AND QUANTUM MONTE CARLO SIMULATION
The Hamiltonian of Eq. (1) in the main text is illustrated by Fig. S1. There are 𝑖, 𝑗 = , · · · , 𝑀 quantum dots, each dot acquires 𝛼, 𝛽 = , · · · , 𝑁 flavors of fermions. Fermions are Yukawa coupled via the random hopping 𝑡 𝑖 𝑗 and anti-symmetric bosonicfield 𝜙 𝛼𝛽 . The system go through a phase transition from non-Fermi liquid to a pairing state when the temperature is below thesuperconducting critical temperature ( 𝑇 𝑠𝑐 ). FIG. S1. Yukawa-SYK model. There are 𝑖, 𝑗 = , · · · , 𝑀 quantum dots, each dot acquires 𝛼, 𝛽 = , · · · , 𝑁 flavors. Fermions are Yukawacoupled via the random hopping 𝑡 𝑖 𝑗 and anti-symmetric bosonic field 𝜙 𝛼𝛽 . The system goes through a phase transition from non-Fermi liquidto a pairing state when the temperature is below the superconducting critical temperature ( 𝑇 𝑠𝑐 ). We use the DQMC method to simulate this Hamiltonian, and the starting point is the partition function of the system 𝑍 = Tr (cid:110) 𝑒 − 𝛽 ˆ 𝐻 (cid:111) = Tr (cid:26)(cid:16) 𝑒 − Δ 𝜏 ˆ 𝐻 (cid:17) 𝐿 𝜏 (cid:27) = ∑︁ { Φ } Tr F (cid:104) Φ | 𝑒 − Δ 𝜏𝐻 | Φ 𝐿 𝜏 (cid:105)(cid:104) Φ 𝐿 𝜏 | 𝑒 − Δ 𝜏𝐻 | Φ 𝐿 𝜏 − (cid:105) · · · (cid:104) Φ | 𝑒 − Δ 𝜏𝐻 | Φ (cid:105) (S1)where we divide the imaginary time axis 𝛽 into 𝐿 𝜏 slices, then we have 𝛽 = 𝐿 𝜏 × Δ 𝜏 . Using Trotter-Suzuki decomposition to theHamiltonian in Eq. (S1), 𝑒 − Δ 𝜏 ˆ 𝐻 ≈ 𝑒 − Δ 𝜏 ˆ 𝐻 𝑓 𝑏 𝑒 − Δ 𝜏 ˆ 𝐻 𝑏 (S2)where 𝐻 𝑓 𝑏 is the boson-fermion term, 𝐻 𝑏 is the boson term in the Hamiltonian.Then the partition function can be written as 𝑍 = ∑︁ { Φ } 𝜔 𝐵 [ 𝜙 ] 𝜔 𝐹 [ 𝜙 ] . (S3)As for the bosonic part of the partition function, 𝜔 𝑏 [ 𝜙 ] = 𝐶 𝐿 𝜏 (cid:169)(cid:173)(cid:171) 𝐿 𝜏 (cid:214) 𝑙 = 𝑁 (cid:214) 𝛼,𝛽 = 𝑒 − Δ 𝜏 𝑚 𝜙 𝛼𝛽,𝑙 (cid:170)(cid:174)(cid:172) · (cid:169)(cid:173)(cid:171) (cid:214) (cid:104) 𝑙,𝑙 (cid:48) (cid:105) 𝑁 (cid:214) 𝛼,𝛽 = 𝑒 − (cid:16) 𝜙𝛼𝛽,𝑙 − 𝜙𝛼𝛽,𝑙 (cid:48) (cid:17) Δ 𝜏 (cid:170)(cid:174)(cid:172) (S4) with (cid:104) 𝑙, 𝑙 (cid:48) (cid:105) stands for the nearest-neighbor interaction in imaginary time direction. As for the fermionic part of the partitionfunction 𝜔 𝐹 [ 𝜙 ] = det [ 𝐼 + B 𝐿 𝜏 B 𝐿 𝜏 − · · · B 𝑙 · · · B B ] , (S5)where 𝐵 𝑙 = 𝑒 − Δ 𝜏𝑉 ( Φ 𝑙 ) (S6)and 𝑉 ( Φ 𝑙 ) = 𝑖 √ 𝑀 𝑁 𝜎 𝑧 × ⊗ (cid:0) 𝑡 𝑖 𝑗 (cid:1) 𝑀 × 𝑀 ⊗ (cid:0) 𝜙 𝛼𝛽,𝑙 (cid:1) 𝑁 × 𝑁 (S7)With these notations prepared, finally the partition function in Eq. (S3) can be written as 𝑍 = ∑︁ { Φ } 𝐿 𝜏 (cid:214) 𝑙 = 𝐶 𝐿 𝜏 (cid:169)(cid:173)(cid:171) 𝐿 𝜏 (cid:214) 𝑙 = 𝑁 (cid:214) 𝛼,𝛽 = 𝑒 − Δ 𝜏 𝑀𝑁 𝑚 𝜙 𝛼𝛽,𝑙 (cid:170)(cid:174)(cid:172) (cid:169)(cid:173)(cid:171) (cid:214) (cid:104) 𝑙,𝑙 (cid:48) (cid:105) 𝑁 (cid:214) 𝛼,𝛽 = 𝑒 − ( 𝜙𝛼𝛽,𝑙 − 𝜙𝛼𝛽,𝑙 (cid:48) ) Δ 𝜏 (cid:170)(cid:174)(cid:172) Det (cid:2) I + B 𝐿 𝜏 B 𝐿 𝜏 − · · · B 𝑙 · · · B B (cid:3) (S8)This partition function is free from the minus-sign problem with any 𝜇 . For the part of boson-fermion term of the Hamiltonian,it is invariant under a time-reversal symmetry operation T = 𝑖𝜎 𝑦 K . The boson-fermion term of the Hamiltonian can be writtenas ˆ 𝐻 𝑓 𝑏 = 𝑁 ∑︁ 𝑖 𝑗 = 𝑁 ∑︁ 𝛼,𝛽 = 𝑖 √ 𝑀 𝑁 𝑡 𝛼𝛽 𝜙 𝑖 𝑗 𝑐 † 𝑖𝛼 ↑ 𝑐 𝑗 𝛼 ↑ − 𝜇𝑐 † 𝑖𝛼 ↑ 𝑐 𝑖𝛼 ↑ − 𝑖 √ 𝑀 𝑁 𝑡 𝛼𝛽 𝜙 𝑖 𝑗 𝑐 † 𝑖𝛼 ↓ 𝑐 𝑗 𝛼 ↓ − 𝜇𝑐 † 𝑖𝛼 ↓ 𝑐 𝑖𝛼 ↓ (S9)Under the transformation of T = 𝑖𝜎 𝑦 K , we have T 𝐻 𝑓 𝑏 T − = 𝑁 ∑︁ 𝑖 𝑗 = 𝑁 ∑︁ 𝛼,𝛽 = − 𝑖 √ 𝑀 𝑁 𝑡 𝛼𝛽 𝜙 𝑖 𝑗 𝑐 † 𝑖𝛼 ↓ 𝑐 𝑗 𝛼 ↓ − 𝜇𝑐 † 𝑖𝛼 ↓ 𝑐 𝑖𝛼 ↓ + 𝑖 √ 𝑀 𝑁 𝑡 𝛼𝛽 𝜙 𝑖 𝑗 𝑐 † 𝑖𝛼 ↑ 𝑐 𝑗 𝛼 ↑ − 𝜇𝑐 † 𝑖𝛼 ↑ 𝑐 𝑖𝛼 ↑ = 𝐻 𝑓 𝑏 (S10)At the same time for the fermion determinantdet [ + 𝐵 ( 𝛽, )] = det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105) det (cid:104) + 𝐵 ↓ ( 𝛽, ) (cid:105) = det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105) det (cid:104) T (cid:16) + 𝐵 ↓ ( 𝛽, ) (cid:17) T − (cid:105) ∗ = det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105) det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105) ∗ (S11) = (cid:12)(cid:12)(cid:12) det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105)(cid:12)(cid:12)(cid:12) (S12)where 𝐵 ( 𝛽, ) = B 𝐿 𝜏 B 𝐿 𝜏 − · · · B 𝑙 · · · B B . The determinant is a positive and real number. Also for the boson part of the weight 𝜔 𝑏 [ 𝜙 ] is positive and real. Therefore we have proved that this Hamiltonian is sign problem free.From a simpler viewpoint, just by looking at the matrix elements of matrices corresponding to different spins: 𝐵 ↑ ( 𝛽, ) and 𝐵 ↓ ( 𝛽, ) , we can see this model doesn’t have sign problem. Every element of 1 + 𝐵 ↓ ( 𝛽, ) is individually complex conjugatewith the corresponding element of 1 + 𝐵 ↑ ( 𝛽, ) , which means :det (cid:104) + 𝐵 ↑ ( 𝛽, ) (cid:105) = det (cid:104) + 𝐵 ↓ ( 𝛽, ) (cid:105) ∗ (S13) II. THE CRITICAL TEMPERATURE OF SUPERCONDUCTING WITH DIFFERENT 𝑚 The influence of the ratio 𝜔 / 𝑚 to the superconducting transition temperature has been studied quantitativly by large- 𝑁 limit calculation. The inverse transition temperature 𝛽 𝑐 from nFL to superconductivity as a function of the ratio 𝜔 / 𝑚 for 𝑁 = 𝑀 and 𝑁 = 𝑀 are discussed in the main text here and in the Ref. [27]. By QMC simulation ,we also obtained that when 𝑚 = , 𝜔 = 𝜇 =
0, the superconducting transition temperature is around 𝑇 𝑠𝑐 ∼ . ( 𝛽 𝑠𝑐 ∼ ) , which is higher than thesuperconducting transition temperature around 𝑇 𝑠𝑐 ∼ . ( 𝛽 𝑠𝑐 ∼ ) at 𝜔 = , 𝑚 = 𝜇 =
0. The results are shown inFig. S2 and consistent with theoretical analysis in Ref. [27]. (a)(b)
FIG. S2. Pair susceptibility for different 𝛽 and data collapse, 𝜔 = 𝑚 = 𝜇 =
0. The superconducting transition temperature is around 𝑇 𝑠𝑐 ∼ . ( 𝛽 𝑐 ∼ ) , which is higher than the superconducting transition temperature around 𝑇 𝑠𝑐 ∼ . ( 𝛽 𝑠𝑐 ∼ ) at 𝜔 = , 𝑚 = 𝜇 ==