Near-degeneracy of several pairing channels in multiorbital models for the Fe-pnictides

Weak-coupling approaches to the pairing problem in the iron pnictide superconductors have predicted a wide variety of superconducting ground states. We argue here that this is due both to the inadequacy of certain approximations to the effective low-energy band structure, and to the natural near-degeneracy of different pairing channels in superconductors with many distinct Fermi surface sheets. In particular, we review attempts to construct two-orbital effective band models, the argument for their fundamental inconsistency with the symmetry of these materials, and the comparison of the dynamical susceptibilities in two- and five-orbital models. We then present results for the magnetic properties, pairing interactions, and pairing instabilities within a five-orbital Random Phase Approximation model. We discuss the robustness of these results for different dopings, interaction strengths, and variations in band structure. Within the parameter space explored, an anisotropic, sign-changing s-wave state and a d_x2-y2 state are nearly degenerate, due to the near nesting of Fermi surface sheets.


I. INTRODUCTION
The undoped Fe-pnictides are semi-metallic materials which exhibit structural and spin density wave (SDW) antiferromagnetic transitions. Photoemission 1 and density functional theory (DFT) calculations 2,3,4,5 find that over an energy range near the Fermi energy the electron bands are made of states from the Fe-pnictide layers ( Fig. 1) of predominantly Fe character. The two dimensional nature of these layers is such that the Fermi surfaces consist of hole cylinders around the Γ point and electron cylinders around the M point of the 2 Fe/cell Brillouin zone (see Fig. 2). In the undoped system, the near nesting of the hole and electron Fermi surfaces can give rise to a colinear antiferromagnetically ordered state within a weak coupling approximation 6 , and this state has been confirmed in neutron experiments 7 . In both the electron and hole doped cases, superconductivity appears in proximity to or coexists with the antiferromagnetic SDW order 7,8,9,10 . It is therefore natural to consider the possibility that spin fluctuations provide the pairing mechanism and various random phase 11,12,13,14,15 (RPA), fluctuating exchange 16,17,18 (FLEX) and renormalization group 19,20 calculations have been reported. These have made use of different approximations to the band structure and obtained a variety of different gap structures.
At the same time, various experimental probes of different Fe-pnictides appear to indicate quite different symmetries of the order parameter 21,22,23,24,25,26,27,28,29,30,31,32 ; as discussed below, there is a real possibility that different symmetries are realized in different materials. It is therefore important to investigate the origin of the discrepancies among the various theories with regard to their predictions for the superconducting ground state. In the case of the cuprates, where a single Cu-O hybridized band predominates at the Fermi surface, different methods and band structures all lead to d-wave pairing. If the situation is qualitatively different in the Fe-pnictides, we need to isolate which of the apparently small differences in the various approaches is responsible for the different results obtained. The reward for such an effort may be considerable, since the evidence suggests that correlation effects in these materials may be modest compared to the cuprates, raising the possibility of a quantitative theory of the superconducting state.
Here we describe the results of an RPA calculation of the pair coupling strength and the momentum dependence of the gap function for a multiorbital tight binding parametrization of a DFT band structure for the LaOFeAs material by Cao et  al. 4 . We begin in Sec. II by discussing the two orbital model for the electronic structure of the pnictides that was introduced by some of the present authors and has been studied by various groups, along with some of the inconsistencies associated with this approach. In Sec. III, we discuss our results for the electronic structure using a five-orbital fit to DFT. Then in Sec. IV we examine the spin and orbital susceptibilities and the effective pairing interaction for both models. In Sec. V, an effective pairing interaction strength λ and gap function g(k) are introduced and in Sec. VI we discuss our results for λ and g(k). In Sec. VII, we examine the spatial and orbital structure of the pairing interaction and the resulting pairs. Our conclusions and comparison with previous work are contained in the final section VIII.

