A real space auxiliary field approach to the BCS-BEC crossover
The BCS to BEC crossover in attractive Fermi systems is a prototype of weak to strong coupling evolution in many body physics. While extensive numerical results are available, and several approximate methods have been developed, most of these schemes are unsuccessful in the presence of spatial inhomogeneity. Such situations call for a real space approach that can handle large spatial scales and retain the crucial thermal fluctuations. With this in mind, we present comprehensive results of a real space auxiliary field approach to the BCS to BEC crossover in the attractive Hubbard model in two dimensions. The scheme reproduces the Hartree-Fock-Bogoliubov ground state, and leads to a scale that agrees with quantum Monte Carlo estimates to within a few percent. We provide results on the , amplitude and phase fluctuations, density of states, and the momentum resolved spectral function over the entire interaction and temperature window. We suggest how the method generalises successfully to the presence of disorder, trapping, and population imbalance.
The BCS to BEC crossover in attractive fermion systems is a topic of enduring interest bcs-bec-rev . With increasing interaction, the ground state of a weak coupling ‘BCS superconductor’ bcs , with pair size much larger than the interparticle separation , evolves smoothly eagles ; leggett ; noz ; micnas ; rand-rev into a ‘Bose-Einstein condensate’ (BEC) of preformed fermion pairs with . , above, is the Fermi wavevector. The ‘high temperature’ normal state changes from a conventional Fermi liquid at weak coupling to a gapped phase at strong coupling. While the zero temperature pairing gap increases with coupling strength, the superconducting in lattice systems reaches a maximum at intermediate coupling and falls thereafter. A striking consequence of the separation of pairing and superconducting scales is the emergence of a (pseudo)gapped normal phase, with preformed fermion pairs but no superconductivity due to strong phase fluctuations.
The early work of Leggett leggett and Nozieres and Schmitt-Rink noz provided the intuitive basis for understanding this problem. It has since been followed up by powerful semi-analytic schemes sa-2psc ; sa-deisz ; sa-kopec ; sa-dupuis ; sa-miyake , extensive quantum Monte Carlo (QMC) work qmc-scal1 ; qmc-moreo1 ; qmc-moreo2 ; qmc-spingap ; qmc-nfl ; qmc-singer ; qmc-trm1 ; qmc-beck ; qmc-paiv1 ; qmc-paiv2 , and most recently dynamical mean field theory (DMFT) dmft-keller ; dmft-castel ; dmft-garg ; dmft-bauer ; dmft-dup ; dmft-koga . The efforts have established the non-monotonic , and the presence of a pseudogap in the single particle spectrum beyond moderate coupling and temperature .
While the success of multiple methods in capturing the crossover is remarkable, most of them depend on translation invariance in the underlying problem. They do not naturally generalise to problems that involve the presence of disorder, or a confining potential, or the emergence of spontaneous modulation. These situations, for example, occur in the context of the disorder driven superconductor-insulator transition sit-rev , trapped fermions in an optical lattice opt-rev , and Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states fflo-rev in population imbalanced systems. Such problems, in general, require a real space approach. This paper presents such a real space implementation based on a static auxiliary field (SAF) scheme.
Our main results are the following: (i) We demonstrate the ability of our approach to quantitatively capture the scale across the BCS-BEC crossover, confirming its usefulness at all interaction strengths. (ii) We quantify the crossover from an amplitude fluctuation dominated regime to a phase fluctuation dominated regime, through the ‘high ’ intermediate coupling window where both are important. (iii) We present the thermal evolution of the single particle density of states with interaction and temperature, and, more importantly, the momentum resolved spectral function . We compare our results to QMC data wherever available.
The paper is organised as follows. In Sec.II we quickly compare the existing analytic and numerical methods used to study the BCS-BEC crossover in lattice models. Sec.III presents our model and describes the method used in detail. Sec.IV shows our results on thermodynamic indicators, the nature of fluctuations, density of states, and . Sec.V discusses the limitations of our method and the scope for further work.
Ii Earlier work
Since the BCS-BEC crossover is a prototype of weak to strong coupling evolution, several methods, of increasing sophistication, have been brought to bear on it. These include mean field theory (MFT) leggett , MFT corrected by gaussian fluctuations noz , the self consistent T-matrix approach (SC-TMA) sa-deisz ; sa-miyake , a two particle self consistent (2PSC) scheme sa-2psc , the mapping to XY models sa-dupuis ; sa-kopec , quantum Monte Carlo (QMC) qmc-scal1 ; qmc-moreo1 ; qmc-moreo2 ; qmc-spingap ; qmc-nfl ; qmc-singer ; qmc-trm1 ; qmc-beck ; qmc-paiv1 ; qmc-paiv2 , and recently dynamical mean field theory (DMFT) dmft-keller ; dmft-castel ; dmft-garg ; dmft-bauer ; dmft-dup ; dmft-koga .
There are detailed descriptions of these methods available in the original literature so we just provide a table that compares the strengths and limitations of these methods in the light of a few crucial indicators. These in our opinion include (i) thermodynamics: , (ii) single particle spectra, (iii) two particle properties, e.g, conductivity, and (iv) handling inhomogeneity, e.g, disorder or trapping.
|Method||Spectra at large||Transport at large||Handling inhomogeneity|
|MFT||Correct only when .||No gap/pseudogap in the normal state.||No access to transport.||Real space MFT is reasonable at , can access size .|
|MFT+fluctuations||No results on the full in the lattice model.||PG above mean field , incapable of capturing BKT physics metzner||No results||Not systematically explored.|
|SC-TMA||Accurate upto intermediate coupling; captures non monotonic behaviour correctly to .sa-deisz .||Shows a PG, but quantitatively inaccurate due to neglect of vertex corrections metzner .||No results||Not generalised.|
|2PSC||Fluctuations drive to zero metzner .||Accurate upto intermediate coupling sa-2psc .||No results||Not generalised.|
|XY models||Captures non monotonic but not quantitatively accurate sa-kopec ; sa-dupuis .||PG inferred from different pairing and scales.||Not explored.||Disordered XY model needs to be derived from the Hubbard model.|
|QMC||Accurate, sets the benchmark qmc-moreo1 ; qmc-paiv2 .||Accurate in principle but involves uncertainties due to analytic continuation from imaginary frequencies qmc-moreo2 ; qmc-singer .||Contains the relevant physics but dynamical properties are difficult to extract due to analytic continuation.||Handles inhomogeneity bouad but sizes limited to . Limited to .|
|DMFT||Captures non monotonic but quantitatively inaccurate when used in the 2D context dmft-keller .||Accurate dmft-dup .||Transport misses bosonic contribution.||Requires ad hoc real space generalisation.|
|SAF||Accurate, matches quantitatively with QMC.||Accurate, compares reasonably with QMC.||Transport misses bosonic contribution.||Handles inhomogeneity, method, readily accesses size trans .|
Iii Model and Method
We study the attractive two dimensional Hubbard model (A2DHM),
For us denotes nearest neighbour tunneling amplitude on a square lattice. is the chemical potential. We set and measure all other energies in terms of it. is the strength of onsite Hubbard attraction. We will focus on the density which is close to half-filling but avoids the density wave features of .
The model is known to have a superconducting ground state for all , while at there is the coexistence of superconducting and density wave (DW) correlations in the ground state. For the ground state evolves from a BCS state at to a BEC of ‘molecular pairs’ at . The pairing amplitude and gap at can be reasonably accessed within mean field theory or a simple variational wavefunction.
Mean field theory, however, assumes that the electrons are subject to a spatially uniform self-consistent pairing amplitude . At small this vanishes when , but at large it vanishes only when . The actual at large is controlled by phase correlation of the local order parameter, rather than finite pairing amplitude, and occurs at , where is a function of the density. The wide temperature window, between the ‘pair formation’ scale and corresponds to equilibrium between unpaired fermions and hardcore bosons (paired fermions).
We use an auxiliary field scheme, based on the Hubbard-Stratonovich (HS) transformation, to retain the thermal amplitude and phase fluctuations explicitly. Formally, the HS transformation (below) allows us to recast the original interacting problem into that of electrons coupled to space-time fluctuating auxiliary fields. Usually the auxiliary fields are either the ‘pairing fields’ or a ‘density’ field . The pairing field based decomposition is employed in the Bogoliubov-de Gennes (BdG) approach, while QMC uses the density based decomposition.
While these two representations lead to the same answer when the resulting problem is treated exactly, any approximation leads to answers that are decomposition dependent. ‘Single channel’ decompositions solms ; dubi , above, have the shortcoming that they do not reproduce the conserving Hartree-Fock-BdG (HFBdG) theory in their static limit, crucial when density inhomogeneities are present.
Following recent suggestions conduit we rewrite the A2DHM in terms of classical auxiliary fields in both pairing and density channels, see below, so that the saddle point reproduces HFBdG theory.
We treat the as time independent, i.e classical fields. The partition function is:
As , it is obvious that will be dominated by the saddle point configuration, i.e, , , that minimise . Thus at : , . These minimum conditions are easily recognised as the self-consistency requirements of HFBdG theory. Minimising is thus the same as solving the Bogoliubov-de Gennes (BdG) equations. While this is good, we want to capture the thermal physics as well.
At finite the BdG approach would promote the consistency condition above to thermal averages, and the electrons would see only the average background at any temperature. It is easy to see the limitation of this approach when we consider the clean problem at large . Since the system is translation invariant, averages like are site independent. Their self-consistent value, and hence superconductivity, vanishes when , way above the correct . To remedy this we need to average the electron motion over disordered thermal configurations rather than solve for electron motion in the disorder averaged configuration.
The Boltzmann weight for the occurence of a particular configuration is
This is related to the electron free energy in a particular background. If the fields are large and random the trace cannot be analytically computed. We generate the equilibrium configurations by a Monte Carlo technique dubi ; dag . This involves diagonalisation of the electron Hamiltonian for every attempted update of the auxiliary fields and the Metropolis condition for acceptance of this move. We use a cluster based implementation of the Monte Carlo tca and can access sizes upto .
The expectation value of any observable can then be computed as follows. For each configuration , the first step is to express the fermionic operators in terms of the quasiparticle operators , with the HFBdG eigenvectors as coefficients. The electronic trace over the resultant expressions in terms of the quasiparticle operators can be easily evaluated. This final expression is then averaged over all thermal configurations . We will describe the indicators separately in detail in their respective sections.
iv.1 Thermodynamic indicators
Fig.1.(a) shows our result for the structure factor corresponding to the growth of superconducting order. We compute the thermally averaged pairing field correlation at . This is like the ‘ferromagnetic’ correlation between the , treating them as two dimensional moments. If the component, , is it implies that the pairing field has a non-zero spatial average and would in turn induce long range order (power law correlation in 2D) in the thermal and quantum averaged correlation . We locate the superconducting transition from the rise in as the system is cooled. The results are not reliable below , since the correlation length becomes comparable to our system size, but compare very well with available QMC data for . We will discuss the interpretation of the results in detail at the end of the paper.
Panel 1.(b) shows the result for showing the clear peak around . We will compare this to the result from QMC, and also discuss the system size dependence, at the end of the paper.
Panel 1.(c) highlights the rapid rise in the ‘gap’ to ratio with increasing interaction. In the weak coupling limit this value is , at it is already , quite beyond BCS, and grows roughly as at large . Needless to say, the gap is not an indicator of the robustness of the superconducting state once we go beyond weak coupling.
The large gap but low leaves its imprint on several physical properties. The phase diagram, Fig.1.(d), highlights this. At weak coupling the vanishing of SC order also means the vanishing of the gap in the density of states. The high regime at small is ungapped (UG). The regime is a ‘renormalised BCS’ window, although the gap to ratio is large. For the phase has a pseudogap (PG), which we show in Fig.2, while for the regime is gapped. Notice that the normal state gap appears before the peak is reached, i.e, on the ‘BCS’ side of the crossover.
iv.2 Background fields
To understand the spatial behaviour of the system and its evolution with and we examine the variation of the background fields and . Fig.2 shows single snapshots of (upper row in each set), normalised by the mean field value, and the phase correlation (lower row in each set), where is the angle at fixed site in the lattice.
Of the three sets in Fig.2, the top set is for , which we will use as typical of ‘weak’ coupling, the middle set is for , typical of intermediate coupling, and the top set is for , strong coupling. The rows are for .
At all the snapshots show almost uniform and perfect phase locking at all . This is just the mean field state. As we move to higher , however, we see a clear difference in the amplitude fluctuations of the three systems. While the plots show an increasing inhomogeneity and a steadily rising value of throughout the system, the case hardly shows any change. The behaviour is intermediate. This shows that with increasing , the system moves smoothly from an amplitude fluctuation dominated regime to one in which amplitudes are effectively constant, the transition being driven by the phase fluctuations.
The phase maps, on the other hand, show how the system breaks up into correlated patches with temperature. The middle column corresponds to , and show large correlated clusters, as expected for a system close to criticality. As is increased, the correlation length decreases, as evident from the right column at .
While it is phase fluctuations that ultimately destroy order at all , the amplitude fluctuations are quantitatively important at weak coupling. To highlight this, we plot the distribution of for the three values at four temperatures , , and for each in Fig.3 (a),(b) and (c). (d) shows the temperature dependence of the mean and its variance for and . We find that the distributions widen for each case with increasing , but the increase is much more pronounced at weak coupling, and decreases systematically with increasing coupling. The distribution is noticeably non gaussian at high temperature in the weak coupling case. The temperature dependence of the mean and width of is shown in Fig.3.(d). A detailed discussion is postponed to the end of the paper.
iv.3 Density of states
Fig.4 highlights the behaviour of the single particle density of states (DOS). We show results at the ‘BCS end’ , near the peak , and in the BEC regime . The results in all cases are described by the canonical DOS, , where is the full gap in the single particle spectrum. There is a gap in the spectrum at all , a ‘coherence peak’ at the gap edges, arising from electron propagation in a perfectly pair correlated background, and a featureless fall at high energies. The oscillatory pattern in the DOS at for is a consequence of finite size, showing up even on a lattice.
While the DOS is just a mean field result, and the dependence at is expected, the dependence at and is not obvious. At , the system has a ‘hard gap’ persisting to . This is reflected in the phase diagram in Fig.1. At both and the transition to SC occurs from a gapped fermion state rather than a Fermi liquid. The DOS indicates that down to values around peak (and even somewhat below) the qualitative physics remains similar to the BEC end.
Another striking feature is the large transfer of spectral weight that occurs on a modest change of temperature. For the case, for example, at there is weight transfer over a scale . The reason is fairly simple: the magnitude in this limit are almost independent, but the phase correlation between them is destroyed at a temperature . As a result, over a small window the system evolves from a state with perfectly ordered , to one where these large amplitudes are randomly oriented. The strong ‘disorder’ in the lead to the broadening of the density of states.
The large size of even in the normal state preserves the gap feature, but the randomness smears the band edges.
iv.4 Spectral functions
Now we turn to the spectral functions. Fig.5 shows intensity plots of the spectral function for (top to bottom) for (left to right). Let us start with the low temperature results, where mean field theory is a good starting point.
Our density involves a non interacting Fermi surface that is almost a square and rotated by with respect to the Brillouin zone. So, the separation between the two branches, , of the mean field dispersion is smallest near (,) and (,). Within MFT the spectral function is given by . In the BCS limit, is either one or zero for or , with a small region around where it crosses from one to the other.
At , for , and for , with significant mixing only near . As a result shows either a lower branch or an upper branch, but not both - except for .
The growing symmetry in the plots, about the horizontal line, with increasing arises from the changing character of and . For , and are both all over the Brillouin zone, since ‘pairing’ is no longer limited to the vicinity of the non interacting Fermi surface (FS). The cases and are already in this regime although some residual asymmetry is visible. The large difference between weak and strong coupling in terms of the pairing amplitude decides the finite state.
At finite , thermal fluctuations broaden the delta functions and the detailed lineshapes for and are shown in Fig.6. As expected, the gap closes for , while it does not for and , though there is a noticeable decrease in the former.
Quantum Monte Carlo work qmc-singer had suggested the presence of non-trivial structure in the spectral function near the zone boundary (). The top row in Fig.6 shows the spectral function at this point for at three temperatures, . The energy is measured in units of the dispersion .
We start with the top row: . At high temperature the case shows only a single broad peak at positive energy, whereas show a second peak at a smaller negative energy value. This two peak structure with a gap around is an indicator of pairing without global coherence. The complete absence at , and the increase in peak height from to bolsters this interpretation.
We do not find any peak near for medium to large . However from downwards, another peak becomes visible at negative energies, at . This peak is indicative of the global coherence setting in below . As is decreased, this peak slowly gains weight while weight in the ‘pairing feature’ becomes smaller, with its maximum shifts to larger negative energies as the gap becomes larger.
Thus, we find some degree of consistency with the QMC work which mainly deals with temperatures larger than , and also find another peak, indicative of global coherence, that starts to develop below .
For and high temperature the spectral functions at and have a gap at , as before, while the result is gapless. In contrast to however we cannot disentangle the lower temperature coherence feature, at from the overall broad band.
Our model has been set up with the explicit constraint that it reproduce the standard mean field (or HFBdG) result at . It ignores quantum fluctuations of the pairing field. The impact of these fluctuations have been discussed using DMFT by Garg et. al. dmft-garg and Bauer and Hewson dmft-bauer . They find that the qualitative results for the order parameter, spectral gap, occupation probability and superfluid stiffness are all given correctly by the mean field method, though it tends to overestimate the spectral gaps and order parameter values at intermediate coupling. For most of the window, however, the mean field results are reasonable. The results from our method should get better at finite temperature as thermal fluctuations become more important than zero point quantum fluctuations. A comparison of our with QMC estimate bears this out.
Accuracy of estimate: Fig. 7(a) compares our with different methods. We find that our results compare well with QMC and sophisticated semi-analytic methods with a slight underestimate at medium to large coupling. Fig.8(b) shows the size dependence of the estimate with data for . We see that while the ‘critical temperature’ decreases noticeably from to , it does not change significantly beyond . Thus, the estimate we obtain at should be a fair approximant to the bulk .
Effective classical functional: The BdG framework involves fermions coupled to the fields and . For simplicity let us focus on the since the do not play a crucial role in a translation invariant system.
The physics of fermions in an arbitrary background is not obvious. It is therefore helpful to have an explicit classical functional involving only the since the minimum and possible fluctuations in are easier to estimate.
If the are small compared to the kinetic energy, as would happen when , the functional, , can be obtained via a standard cumulant expansion:
The superscript in is to indicate character. , being the non-local pairing susceptibility of the free Fermi system, and can be computed from a convolution of four free Fermi Greens functions.
If we had then would have to be expanded to higher order in . In that situation it actually helps to extract the functional by expanding in powers of , leading to the strong coupling limit:
can be obtained from the atomic problem. The leading intersite term is calculated perturbatively and connects only nearest neighbour sites.
While these BCS and BEC limits are easy, obtaining an usable functional at arbitrary does not seem possible. We have therefore tried a parametrisation of the (local) amplitude fluctuation spectrum and the phase correlations in terms of the following phenomenological model. It is valid at all and over the temperature window of interest.
The first term defines an effective XY model involving only the phases, but, as we will see, the needs to be temperature dependent to incorporate the effect of amplitude fluctuations. The amplitude part of is purely local, and to that extent misses out on spatial correlation between amplitude fluctuations.
The parameters and are extracted from a fit to the that we obtain from the full MC, see Fig.3. With the moments of fixed by and , the is obtained by imposing the following equality:
The left hand side is the MC based order parameter, Fig.1(a). The right hand side computes the same quantity within the phenomenological model (in which the and averages factorise) by using a Monte Carlo estimate of in the XY model.
Fig.8.(a)-(b) shows the dependence of and for , and . Because of the large difference in scales between the weak and strong coupling we have normalized the parameters by their values.
The parameters show large change with . The normalized quickly increases and becomes positive, while rapidly decreases from its . Both parameters tend to saturate for . The 6th order term in the expansion would be necessary to describe the case accurately.
With increasing , the thermal change of the parameters slows down, and even at the parameters show a much weaker dependence on . By , they are essentially constant at their values, indicating that only phase fluctuations are relevant in this regime.
To gain more perspective, we consider an expansion of the distribution about its mean value, , where the first term represents the gaussian stiffness of the distribution and the other terms represent non gaussian contributions. At strong coupling, the physics is driven completely by the phase fluctuation term, and the magnitude of the amplitudes is almost fixed. Thus, the stiffness coefficient is very large. As the coupling decreases, amplitude fluctuations increase, hence signalling a decrease in the stiffness. Apart from the increase in amplitude fluctuations, the mean value of also shows a remarkable increase with at weak coupling, signalling the importance of the non gaussian terms in the expansion. We next turn to examine the phase stiffness which is the crucial coupling at large .
Fig.8.(c) shows the dependence of for the three couplings, again normalized by their values, . The case shows a pronounced decrease with , while the other two are effectively constant. An XY description with a independent coupling is reasonable for and but inadequate at .
Role of the ‘density’ field: An important addition in our model is the field , coupling to the density operator. It serves a twofold purpose: first, it is indispensable in a disordered system since it provides a site dependent background field that renormalises the total disorder, and is crucial to get the correct scales; and second, it incorporates fluctuations in the charge sector, which play an important role for .
At the negative Hubbard model can be mapped to its positive counterpart, with the components of the magnetization field, of the positive model corresponding to the and . The symmetry the model is increased from to , and hence, there can be no superconducting order at finite temperature in 2D. At , the superconducting state is degenerate with the charge density wave state. This degeneracy is built into the structure of our model, and simulations at do actually show both superconductivity and charge density wave order at low . The two field decomposition captures the correct ground state and relevant fluctuations in the model.
Handling inhomogeneity: As remarked earlier, this method is particularly well suited to dealing with inhomogeneous systems, including disordered systems and systems in a trap. We have extensively studied both of these, the former in the context of the disorder induced superconductor-insulator transition trans ; spat , and the latter in the context of superconductors in a harmonic trap catom . In such inhomogeneous systems, the Hartree feedback plays a crucial role in modifying the effective potential that the electronic system sees. In the former, this is crucial in determining the correct critical disorder and moderate disorder charge transport properties trans , and plays a major role in the spatial fragmentation of the system spat . In a harmonic trap, the resultant inhomogeneous density profile can drastically alter the spectral properties of the system catom , compared to a flat one. Similarly for FFLO phases in imbalanced Fermi systems the real space treatment on large lattices allow access to a wealth of non trivial modulated phases.
Quantum fluctuations: The major approximation in our model is the neglect of temporal fluctuations in the auxiliary fields. The primary effect is the absence of low energy ‘bosonic’ modes (due to preformed pairs) at strong coupling. This does not affect the thermodynamics and single particle spectrum significantly. Two particle correlations like conductivity also give accurate results as long as we are at moderate . However, as we increase the system develops a pseudo-gap (or a gap) above , and we get a resistivity with steadily increasing insulating character. In the complete treatment the bosonic modes allow a parallel channel of conduction. A purely static approximation misses this contribution at large , as does DMFT.
We have presented results on BCS-BEC crossover in an attractive Fermi system in the context of the two dimensional Hubbard model. We use an auxiliary field decomposition, treat these fields as classical, and solve the resulting problem through a real space Monte Carlo technique. The inclusion of all spatial thermal fluctuations allows us to capture the correct all the way from the BCS to the BEC end. It allows conceptual clarity about the amplitude and phase fluctuation dominated asymptotes and the crucial intermediate coupling window where both these fluctuations are relevant. We provide a detailed characterisation of the auxiliary field behaviour that dictates fermion physics and access results on the density of states and angle resolved spectral features without any need for analytic continuation. We lay the groundwork for the study of disordered superconductors, trapping effects in superfluids, and spontaneous inhomogeneity in imbalanced systems.
Acknowledgments: We acknowledge use of the High Performance Computing Cluster at HRI. PM acknowledges support from a DAE-SRC Outstanding Research Investigator Award.
- (1) For a recent review, see Q. Chen, J. Stajic, S. Tan and K. Levin, Phys. Repts. 412, 1 (2005).
- (2) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
- (3) D. M. Eagles, Phys. Rev. 186, 456 (1969).
- (4) A. J. Leggett in Modern Trends in the Theory of Condensed Matter, Springer-Verlag, Berlin.
- (5) P. Nozieres and S. Schmitt-Rink, J. Low. Temp. Phys. 59, 195 (1985).
- (6) Daniel Rohe and Walter Metzner, Phys. Rev. B 63, 224509 (2001).
- (7) For an early review, see R. Micnas, J. Ranninger and S. Robaszkiewicz, Rev. Mod. Phys. 62, 113 (1990).
- (8) M. Randeria in Bose-Einstein Condensation, Cambridge University Press (1995).
- (9) J. J. Deisz, D. W. Hess, and J. W. Serene, Phys. Rev. B66, 014539 (2002).
- (10) H. Tamaki, Y. Ohashi and K. Miyake, Phys. Rev. A77, 063616 (2008).
- (11) B. Kyung, S. Allen, and A.-M. S. Tremblay, Phys. Rev. B64, 075116 (2001).
- (12) N. Dupuis, Phys. Rev. B 70, 134502 (2004).
- (13) T. K. Kopec, Phys. Rev. B 65, 054509 (2002).
- (14) R. T. Scalettar, E. Y. Loh, J. E. Gubernatis, A. Moreo, S. R. White, D. J. Scalapino, R. L. Sugar and E. Dagotto, Phys. Rev. Lett. 62, 1407 (1989).
- (15) A. Moreo and D. J. Scalapino, Phys. Rev. Lett. 66, 946 (1991).
- (16) A. Moreo, D. J. Scalapino and S. R. White, Phys. Rev. B45, 7544 (1992).
- (17) M. Randeria, N. Trivedi, A. Moreo, and R. T. Scalettar, Phys. Rev. Lett. 69, 2001 (1992).
- (18) N. Trivedi and M. Randeria, Phys. Rev. Lett. 75, 312 (1995).
- (19) S. Allen, H. Touchette, S. Moukouri, Y. M. Vilk and A. M. S. Tremblay. Phys. Rev. Lett. 83, 4128 (1999).
- (20) J. M. Singer, T. Schneider and P. F. Meier, EPJB 7, 37 (1999).
- (21) A. Sewer, X. Zotos and H. Beck, Phys. Rev. B 66, 140504 (2002).
- (22) T. Paiva, R. R. dos Santos, R. T. Scalettar and P. J. H. Denteneer, Phys. Rev. B. 69, 184501 (2004).
- (23) T. Paiva, R. Scalettar, M. Randeria and N. Trivedi, Phys. Rev. Lett. 104, 066406 (2010).
- (24) K. Bouadim, et. al., Nat. Phys. 7, 884 (2011).
- (25) M. Keller, W. Metzner and U. Schwollwock, Phys. Rev. Lett. 86, 4612 (2001).
- (26) M. Capone, C. Castellani and M. Grilli, Phys. Rev. Lett. 88, 126403 (2002).
- (27) A. Garg, H. R. Krishnamurthy and M. Randeria, Phys. Rev. B 72, 024517 (2005).
- (28) J. Bauer and A. C. Hewson, Europhys. Lett. 85 27001 (2009).
- (29) J. Bauer, A. C. Hewson and N. Dupuis, Phys. Rev. B 79, 214518 (2009).
- (30) A. Koga and P. Werner, Phys. Rev. A 84, 023638 (2011).
- (31) V. F. Gantmakher and V. T. Dolgopolov, Phys. Usp. 53, 3 (2010).
- (32) I. Bloch, Nat. Phys. 1, 23 (2005), I. Bloch, et. al., Rev. Mod. Phys. 80, 885 (2008).
- (33) S. Giorgini, et. al., Rev. Mod. Phys. 80, 80, 1215 (2008), R. Casalbouni, et. al., Rev. Mod. Phys. 76, 263 (2004), X-Wen Guan, et. al., Rev. Mod. Phys. 85, 1633 (2013)
- (34) F. Solms, H. G. Miller and R. M. Quick, Phys. Rev. B49, 15945 (1994).
- (35) Look at supplementary material, in K. Bouadim, Y. L. Loh, M. Randeria and N. Trivedi, arxiv: 1011.3275
- (36) G. J. Conduit and Y. Meir, Phys. Rev. B 84, 064513 (2011).
- (37) Y. Dubi, et al., Nature, 449, 876 (2007).
- (38) M. Mayr, G. Alvarez, C. Şen, and E. Dagotto, Phys. Rev. Lett. 94, 217001 (2005).
- (39) P. G. de Gennes, Superconductivity of metals and alloys, Addison Wesley (1989).
- (40) S. Kumar and P. Majumdar, Eur. Phys. J. B, 50, 571 (2006).
- (41) S. Tarat and P. Majumdar, arxiv:1311.6951
- (42) S. Tarat and P. Majumdar, to be published
- (43) Sanjoy Datta, Viveka Nand Singh and Pinaki Majumdar, arxiv:1312.5761