Superconductor-Insulator Transition and Fermi-Bose Crossovers

The direct transition from an insulator to a superconductor (SC) in Fermi systems is a problem of long-standing interest, which necessarily goes beyond the standard BCS paradigm of superconductivity as a Fermi surface instability. We introduce here a simple, translationally-invariant lattice fermion model that undergoes a SC-insulator transition (SIT) and elucidate its properties using analytical methods and quantum Monte Carlo simulations. We show that there is a fermionic band insulator to bosonic insulator crossover in the insulating phase and a BCS-to-BEC crossover in the SC. The SIT is always found to be from a bosonic insulator to a BEC-like SC, with an energy gap for fermions that remains finite across the SIT. The energy scales that go critical at the SIT are the gap to pair excitations in the insulator and the superfluid stiffness in the SC. In addition to giving insights into important questions about the SIT in solid state systems, our model should be experimentally realizable using ultracold fermions in optical lattices.

Understanding superconductor-insulator transitions has long been an important challenge in condensed matter physics. The Bose Hubbard model and the Josephson junction array have led to key insights into the SIT in bosonic systems [1][2][3][4]. In contrast, there has been less progress in understanding the SIT in Fermi systems, despite the existence of many electronic systems that exhibit a direct transition from an insulator to a SC. These include, e.g., the disorder-driven SIT in thin films [5][6][7], superconductivity in doped band insulators like SrTiO 3 [8] and at oxide interfaces [9], and the SIT induced in Mott insulators by doping (high-T c cuprates) [10,11] and pressure (organics) [12].
Here we introduce and analyze a two-dimensional (2D) translationally-invariant lattice fermion model with local interactions that exhibits a direct quantum phase transition from an insulator to a SC. Our goal is to gain insights into a number of key issues in the field of SIT through a simple (disorder-free) model that can be analyzed in great detail, rather than to describe a specific experimental condensed matter system. We also note that our model can be realized experimentally using ultracold Fermi atoms in optical lattices.
In the standard Bardeen-Cooper-Schrieffer (BCS) paradigm, superconductivity is a Fermi surface instability of a normal metal. The challenge here is to understand how superconductivity arises out of an insulator that has no Fermi surface. The key insight from our work is that there are Fermi-Bose crossovers in both the insulating and in the superconducting phases, and the SIT is always between a bosonic insulator and a SC in a Bose-Einstein condensation (BEC) regime; see Fig. 1. In the weak coupling limit, though, the bosonic regimes on either side of the SIT may be narrow.
Our main results are: 1) The pairing susceptibility χ in the insulator diverges and the gap to pair excitations ω pair in the insulator van-ishes upon approaching the SIT.
2) The single-particle energy gap E g for fermions remains finite in both phases across the SIT.
3) The SC state is characterized by a pairing amplitude ∆ and a superfluid stiffness D s , both of which vanish approaching the SIT. 4) The insulating state near the SIT is bosonic in the sense that the gap to pair excitations is much smaller than the fermionic gap. 5) The SC in the vicinity of the SIT is a BEC-like regime with several unusual properties. Its fermionic energy gap E g = [(E 0 g ) 2 + ∆ 2 ] 1/2 depends upon both pairing ∆ and the insulating gap E 0 g , and its superfluid stiffness D s is much smaller than the energy gap. Its Bogoliubov excitation spectrum exhibits a minimum gap at k = 0 rather than on an underlying Fermi surface, giving rise to a gap edge singularity that is qualitatively different from the BCS result.
We work with a half-filled attractive Hubbard model on a triangular lattice bilayer; the reasons for this particular choice of lattice are explained in detail below. Our results are based on a variety of analytical approaches, including a strong coupling analysis about the atomic limit, a weakcoupling analysis of the pairing instability in an insulator and mean field theory (MFT). We also present numerical results from determinant quantum Monte Carlo (DQMC) simulations that are free of the fermion sign-problem.
Before describing our work in detail, we comment on its relationship with the classic paper by Nozieres and Pistolesi on "pairing across a semiconducting gap" [13]. They used MFT and estimates of phase fluctuations to analyze superconductivity in a system with a band gap that separates two bands, each with a constant density of states (DOS). Building on their ideas, our work goes beyond their analysis in terms of what we calculate and the methodology used, and this leads to new insights into the problem. Our explicit lattice Hamiltonian permits us arXiv:1507.05641v1 [cond-mat.supr-con] 20 Jul 2015 to use DQMC and is of a form that can be realized in cold atom experiments. We conclude with a comment on the connection between our results and other problems -such as the disorder-tuned SIT, the superfluid-Mott transition for bosons and the BCS-BEC crossover in multi-band systems.
Model: We begin with the constraints on a fermion Hamiltonian that realizes a band insulator-SC transition. First, we need at least two sites (or orbitals) per unit cell to describe a band insulator. Second, we must ensure that the attraction needed for SC does not lead to other broken symmetries. The attractive Hubbard model on a bipartite lattice has an SU(2) symmetry at half-filling, with a degeneracy between SC and charge density waves (CDW) that leads to T c = 0 in 2D. To avoid this, we choose a non-bipartite lattice. Finally, we want to tune the SIT at a fixed commensurate filling. Away from this filling the band insulator becomes a metal, and we do not get an insulator to SC transition. A simple model that meets these criteria is the attractive Hubbard model on two coupled triangular lattices (inset of Fig. 2) with the Hamiltonian The spin σ =↑, ↓ fermion operators at site i are c † iσ and c iσ , with hopping t between in-plane neighbors ij and t ⊥ between interlayer neighbors ij ⊥ . The chemical potential is µ, the local attraction is |U |, and n iσ = c † iσ c iσ . Recently the SIT has been studied [14] in an attractive Hubbard model on a square lattice with near-and next-near-neighbor hopping and a staggered ("ionic") potential to double the unit cell. This model differs from ours in that it has one additional parameter and exhibits CDW order in a limiting case. More importantly, the questions we address and the methodology we use are quite different.
Non-interacting and Atomic Limits: We begin with exactly solvable limits in the phase diagram in Fig. 1. First, consider the noninteracting (U = 0) system with dispersion Here (k x , k y ) lie in the triangular lattice Brillouin zone and k z takes on values 0 (and π) for the bonding (and antibonding) band. For fixed n = 1, this implies a transition from a band insulator (gap E 0 g = 2t ⊥ − 9t) to a metal at t = 2t ⊥ /9; see Fig. 1(a). At U = 0 the (t, µ) phase diagram looks qualitatively similar to Fig. 1(b), except that the insulating "lobes" are triangular in shape with µ c /t ⊥ = ±1 and the SC is replaced by a metal. (See Fig. 6a in Appendix A).
Next consider the "atomic" limit t = 0, for which the lattice breaks up into disconnected (vertical) rungs. We solve in Appendix B the two-site Hubbard model on a rung for arbitrary t ⊥ and |U |. For n = 1 there is a crossover from a fermionic insulator to a bosonic insulator at |U |/t ⊥ 2. For |U |/t ⊥ < 2, the lowest energy excitation is a single fermion (particle or hole), and hence similar to the fermionic band insulator, which is ground state for U = 0. On the other hand, for |U |/t ⊥ > 2 the lowest energy excitation is a pair of fermions. In the large U limit, the ground state is a Mott insulator of bosons, with one boson per rung.
At t = 0 we also find that the extent of the n = 1 insulating phase, (−µ c , µ c ) in Fig. 1(b), is reduced with increasing |U |/t ⊥ . Thus the n = 0 and 2 lobes grow in size relative to n = 1 as |U |/t ⊥ increases. We note that the description of the atomic limit phases as 'insulators' is justified given that the gap is robust to turning on a small t = 0 hopping.
Pairing Instability in the Insulator: We next describe a weak-coupling theory for the dominant instability in the band insulator as we turn on an attraction |U |. The q = 0 dynamical pairing susceptibility in the ladder approximation is given by Here N is the number of lattice sites, f k is the Fermi function, and ε k = ε 0 k − µ H includes the Hartree shift µ H = |U |(n − 1)/2. We analyze the problem for n = 1, choosing µ in the insulator so that we take a trajectory in Fig. 1(a) that goes through the tip of the lobe. (This choice of chemical potential is described in Appendix D).
The divergence of χ ≡ χ(ω = 0) at the SIT is shown in Fig. 2. We tune through the SIT by varying t/t ⊥ , which controls the gap in the band structure, keeping |U |/t ⊥ fixed. It is energetically favorable to create pairs of particles and of holes when the gain in pair binding energy exceeds the band gap. This triggers the SIT, a particle-particle channel analog of exciton condensation in semiconductors [15,16].
The dynamical pair susceptibility χ(ω) exhibits a pole at the two-particle gap ω pair . We see in Fig. 2 how ω pair , the energy to insert a pair into the insulator, goes soft and vanishes at the SIT. In the SC, where pairs can be inserted into the condensate at no cost, ω pair ≡ 0. We show below that the single-particle gap E g remains finite across the SIT.
Mean Field Theory of SC state: The pairing amplitude ∆ = |U | i c i↑ c i↓ /N characterizes the SC state. We find ∆ and µ from the mean field equations from which we calculate the energy gap E g = min E k in the single-particle DOS N (ω).
The evolution of the T = 0 energy gap E g across the SIT is shown in Fig. 2. In the insulator we call the gap E g = E 0 g , but once SC sets in, we find that the SC and insulating gaps add in quadrature E g = [(E 0 g ) 2 + ∆ 2 ] 1/2 . For large t/t ⊥ , the two bands merge, the insulating gap E 0 g collapses and E g = ∆, as in BCS theory. A single-particle gap E g that remains finite and a pairgap ω pair that vanishes at the SIT implies that the insulating state close to the SIT is a boson insulator. As discussed in the Conclusion, these results are very similar to those in simple models of the disorder-tuned SIT [17][18][19].
Given that the gap E g remains finite, what is the critical energy scale as the SIT is approached from the SC? We show in Fig. 2 that the superfluid stiffness D s goes soft at the SIT. D s is obtained from the q → 0 limit of the transverse current-current correlation function [20].
BEC-BCS crossover: We next argue that the SC state near the SIT is more akin to the BEC regime than to the BCS regime [21]. This is best seen from the excitation spectrum in Fig. 3. For large t/t ⊥ , we are in a BCS regime, with the usual Bogoliubov dispersion ±E k with weights u 2 k and v 2 k . The minimum gap ∆ occurs at a finite wavevector (k F in weak coupling), which identifies an "underlying" Fermi surface (FS) [22]. This gap minimum located on a 1D FS contour in 2D k-space leads to the well known (E − E g ) −1/2 singularity in the DOS (in addition to van Hove singularities in the band structure).
The BEC regime near the SIT differs from BCS in a variety of ways. The energy gap E g = [(E 0 g ) 2 +∆ 2 ] 1/2 is located at k = 0. The fact that min E k occurs at a point (not a contour) leads to a jump discontinuity in 2D at the gap-edge (not a square-root singularity). This qualitative difference in the DOS singularity in the BEC and BCS regimes seems not to have been recognized earlier.
A comparison of the superfluid stiffness D s and the pairing ∆ is also illuminating; see Fig. 2. For large t/t ⊥ we find D s ∆ as in BCS theory. Close to the SIT, however, we find D s ∆, a BEC-like regime with wellformed pairs but a small phase stiffness.
Limitations of MFT: The MFT results described above give many important insights into the SIT as a function of t/t ⊥ for fixed |U |/t ⊥ 2, but the approximations involved give misleading results for large |U |. MFT overemphasizes order and predicts a SC state at all t/t ⊥ for |U |/t ⊥ > 2. (See Appendix D for details.) It thus fails to describe the Fermi to Bose insulator crossover, shown in Fig. 1(a), as well as the SIT that occurs at t/t ⊥ of order unity for |U |/t ⊥ → ∞.
DQMC results: To investigate the role of quantum  We cannot obtain ∆ directly from DQMC, so we compute P s = 1/N i,j c i↑ c i↓ c † j↓ c † j↑ , the q = 0 pairing structure factor, which equals |∆| 2 in the infinite-size limit. We also compute the superfluid stiffness D s using the standard Kubo formula [20]. We see the SIT in Fig. 4 from the onset of both |∆| 2 and D s as a function of t/t ⊥ at fixed |U |/t ⊥ = 4. We found similar results at |U |/t ⊥ = 3 (not shown).
We see that D s increases monotonically with t/t ⊥ but P s exhibits non-monotonic behavior, similar to the MFT results for ∆, which we can understand as follows. In the BEC regime, close to the SIT, the order parameter ∆ increases with t/t ⊥ , however, eventually the increase in bandwidth leads to a smaller normal-state DOS, and the BCS ∆ decreases.
Finally, we show in Fig. 5 that the persistence of a finite single-particle gap E g from the insulator to the SC is not an artifact of MFT, and is clearly seen in the DQMC  results for the DOS. The DOS N (ω) was obtained from analytic continuation of DQMC data using the maximum entropy method [27]. A detailed analysis of our DQMC results will be presented elsewhere. Conclusions: Our main results are summarized in the Introduction. We conclude by discussing the parallels between our results and several problems of current interest. First, we comment on the very different problem of the disorder-driven SIT in the attractive Hubbard model, analyzed using a spatially inhomogeneous Bogoliubovde Gennes MFT [17,18] and DQMC [19,28,29]. The predicted [17,18] persistence of the single-particle gap across the SIT has been seen in STM experiments [30][31][32]. The mechanism leading to the insulating gap crucially involves inhomogeneity and an unusual insulator with localized pairs, and is quite different from the gap that persists across the SIT analyzed here.
The collapse of the two-particle gap ω pair on the insulating side of the SIT and that of the stiffness D s on the SC side also have interesting parallels between the disorder-driven SIT [19] and the SIT studied here. However, the disorder-driven SIT is hard to analyze analyti-cally in the manner presented here given the complexity of treating both disorder and interactions.
Next, we note some similarities of our results which are on a fermionic system with the superfluid-Mott transition in the repulsive Bose-Hubbard model (BHM) [1]. Of course, the critical behavior at the SIT in our model is expected to be in the same universality class as that in the BHM. What is much more interesting, however, is that the phase diagrams of the two models in the (µ, t)plane are so similar -see Fig. 1(b) -even though, for small |U |, our fermionic system cannot be mapped onto bosons.
Finally, we comment that the BCS-BEC crossover problem [21], investigated in great detail for ultracold Fermi gases, has almost exclusively focussed on the case of a single band. The BCS-BEC crossover described here has interesting new features that relate to the two bands and the underlying insulating band gap. The study of the BCS-BEC crossover in multi-band systems is in its infancy and there are indications that these ideas may be relevant to SC in materials such as FeSeTe [33].

