Two-stage superconductivity in the Hatsugai-Kohmoto-BCS model

Superconductivity in strongly correlated electrons can emerge out from a normal state that is beyond the Landau's Fermi liquid paradigm, often dubbed as"non-Fermi liquid". While the theory for non-Fermi liquid is still not yet conclusive, a recent study on the exactly-solvable Hatsugai-Kohmoto (HK) model has suggested a non-Fermi liquid ground state whose Green's function resembles the Yang-Rice-Zhang ansatz for cuprates [P. W. Phillips, L. Yeo and E. W. Huang, Nat. Phys. $\bf{16}$, 1175 (2020)]. Similar to the effect of on-site Coulomb repulsion in the Hubbard model, the repulsive interaction in the HK model divides the momentum space into three parts: empty, single-occupied and double-occupied regions, that are separated from each other by two distinct Fermi surfaces. In the presence of an additional Bardeen-Cooper-Schrieffer (BCS)-type pairing interaction of a moderate strength, we show that the system exhibits a"two-stage superconductivity"feature as temperature decreases: a first-order superconducting transition occurs at a temperature $T_{\rm c}$ that is followed by a sudden increase of the superconducting order parameter at a lower temperature $T_{\rm c}^{\prime}<T_{\rm c}$. At the first stage, $T_{\rm c}^{\prime}<T<T_{\rm c}$, the pairing function arises and the entropy is released only in the vicinity of the two Fermi surfaces; while at the second stage, $T<T_{\rm c}^{\prime}$, the pairing function becomes significant and the entropy is further released in deep (single-occupied) region in the Fermi sea. The phase transitions are analyzed within the Ginzburg-Landau theory. Our work sheds new light on unconventional superconductivity in strongly correlated electrons.