II. TWO-ORBITAL MODEL
Effective models of the band structure based entirely on Fe orbitals should be possible because DFT tells us that the states due to the As 4p orbitals are located approximately 2 eV below the Fermi level. The As orbitals allow for hybridization with the Fe 3d states, however, and therefore an effective Fe-Fe hopping Hamiltonian can be constructed, provided the symmetries of the entire FeAs layer are respected. It is conceptually simplest to work therefore with a square lattice with sites corresponding to Fe atoms and introduce a set of effective hoppings between these sites. If the primitive unit cell is taken to be a square containing a single Fe atom, the effective Brillouin zone is the square shown in Fig. 3a. Since the true primitive unit cell contains two Fe atoms, the true Brillouin zone is smaller by a factor of two, as shown in Fig. 3b. The Fermi surface in the correct small zone should be obtained from the effective large zone by folding the effective zone about the black dashed lines in Fig. 3a. Sheets around the M point in Fig. 3a thus correspond to sheets around the Γ point in Fig. 3b.
In Qi et al. 12 , a 2-orbital tight binding model was used to carry out an RPA calculation of the spin and orbital fluctuation pairing interaction. Here we briefly describe this model and discuss why it is insufficient to approximate the full band structure, especially with regard to the correct orbital weights along the Fermi surface sheets. The Hamiltonian for the two band model can be written as is a two-component field which describes the two "d xz " and "d yz " orbitals. The energies ǫ ± (k) and ǫ xy (k) are parameterized in terms of four hopping parameters t i : Here, γ νσ (k) destroys an electron in the ν = ± band with spin σ and With the addition of a chemical potential the Hamiltonian takes the following form where the band energies are E ± (k) = ǫ + (k) ± ǫ 2 − (k) + ǫ 2 xy (k). Measuring energy in units of |t 1 |, for t 1 = −1, t 2 = 1.3, t 3 = t 4 = −0.85 and µ = 1.5 we exhibit the Fermi surface in Fig. 4, displayed in the large effective zone. To obtain the corresponding Fermi surface in the small zone, a folding across the line joining (π, 0) and (0, π) is required, as discussed above. The Fermi surfaces α 1 around (0, 0) and α 2 around (π, π) are hole pockets associated with E − (k f ) = 0 and the β 1 and β 2 Fermi surfaces around (π, 0) and (0, π) are electron pockets from E + (k f ) = 0. We note that the displacement of the α 2 Fermi surface from the Γ-point to (π, π) is an artifact of the 2-orbital approximation. As shown by various DFT calculations 3,4,5 and noted by Lee and Wen 33 , the orbital states that have significant weight near the β Fermi surfaces include, in addition to the d xz and d yz orbitals, a d xy orbital. In addition, while the Fermi surface sheets shown in Fig. 4 fold down to give two hole Fermi surfaces around the Γ point of the 2Fe/cell Brillouin zone, there should in fact be two hole Fermi surfaces around the Γ point of the large, effective Brillouin zone. This is known from the wave functions found in the band structure calculations. Finally, the 2-orbital model lacks the flexibility to fit the Fermi velocities found in the band structure calculations, giving Fermi velocities on the electron β Fermi surfaces that are anomalously small compared to the velocities on the hole Fermi surfaces.

