Coupled multiple-mode theory for s± pairing mechanism in iron based superconductors

We investigate the interplay between the magnetic and the superconducting degrees of freedom in unconventional multi-band superconductors such as iron pnictides. For this purpose a dynamical mode-mode coupling theory is developed based on the coupled Bethe-Salpeter equations. In order to investigate the region of the phase diagram not too far from the tetracritical point where the magnetic spin density wave, (SDW) and superconducting (SC) transition temperatures coincide, we also construct a Ginzburg-Landau functional including both SC and SDW fluctuations in a critical region above the transition temperatures. The fluctuation corrections tend to suppress the magnetic transition, but in the superconducting channel the intraband and interband contribution of the fluctuations nearly compensate each other.

interband scattering due to non-magnetic impurities is destructive for s ± superconductivity, which is considered to be realized in most of FeSC [18][19][20][21][22] and magnetic interband scattering on the contrary stabilizes the s ± coupling mechanism 18,23 . The kinematics of this stabilization mechanism is similar to that in so called π-junctions in heterostructures superconductor/insulator/superconductor, where the Josephson tunnelling through a barrier with paramagnetic impurities is accompanied with the sign change of the superconducting order parameter Δ 24 . Another remarkable response to doping occurs when nominally nonmagnetic defects such as As vacancies [25][26][27] or Ru substitutions 28 are introduced in 1111 materials. These impurities trigger a strong paramagnetic response (see ref. 29 for explanation of this effect), and their observed modification of T c does not fit with available theories.
To this end we introduce and develop a multi-mode theory that takes into account the fact that the superconductivity in these materials arises as a result of competition of at least two coupling mechanisms: superconducting pairing with the scalar order parameter Δ and the spin density wave ordering with vector order parameter m (the magnetic moment). These modes are coupled by the interband electron-electron interaction, and we consider the effect of impurity scattering on this coupling mechanism.
The paper is written as follows. We start with the derivation of a system of Bethe-Salpeter (BS) equations containing coupled SDW and Cooper channels for temperatures above the magnetic and superconducting transition temperatures in the next section. Based on BS equations we show that dynamical fluctuations play an important role close to the tetracritical point and must be taken into account for the correct description of the tetracritical region.
Furthermore we derive the Landau-Ginzburg (LG) functional and calculate fluctuation corrections to the phase transition temperatures to superconductor and SDW phases (T c and T s , respectively). Detailed derivation of BS and LG equations is given in the Section "Methods".
Where needed, we will be specific and concentrate on the multi-valley semi-metallic FeSCs from the 1111 (ReFeAsO) and 122 (AeFe 2 As 2 ) families under electron and hole doping (Re and Ae stay for the rare-earth and alkaline-earth elements) respectively.