Appendix A: Non-interacting limit
First consider the non-interacting tight-binding model on a triangular bilayer (Eq. (1) with |U | = 0), whose dispersion relation is given by Eq. (2). The density of states (DOS) is where g tri (ε) = − 1 π Im G tri (ε+i0 + ) is the density of states of a single triangular lattice, and G tri is the triangular lattice Green function [34], G tri (ε) = π 0 dp dq π 2 1 ε − 2(cos p + cos q + cos p cos q) where K is the complete elliptic integral of the first kind as implemented by the Mathematica function EllipticK.
For small t, this model possesses a valence band with energies −t ⊥ − 6t ≤ E ≤ −t ⊥ + 3t and a conduction band with energies +t ⊥ − 6t ≤ E ≤ +t ⊥ + 3t. When t > 2t ⊥ /9 the bands overlap to form a single band. The model behaves as a band insulator or metal depending on where the Fermi energy µ lies within the bandstructure. This leads to the phase diagram in Fig. 6a.

Appendix B: Atomic Limit
Now consider the triangular bilayer Hubbard attractive model in the "atomic limit" (t = 0, |U | > 0). The system now consists of independent two-site Hubbard models on disconnected rungs of the bilayer, each rung described by where i = A, B distinguishes the two sites. Each two-site system can be occupied by N = 0, 1, 2, 3, or 4 fermions. Exact diagonalization in the basis of 16 Fock states shows that the lowest-energy state in each subspace of fermion number N is Fig. 7 shows the energy levels as a function of U (for µ = 0). The ground state is always at half-filling (N = 2, black curve). For |U | /t ⊥ < 2, the lowest energy excitation out of the N = 2 subspace is a fermionic excitation, i.e., adding or removing one fermion to get to N = 1 or N = 3 (red line). In this regime the model might be said to be a "Fermi insulator," with a single-particle gap For |U | /t ⊥ > 2, the lowest energy excitation out of the N = 2 subspace is a bosonic excitation, i.e., adding or removing a pair of fermions to get to the N = 0 or N = 4 subspace (blue line). In this regime the model is a "Bose insulator," with a two-particle gap Appendix C: Atomic Limit plus Hopping We now consider the effects of a small in-plane hopping t t ⊥ on the atomic limit results.
Corrections to ground state: For t = 0, the halffilled ground state of the system can be schematically written as |. . . 2222 . . . . Here each "2" represents the fact that there are two fermions on every rung in the bilayer lattice when the density per site is n = 1. In fact, for each rung, "2" is a specific linear combination of six states |↑↓, 0 , |0, ↑↓ , |σ, σ with two fermions on one rung, with definite amplitudes that depend on |U |/t ⊥ . In the presence of finite in-plane hopping t > 0, the ground state develops an admixture of states such as |. . . 2312 . . . and |. . . 2402 . . . . We will ignore these corrections.
Corrections to excited state with one extra particle: If t = 0, the lowest single-particle excited state is of the form |. . . 2322 . . . , and the single-particle gap is If t > 0, the lowest single-particle excited state is mainly a superposition of states such as |. . . 2322 . . . and |. . . 2232 . . . , connected by hopping t. Thus we might expect that the new single-particle gap is where α is a constant of order unity. In the strong-coupling limit |U | t ⊥ , E g ≈ |U | /2, so we expect the single-particle gap to collapse to zero at about t ∼ |U |. However, this is not an insulator-metal transition (as we see in the next section). Rather, the collapse of the "nominal" single-particle gap represents a crossover from a BEC superfluid to a BCS superfluid regime. This behavior, t BEC-BCS ∝ |U |, is shown as a dashed line in Fig. 1(a) of the text.
In the weak-coupling limit |U | t ⊥ , the above approximation suggests a transition at t ∼ t ⊥ . This is a crude estimate; we actually know that the transition actually occurs at t ≈ 2t ⊥ /9 as |U | vanishes, as shown earlier.
Corrections to excited state with two extra particles: If t = 0, the lowest two-particle excited state is of the form |. . . 2422 . . . , and the two-particle gap is If t > 0, a local two-particle excitation can hop from one rung to an adjacent rung in a two-step process. The effective boson hopping is of the order of t boson = 4t 2 / |U |. Hence the two-particle gap is reduced to where α is a constant of order unity. In the limit |U | t ⊥ , we have ω pair (0) ≈ 4t ⊥ 2 / |U |. This suggests that the two-particle gap ω pair falls to zero when 4t ⊥ 2 / |U | − α × 4t 2 / |U | = 0, i.e., when t ∼ t ⊥ . Then the system becomes overrun by two-particle excitations and makes a transition from BI to BEC. This estimate is crude, as it neglects a large number of corrections to the initial and final state. Nevertheless, we expect that the general idea is correct, i.e., for |U | t ⊥ the BI-BEC phase boundary in the (t, U ) plane has a vertical asymptote at some finite value of t/t ⊥ . This behavior is illustrated in Fig. 1(a) of the text.

