Mechanism of the insulator-to-metal transition and superconductivity in the spin liquid candidate NaYbSe2 under pressure

The quantum spin liquid candidate NaYbSe2 was recently reported to exhibit a Mott transition under pressure. Superconductivity was observed in the high-pressure metallic phase, raising the question concerning its relation with the low-pressure quantum spin liquid ground state. Here we combine the density functional theory and the dynamical mean-field theory to explore the underlying mechanism of the insulator-to-metal transition and superconductivity and establish an overall picture of its electronic phases under pressure. Our results suggest that NaYbSe2 is a charge-transfer insulator at ambient pressure. Upon increasing pressure, however, the system first enters a semi-metallic state with incoherent Kondo scattering against coexisting localized Yb-4f moments, and then turns into a heavy-fermion metal. In between, there may exist a delocalization quantum critical point responsible for the observed non-Fermi liquid region with linear-in-T resistivity. The insulator-to-metal transition is therefore a two-stage process. Superconductivity emerges in the heavy-fermion phase with well-nested Yb-4f Fermi surfaces, suggesting that spin fluctuations may play a role in the Cooper pairing. NaYbSe2 might therefore be the 3rd Yb-based heavy-fermion superconductor with a very “high” Tc than most heavy-fermion superconductors.

Recently, superconductivity has been reported by pressurizing the QSL candidate NaYbSe 2 19,20 . In contrast to YbMgGaO 4 which has the same structure (space group R-3mH) but with Mg/Ga disorders [21][22][23][24][25][26][27][28][29][30] , NaYbSe 2 contains no intrinsic site disorder and has arguably the simplest crystal structure and chemical formula among all existing QSL candidates 31 . At ambient pressure, analyses of transport and neutron data found an insulating charge gap of 1.9 eV 32 and a crystalline excitation energy of about 15.7 meV [33][34][35] . No magnetic or structural transition was observed down to 50 mK [32][33][34] . Neutron scattering measurement has revealed the signature of magnetic excitation continuum and supported a QSL ground state 34 . The low-energy physics may thus be described by an effective spin- 1 2 model on a perfect triangular lattice formed by Yb ions. At 11 GPa, a structural transition was reported 19 . The high-pressure structure has a lower symmetry (space group P-3m1) with two inequivalent Yb ions of different Yb-Se distances 19 . Upon further increasing pressure, an insulator-tometal transition appears at about 58.9 GPa and superconductivity emerges over a wide pressure range above 103.4 GPa with a maximum T c of about 8 K 19 . The insulator-to-metal transition was claimed to be of Mott type, raising the question concerning the relationship between the low-pressure QSL and the high-pressure superconductivity.
To clarify the nature of the "Mott" transition and the superconducting pairing mechanism, we investigate here the electronic structure of NaYbSe 2 and its pressure evolution using the density functional theory (DFT) plus the dynamical mean-field theory (DMFT) calculations [36][37][38][39] . We find that NaYbSe 2 is a charge-transfer insulator at ambient pressure, with the charge gap formed between Se-3p bands and Yb-4f upper-Hubbard bands. The charge gap is closed near about 60 GPa, giving rise to an intermediate semi-metallic phase with small electron and hole pockets from Yb-5d and Se-3p bands, respectively. The Yb-4f electrons remain localized at this pressure and the system is a typical low carrier density Kondo system, which explains the logarithmic temperature dependence of the resistivity observed in transport measurements. Above 73 GPa, the Yb-4f moments are fully screened on one (Yb1) of the two inequivalent Yb-ions, giving rise to flat heavy electron bands near the Fermi energy with sharp quasiparticle peaks in the density of states. The system turns into a heavy-fermion metal. But the 4f electrons on the other Yb-ions (Yb2) remain localized and only weakly hybridized with conduction electrons. Superconductivity emerges in the heavy-fermion phase when the Yb1-4f Fermi surfaces develop a well-nested quasi-two-dimensional structure. The nesting wave vector coincides with the antiferromagnetic wave vector of the localized spins, indicating possible involvement of heavy electrons and antiferromagnetic spin fluctuations. NaYbSe 2 might therefore be the 3rd Yb-based heavy-fermion superconductor with a rather "high" transition temperature (T c~8 K) among all heavy-fermion superconductors.