Multimode approach to SDW and BCS instabilities
In this Section we formulate a mode coupling approach for multiband ferropnictide superconductors with nearly nested hole and electron pockets of Fermi surface. Although a complete consensus is lacking on the type of superconductor pairing in these system, it appears that the experimental arguments in favor of s ± mechanism are rather solid. This type of ordering was proposed originally for superconductor-excitonic instability in semimetals 30 and has been reformulated for superconductor pairing in FeSCs in refs 31, 32 (for a review see ref. 5). It is important to note that mode-mode coupling is an inherent constituent of s ± pairing [33][34][35] : in the two-dimensional systems with nesting conditions between electron and hole pockets an interplay between the SDW and Cooper modes results in additional coupling within each channel induced by all other channels. This theory is based on a renormalization group (RG) approach, which implies logarithmic renormalization of the relevant vertices. It was shown that the Cooper and SDW decouple at the scale of Fermi energy and further flow goes independently in both channels. The natural limitation of such approach is the demand of perfect nesting in SDW channel. In real FeSC systems the nesting conditions are satisfied only approximately. It leads to the suppression of the SDW channel at some energy.
The general approach to deal with mode-mode coupling that we introduce here is based on the mode coupling theory of critical phenomena. We start from the "high temperature" region, where the critical fluctuations are already well developed, while long-range SC and SDW order are not yet established. The corresponding vertex parts Γ sc (q, Ω, T) and Γ sdw (Q + q, Ω, T), are represented at this temperature by the fluctuating modes D sc (q, Ω) and D sdw (q, Ω) respectively 36,37 : c c c s s s and γ sc/sdw = 8T/π. Here we consider the limit of small momenta c Δ/m q 2 ≪ 1, ν is a density of states at the Fermi level, ξ = ∆ c s 2 and ξ = c m m 2 are superconducting and magnetic coherence lengths respectively, see Section "Ginzburg-Landau approach".
As the pole iΩ sc/sdw = γ sc/sdw (τ c/s + c Δ/m q 2 ) of the vertex part Γ sc or Γ sdw at q → 0 tends to zero, the corresponding instability results in the phase transition. When approaching the transition temperature T c or T s , the corresponding vertex part diverges in the static limit: The mode coupling approach is valid in the vicinity of the tetracritical point on the phase diagram shown in Fig. 1 where the two temperatures T c and T s are comparable. Then the system of equations for the vertex parts at small Ω and q can be derived and solved. These are the Bethe-Salpeter equations for the most relevant vertices. As Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 will be shown below, this system contains at least four bare vertices which should be taken into account. The divergences (3) provide us with a criterion of selecting the corresponding diagrams. Namely, the terms containing polarization operators in the electron-hole channel with transmitted momentum close to Q, and the electron-electron polarization operators of Cooper type should be collected in all orders. As a result, the system of Bethe -Salpeter equations acquires a parquet-like structure, and the solution of the mode-coupling equations obtained at T > T c may be mapped on the mean-field solutions of renomalization group (RG) equations below T c 33,34 . We assume that the superconducting transition temperature T c is imposed on the system by a dominant s ± pairing. In this Ansatz the interband electron-hole Coulomb interaction is the main source of Cooper pairing, so that a weak BCS attractive interaction in the electron and hole pockets would result in SC instabilities with Bethe-Salpeter equations. Following the Ansatz formulated above (see details in the Section "Methods"), we derive the system of BS equations for the relevant vertex parts Γ i (see Fig. 2). The starting Hamiltonian is: where f σk , c σk stand for annihilation operator in electron and hole bands; ε e h k / and μ are dispersions and chemical potential, respectively.
Based on the results of RG calculations 33,34 , we choose four vertices relevant to the s ± coupling (the vertices u 1 , u 3 , u 4 = u s , u 5 = u s ). The vertex u 2 which has been shown to be irrelevant 34 is excluded.
The vertices Γ 1i with i = 1, 2 (see Fig. 2) describe the interactions in the density-wave block, the vertices Γ 4,5 describe the singlet Cooper pairing in the electron and hole pockets, respectively, and the vertex part Γ 3 includes the interactions responsible for the interband (e-h) Cooper pairing. The system of Bethe--Salpeter equations may be written in the symmetric form using a matrix Λ : The matrix − Λ (1 ) is the secular matrix and  is − Λ = det(1 ) 0. In these notations both coupling constants u i = νU i and the polarization loops Π i (q = 0, ω = 0) = ln(W/T) are dimensionless (here W is the bandwidth). We neglect below the difference between the band dispersion in the electron and hole pockets, when calculating the polarization operators Π s = Π se = Π sh . One should note that his approximation does not imply that we rely on perfect nesting. In electron doped materials the strict nesting conditions are of course violated, but the small difference between the electron and hole polarization operators is insufficient for our theory because two contributions are summed in the critical mode D 32 which is central to our scenario. The criticality in this mode is dominated by the term u 1 Π 1 . We use the sign convention that both Π 1 and Π s are positive. Under this convention u s > 0 is attractive in the Cooper channel and u 1 > 0 facilitates the instability in the SDW channel.
In the explicit form the secular equation reads: One can see that in this approximation the  ±1 triplet Γ 31 -channel described by the second row of the secular equation decouples from the rest, which corresponds to the 0  -SDW and superconducting channels. The D 32 -mode plays a very particular role. On the one hand it contributes to the density wave channels (both CDW and SDW). On the other hand it strongly affects the superconducting channel. In the approximation used in the present paper the SDW and CDW are degenerate. The interband interaction u 2 > 0 with momentum transfer Q lifts out the degeneracy favouring the SDW transition.
We will consider the part of the phase diagram concentration-temperature (c − T) close to the point of the degeneracy of the s ± and SDW channels (i.e. the tetracritical point) shown in Fig. 1. In this region the SC instability takes place in the presence of critical SDW fluctuations. Above T s the fluctuation modes arise at momenta p = Q + q close to the nesting vector Q connecting the Γ and X points in the Brillouin zone and at small ω → 0.
We assume that the two Cooper propagators D 4 and D 5 are far from any divergence, namely Π  u 1 s s at these temperatures (in case of a dominant interband pairing mechanism as in the pronounced s ± case adopted here, the purely intraband Cooper instabilities develop at temperatures much less than the actual T c ). In this approximation the vertices are real, and the channels 4, 5 are represented by a single row and column.
Away from the tetracritical point. We consider first the case when the temperature is less than the Fermi-energy T < ε F and the doping c is away from the tetracritical point. Then the divergence of the vertices is strongly peaked at particular momenta. For instance, putting u 3 → 0 in the integral equation for the Γ 11 (p 1 , p 2 , p 3 , p 4 ) (see Fig. 2a), one immediately gets that this vertex is divergent only for p 1 + p 2 = Q. A similar analysis shows that Γ 32 diverges for p 1 + p 2 → 0 and p 1 − p 3 → Q. The splitting of momenta, at which the vertices diverge, decouples the matrix Eq. (5) into the density-wave and Cooper channels. To study the properties of the vertex functions in the vicinity of their singularities we introduce Γ ≈ Γ have poles at k, p → 0 correspondingly. For the vertices Γ 32 and Γ 4 we use the s ame de comp osit ion: . To identify the corresponding transition temperatures we consider the Bethe-Salpeter equations in the static limit. For small momenta  p k q , , the system Eq. (5) splits up in the logarithmic approximation (see Section "Methods") into cooper channel: and density-wave channel With these vertices we now calculate the susceptibility in the cooper and density-wave channels. In the particle-particle s ± and s ++ channel drawn in Fig. 3 we obtain Substituting Eq. (9) into Eq. (11), we get: This equation has the typical pole structure for the superconducting susceptibility in the vicinity to the transition temperature. The susceptibility χ = increases with decreasing temperature and diverges at the critical transition temperature. This happens under condition 1 − (u s ± u 3 )Π s (k = 0, T = T c ) = 0. The two solutions correspond to s ++ and s ± superconducting order parameter. For CDW and SDW channels we get: The SDW transition is determined by 1 − (u 1 + u 3 )Π 1 (q = 0, T c ) = 0. One can define a tetracritical point, when both susceptibilities diverge. It happens under the condition (u s + u 3 )Π s (k = 0, T c ) = (u 1 + u 3 )Π 1 (k, T c ) → 1. In this case and in the vicinity of the tetracritical point one should find the divergence of the susceptibility taking into account the full matrix Eq. (5).
Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 Close to the tetracritical point: dynamical multi-mode coupling theory. Now we consider the interplay between magnetic and superconducting degrees of freedom in the vicinity of the tetracritical point. In this case the approach based on the effective hydrodynamic action is much easier than the direct solution of the BS equations. Besides, this approach which also accounts for the multi-point correlation functions goes beyond the conventional BS paradigm limited by two-and four-point Green's function 16,17 . The derivation of the effective hydrodynamic action in terms of magnetic → → m x t ( , ) and superconducting ∆ → x t ( , ) i (i = e, h) dynamical fluctuations is done by integrating out the BS equations with respect to the "fast" (with the energy of the order of the bandwidth) degrees of freedom and is presented in Section "Methods".
Within our analysis we start discussing the limiting case  u u s 3 and Πũ 1 s 3 corresponding to the interplay between two fluctuating modes: one is the s ± superconducting and another one is the SDW magnetic. This regime is believed to be present in most of FeSC [e.g., in broadly studies doped "A-122" systems (where A = Ba, Sr, Ca)] The Lagrangian of the two-mode coupling theory takes the form: Here the Fourier transforms of superconducting and SDW fluctuators computed on Matsubara frequencies iΩ m read:  Here ψ(x) and ψ′ ′ (x) are the di-gamma function and its second derivative respectively, 〈 ...〉 FS denotes the averaging over the Fermi surface (here a parabolic dispersion and equal mass were assumed for electron and hole pockets). The analytic continuation of the superconducting and SDW fluctuators to the upper half-plane iΩ m → Ω + i0 + for Ω  T T ( , ) c s is given by Notice that we normalize the fluctuators by the density of states to make them dimensionless. The equations defining the coefficients A, B, C 1 , C 2 and c Δ/m are derived in Section "Ginzburg-Landau approach". The static q -independent part of the Lagrangian corresponds to the Landau expansion of the free energy. Inclusion of the gradient terms generalizes the Landau theory to the Ginzburg-Landau functional (see Section "Ginzburg-Landau Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 approach"). The effective Lagrangian describes effective non-linear theory of interacting mode and therefore goes beyond linear Bethe -Salpeter approach. Hidden non-linearity of Bethe-Salpeter equations is associated with curved ("non-flat") phase space: Γ 32 enters both SDW and superconducting equations being singular both at small total moment/energy and small deviating from Q momentum transfer. The Lagrangian (14) is U(1) × SU (2) symmetric.
Although the situation in 1111 and 122 pnictides is described in the framework above, we consider for the sake of generality also the opposite limit u 3 ≪ u s and u s Π s ~ 1. In this case one deals with the three-mode coupling theory: two superconducting fluctuating modes Δ e and Δ h and one SDW → m magnetic mode:  As in the static case, we neglect the difference between the contribution of electron and hole pockets in the fluctuation modes. In principle, the SDW and CDW modes are degenerate when Π  u 1 3 1 and one needs to consider a four-mode coupling theory. We however, assume that the SDW-CDW degeneracy is lifted out by additional inter-band processes u 2 (see Eq. (4) and further details in refs 33, 34), omitted through the derivation of the BS equations, T CDW < T s and restrict ourselves by the three-mode coupling theory.
The Lagrangian (17)  Fluctuation corrections to T c and T s . In order to find the fluctuation correction to the superconducting transition temperature given by the mean-field analysis in the vicinity of the tetracritical point, we integrate out the remaining slow magnetic fluctuations in (14) and finally obtain the effective Lagrangian describing the superconducting system:  Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 where T c 0 is the transition temperature defined by the static solution of BS equations (mean-field theory) and  are bare Green's Functions for electron/hole bands respectively. We observe that there are two types of competing processes: the one with + C 1 leads to K 1 > 0 and results in suppression of T c while the second one with − C 2 corresponds to K 2 < 0 and therefore enhances T c .
Similarly, for the second case described by the Lagrangian (17) characterized by magnetic fluctuation broken U(1) × U(1) we get the following effective Lagrangian  Fig. 5(a,b) while the off-diagonal Δ e − Δ h coupling K 2 is defined by the diagrams Fig. 5(c,d). The corresponding equations for the d-dimensional cases (d = 2, 3) are given by:   While the diagrams Fig. 5(a,b) always reduce the effective temperature of the superconducting transition, the diagrams Fig. 5(c,d) lift the degeneracy between e-h transition temperatures. As a result, assuming that K 1e = K 1h = K 1 , K 2eh = K 2 we get the fluctuation correction to the critical temperature T c : Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 To construct the effective field theory for the influence of superconducting fluctuations on the SDW dynamics, we integrate out the slow superconducting fluctuations in (17) and finally obtain the effective Lagrangian for paramagnetic SDW fluctuations: where for two-mode coupling theory     The diagrams defining the fluctuation corrections to T s are shown on Fig. 5(e,f). Computing the correcting terms for the three-mode coupling theory we obtain that the fluctuation correction to the SDW transition temperature always leads to its reduction: . Due to the assumed particle-hole symmetry we get = ∼ K K 1 1 . Finally, combining (28) and (35) we get and therefore the new tetracritical point corresponds to lower values of the concentration c < c 0 (see Fig. 1). Let us illustrate, as an example, the calculation of the fluctuation corrections diagrams K 1,2 containing the superconducting fluctuator L e/h (q, Ω) (the evaluation of the diagrams containing the SDW fluctuator can be performed in a similar fashion). The evaluation of the Matsubara sums and integrals can be done in several steps: (i) First, the main contribution to the sum over Ω m is given by the m = 0 term. Therefore we can replace ii Second, we notice that the integral over q is determined by the small q: 2 for d-spatial dimensions d = 2, 3. We can therefore neglect the q-dependence in the e/h Green's functions: e h n eh n / 0 / 0 (iii) Next, the sum over fermionic Matsubara frequency ε n and summation over k in K 1 and K 2 is performed in the same way as explained below in the Section "Ginzburg-Landau approach" (see Section "Methods"). As a result we obtain: Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 Notice, that after these approximations there is no difference between K 1 and K 2 .
(iv) The remaining integral over momentum q is performed using the static Ornstein-Zernike correlator L e/h (q, 0) (see ref. 36) Performing similar calculations for the renormalization of the SDW critical temperature T s we obtain: We emphasize that the constants K 1 and K 2 characterizing the fluctuational shift of the critical temperature are found to depend on the Ginzburg number Gi (d) and the dimensionality of the system d only.
Eq. (42) represent the central result of the paper: fluctuation corrections are responsible for the competition between the spin density wave and superconducting critical modes in the vicinity of the tetracritical point. We have obtained this results by calculating the fluctuation corrections to the GL functional, but the same result may be obtained in the framework of renormalization group approach. Since the U(1) × U(1) symmetry is a priori broken in the superconducting sector of the model, the influence of the magnetic fluctuations onto the superconducting transition is two-fold: first, the intra-band corrections tempt to reduce the transition temperature; secondthe inter-band corrections do exactly the opposite -split the two transition temperatures and effectively facilitate the superconducting transition. On the other hand, the contribution of superconducting fluctuations in the SU (2) magnetic sector of the model does result only in a suppression of the SDW transition. The contribution of "eigen" fluctuations, namely the superconducting-superconducting and magnetic-magnetic is alike and therefore it is dropped off from the difference of the critical temperatures. The asymmetric character of the fluctuation corrections results in the shift of tetracritical point towards lower values of the carriers concentration. For the "dirty" limit of the multi-mode coupling theory the Ginzburg numbers should be updated accordingly 36 .
The strong renormalizations of T c and T s are expected in iron based superconductors, which are rather two-dimensional and many of them have rather small Fermi energy. In these materials the value of the Ginzburg number may be quite large. For instance, the characteristic value of the Fermi energy in 1111 and 11 compounds is 100 meV [38][39][40] , while the superconducting critical temperature is about 20-50 K. With this numbers the one finds large Ginzburg number Gi (2) ~ T c /ε F ~ 0.02-0.05 and considerable renormalizations of the critical temperatures (T c − T s )/T c ~ 0.1.