Appendix D: Pairing Susceptibility and Mean-Field Theory
We next present a mean-field approach, which works best in the |U | t regime, and thus complements the atomic limit analysis presented in the previous two appendices.
Pairing instability: As discussed in the text, we first calculate the bare particle-particle channel susceptibility, with center-of-mass momentum q = 0, given by Here ε k = ε 0 k − µ H , takes into account the Hartree shift µ H = |U |(n − 1)/2 with the bare dispersion ε 0 k given by Eq. (2). n is the density, N the number of lattice sites and f k the Fermi function. The k-sum is over the two bands k z = 0, π and over the 2D Brillouin zone (k x , k y ).
In the metallic phase χ 0 has the well-known ln ω divergence of BCS theory. However, we concentrate here on the pairing instability in the insulating phase, where χ 0 is finite. The T = 0 result for the n = 1 insulator, which corresponds to a choice of chemical potential −t ⊥ + 3t < µ < t ⊥ − 6t, can be shown to be We are able to write the susceptibility χ 0 in terms of a single-particle Green's function because we are looking only at pairs with total momentum q = 0, built up from two single-particle excitations with momenta ±k, and ε 0 k = ε 0 −k . We also find that similar expressions hold for the n = 0 and n = 2 insulating states.
We then calculate the pairing susceptibility in the ladder approximation, We find the critical in-plane hopping t c in MFT by locating the value of t at which χ 0 (ω = 0) = 1/|U |, so that χ(ω) diverges. Figure 6c shows the MFT phase diagram in the (t, µ) plane, for a fixed coupling |U | /t ⊥ = 1.8. In the insulating phase there is a finite pairing susceptibility χ 0 , which can lead to a pairing instability in the presence of attraction |U |. As |U | increases, more and more of the insulating lobe is converted into a BEC superfluid. Figure 8 shows the MFT phase diagram in the (t, U ) plane at half-filling. As |U | increases, the system undergoes a transition from an insulator (INS) to a superfluid (BEC), due to the pairing instability. At the end of this Appendix we will discuss the limitations of MFT and why the MFT results in Fig. 8 look so different from the more correct results shown in Fig. 1(a).
Ordered state: In the superconducting state, there is a finite pairing amplitude ∆ = (|U |/N ) i c i↑ c i↓ . We solve the mean-field theory (MFT) equations to find ∆ and µ, where the Hartree-shifted dispersion ε k = ε 0 k − µ H determines the Bogoliubov quasiparticle spectrum E k = √ ε k 2 + ∆ 2 . Using this we can calculate the energy gap E g = min E k in the DOS N (ω) for single-particle excitations. (The k-integrals over the 2D Brillouin of the triangular lattice are performed by reducing them to 1D integrals using the exact density of states (DOS) in terms of the triangular lattice Green function.) The BEC-BCS boundary (dotted black curves in Figs. 6b and 6c) is determined by the criterion that the dispersion minimum occurs at k = 0 in the BEC regime and at k = 0 in the BCS regime.
Choice of chemical potential: We are interested in studying the superconductor-insulator transition tuned by hopping t, rather than by fermion density n. Thus, in the superconducting phase we choose the chemical potential µ (according to Eq. (A17)) such that the average density corresponds to half-filling, n = 1. In the insulating phase, we choose µ such that the two-particle excitation gap is particle-hole-symmetric. In other words, the energy cost of adding a pair of fermions, ω e2 , is equal to the energy cost of removing a pair of fermions, ω h2 . This choice of µ corresponds to the dashed blue curve in Figs. 6b and 6c, which bisects the n = 1 insulator lobe, and passes through the tip of the lobe. The quantities in Fig. 2 of the text are plotted for this choice of µ.
Note that the single-particle excitation gap, E g , is not particle-hole-symmetric along this trajectory. The value of E g plotted in Fig. 2 of the text corresponds to the smaller of the particle gap and the hole gap, min(E e , E h ).
If one attempts to choose µ to make E e = E h , one finds that ω e2 = ω h2 . This µ(t) trajectory exits through the side of the lobe rather than the tip of the lobe, representing a density-tuned transition rather than a hoppingtuned one.
Strengths and Limitations of MFT: MFT allows us to track the behavior of many quantities (see Fig. 2 in the text) that are intuitively meaningful, but not accessible in DQMC or other methods. For fixed U and varying t, MFT gives a good idea of the general behavior of these quantities.
In MFT the interaction is decoupled exclusively in terms of the order parameter. Hence MFT fails to capture important correlations in the insulating state. In the atomic limit (t = 0), as |U | increases, the ground state wavefunction adjusts itself to take advantage of the attraction to produce a binding energy. There is no superfluid state at t = 0 (as would have been predicted by MFT). Compare Fig. 8 and Fig. 1(a). MFT gives bad results as a function of U .
In MFT, the insulator-superfluid boundary |U | (t) is a decreasing function, whereas atomic limit and DQMC results suggest that it is an increasing function as shown in Fig. 1(a). Thus we can no longer think of the insulator-SC transition as due to a pairing instability triggered by increasing |U |. Nevertheless, we can still understand the insulator-SC transition at fixed U , where a change in bandstructure (e.g., reduction in band gap) that causes χ 0 to increase beyond 1/ |U | and precipitates a pairing instability.