RESULTS AND DISCUSSION Optimization of crystal structures
The crystal structures of NaYbSe 2 at low and high pressures are compared in Fig. 1a, both composed of edge-shared YbSe 6 and NaSe 6 octahedra. The Yb ions form a layer structure of flat triangular lattices with the shortest Yb-Yb distance in the ab plane for both structures 19 . The lattice parameters have been measured below 32 GPa in experiment 19,20 . For numerical calculations, we have used DFT+U to optimize the lattice parameters for both structures [40][41][42] and the results are compared with existing experimental data in Fig. 1b. Their good agreement over the whole low-pressure range validates the starting point of our numerical calculations.

Electronic structures at low pressures
The electronic structures of NaYbSe 2 were obtained using DFT +DMFT with the hybridization expansion continuous-time quantum Monte Carlo (CT-HYB) as the impurity solver 43,44 , taking into account both the spin-orbit coupling and electronic correlations. A Coulomb interaction U = 6 eV was chosen to reproduce the correct charge gap of about 1.9 eV at ambient pressure estimated from absorption spectra measurements 32 . Figure 1c plots the calculated density of states and spectral function at T = 10 K and ambient pressure. The Yb-4f Hubbard bands are located at about −3.0 and 2.0 eV, respectively, with a clear charge gap between the Se-3p bands and the Yb-4f upper-Hubbard bands. Hence NaYbSe 2 should be identified as a charge-transfer insulator at ambient pressure. Without DMFT, DFT+U calculations alone yield a much larger Yb occupation number and predict incorrectly a metallic ground state, in disagreement with experiment. The band structures at 20 GPa are compared in Fig. 1d. We have now two inequivalent Yb-ions. Their 4f Hubbard bands become blurred, and the Yb-5d bands move down closer to the Fermi energy to form a smaller normal band insulating gap with the slightly upshifted Se-3p bands.

The intermediate semi-metallic state
Upon further increasing pressure, our calculations reveal a sequence of electronic phases, in good correspondence with the experimental observation. This is summarized in the upmost plot of Fig. 2. The first transition occurs near about 60 GPa, where the charge gap is closed and the spectral function reveals a semimetallic state with a small electron pocket at the M point from Yb-5d bands and a hole pocket at Γ from Se-3p bands, as shown in Fig. 2a. However, the Yb-4f electrons are still localized at this pressure. In Fig. 2d, we see a tiny pseudogap structure near the Fermi energy in the 4f partial density of states, indicating a weak hybridization with conduction bands. Thus NaYbSe 2 should be better identified as a low carrier density Kondo system at this  19,20 . c, d Orbital-resolved density of states and the spectral function at ambient pressure and 20 GPa obtained by DFT+DMFT at 10 K.
pressure. Consequently, there are not enough conduction electrons to screen the Yb-4f local moments so that the resistivity should be dominated by incoherent Kondo scattering and exhibit typical logarithmic temperature dependence 45,46 . This has indeed been observed in experiment between 58.9 and 74.8 GPa, but was initially attributed to weak localization 19 . We argue that similar properties, including the logarithmic-in-T resistivity and a pseudogap in the quasiparticle density of states, have previously been discussed in other low carrier density Kondo systems [45][46][47][48] . Since no magnetic long-range order has been observed in experiment under pressure 19 , the weak hybridization in this regime suggests that the 4f spins on both Yb ions should still be in their respective spin liquid states as reported by neutron scattering at ambient pressure 34 .
Heavy-fermion metal at high pressure The second transition occurs close to 75 GPa beyond which flat heavy electron bands are seen to develop near the Fermi energy in Fig. 2b, c. Correspondingly, a sharp quasiparticle peak appears at the Fermi energy in the Yb1-4f density of states shown in Fig.  2e, f due to the many-body Kondo effect. The system now turns into a good heavy-fermion metal. To determine the exact transition point, we plot in the inset of Fig. 2e the full width at half maximum (FWHM) of the Yb1-4f quasiparticle peak. Its variation with pressure reflects tentatively the strength of the Kondo hybridization 49 . A linear extrapolation to zero suggests that the transition occurs at about 73 GPa, which separates the low carrier density Kondo phase and the heavy-fermion metal. Thus, it marks a delocalization transition of the Yb1-4f electrons. Interestingly, experiment did observe a non-Fermi liquid state with linearin-T resistivity near 74.8 GPa and beyond that a gradual crossover to the Fermi liquid at 126 GPa 19 . This agrees with our theory and strongly suggests the possible existence of a delocalization quantum critical point at zero temperature. The nature of the non-Fermi liquid phase is still under debate 50,51 . However, the broad crossover region is quite unusual and might be associated with the special layer structure of NaYbSe 2 to be discussed later. Nevertheless, the insulator-to-metal transition is not a simple Mott transition, but a two-stage process first from a normal band insulator to a semimetal with incoherent Kondo scattering and then to a heavy-fermion metal. Similar transition and non-Fermi liquid property have also been reported in the low carrier density Kondo system CeNi 2−δ (As 1−x P x ) 2 46 . It may be interesting to mention that pressure seems to favor the heavy-fermion phase in NaYbSe 2 , while in many other Ybbased heavy-fermion compounds, the Kondo effect is typically weakened under pressure and gives in to long-range magnetic orders because of the upshift of the Yb-4f hole (upper-Hubbard) bands further away from the Fermi energy, suppressing possible charge fluctuations and making Yb-4f electrons more localized. In NaYbSe 2 , however, our results show that the localized Yb-4f levels are already away from the Fermi energy with a stable Yb 3+ configuration and the major effect of pressure is to increase the carrier density and enhance the hybridization as will be discussed below. Such anomalous increase of the Kondo effect with pressure has previously been observed in Yb 2 Ni 2 Al and YbCu 2 Si 2 52,53 but attributed to the electron-lattice coupling 54 . Here NaYbSe 2 provides an interesting alternative for further investigation.