Discussion
We have developed a high-temperature approach to the problem of interplay between magnetic and superconducting ordering in multi-band systems. Both static and dynamical (fluctuation related) contribution to the mode-mode coupling are discussed. It is shown that the fluctuation corrections in the vicinity of the tetracritical point reduce the magnetic transition in accordance with Eqs (32), (36). On the superconductor side of the phase diagram the situation is different. The intraband and interband contributions of spin fluctuations to the Ginzburg-Landau functional nearly compensate each other.
Scientific RepoRts | 6:37508 | DOI: 10.1038/srep37508 The approach that we have introduced is more general than for instance a mean-field description of competing phases 4,5 as it accounts for dynamical fluctuations and goes beyond the static scaling paradigm. Apart from this, the framework presented here has the advantage that it can be generalized to other multi-mode regimes in a straightforward manner, which can include e.g. the presence of a nematic instability, the competition between s ± and s ++ and/ or singlet/triplet pairings. Besides, it is highly appealing to discuss in a framework of this approach the influence of conventional magnetic defects (transition metal ions substituting for Fe) and weak magnetic defects produced by Zn impurities and As-vacancies on the superconducting transition temperature. This work is now in progress.

Methods
Integral Bethe-Salpeter equations. The integral Bethe-Salpeter equations in d = 2, 3 dimensions take the form   where we used the short-hand notations n , = → Ω q q i ( , ) m with fermionic Matsubara frequencies ω n = ε n = 2πT(n + 1/2) and bosonic Matsubara frequency Ω m = 2πmT and the irreducible vertices ℑ j0 include all irreducible diagrams in the channel j. The spin, momentum and energy are conserved in each vertex: p 4 = p 1 + p 2 − p 3 . The simplified matrix BS equations are obtained by a replacement , , s s 10 1 3 0 3 0 Due to momentum conservation there are only three independent momenta. It worth to introduce new variables p = p 1 + p 2 = p 3 + p 4 , q = p 3 − p 2 and t = p 3 − p 1 . In the new variables Γ =Γ  p q t p p p p ( , , ) ( , , , 4 the first integral in Eq. (43) has the form:   The eigen modes of superconducting fluctuators are given by: ii we insert vector bosonic fields describing magnetic fluctuations by means of δ-function iii implement the integral representation for the functional δ-function The next steps are similar to derivation of superconducting fluctuating part: integrating over Grassmann variables and integrating over → n by means of the saddle point approximation. As a result, the magnetic fluctuators are described by the Lagrangian: Performing a loop expansion (see Fig. 4) we get where C 1 and C 2 are defined in the next section.
Ginzburg-Landau approach. In order to establish a correspondence between the microscopic fluctuation theory of competing modes and the macroscopic thermodynamic description, we can derive an effective Ginzburg-Landau (GL) functional. The GL functional for = − d ord n orm    should be written in terms of order parameters corresponding to competing modes (see, e.g. refs [41][42][43][44][45].
We present the GL functional for two important cases discussed in the paper. The coefficients α Δ = ln(T/T c ) and α m = ln(T/T s ) change sign at the corresponding transition temperatures (in the case of the transition between the ordered states SDW or SC and the coexistent states, the closed loops are constructed from corresponding anomalous Green functions). The coefficients ξ = ∆ c s 2 and ξ = c m m 2 in front of the gradient terms are obtained in a standard way from the small momentum expansion of the polarization loops Π s and Π 1 , ξ s and ξ m are superconducting and magnetic coherence lengths respectively (in the "dirty" limit the coherence length ξ s/m should be replaced by ξ ⋅ l s m / where l is the mean-free path. This replacement changes the Ginzburg number accordingly), see details in ref. 16, 17, 36. The fourth order terms are given by the loops containing four Green's functions (see Fig. 4