Introduction. The pairing mechanism of unconventional superconductivity remains one of the central issues in condensed matter physics. Conventional superconductivity has been well captured by the classic Bardeen-Cooper-Schrieffer (BCS) theory [1], in which a second-order superconducting phase transition occurs as a result of the Cooper pairing instability of the Fermi liquid normal state [2]. However, such a Fermi liquid normal state is absent in many, if not most, unconventional superconductors. Instead, the corresponding normal state is often referred as a "non-Fermi liquid" (NFL) or "unconventional metal" state [3][4][5]. In contrast to Fermi liquids that can be adiabatically connected to a gas of non-interacting fermions and be well depicted by interactions between quasi-particles [6,7], a generic paradigm for NFLs has not yet been established so far [8,9]. However, some experimental criteria for NFLs are commonly accepted. For instance, electric resistivity deviates from the ρ(T ) ∝ T 2 temperature dependence, and specific heat C V (T ) is no longer linearly temperature-dependent [4,[10][11][12]. Moreover, a variety of realistic materials exhibit NFL behaviors, which include but are not limited to cuprates [11,12], ironpnictides and chalcogenides [13,14], and heavy-fermion compounds [10,15]. The superconducting phase emerges from a NFL normal state in these materials [10,13,16]. It is illuminating to understand their pairing mechanisms from studying the pairing from a unconventional metal which beyond Landau's Fermi liquid theory.
On the theoretical side, despite the lack of a well recognized paradigm for NFLs [8,[17][18][19][20], the mechanisms and their superconducting instabilities have been extensively explored in quantum critical models from several different approaches in recent years, such as: coupling of the Fermi sea and the bosonic fluctuations [21][22][23][24], and the system of of fermions with strong random interactions [25][26][27][28][29][30][31], and the phenomenological fermion propagators with anomalous retardations [32][33][34][35][36], etc.. Among them, several exactly-solvable models are of particular interest that include the Hatsugai-Kohmoto (HK) model [37,38]. The interacting part in this model can be viewed as a momentum-space counterpart to the on-site Hubbard interaction, while the non-interacting part is the same as the Hubbard model. The HK model can host a NFL state with non-Landau's quasi-particle excitations [39,40], such that it violates the Luttinger's theorem and gives rise to a Green's function that resembles the Yang-Rice-Zhang (YRZ) ansatz for cuprates [41,42]. Indeed, the zeros of the YRZ-like Green's function G(k, ω = 0) enclose a Luttinger surface instead of a usual Fermi surface [43][44][45][46], indicating the Mottness in the strong-coupling limit and an unconventional metal or NFL in the region of weak or intermediate-coupling [39,47]. The possible Cooper pairing instability and associated dynamic spectral weight transfer were also investigated [39,48,49]. More interestingly, it was demonstrated that Fermi arcs and a pseudo gap will show up in such an unconventional metal, when the "on-site" interaction becomes k-dependent and changes sign in momentumspace [50]. Very recently, taking account of additional BCS pairing terms, J. Zhao et al. studied the thermodynamics of the HK-BCS model in the strong pairing limit, and revealed a first-order superconducting transition instead of the continuous phase transition in the BCS theory [51]. (1) At T > Tc, ∆ = 0, the normal state is an unconventional metal (NFL), on which the momentum space is divided into three regions: empty (Ω0), single-occupied (Ω1) and double-occupied region (Ω2). These regions are separated by two Fermi surfaces (brown dash lines).
(2) As temperature decreases, a first-order superconducting transition occurs at T = Tc. For T c < T < Tc, the superconductivity comes from the electron pairing in the vicinity of the two Fermi surfaces and leads to the superconducting phase in regime 1 (SC1). (3) At a lower temperature T c (T c < Tc), a second jump of ∆(T ) takes place, resulting in the superconducting phase in regime 2 (SC2), where Cooper pairs come into being inside the singly-occupied region (Ω1) and play a significant role.
ing limit studied in Ref. [51]. We calculate the binding energy of a Cooper pair, and study the phase diagram. Unexpectedly, we find that the system undergoes a "two-stage" process as temperature decreases. As illustrated in Fig. 1, in addition to a first-order superconducting transition at T c , the superconducting order parameter ∆(T ) has a jump to a larger value at a lower temperature T c (< T c ), accompanying with a sudden drop in entropy. The underlying physics is interpreted in accordance with the pairing function and the entropy release in momentum space, and the nature of discontinuity in the SC order parameter and entropy as a function of temperature is analyzed in the Ginzburg-Landau theory. HK model revisit. The Hatsugai-Kohmoto model [37] describes strongly correlated electrons with momentum-space on-site interaction. The Hamiltonian takes a form of where c † k,σ (c k,σ ) creates (annihilates) a fermion at momentum k with spin σ =↑, ↓, and n kσ = c † k,σ c k,σ relates to its density distribution. k = −2t(cos k x + cos k y ) is the singleparticle energy dispersion and µ is the chemical potential, in which t is the hopping integral. Without loss of generality, we set t = 1 as the energy unit hereafter. U > 0 represents an onsite repulsion in the momentum space. The locality of the interaction allows us to factorize the huge Hilbert space into the direct product of the k-subspace that is spanned by the basis Ground states of the HK model can be obtained from the fermion occupation in the momentum space, as illustrated in Fig. 2 (a). In the presence of a positive U , the momentum space will be divided into three regions in a ground state: empty region (Ω 0 ), single-occupied region (Ω 1 ), and double-occupied region (Ω 2 ). This gives rise to two distinct Fermi surfaces [52] and two corresponding Fermi levels at µ and µ − U respectively. Here the chemical potential µ is determined by the filling number n using the relation where Θ is the Heaviside function and V 0 is the volume of the Brillouin zone. The two Fermi levels can be viewed from the distribution function n(ω, T ) ≡ n kσ ( k = ω, T ) as well [53], where two sudden jumps occur at µ and µ − U , as shown in Fig. 2 (b). When U = 0, the region Ω 1 vanishes and the two Fermi surfaces merged into a single one as in the free-fermion model. Note that the single-occupied region Ω 1 always exists as long as U > 0, while the double-occupied region Ω 2 may vanish if a filling number is chosen such that µ − U exceeds the bottom of the energy band.
The retarded Green's function for this exactly solvable model reads, where ξ k = k − µ andσ is the opposite spin index to σ. Note that G σ (k, ω) does not depend on σ, and can be abbreviated as G (k, ω). The electron spectral function A(k, ω) = −2ImG (k, ω) has been found at zero temperature and plotted in Fig. 2 (c). It displays two "truncated" bands separated by U , which originate from a double-occupied to single-occupied excitation and a single-occupied to empty excitation, respectively. Residual entropy. It is worth noting that the positive U imposes the single occupancy constraint at each k-point in the single-occupied region (Ω 1 ) that gives rise to a huge ground state degeneracy and a finite entropy density at zero temperature, which is proportional to the volume of Ω 1 . This violates the third law of thermodynamics, and resembles the residual entropy in classical spin liquids on geometrically frustrated lattices [54]. As will be discussed later, an extra pairing interaction will lift the huge ground state degeneracy and release the entropy, resulting in a two-stage superconductivity.
Cooper pair problem. As investigated in Ref. [39], an infinitesimal pairing interaction will cause superconducting pairing instability in the HK model. The bound-state energy E C for the formation of a Cooper pair on top of the Fermi sea has been estimated, where the spin polarization in the singleoccupied region was assumed [39]. However, there is a huge spin degeneracy in the Ω 1 region, and spin polarization configuration is not favorable for the Cooper pairing. Here we revisit the Cooper pair problem without assuming the spin polarization in the Ω 1 region [53]. Consider a generic situation when both Fermi levels locate within the bandwidth W = 8t, thereby −W/2 < µ < W/2 and U < W , we find that the bound-state energy E C can be determined as follows [53], In the limit of U → 0, it yields E C ≈ −(W − 2µ)e − 2W V , which restores the BCS solution. In the presence of a weak or intermediate pairing interaction V and a relative large U , namely, when V W and |E C | < U W , we find an asymptotic solution to Eq. (3), which deviates from the BCS solution apparently. For a fixed electron filling number n, numerical solutions to Eq. (3) can be found self-consistently. As plotted in Fig. 3, the binding energy |E C | = −E C increases as U and/or V increases, suggesting the enhancement of Cooper instability by the repulsive U . We should note that, our results are different from Ref. [39], in which the single-occupied region plays no role to the bind-energy, and the pairing instability is underestimated since the instability of the Fermi surface on µ−U is neglected. HK-BCS model. For further studying the superconductivity in the HK model, we introduce a BCS pairing interaction in the mean-field level that gives rise to the HK-BCS model as follows, where ∆ is the superconducting pairing gap and V > 0 refers to an attractive pairing strength. This mean-field Hamiltonian can be exactly diagonalized at each k point, and similar to Ref. [51], the superconducting order parameter ∆ ≡ −V k c −k↓ c k↑ can be found through searching the global minimum of the free energy, with the help of ∂FS[∆] ∂∆ = 0. Where Z k ≡ n e −βE n,k , in which {|n, k } and {E n,k } are the eigenstates and eigenspectra obtained from diagonalizing H within the tensor-product In this treatment, the HK interaction remains intact and plays a crucial role in the unconventional superconductivity.
Two-stage transitions. With fixed electron filling n at T = 0 in the normal limit, and we neglect the temperature evolution of µ for simplicity, the superconducting gap ∆(T ) can be found out through the minimization of the free energy for given U and V . As observed in Ref. [51], there occurs a first-order superconducting phase transition as long as U > 0, in contrast to a continuous phase transition in the U = 0 BCS limit. Namely, as temperature is lowering, ∆(T ) jumps from zero to a finite value at T c abruptly.
Surprisingly, in addition to the first-order transition at T c , we find that there exists an extra ∆(T ) jump at a lower temperature T c (< T c ) when U is sufficiently large in comparison with V . As demonstrated in Fig. 4 (a), for n (T = 0) = 0.4 and V = 1, when 0 < U < U c ∼ 2.7, there is only one firstorder transition at T c ; while when U > U c , there emerges a sudden jump at T c (< T c ). Free energy. These first-order phase transitions can be understood from tracking the temperature evolution of the minimum of the free energy as a function of ∆. We compute the free energy difference between the superconducting state and the normal state, F S − F N , and plot it in Fig. 4 (b) as a function of ∆/∆ 0 , where ∆ 0 = ∆(T = 0). Here F N is the normal state free energy calculated at ∆ = 0.
As shown in Figs. 4(b) and 4(c): (1) When T > T c , the minimum of the free energy locates at ∆ = 0; (2) When T goes across T c , the free energy minimum switches from ∆ = 0 to a finite value ∆ = ∆ min , suggesting a first-order superconducting phase transition at T c [see Fig. 4 (a)]; (3) As temperature is lowering, in the region of T c < T < T c , there develops an extra local minimum at a larger value, ∆ = ∆ min (> ∆ min ), while the global minimum (i.e., the one associated with the lowest free energy) evolves from the one (∆ = ∆ min ) arising at T + c continuously; (4) When T decreases further and goes across T c < T c , the global free energy minimum switches from the smaller one ∆ = ∆ min to the larger one ∆ = ∆ min [see Fig. 4 (c)], and then the global minimum ∆ min approaches to ∆ = ∆ 0 as T → 0.
Phase transitions. Near T c , Zhao et al. [51] has explained the first-order nature of the SC transition via the analysis of the Ginzburg-Landau (GL) approach up to sixth-order terms of the free energy. While, in here, it is remarkable that such the two-minimum feature in free energy around T c requires eighth-order terms in the GL free energy functional, which takes the form of where α, β, γ and η are the expansion coefficients and depend on temperature T , and η > 0, or η = 0 and γ > 0, ensures the stability of the system. To study phase transitions for the HK-BCS model, we consider critical regions: T ≈ T c and T ≈ T c .
(1) T ≈ T c : It turns out that the occurrence of a first-order transition at T = T c impose constraints for expansion coefficients at this critical point as follows [53], And the superconducting gap at T c reads In the limit of η = 0, it becomes which restores the result in Ref. [51].
(2) T ≈ T c : In the presence of the first-order-like jump at T = T c , the sign of expansion coefficients can be determined in the critical region as follows [53], α < 0, β > 0, γ < 0 and η > 0.
The critical condition at T = T c is given by and the temperature regions T > (<)T c are separated from each other in accordance with the inequality as follows, The superconducting order parameters at T ± c read To study the temperature dependence ∆(T ) around T c , we introduce the dimensionless parameter t = (T − T c ) /T c , and find for small t , where b ± > 0 are two positive parameters that can be determined from experimental data or microscopic theory [53]. We would like to remark that the first-order-like jump at T = T c will be rounded and become a crossover when ∆ min = ∆ min , which leads to an extra condition for the crossover, in addition to Eq. (12) Entropy release. As mentioned earlier in this paper, the HK model has a huge residual entropy at zero temperature, which is proportional to the volume of Ω 1 . This residual entropy will be released by the superconducting pairing. We find that major entropy release will take place below T c , as long as there exist two-stage process; while there is a minor entropy release at T c . On the contrary, when there is only one firstorder transition, or the two-stage process merge to a single one, there will be a significant entropy release at T c . Typical examples for entropy release have been demonstrated in Fig. 5. Discussions. (i) On the sudden changes at T c : Though the first order derivative of the free energy ∂F/∂T are discontinuous at T c , we didn't mark it as a real phase transition since there is only one order parameter in here. This first-order-like changes can be understood by using the Ginzburg-Landau theory as well, which can be circumvented around some critical point in the phase diagram, resembling the liquid-gas phase transition [55]. The end of such a first-order-like change in the parameter space is indicated by the dimension reduction of the critical hyper-surface, i.e., the extra constraint in Eq. (16) reduce the dimensionality of the critical hyper-surface [given by Eq. (12)] by one.
(ii) The existence of the two-stage process can be attributed to a sufficiently large U . Note that U is the energy separation between the two Fermi levels, µ and µ − U , i.e., the energy width of the single-occupied region Ω 1 . When U is not sufficiently large, the pairing interaction V will pair up all the k points in the small Ω 1 at T = T − c . Otherwise, at the fist stage, say, T c < T < T c , V will pair up the states in the vicinity of Fermi surfaces only, while leave those deep inside Ω 1 unpaired.
(iii) Microscopically, the two-stage superconductivity can be visualized by the momentum distribution of the Cooper pairing function b k = c −k↓ c k↑ and the entropy s k = − 1 2 n ρ n,k ln ρ n,k in the Brillouin zone, where ρ nk = e −E n,k /T /Z k and the factor 1 2 come from the folding of the Brillouin zone. As demonstrated in Fig. 6 (a) and (b): (1) when T > T c , the pairing function b k = 0, i.e., Cooper pairs are absent, and the entropy is dominated by the singleoccupied region Ω 1 ; (2) when T c < T < T c , Cooper pairs come into being in the vicinity of the two Fermi surfaces, associated with a weak and in-situ entropy release, while the entropy inside Ω 1 region remains to be ln 2 [see Fig. 6

(d)];
(3) when temperature is below T c , the distribution of Cooper pairs starts to extend to the whole Brillouin zone, in particular, single-occupiled Ω 1 region and double occupied Ω 2 region, and the residual entropy in Ω 1 region are released entirely [see Fig. 6 (c) and (d)].
Summary. We have studied the HK-BCS model and revealed that, in addition to a first-order superconducting transition occurs at T c , there allows an extra first-order-like changes at a lower temperature T c < T c , as long as the momentum space on-site repulsion U dominates over the superconducting pairing strength V . This type of two-stage process have been formulated within a Ginzburg-Landau theory consisting of eighth-order terms. The underlying physics for the formation of two-stage superconductivity has been discussed. Our results provide basic qualitative understanding of pairing in NFL systems, such as, lightly doped cuprates, heavy fermions etc. These results are also useful in understanding superconductivity in the material-specific large-scale computational methods for strongly correlated materials.
-SUPPLEMENTAL MATERIALS-Appendix A. The Green's function of the HK model Consider the the single-particle Green's function in the imaginary time space, whereT is the time-ordering operator. By use of the equation of motion, 3) Then, consider the equation of motion on Q σ (k, τ ), it gives Perform the Fourier transformation by using the relation Q σ (k, τ ) = 1 β ωn e −iωnτ Q σ (k, iω n ), where ω n = (2n + 1)πT (n ∈ Z) is the fermionic Matsubara frequency and T is the temperature, one can obtain the analytical expression on Q σ (k, iω n ), and substitute it into the equation of motion of G σ , one can obtain the Green's function The ground state of the HK model can be written as where φ k is an arbitrary phase which has no effect on the physical quantities and we take φ k = 0 for brevity, and |0 is the vacuum state. Due the spin uncertainty in the Ω 1 region, there are large spin degeneracies are persist in the low temperature.
Since H HK (k) are composed of 4 states for each k point: , with their energies are 0, ξ k , ξ k and 2ξ k + U , respectively. One can easily obtain the partition function as where β = 1/k B T , with k B is the Boltzmann constant, which is taken as 1 for brevity. The particle density distribution can be expressed as n kσ (T ) = 1 Z HK Tr n kσ e −βHHK = e −βξ k + e −β(2ξ k +U ) 1 + 2e −βξ k + e −β(2ξ k +U ) .

Appendix B: Revisit the Cooper instability
Before the study of the superconductivity, we investigate the microscopic picture of Cooper instability from the unconventional metal of the HK model.
With adding a BCS pairing interaction The average energy of HK-BCS model can be evaluated as Along with the Cooper's approach [39,56], we construct the Cooper-pair wavefunction, where b † k = c † k↑ c † k↓ , and the normalization condition for |ψ gives the relations for the coefficients α k , β k , The average energy of the |ψ can be evaluated as Then, the energy change with adding two electrons is E C = E − E 0 . Including the normalization condition by introduce a Lagrange multipler λ, we define the function and using the variational conditions ∂Q ∂α * k = 0 and ∂Q ∂β * k = 0, the equations for α k and β k can be obtained as Summing α k and β k over Ω 0 and Ω 1 regions respectively, we can obtain the self-consistent equation For brevity, if we consider ρ (ω) = k δ (ω − ε k ) = 1 W for −W/2 < ω < W/2, and assume U < W and W > 2µ, then, By using the relation k∈Ωi (...) = ρ i (ω)(...)dω, the self consistent equation eq.(B.8) can then be rearranged as Define the dimensionless quantities u = U/W , υ = V /W , ε = E C /W ,μ = µ/W , the above equation reduces as In BCS case, u = 0, the equation gives 1 = υ 2 ln 1−2μ−ε −ε , which yielding the solution: ε = − 1−2μ . One can see that, an infinitesimal attraction can induce a stable Cooper-pair bound state (i.e., ε < 0).
For u = 0, as υ << 1, and |ε| << u, the negative solution for ε is given by the equation −e − 4 υ ≈ ε 3 (1−2μ) 2 u , which yielding to ε = − (1 − 2μ) 2/3 u 1/3 e − 4 3υ . As is shown in here and in main text, the Cooper instability is still exact for finite U case. [ξ k (n k,↑ + n −k,↑ + n k,↓ + n −k,↓ ) + U (n k↑ n k↓ + n −k↑ n −k↓ ) which can be exactly diagonalized in the Fock space spanned by the 4-fermion occupation states {|n k,↑ , n k,↓ , n −k,↑ , n −k,↓ } for each k point. By performing the exact diagonalization, we can obtain the eigen states {|n, k } and the eigenspectra {E n,k }. Then, the free energy can be computed as, where Z k = n e −E n,k /T . And the superconducting gap ∆ ≡ −V k c −k↓ c k↑ can be found by searching the global minimum of the free energy, with the help of its minimization, S1 show the numerical data of the evolution of the changes of the free energy F = F S − F N for U = 4, n(T = 0) = 0.4 and V = 1, where the F N refers to the normal-state free energy with taken ∆ = 0.

Appendix D: Ginzburg-Landau analysis
Define the Ginzburg-Landau free energy functional as, where η > 0, or η = 0 and γ > 0, ensures the stability of the system, and α, β, γ are the expansion coefficients. For brevity, we introduce x ≡ ∆ 2 , and define a function f (x) as Then, the minimization of the free energy is equivalent to find the minimum of f (x) for x ≥ 0. By requiring f (x) = ∂f (x) ∂x = 0, the extreme of f (x) is determined by the cubic equation, To study phase transitions for the HK-BCS model, we consider critical regions: T ≈ T c and T ≈ T c . (1) T ≈ T c : to make the first-order transition happen, the superconducting order parameter will be switched to a finite value ∆ = ∆ min from ∆ = 0. There are at least one maximal x = x 1 and one minimum x = x 2 occur (x 2 > x 1 > 0), as shown in Fig. S2(a), and the critical condition is determined by the solution of the equations f (x 2 ) = f (x 2 ) = 0, i.e., Combining the above two equations to eliminate α or η, we can obtain the two equations These two quadratic equations should share a same positive root x 2 , i.e., which restores the result in Ref. [51].