Alternating heavy-fermion and QSL layers
To have a better understanding of the peculiarity of NaYbSe 2 , we compare in Fig. 3a, b the imaginary part of the Yb-4f hybridization function at four typical pressures. Quite remarkably, this reveals an unexpected but substantial difference between the two inequivalent Yb ions. While no hybridization is observed near the Fermi energy on both ions at low pressure (20 GPa), it begins to develop on Yb1 at 60 GPa and becomes large above 75 GPa, but on Yb2 the hybridization remains weak over the whole pressure range and shows no significant enhancement. These suggest that only Yb1 ions are responsible for the flat bands and heavy-fermion properties at high pressure, while Yb2-4f electrons remain localized and weakly interacting with conduction electrons. These observations are further confirmed in the local spin susceptibility, χ zz ðωÞ ¼ hS z ðτÞS z ð0Þi ω , calculated using the CT-HYB impurity solver. Here S z is the local spin operator of the Yb-4f electrons and τ is the imaginary time. As shown in Fig. 3c, the static susceptibility of Yb1 deviates from the Curie-Weiss law at low temperature above 75 GPa due to the Kondo hybridization, but that of Yb2 always follows the Curie-Weiss law and keeps the local moment behavior. This explains the absence of quasiparticle resonance in the Yb2-4f density of states plotted earlier in Fig. 2f. To understand these huge differences between Yb1 and Yb2 ions, we compare in Fig. 3d their local environment, namely their respective Yb-Se bond distances. With increasing pressure from 20 to 100 GPa, the Yb1-Se distance is seen to decrease from 2.67 to 2.37 Å, while the Yb2-Se distance decreases from 3.0 to 2.66 Å. Hence the Yb2-Se distance at 100 GPa is roughly equal to the Yb1-Se distance at 20 GPa, explaining why the Yb2-4f electrons are still localized at such high pressure.
Note that it is not our purpose here to predict either magnetic ordering or spin liquidity based purely on numerical calculations. Our strategy is to connect the calculated electronic structures with experimental observations of these physical properties. Since in experiment there is no magnetic ordering observed under pressure 19 , the Yb2-4f spins in this pressure region should remain in a spin liquid state 34 . We thus conclude that the high-pressure structure of NaYbSe 2 contains alternating Yb1 heavy-fermion layers and Yb2 QSL layers, a very peculiar structure for studying the interplay between the two exotic quantum many-body states. It is wondering if this might be responsible for the broad crossover to the Fermi liquid behavior observed in the resistivity 19 .