III. 5-ORBITAL MODEL
In principle, one could capture the correct behavior near the α and β Fermi surface sheets by treating a 3-orbital d xz , d yz , d xy model. However, with short range hoppings, this leads to the appearance of an extra unphysical Fermi surface and a fourth orbital is required to remove it 33 . In addition, recent theoretical calculations using 5 Wannier d-orbits per Fe site find that one can obtain an excellent representation of the electronic structure within a ±2eV window of the Fermi energy 11 . Furthermore, the values of the Coulomb interaction parameters obtained for the Wannier basis are such that the average Coulomb interaction is small compared to the bandwidth, implying that one is dealing with a weakly coupled system 34 .
At this point it therefore seems best to use all five d orbitals in developing a tight-binding model 11 . Here, using a Slater-Koster based parametrization which respects the symmetry of the FeAs layers, we fit a five-orbital (d xz , d yz , d xy , d x 2 −y 2 , and d 3z 2 −r 2 ) tight binding model to the DFT band structure determined by Cao et al. 4 . We will use an orbital basis that is aligned parallel to the nearest neighbor Fe-Fe direction rather than the Fe-As direction. With this choice we avoid the necessity of introducing a second, rotated coordinate system in addition to the one that is used to describe the single Fe unit cell.
The Hamiltonian for the 5 band model takes the following form Here d † m,σ (k) creates a particle with momentum k and spin σ in the orbital m. The kinetic energy terms ξ mn (k) together with the parameters for a 5 band tight binding fit of the DFT band structure by Cao et al. are listed in the appendix. A diagonalization of this Hamiltonian yields the eigenenergies and the matrix elements analogously to the two band case discussed above.
In Fig. 5 (a), we have plotted the resulting band structure in the backfolded "small" Brillouin zone while the Fermi surface sheets for zero doping are shown in Fig. 5 (b). The colors correspond to the dominant orbital weight of each band in momentum space. The gray lines represent the DFT band structure by Cao et al. and the comparison shows, that the 5 band fit approximately reproduces the DFT bands, especially in the vicinity of the Fermi level. It is obvious that the d xy contribution plays an important role in building up the electron-like Fermi surfaces (β sheets) around the M point of the "small" Brillouin zone. In the unfolded Brillouin zone ( Fig. 5 (b) the d xy orbital and d yz (d xz ) orbital contribute the dominant weights to the band states at the β 1 (β 2 ) Fermi surface. To confirm this we also show the orbital weights as a function of the winding angle on the different Fermi surfaces in Fig. 6.

A. Noninteracting susceptibilities
In this section we examine the RPA enhanced spin and charge susceptibilities for the multiorbital model that we introduced in the previous section. The spin operator for an orbital s is defined as where α and β are spin indices. The spin susceptibility can then be calculated from the Matsubara spin-spin correlation

function
with τ the imaginary time and ω a Matsubara frequency. In the same way we can define the Fourier component of the charge density for the orbital s as and we can calculate the charge susceptibility from In a more general formulation the susceptibilities are functions of four orbital indices. For the non-interacting case (χ 0 ) pq st and (χ 1 ) pq st are equivalent and can be written Now we can derive an explicit expression for the noninteracting susceptibilities from G sp (k, iω n )G qt (k + q, iω n + iω) (12) where N is the number of Fe lattice sites, β = 1/T , and the spectral representation of the Green's function is given as Here the matrix elements a s µ (k) = s|µk connect the orbital and the band space and are the components of the eigenvectors resulting from the diagonalization of the initial Hamiltonian. Performing the Matsubara frequency summation and setting iω n → ω + i0 + we find the retarded susceptibility B. Random phase approximation for multiorbital system Now we consider Coulomb interactions of the electrons on the same Fe atom in an RPA framework. We distinguish between an intraorbital interaction U of electrons in the same orbital and an interorbital interaction of electrons in different orbitals V . We can also take the Hund's rule coupling into account that favors the parallel alignment of electron spins on the same ion and is described by an energy J > 0 and the pair hopping energy denoted by J ′ . These interactions are generated automatically in multiorbital models with general twobody interactions using a Hubbard-type approach restricted to intrasite processes 35,36,37 . Thus in general one can write where n is = n i,s↑ + n is↓ . We have separated the intraorbital exchange J and "pair hopping" term J ′ for generality, but note that if they are generated from a single two-body term they are related by J ′ = J/2. If we now split the interaction Hamiltonian into singlet (H and where σ is = 2 S is , as well as The RPA susceptibilities are obtained in the form of Dysontype equations as and where repeated indices are summed over. Here the non-zero components of the matrices U c and U s are given as where a = b. Our notation of the interaction parameters U , V , J and J ′ can be compared to the notation in Kubo 37 as whereŨ ,Ṽ ,J andJ ′ denote the interaction parameters introduced in Ref. 37. In Fig. 7 we show results for the static, homogeneous bare spin susceptibility as a function of q in the first quadrant of the effective "large" Brillouin zone as well as cuts along the main symmetry directions. We compare calculations that have been performed correctly including the matrix elements to calculations where the matrix elements have been considered to be constant a s µ (k)a p * µ (k) = 1/5, as has been reported in certain ab initio electronic structure calculations. Here we used the five band fit discussed in the previous section and chose a chemical potential that corresponds to the undoped compound (x = 0). As can be seen from Fig. 7 the matrix elements play a very important role by "filtering" special structures of χ S (q) in momentum space, e.g. suppressing the weight for small q. The matrix elements for the d x 2 −y 2 and d 3z 2 −r 2 are very small on the Fermi surfaces ( Fig. 6) so that the actual χ S (q) is reduced when these matrix elements are properly taken into account. Neglecting the matrix elements will lead to a wrong result and generally to a higher value and more homogeneous distribution of the susceptibility in momentum space. Important features like the Q = (π, 0)-peak, that is responsible for the antiferromagnetic SDW instability, can be under-or overestimated.
For the single band susceptibility the inclusion of interactions within an RPA approach is known to enhance existing structures in the bare susceptibility as U χ ′ 0 (q) approaches 1. In the case of a multiorbital susceptibility with onsite intraorbital repulsion U , onsite interorbital repulsion V and Hund's rule coupling J it is not obvious how the different structures in the spin and in the charge susceptibility are changed by the variation of these three parameters. For a simplified and more transparent discussion we will first study the susceptibility as a function of U and we will choose V = U and J = J ′ = 0. In Fig. 8 we compare the spin and charge susceptibilities of the electron-doped compound with x = 0.125 for a value of U = 1.65 that is chosen such that the main peaks in the spin susceptibility are considerably enhanced. Here, as elsewhere in this paper, energy units are in eV. For this value of U we find that the charge susceptibility is more than one order of magnitude smaller than the spin susceptibility. In addition the charge susceptibility has no pronounced structures in momentum space whereas the spin susceptibility shows two distinct peaks between (π, 0) and (π, π). While in the undoped compound (see Fig. 7) we find the main peak in the spin susceptibility at the antiferromagnetic wave vector Q = (π, 0) (corresponding to (π, π) in the "small" BZ) we observe a shift of this peak towards an incommensurate wave vector Q * = (π, 0.16π) in the electron-doped com- pound (x = 0.125). This shift can be attributed to an imperfect nesting of the α and β FS sheets due to the opposite growth of the electron and hole sheets with doping.
When we finally compare in Fig. 9 the change of the spin susceptibility as a function of U for the antiferromagnetic wave vector Q = (π, 0) to the one for the incommensurate wave vector Q * = (π, 0.16π) we find a moderate and nearly The temperature dependence of the magnetic susceptibility and NMR spin-lattice relaxation rate in the paramagnetic phase are immediate tests of our RPA calculation of χ S (q, ω). The susceptibility of the F-doped LaOFeAs 1111 material which we study here has been measured both by direct magnetization and NMR methods 38,39,40 and found to increase quasi-linearly up to several hundred degrees K. In addition, Zhang et al. have pointed out that this behavior appears to hold for the 122 class of materials as well, and that the slope of the linear-T behavior appears to be roughly dopingindependent 41 . Within our theory, we can calculate the homogeneous, static bare spin susceptibility by summing over spins and using Eq. 22. At T = 0 it reduces to which is of course just the Pauli susceptibility proportional to the total density of states at the Fermi level (N ν (0) is the single-spin density of states at the Fermi level in band ν.). At finite temperatures one might assume that the susceptibility of an itinerant electron system should decrease. This is the usual case for a single parabolic band, but in the presence of band structure effects the susceptibility may first increase. For example, in a single band system the band structure enters in a simple way as density of states at the Fermi level. Our situation is also somewhat unusual as the chemical potential sits in a region of rapidly varying density of states, so one may expect the T dependence of the susceptibility to differ from the textbook results.
We now plot in Fig. 10 the temperature dependence of both χ 0 and the RPA enhanced χ for the undoped and doped cases, for a particular choice of interaction parameters. It is seen that in both the doped and undoped cases, the susceptibility indeed increases quasilinearly over the experimentally interesting range of 100 to 500K. This qualitative result appears to be relatively independent of the choice of interaction parameters since it arises from the band structure. The RPA enhancement for the cases shown is of order 30%. While this factor increases as one approaches the critical U at which the Stoner instability is reached, it cannot increase dramatically, because for the electronic structure of the current model, the instability at large q occurs first. Thus a large enhancement (Ref. 42) of the q = (0, 0) susceptibility does not occur in this model. We conclude that the observed susceptibility contains a significant temperature-independent interband or van Vleck component which is not included in the χ S calculated here. Similar conclusions were reached by the authors of Ref. 39.
We also calculated the NMR relaxation rate 1/(T 1 T ) that is proportional to to the imaginary part of the local (q integrated) spin susceptibility where A(q) is a geometrical structure factor and ω L is the NMR frequency. In Fig 11 we show 1/(T 1 T ) as a function of temperature for the same set of parameters that we have used for the q = 0 susceptibility. Here we find for the interacting system an initial decrease of the relaxation rate as a function of temperature that is most significant for the undoped compound with finite values of J. In the latter case we find a very distinct upturn around 300 K followed by a slight and nearly linear increase. For the other cases we find a similar tendency, although not as pronounced. In panels 11 (a) and (b) we take the hyperfine form factor of the 19 F nucleus used e.g. in Ref. 40 to be constant, since the F is located directly above the Fe in the LOFFA material. In panels (c) and (d), we have modelled the form factor of the 75 As nucleus as A(q) = cos q x /2 cos q y /2, which suppresses the (π, 0) fluctuations, since the As atom is in the center of a square containing four Fe.  38,39,40 , although the increase with increasing temperature for both quantities found here is somewhat weaker than observed. We have not attempted to fit experiments in detail.

V. THE EFFECTIVE INTERACTION AND THE GAP FUNCTION
Assuming that the pairing interaction responsible for the occurrence of superconductivity in the iron pnictides arises from the exchange of spin and charge fluctuations, we can calculate the pairing vertex using the fluctuation exchange approximation 43 . For the multi-orbital case 36 , the singlet pairing vertex is given by The χ RP A 1 term describes the spin fluctuation contribution and the χ RP A 0 term the orbital (charge) fluctuation contribution. Fig. 12 shows the Matsubara frequency dependence of these two terms for a momentum transfer (π, 0) and a typical interaction strength. Here one sees the that the spin fluctuation contribution is dominant and falls off on a frequency scale which is small compared with the bandwidth. Thus while the gap equation depends upon the kernel ImΓ pq st (k, k ′ , ω)), the important k and k ′ values are restricted by this frequency cutoff to remain near the Fermi surfaces. Then, just as for the electron-phonon case, the strength of the pairing interaction is characterized by a frequency integral of this kernel weighted by ω −1 . Making use of the Kramers-Kronig relation we can proceed further by considering only the real part of the ω = 0 pairing interaction. If we now confine our considerations to the vicinity of the Fermi surfaces we can determine the scattering of a Cooper pair from the state (k, −k) on the Fermi surface C i to the state (k ′ , −k ′ ) on the Fermi surface C j from the projected interaction vertex where the momenta k and k ′ are restricted to the different Fermi surface sheets with k ∈ C i and k ′ ∈ C j . If we decompose the superconducting gap into an amplitude ∆ and a normalized symmetry function g(k) we can define a dimensionless pairing strength functional 44 as Here Γ ij is only the symmetric part of the full interaction, which gives identical results within the spin singlet subspace. The Fermi velocity is defined to be v F (k) = |∇ k E ν (k)| for k on the given Fermi surface.
From the stationary condition we find the following eigenvalue problem The kernel Γ ij (k, k ′ ) is evaluated at temperatures well below the characteristic temperature at which the spin fluctuation spectrum has formed. In this temperature region, the interaction is independent of temperature. From the above eigenvalue problem we will determine the leading eigenfunction g α (k) and eigenvalue λ α for a given interaction vertex Γ ij (k, k ′ ). The largest eigenvalue will lead to the highest transition temperature and its eigenfunction determines the symmetry of the gap. The next leading eigenvalues and eigenfunctions further characterize the pairing interaction and can indicate the structure of possible collective states. With the knowledge of the pairing function g α (k) we can also determine the individual contributions of the different intra-and interorbital scattering processes to the total pairing strength λ α .

A. Eigenvalue problem
In this section, we present and discuss results for the pairing strength λ and the symmetry function g(k) obtained within our weak-coupling approach.
Solving the eigenvalue problem given in Eq. 30, we find a set of eigenfunctions g α (k) defining the k dependence of the gap on the FS sheets together with a set of corresponding eigenvalues λ α denoting the dimensionless pairing strengths associated with a given g α (k). We first have to classify the different eigenfunctions according to basic symmetry operations. In the following we will speak of an s wave state if while we have a d x 2 −y 2 wave state if g(−k x , k y ) = g(k x , −k y ) = g(k x , k y ) g(k y , k x ) = −g(k x , k y ) Since there are two different d x 2 −y 2 wave states among the first few eigenfunctions, we have to distinguish them further, e.g. by comparing the sign of g ν (k) for ν = α 1 and ν = α 2 in the same direction in momentum space. Here we will label the state that changes sign between α 1 and α 2 with d x 2 −y 2 (1), the one without sign change with d x 2 −y 2 (2). Furthermore we can distinguish a d xy wave state that is given by and a g wave state with First we consider a case where the Hund's rule coupling and the pair hopping energy are negligible compared to the intraand interorbital Coulomb interactions, so we set J = J ′ = 0 and V = U . In Fig. 13 (a) we show the pairing strength eigenvalues λ α for the four leading eigenvalues as a function of U for the electron-doped compound (x = 0.125). Approaching the critical value of U where the eigenvalues start to diverge we find a clear separation of the two leading eigenvalues from the next two eigenvalues. However the two leading eigenvalues, corresponding to the s and the d x 2 −y 2 symmetry remain very similar in size and we find a crossover from the d x 2 −y 2 to the s symmetry around U = 1.65. We also show the s wave (d) and the d x 2 −y 2 wave (e) pairing functions on the four FS sheets for U = 1.73 close to its critical value 45 . The extended s wave state (which we have labelled s) is characterized by i) a sign change on the β FS sheets with nodal points displaced from the generic (0, π)-(π, 0) direction that would result from a pure cos k x +cos k y state, and ii) a change of average sign on the β sheet relative to that on the α sheets. On the α FS sheets we find a nodeless but anisotropic gap distribution with higher weight on the small α 1 compared to the larger α 2 sheet. The d x 2 −y 2 state features an anisotropic gap distribution on the β FS sheets and a rather conventional d wave gap distribution on the α FS surfaces, with a sign change between the α 1 and α 2 sheet. In Fig. 14 we show the corresponding results for the undoped (x = 0) compound. Here we find similar results, although the eigenvalues close to the instability are better separated and the crossover from d x 2 −y 2 to s appears already for a rather small value of U . If we compare the s wave symmetry function for x = 0 close to the instability (see Fig. 14 (d)) to the corresponding symmetry function for the electron-doped compound, we find that the nodal points on the β FS have moved even closer to the tips of the sheets and the negative dip between them has decreased considerably. For the d x 2 −y 2 symmetry (see Fig. 14 (e) we find that the weight on the α 2 FS sheet has nearly vanished.
In Figs. 13 and 14, as well as in the following Figs. 16 and 15 we compare in panels (b) and (c) the contributions of the different intra-and interband processes. Here λ αα sums up all contributions of λ ij resulting from scattering within each α sheets and in between the two α sheets as λ αα = i,j=1,2 λ ij  where i = 1 refers to the α 1 , i = 2, to the α 2 , i = 3 to the β 1 and i = 4 to the β 2 Fermi surface. λ intra ββ = λ 33 + λ 44 denotes the intraband contributions of the β FS sheets, while λ inter ββ = λ 43 + λ 34 is the corresponding interband contribution. Finally λ αβ sums up all the remaining contributions, resulting from scattering between the α and β Fermi surfaces. It is obvious that in the case of J = J ′ = 0 for both of the main pairing symmetries the intersheet λ αβ contribution is responsible for the rapid increase of the pairing strength λ with U , reflecting the nesting between the α and β FS sheets with nesting vector Q * . All other contributions are small and mainly negative. For finite values of J and J ′ , which we will next discuss we find that for the d x 2 −y 2 pairing symmetry the interorbital contributions between the two β sheets become also important and for the electron-doped compound even dominant. Now we consider a case of finite J and J ′ . Here we will choose J ′ = J/2 and we will fix V = U − 3/4J − J ′ . These choices are consistent with the generation of intrasite couplings from a Hubbard argument 35,37 , as mentioned above, as well as with the range of values for J/U found by Anisimov et al. in an ab initio calculation 34 . In Fig. 15 (a) we show for the undoped compound (x = 0) and for J ′ = J/2 = U/4 the same four eigenvalues as a function of U as in the previous figures. Here we find again that the extended s wave state (d) is the most stable pairing configuration followed by the d x 2 −y 2 wave state (e). For the d x 2 −y 2 state the main contributions to the pairing function are along the β FS sheets, while the α FS sheets contribute less significantly. Since this state has no nodes along the β sheets, it can be considered as a mainly nodeless state, in contrast to the states found for J = J ′ = 0. The same is true for the electron-doped case with x = 0.125 (Fig. 16). Here the d x 2 −y 2 state is the most Although the pairing is constrained to the close vicinity of the Fermi surfaces, one can try to find an approximation of the g ν (k) on the different FS sheets that extends further into the Brillouin zone. Here we want to restrict ourselves to the smallest number of harmonics necessary to find a reasonable fit of the pairing function on the different FS sheets. For the s wave symmetry we can write g ν (k) = 2A ν [cos k x + cos k y + 2w ν,xy cos k x cos k y +2w ν,4x4y cos 4k x cos 4k y ] where we find for U = 1.5 the following parameters: A α1 = 0.051, w α1,xy = −0.35, and w α1,4x4y = −1.35 for the α 1 sheet, A α2 = 0.02, w α2,xy = −0.1, and w α2,4x4y = 1 for the α 2 sheet, and A β = 0.15, w β,xy = 0.17, and w β,4x4y = −0.1 for the β sheets. In the same way, we can find an approxima- tion to the d wave symmetry function as with A α1 = −1.2 and w α1,2x = 0 for the α 1 sheet, A α2 = 0.12 and w α2,2x = 0, for the α 2 sheet, and A β = −0.018 and w β,2x = −1.12 for the β sheets. The results of the fitting are shown in Figs. 17 and 18.

D. Role of nesting
To study the influence of nesting on the spin susceptibility, especially on the (π, 0) peak, we slightly modify the hopping parameters, creating a toy model with perfectly nested FS sheets. This means that we try to find an approximation where the α 2 and the β FS sheets are of approximately the same size and shape. The band structure used within this toy model is very similar to the band structure found formerly by   a fit to the DFT bands, to assure that all basic properties of the model are reasonable, i.e. that the matrix elements are similar to the correct matrix elements found from the DFT fit. In Fig. 19 we show the near circular FS sheets of this model. The bare spin susceptibility χ 0 (q) for this model in the undoped case is exhibited in Fig. 20, and may be compared to the results of the original model in Fig. 7. We see that the (π, 0) peak remains but is not particularly enhanced by the enhanced nesting, implying that the original model was already quite close to a nesting condition. On the other hand, some subtle differences which are quite interesting appear when we examine the pairing functions. If we consider an even more Gedanken-type model where the two α sheets have degenerate radius a, which is also the same as the two β sheets, we see that the Fermi surface in the 1st effective Brillouin zone is invariant under a translation by (π, 0). Under this transformation, simple extended-s and d x 2 −y 2 functions cos k x + cos k y and cos k x − cos k y map into one another identically. We therefore expect that the pairing eigenvalues for s and d x 2 −y 2 will become degenerate. The toy model band structure considered here and shown in Fig. 19 is nearly the same as the Gedanken model, but has two slightly nondegenerate α sheets. We see nonetheless in Fig. 21 (b,c) and (d,e) that the two com- peting states are indeed nearly exactly degenerate. Fig. 21 (a) exhibits the dependence of the leading eigenvalues on U .

VII. THE SPATIAL AND ORBITAL STRUCTURE
The gap function g ν (k) contains information on the spatial and orbital structure of the pairs. In the previous section, we have calculated g ν (k) and found that its behavior on the Fermi surfaces can be fit by low order harmonics for both the s and d-wave eigenstates. Here we use these results to determine simple pictures of the internal pair structure. A gap operator can be written as with γ νσ (k) the destruction operator for an electron in the ν th band with wave vector k and spin σ. Using the band-orbital matrix elements to relate the band operator to the orbital operator γ νσ (k) = n a n ν (k)d nσ (k) (38) we have with a nm (l) = 1 N kν g ν (k)a n ν (k)a m ν (−k)e −ikl (40) where l = l 1 −l 2 . The amplitude a nm (l) describes the internal spatial and orbital structure of the pair. Using the harmonic approximations for g ν (k) for the s and d x 2 −y 2 gaps discussed in the previous section, we have calculated a nm (l). In principle the k sum should be cut-off when k − k F exceeds the inverse coherence length. Here we have used a Gaussian cutoff of the k sum which decays when |k − k F | exceeds 2π/ξ 0 , where ξ 0 is taken to be of order 3 times the Fe-Fe spacing a. This provides a local picture of the internal orbital structure of a pair. This basic structure continues out to a radius set by the coherence length ξ 0 .
We find that the off-diagonal amplitudes a nm are negligible compared to the diagonal ones, and that the orbitals that contribute are the d xz , d yz , and d xy orbitals which have weight near the Fermi surfaces. The amplitudes a nn (l) for the s and the d x 2 −y 2 gaps are shown in Figs. 22 and 23, respectively. For the s case, one sees that the internal structure of a pair consists of a superposition of (xz ↑, xz ↓) singlets which are formed between a central site and sites displaced by an odd number of lattice spacings in the y-direction, (yz ↑, yz ↓) singlets between the central site and odd numbered sites in the xdirection, and weaker (xy ↑, xy ↓) singlets with a more intricate s-wave arrangement. The internal structure of the d x 2 −y 2 state consists of a similar superposition but with a negative phase difference between the (xz ↑, xz ↓) and (yz ↑, yz ↓) singlets. In addition, there is a significant d x 2 −y 2 contribution from the (xy ↑, xy ↓) singlets confined primarily to the nearest neighbor sites.

VIII. SUMMARY AND CONCLUSIONS
We have studied a weak-coupling spin and charge fluctuation model of the pairing interaction using a tight-binding parametrization of the electronic band structure obtained from DFT calculations for the Fe-pnictides. We initially reviewed criticism of the 2-orbital models used in early studies of this kind, and compared results with more accurate 5-orbital models. Within the 5-orbital framework, we have calculated the multiorbital susceptibility. We calculated the static homogeneous spin susceptibility within this framework, and showed that it qualitatively agrees with experimental measurements of the T dependence, but that an interband susceptibility component was necessary to understand the magnitude of the measured low-temperature limiting susceptibility. The Tdependence of both χ(T ) and (T 1 T ) −1 is qualitatively similar to experiments, and may be consistent with them for some values of the parameters. A detailed fit was not attempted here.
We then constructed the pairing vertex from the generalized RPA susceptibilities, and calculated the pairing eigenvalues for dopings corresponding to undoped and electron-doped materials, for a variety of interaction parameters corresponding to intra-and inter-orbital Coulomb matrix elements and Hund's rule couplings. We found that in the parameter region where there is a significant coupling strength, the pairing interaction for our model of the Fe-pnictides arises from the exchange of spin fluctuations. These give rise to both intra-and inter-Fermi surface scattering processes. In a number of cases the dominant pairing contribution came from particle-particle scattering processes from the hole Fermi surface around the Γ point to the electron Fermi surfaces around the (π, 0) and (0, π) points in the unfolded Brillouin zone. These scatterings involve momentum transfer of order (π, 0). We find that there are two leading pairing channels, one with "sign-changing" swave (A 1g ) symmetry and one with d x 2 −y 2 (B 1g ) symmetry, and that the gap functions corresponding to both have nodes on the Fermi surface. For values of the interaction parameters corresponding to significant pairing strength, these two eigenvalues can become very close. The s-wave gap that exhibits a sign change between its average values on the electron and hole Fermi surfaces nevertheless also displays nodes on the electron sheet, whereas the d-wave state has its nodes on the hole sheet. The s-wave gap nodes are not required by symmetry, and could be absent for a different choice of parameters, but appear to be robust within the manifold of apparently re- alistic interaction parameters we have studied.
An obvious question which arises in such studies is the sensitivity of results to the particular choice of band structure. We have relied upon the Fermi surface for the paramagnetic system determined in the Generalized Gradient Approximation (GGA) calculations of Cao et al. 4 , but other band calculations find subtle differences in Fermi surfaces for the La-1111 material. We have attempted to address this question by considering variations of our effective 5-band model hopping parameters chosen originally to fit the Cao et al. band structure. In particular, we considered variations which led to a perfect nesting of the outer α and β Fermi surface sheets, to try to maximize the (π, 0) contribution to the pairing vertex. We found that our results for the structure of the susceptibility and pairing vertices changed in fact very little. On the other hand, the subtle changes led to a nearly exact degeneracy of the extended-s and d x 2 −y 2 eigenvalues of the linearized pairing problem. We argued that this degeneracy becomes exact for a situation where all sheets have the same radius, and propose this as a simple explanation of our finding that these two pairing channels appear to compete very strongly with one another, even in the realistic cases.
Our results appear to be qualitatively different from several works 15,18,42,46,47 which find, on heuristic grounds or within microscopic calculations, that the A 1g state is closer to the form cos k x × cos k y (in the unfolded zone), in contrast to our result, which is closer to the form cos k x + cos k y in the sense described above. Although both states have identical A 1g symmetry, the difference between them is important given the structure of the Fermi surface; for example, the cos k x ×cos k y form does not have nodes on either of the β sheets shown in Fig. 4. We believe that the strong pair scattering between the β 1 and β 2 sheets found in our calculations, and neglected in some simpler models, is crucial to stabilize the gap structure similar to cos kx + cos ky.
We further studied the spatial structure of the states corresponding to the leading pairing eigenvalues. In each case, the internal structure consists of a superposition of singlets made up from electrons which occupy the same d xz , d yz , or d xy orbitals on different sites. The distribution and relative signs of these superpositions are illustrated in Figs. 22 and 23. The fact that the s-wave and d-wave solutions have very similar coupling constants opens the possibility that different members of the Fe-pnictide superconducting family may have different gap symmetries. Furthermore, it suggests that these systems may have low-lying collective particle-particle modes which would be s-wave like in a d-wave superconductor, and d-wave like in an s-wave superconductor 48 .
There have been several earlier weak-coupling calculations of pairing in the Fe-pnictide superconductors. An early suggestion of an extended s-wave sign reversed gap made by Mazin et al. 42 was based upon the Fermi surface structure found in a DFT calculation. There it was proposed that (π, 0) antiferromagnetic spin fluctuations could lead to a superconducting state with an isotropic s-wave gap which had opposite signs on the electron and hole pockets.
Using a 5-orbital parameterization of the DFT bandstructure, Kuroki et al. 11 have carried out a calculation in which an RPA form for the spin and orbital contributions to the interaction was used along with bare single particle Green's functions to construct a linearized gap equation. In a parameter range similar to ours, they found that the leading pairing instability had s-wave symmetry with nodes on the electron-like Fermi surface. They also noted that the next leading channel had d x 2 −y 2 symmetry and discussed conditions under which this could become the leading pairing channel. The behavior of their s-wave gap differs from what we find by a phase factor of -1 on the electron Fermi surfaces. That is, if the gap on the inner hole Fermi surface around the Γ point is taken as positive, then Kuroki et al. find that the gap on the electron Fermi surface at (π, 0) reaches its largest negative value at the point closest to the origin Γ. With the same convention, we find that the gap on the electron Fermi surface takes on its largest positive value at this point. This difference may reflect a difference in the nesting of the hole and electron Fermi surfaces due to fits to slightly different band structures or from the choice of interaction parameters 49 .
Wang et al. 19 have studied the pairing problem using a 5-orbital effective band structure together with the functional renormalization group approach. Using the same tightbinding parameterization as Kuroki et al., they also find that the leading pairing channel has s-wave symmetry and that the next leading channel has d x 2 −y 2 symmetry. For their interaction parameters, they find that there are no nodes on the Fermi surface, but there is a significant variation in the magnitude of the gap. With the sign convention where the gap on the inner hole Fermi surface is positive, their gap function reaches its smallest negative value at the point on the electron Fermi surface which is closest to Γ. They note that for other parameter choices, this s-wave gap function on the electron Fermi surface may have nodes. This would lead to a gap which becomes positive on a region of the electron Fermi surface closest to the Γ point, in agreement with what we have found.
Both of these studies also conclude that the pairing interaction arises from spin fluctuation scattering with the dominant contribution associated with Q ∼ (π, 0) scattering of a pair from the inner Fermi hole surface around the Γ point to the electron Fermi surface around (π, 0) and (0, π). They also find that the magnitude of the gap on the central hole Fermi surface is smaller than the maximum magnitudes of the gap on the inner hole and the electron Fermi surfaces. Thus the five-orbital weak coupling calculations all find singlet pairing in the extended s-wave channel with a nearby d x 2 −y 2 wave state. The question of whether there are gap nodes for the swave, and indeed whether or not the s-wave state is ultimately stable with respect to the d-wave state, would both appear to depend on parameters. This sensitivity appears to us to be a natural consequence of the importance of several orbitals near the Fermi surface. In contrast to the cuprates, where calculations of this kind show a clear d x 2 −y 2 state well separated from other pairing states, it seems possible here that variations in band structure or interaction parameters found in different materials might possibly lead to different symmetry ground states.
We close by pointing out that at this writing, experiments have not conclusively answered either the question of order parameter symmetry or even whether gap nodes are present. Early indications from specific heat 21 , Andreev point contact spectroscopy 28 , and NMR 50 on the 1111 materials suggested nodes in the superconducting order parameter. On the other hand, some penetration depth experiments 31 and ARPES experiments 29,32 appear to find nearly isotropic gaps in the 122 materials. Our calculations find nearly degenerate s and d x 2 −y 2 states at T c , both of which have nodes on the Fermi surface. If the experiments indicating a lack of nodes are correct, there are two ways to imagine reconciling our conclusions with experiments on the 122 materials. The first is to consider the effects of disorder present at fairly high levels in current crystals, which should average the order parameters to a finite quasiisotropic value in the s case. The second possibility follows from the observation that if the d x 2 −y 2 state is the leading eigenvalue at T c , the thermodynamic ground state may be unstable towards an admixture of the s state 51 . A similar phenomenon occurs in the phase diagram of ordinary s and d symmetries in one-band superconductors with these pairing channels competing 52 . Separating these possibilities may be difficult, and will probably require phase sensitive probes as were discussed in the cuprate context 53 . ξ 11/22 = 2t 11 x/y cos k x + 2t 11 y/x cos k y + 4t 11 xy cos k x cos k y ± 2t 11 xx (cos 2k x − cos 2k y ) + 4t 11 xxy/xyy cos 2k x cos k y +4t 11 xyy/xxy cos 2k y cos k x + 4t 11 xxyy cos(2k x ) cos(2k y ) ξ 33 = 2t 33 x (cos k x + cos k y ) + 4t 33 xy cos k x cos k y + 2t 33 xx (cos 2k x + cos 2k y ) ξ 44 = 2t 44 x (cos k x + cos k y ) + 4t 44 xy cos k x cos k y + 2t 44 xx (cos 2k x + cos 2k y ) +4t 44 xxy (cos 2k x cos k y + cos 2k y cos k x ) + 4t 44 xxyy cos 2k x cos 2k y ξ 55 = 2t 55 x (cos k x + cos k y ) + 2t 55 xx (cos 2k x + cos 2k y ) +4t 55 xxy (cos 2k x cos k y + cos 2k y cos k x ) + 4t 55 xxyy cos 2k x cos 2k y ξ 12 = −4t 12 xy sin k x sin k y − 4t 12 xxy (sin 2k x sin k y + sin 2k y sin k x ) − 4t 12 xxyy sin 2k x sin 2k y ξ 13/23 = ±2it 13 x sin k y/x ± 4it 13 xy sin k y/x cos k x/y ∓ 4it 13 xxy (sin 2k y/x cos k x/y − cos 2k x/y sin k y/x ) ξ 14/24 = 2it 14 x sin k x/y + 4it 14 xy cos k y/x sin k x/y + 4it 14 xxy sin 2k x/y cos k y/x ξ 15/25 = 2it 15 x sin k y/x − 4it 15 xy sin k y/x cos k x/y − 4it 15 xxyy sin 2k y/x cos 2k x/y ξ 34 = 4t 34 xxy (sin 2k y sin k x − sin 2k x sin k y ) ξ 35 = 2t 35 x (cos k x − cos k y ) + 4t 35 xxy (cos 2k x cos k y − cos 2k y cos k x ) ξ 45 = 4t 45 xy sin k x sin k y + 4t 45 xxyy sin 2k x sin 2k y