{'fraction_non_alphanumeric': 0.03033488511788557, 'fraction_numerical': 0.009535966361315512, 'mean_word_length': 4.368399838774687, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disordered Hubbard Model with Attraction: Coupling Energy of Cooper Pairs in Small Clusters

We generalize the Cooper problem to the case of many interacting particles in the vicinity of the Fermi level in the presence of disorder. On the basis of this approach we study numerically the variation of the pair coupling energy in small clusters as a function of disorder. We show that the Cooper pair energy is strongly enhanced by disorder, which at the same time leads to the localization of pairs.


I. INTRODUCTION
The superconductor-insulator transition (SIT) in disordered films of metals has attracted widespread interest in recent years 1 . The SIT, driven by adjusting some tuning parameter such as film thickness or magnetic field strength, is particularly interesting in two dimensions (2D), where both superconductivity and metallic behavior are marginal.
In the composite bosons model, Cooper pairs are treated as point-like charge 2e bosons. In this picture, the superconducting state displays a quantum phase transition to an insulating state characterized by a quenching of the condensate of composite bosons 2 . In this scenario, the SIT is caused by the loss of phase coherence between the pairs in different parts of the sample, while the magnitude of the pairing gap remains finite. Numerical studies support this scenario 3 .
However, a few relevant experimental aspects of the SIT 1 seem to be beyond the scope of the composite fermions theory. It is therefore highly desirable to study models where the fermionic nature of charge carriers is not rubbed out from the beginning. The Quantum Monte Carlo studies of the disordered attractive Hubbard model in two dimensions have supported the possibility of a disorder driven superconductor to insulator quantum phase transition 4 . At the same time the mean field approach within the Bogoliubov-de Gennes framework has shown that also space fluctuations of the pairing amplitude should be taken into account in order to give a full picture of the SIT 5 .
In parallel a growing interest has been devoted to the question of what is the coupling energy of pairs placed in small superconducting grains, with the average level spacing of the same order as the superconducting gap [6][7][8][9] . Also the pair properties in small size samples may be related to their properties in the localized phase, where the pair motion is bounded inside the localization domain.
In this paper we study numerically the properties of Cooper pairs in small two-dimensional clusters with disorder. We take the Hubbard attractive interaction between fermionic particles with spin 1/2 which move in a two-dimensional Anderson lattice. Following the approach introduced by Cooper 10 , we consider some part of the particles below the Fermi sea as frozen, while the remaining particles, in the direct vicinity of the Fermi level, can move and interact in the presence of disorder. Recently, such generalized Cooper problem has been considered 11 for the case of two particles in a disordered potential. Here, we further develop this approach for the case of many interacting Cooper pairs, that allows us to study the case of finite particle density. Our numerical studies allow us to determine the dependence of the Cooper pair coupling energy on the strength of disorder. They show that this coupling can be strongly increased by disorder, which however leads to localization of pairs. In the regime of weak disorder the pairs are delocalized but their coupling energy is significantly reduced compared to the localized phase.
The paper is organized as follows: In Section II we introduce the attractive Hubbard model and discuss the numerical method used to study the case of a finite particle density. In Section III we determine the Cooper pair coupling energy (pairing gap) and investigate its dependence on the strength of the disorder. In Section IV we study the disorder-induced pair localization and compare the results with the case of noninteracting particles. In Section V we study the behavior of the superconducting order parameter, obtained from the pairing correlation function. The conclusions are presented in Section VI.

II. MODEL AND NUMERICAL METHOD
We study a disordered square lattice with N fermions on L 2 sites. The Hamiltonian is defined by where c † iσ (c iσ ) creates (destroys) an electron at site i with spin σ, n iσ = c † iσ c iσ is the corresponding occupation number, the hopping term V between nearest neighbors lattice sites characterizes the kinetic energy, the site energies ǫ i are taken from a box distribution over [−W/2, W/2], U measures the strength of the Hubbard attraction (U < 0), and periodic boundary conditions are taken in both directions. We restrict our numerical investigations to the subspace with S z = 0 for even N (N/2 spins up and N/2 spins down) and S z = 1/2 for odd N .
The model (1) at U = 0 reduces to the one body Anderson model, giving localized states in two dimensions at the thermodynamic limit 12 . At W = 0 one gets the clean attractive Hubbard model, which in 2D shows a Kosterlitz-Thouless transition to a superconducting state with power-law decay of the pairing correlations 13 .
We study numerically the model (1) for a finite density of interacting quasiparticles above the frozen Fermi sea: (i) Single particle eigenvalues E α and eigenstates (orbitals) φ α (i) (α = 1, ..., L 2 ) at U = 0 are obtained via numerical diagonalization of the Anderson Hamiltonian.
(ii) The Hamiltonian (1) is written in the orbital basis:

and transition matrix elements
(iii) The Fermi sea is introduced by restricting the sums in (2) to orbitals with energies above the Fermi energy E mF : α, β, γ, δ > m F . We consider a filling factor ν = m F /L 2 = 1/4 (corresponding to 2m F frozen electrons due to spin degeneracy) and a finite density of N interacting quasiparticles above the Fermi level.
(iv) The Slater determinants basis, built from the single particle orbitals φ α , is energetically cut-off by means of the condition gives an effective phonon frequency ω D ∝ M/L 2 .
(v) The ground state of this truncated Hamiltonian is found via the Lanczos algorithm 14 .
In our numerical simulations we considered up to N = 8 interacting quasiparticles in up to M = 20 orbitals. We checked that results are qualitatively similar under variation of the cut-off orbital M . In the following Sections we present data for U = −4V , averaged over N R = 100 disorder realizations.

III. PAIRING GAP
In order to compute the pairing energy, we first compute the total N -body pairing energy as with E g (U ) many body ground state for an attractive interaction U .  In Fig. 1 we show the N -body pairing E p (N ) as a function of the number N of interacting quasiparticles above the Fermi sea, at different disorder strengths. This figure shows a clear even-odd effect, with a much larger increase of the pairing energy when N is pair. This fact has a clear meaning: for N even, it is possible to build a new pair, reducing the ground state energy due to the negative coupling U < 0. For N odd, in an ideal BCS superconductor the additional particle cannot be paired and remains as a quasiparticle excitation. However, a small ground state energy reduction is still present in our numerical simulations, since the unpaired particle weakly interacts with the superconducting pairs.
The jump in the pairing energy from odd to even number of particles, with N even, can therefore be interpreted as the energy necessary to break a superconducting pair; in the ideal BCS case, this would give the superconducting energy gap. We note that the superconducting gap is extracted in a similar way in experiments with single Cooper pair tunneling inside superconducting islands 15,16 .
In Fig. 2 we show the pairing gap ∆(N ) as a function of the disorder strength W , for N = 2, 4, 6, 8. We see that, with the exception of the first jump (N = 2), the other jumps are rather similar. It is clear that ∆ grows significantly with the disorder strength W . We attribute this effect to the fact that at strong disorder particles are trapped in the deepest minima of the random potential. Therefore the pair size becomes smaller that enhances the interaction between coupled particles, hence ∆. In Fig. 3 we show the dependence of ∆ on the system size 6 ≤ L ≤ 10. Since the Debye cut-off frequency should be independent of the system size, we keep constant the ratio M/L 2 ≈ 0.2. The coupling energy ∆ becomes independent of the system size at large W , while it is not yet saturated at small W . As the pair size is determined by 1/∆, this means that at small W the size of the pair becomes comparable with the system size.

IV. PAIR LOCALIZATION
In order to study the localization properties of the system, we consider the fraction ξ of the sample occupied by the N -body wavefunction |Ψ g : where is the charge density of the ground state at the site i.
With the definition (6), N/2L 2 ≤ ξ ≤ 1, the lower limit corresponding to pairs localized in a single site, the upper limit to complete charge delocalization. In Fig.4 we show the fraction of occupied sites ξ as a function of disorder, at different system sizes 6 ≤ L ≤ 10, for N = 6 particles and a fixed Debye frequency (M/L 2 ≈ 0.2). This figure gives a clear indication of the presence of two regimes: at small disorder the wavefunction fills a large fraction of the sample (superconducting regime), while at large disorder ξ decreases with the system size (ξ ∝ 1/L 2 , localized regime). In the inset of Fig. 4 we show the parameter ξ at different system sizes L = 8, 10, and for a constant electronic density of mobile fermions N/L 2 ≈ 0.06. The two curves puts on top of each other, suggesting the existence of a size-independent function ξ(W ) in the thermodynamic limit. The drop of ξ(W ) with W demonstrates that disorder gives localization of Cooper pairs.
Finally we examine the question of to what extent the wavefunction localization is a many body effect instead of a single particle Anderson localization phenomenon. Therefore in Fig. 4 we also show the parameter ξ at U = 0. The comparison between the interacting (U = −4) and the noninteracting (U = 0) case suggests that the interaction makes localization stronger, in agreement with results for two particles in a three dimensional random potential 11 . This effect can be explained qualitatively with the following argument 11 : attractive interaction creates pairs of total mass m p twice the electronic mass. This halves the effective hopping term V eff ∝ 1/m p , thus doubling the ratio W/V ef f . We remark that this rough argument fails in the delocalized regime at small disorder, where the tendency seems to be reversed. Actually, due to localization of single particles states in 2D 12 , for a given number N of particles one should get ξ(U = 0) → 0 when L → ∞.

V. PAIR CORRELATION
The superconducting state can be characterized by the s-wave pair correlation function, where creates a pair at site i. For an s-wave superconducting state, is the order parameter of the superconductor-insulator transition 4 . In Fig. 5 we show the dependence of the order parameter ∆ op on the number of interacting quasiparticles above the Fermi sea. The order parameter is strongly suppressed by disorder (see also Ref. 4 ), an effect which becomes more evident with the addition of particles. We remark that ∆ op shows an approximate linear increase with the number of pairs. This is quite natural if the many body ground state wavefunction is in the BCS form 17 , built from single particle eigenfunctions including disorder 18 : with b † γ = d † γ↑ d † γ↓ , and g γ variational parameters. Using the relation c † iσ = α φ α (i)d † ασ , after lattice and disorder averaging the dominant contributions in P s (r) is proportional to the number of pairs, ∆ op (N ) ∝ α N 2 (here we have taken into account that g γ in the BCS theory changes smoothly with α around the Fermi level 17 and here we considered N ≪ 2n F ). Here α is a parameter which determines the slope of ∆ op variation with N . The order parameter decreases when an unpaired particle is added. In our opinion, this is due to the fact that this extra electron weakly interacts with the paired particles, reducing the pair correlation function.
In Fig.6 we show the slope α of the linear fit of the order parameter ∆ op as a function of the number N of quasiparticles. The suppression of this quantity with disorder is evident, indicating a rather sharp crossover from a superconducting to an insulating behavior in our finite size lattice.

VI. CONCLUSIONS
In this paper we have investigated the localization of Cooper pairs for small clusters in a two-dimensional disordered substrate. We have shown that the Cooper pair coupling energy displays an even-odd asymmetry: this parity effect survives also in the presence of disorder. The pairing gap is strongly enhanced by disorder, which at the same time leads to localization of Cooper pairs (gapped insulator). Therefore, in the insulating regime, the breaking of Cooper pairs should enhance transport. This is consistent with the resistivity drop observed in experiments with an applied magnetic field 19,20 , which might signal the crossover from a Cooper pair insulator to an electronic insulator.