Superconductivity
Experimentally, superconductivity in NaYbSe 2 emerges above 103 GPa in the heavy Fermi liquid phase 19 . To explore possible connection between the pairing and electronic band structures, we compare in Fig. 4a-c the Fermi surface topology at 60, 75, and 100 GPa, respectively. At 60 GPa, the Fermi surfaces only consist of small electron and hole pockets centered at M and Γ points. The small Fermi volume reflects a low carrier density of the semimetallic phase for incoherent Kondo scattering. At 75 GPa, the Fermi surfaces become larger with three cylinders around Γ and a small electron pocket around the M point. At 100 GPa, the cylindrical Fermi surfaces continue to grow but the electron pocket at M is gapped out and diminishes due to strong hybridization. We find superconductivity emerges as the Fermi surfaces become well nested and quasi-two-dimensional. This is seen in Fig. 4d-f, where we plot the real part of the Lindhard susceptibility derived from DFT+DMFT Fermi surfaces. Only beyond 100 GPa, we find a well-defined nesting wave vector. This suggests an intimate connection between superconductivity and the Fermi surface topology and hints that the Fermi surface nesting and its resulting spin density wave fluctuations may play a major role in the Cooper pairing. Moreover, the nesting pattern is in close resemblance to that of the spin excitation spectra of the Heisenberg model on the same lattice as revealed by neutron scattering experiment at ambient pressure 34 . Considering that the Yb2 layers remain in a QSL state, it may be natural to wonder if their associated magnetic fluctuations may also be involved in the electron pairing. At this stage, phonons also cannot be excluded. But in any case, since NaYbSe 2 is in a heavy-fermion state at this pressure, the itinerant heavy Yb1-4f quasiparticles must participate in the superconductivity, making it a candidate heavyfermion superconductor. If this is the case, it would be the 3rd Ybbased heavy-fermion superconductor, in addition to YbRh 2 Si 2 of T c~2 mK 55 and β-YbAlB 4 of T c~8 0 mK 56 . The transition temperature of about 8 K in NaYbSe 2 is strikingly high among all heavyfermion superconductors except Pu-115, probably suggesting a potential cooperation of the Fermi surface nesting, spin fluctuations, and phonons. For such a Fermi surface topology with the antiferromagnetic nesting wave vector, group-theoretical analyses favor a spin-singlet pairing state with gap symmetries listed in Table 1. We have either a nodeless s-wave pairing (A 1g ) or a twocomponent d-wave pairing (E g ), with possible accidental line nodes for appropriate combination of the basis functions, or a less probable g-wave gap (A 2g ) with loop nodes in the k z = 0 plane. More elaborate experiments are needed for better clarification but would be extremely difficult under such high pressures. It will be important if superconductivity can be found at lower pressures in other members of the same family through chemical substitution.

Summary
Our calculations reveal a sequence of electronic phases in pressurized QSL candidate NaYbSe 2 . We find that it is a chargetransfer insulator at ambient pressure and turns into a low carrier density Kondo system with semi-metallic band structures close to 60 GPa and a good heavy-fermion metal above 73 GPa. In between, we predict a delocalization quantum critical point of the Yb1-4f electrons that may explain the observed non-Fermi liquid with linear-in-T resistivity, while the Yb2-4f electrons remain localized over the whole pressure range. The lattice therefore has a special structure constructed with alternating Yb1 heavy-fermion layers and Yb2 QSL layers. Superconductivity emerges in the heavy-fermion metal with well nested quasi-two-dimensional Fermi surfaces, making NaYbSe 2 possibly the 3rd Yb-based heavy-fermion superconductor with a record-high T c . This suggests the possible involvement of both spin fluctuations and phonons as the pairing glues. Our work will be a useful guide for future experimental study of NaYbSe 2 and its other family members.  A 2g k x k z ðk 2 x À 3k 2 y Þ; k y k z ðk 2 y À 3k 2 x Þ line E g ðk 2 x À k 2 y ; k x k y Þ; ðk x k z ; k y k z Þ line * The line * means possible accidental line nodes for certain parameters of the basis functions and otherwise fully gapped.

METHODS DFT calculations
DFT calculations were performed using the full-potential augmented plane-wave plus local orbital method as implemented in the WIEN2k package 40 and the Perdew-Burke-Ernzerhof exchange-correlation potential 42 . The lattice parameters were optimized using DFT+U with an effective Coulomb interaction of 6 eV 41 . The plane-wave cutoff K max was chosen to make R MT K max ¼ 8:0. All calculations were performed on a grid of 1000 k-points in the Brillouin zone. The spin-orbit coupling was included explicitly.

DFT+DMFT calculations
DFT+DMFT calculations were employed to treat the electronic correlations of Yb-4f orbitals. The hybridization expansion continuous-time quantum Monte Carlo (CT-HYB) was used as the impurity solver 43,44 . Yb-4f local orbitals were constructed using projectors with an energy window from −10 to 10 eV relative to the Fermi energy. The Coulomb interaction was chosen to be U = 6 eV and J H = 0.7 eV. A nominal double-counting scheme with n f = 13.0 was used. The self-energy in real frequency was obtained by analytic continuation based on the maximum entropy 38,57 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
We have used the WIEN2k 40 and DFT+eDMFT codes 38 to generate the data. WIEN2k is a commercially available software package and DFT+eDMFT can be obtained from http://hauleweb.rutgers.edu.