Mixed-state effect on quasiparticle interference in iron-based superconductors

The discovery of iron-based superconductors [1] has triggered efforts to elucidate the superconducting (SC) pairing mechanism in these materials. One hotly debated issue is the symmetry and structure of the SC gap. Theoretically, it was initially suggested that the pairing may be established via spin fluctuations, leading to the so-called s_± pairing symmetry ($\Delta _{\vec{t} {k}} \sim \cos k_{x}+\cos k_{y}$ defined in the 2Fe/cell Brillouin zone (BZ)). In this case the SC order parameter (OP) exhibits a sign reversal between the hole pockets (around the Γ point) and electron pockets (around the M point) [2–4]. Later, s-wave pairing symmetry without sign reversal was also proposed to be a possible candidate which is induced by orbital fluctuations due to the electron-phonon interaction [5]. Experimentally, the results about the pairing symmetry remain highly controversial as well. For example, in Ba_0.6K_0.4Fe₂As₂, an optimally hole-doped iron-based superconductor [6], the SC gaps measured by the angle-resolved photoemission spectroscopy (ARPES) can be approximately fitted by $|\Delta _{\vec{t} {k}}| \sim |\cos k_{x}+\cos k_{y}|$, with almost isotropic gaps on all the Fermi surfaces (FS) [7, 8], indicating the possible pairing symmetry to be either s_± or s-wave. The situation is similar for optimally doped FeTe_1−xSe_x with x ≈ 0.45 [9]. Neutron scattering experiments observed a resonance peak at $\vec{t} {Q}=(\pi ,\pi )$ below the SC transition temperature [10–13] as predicted by some theoretical works assuming s_± symmetry [14, 15], thus at first glance they seemed to support each other. However, later theoretical works suggested that the experimentally observed resonance peak can also be reproduced by assuming s-wave pairing [16, 17].

Recently, in order to clarify the pairing symmetry, Hanaguri et al. used scanning tunneling microscopy (STM) to image the quasiparticle interference (QPI) patterns in the SC state [18]. They proposed that the relative sign of the SC OP can be determined from the magnetic-field dependence of quasiparticle scattering amplitudes and claimed that their experimental data were only consistent with the s_± scenario but not the s-wave one. Soon after, the experimental results were put into question and argued instead to arise from the Bragg scattering but not due to the QPI because the observed peaks are too sharp [19]. On the other hand, the magnetic field will induce vortices into the system and lead to the inhomogeneity of the pairing OP in real space, thus affecting the QPI patterns. Theoretical analyzes performed previously have investigated the mixed-state effect on the QPI [20–22]. However, in ref. [20], only amplitude suppression of the OP near the vortex core was considered, without taking into account the phase variation. In ref. [21], the mixed-state effect was treated by using quasiclassical approximation and the QPI derived in this way was rather broad compared to the experimental observation. Furthermore, in ref. [22], only the effect of the Zeeman splitting was considered, which should be negligible with respect to the mixed-state effect. Thus there still lacks direct theoretical confirmation of the mixed-state effect on the QPI.

Therefore in this work we adopt a phenomenological model with s_± pairing symmetry to study the influence of vortices on the QPI by directly solving large-scale Bogoliubov-de Gennes (BdG) equations in real space based on the Chebyshev polynomial expansion. For comparison, the problem is also studied for s-wave pairing. In this way the mixed-state effect on the QPI can be rigorously investigated and the results unambiguously support s_± pairing symmetry in iron-based superconductors.

We begin with an effective two-orbital model on a two-dimensional lattice [23], with a phenomenological form for the intraorbital pairing terms. The Hamiltonian can be written as

$$\begin{eqnarray} H=&-&\sum_{ij,\alpha\beta,\sigma}[t^\prime_{ij,\alpha\beta}+(\mu-V_{imp}\delta_{ii_{imp}})\delta_{ij}\delta_{\alpha\beta}] c^{\dag}_{i\alpha\sigma}c_{j\beta\sigma}\nonumber &&+\sum_{ij,\alpha\beta}(\Delta_{ij,\alpha\beta}c^{\dag}_{i\alpha\uparrow}c^{\dag}_{j\beta\downarrow}+{\rm H.c.}). \end{eqnarray} \tag{ 1 }$$

Here i and j are the site indices while α,β = 1,2 are the orbital ones. σ represents the spin and μ is the chemical potential. Then we consider potential scattering by nonmagnetic impurities through V_imp with i_imp being the locations of the impurities. $\Delta _{ij,\alpha \beta }=\frac {V_{ij}\delta _{\alpha \beta }}{2}(\langle c_{j\beta \downarrow }c_{i\alpha \uparrow }\rangle -\langle c_{j\beta \uparrow }c_{i\alpha \downarrow }\rangle )$ is the intraorbital spin singlet bond OP, where V_ij is the onsite (i = j) or next-nearest-neighbor ($i=j\pm (\hat {x}\pm \hat {y})$) attraction which we choose to achieve the s-wave or s_± pairing symmetry, respectively. Here we need to emphasize that the purpose of the present work is to study the mixed-state effect on the QPI by assuming s_± or s-wave pairing symmetry, but not the mechanism of superconductivity in this type of compound. Thus we choose such an attractive interaction V_ij only to reproduce the desired s_± or s-wave pairing symmetry. Interestingly, based on this phenomenological model, the phase diagram [24] and vortex states [25] we derived are both in agreement with experimental observations, thus justifying the validity of the present approach. In the presence of a magnetic field B perpendicular to the plane, the hopping integral can be expressed as t'_ij,αβ = t_ij,αβexp $[i\frac {\pi }{\Phi _{0}}\int _{j}^{i}\vec{t} {A}(\vec{t} {r})\cdot \mathrm {d}\vec{t} {r}]$, where Φ₀ = hc/2e is the SC flux quantum, and $\vec{t} {A}=(-By,0,0)$ is the vector potential in the Landau gauge. The hopping integrals are

$$\begin{equation} t_{ij,\alpha\beta}=\left\{\begin{array}{@{}l@{\quad}l@{}} t_{1}, &{\alpha=\beta,i=j\pm\hat{x}(\hat{y})}, \\ \frac{1+(-1)^{j}}{2}t_{2}+\frac{1-(-1)^{j}}{2}t_{3}, & {\alpha=\beta,i=j\pm(\hat{x}+\hat{y})}, \\ \frac{1+(-1)^{j}}{2}t_{3}+\frac{1-(-1)^{j}}{2}t_{2}, & {\alpha=\beta,i=j\pm(\hat{x}-\hat{y})}, \\ t_{4}, & {\alpha\neq\beta,i=j\pm(\hat{x}\pm\hat{y})}, \\ 0, & \mathrm{otherwise}. \end{array}\right. \end{equation} \tag{ 2 }$$

The Chebyshev polynomials can be written as $\phi _{k}(x)=\cos [k\arccos x]$ and statisfy [26–29]

$$\begin{eqnarray} &&\phi_{0}(x)=1,\phi_{1}(x)=x,\phi_{k}(x)=2x\phi_{k-1}(x)-\phi_{k-2}(x),\nonumber\\ &&\sum_{k=0}^{\infty}\frac{W(x)}{\nu_{k}}\phi_{k}(x)\phi_{k}(x')=\delta(x-x'), \end{eqnarray} \tag{ 3 }$$

where $W(x)=1/\sqrt {1-x^{2}}$, ν_k = π(1 + δ_k0)/2 and x∈[ − 1,1]. Next we define the Green's function matrix:

$$\begin{equation} G(\tau)=-\langle T\tau C(\tau)C^{\dag}(0)\rangle, \end{equation} \tag{ 4 }$$

with C^† = (··· ,c^†_j1↑,c^†_j2↑,··· ,c_j1↓,c_j2↓,··· ). Equation (1) can be diagonalized as

$$\begin{equation} H=C^{\dag}MC=C^{\dag}QQ^{\dag}MQQ^{\dag}C=\Phi^{\dag}D\Phi. \end{equation} \tag{ 5 }$$

Here Q is a unitary matrix that satisfies (Q^†MQ)_rs = D_rs = δ_rsE_s and Φ = Q^†C. The spectral function can be expressed as [26, 27]

$$\begin{eqnarray} d_{rs}(\omega)&=&-\frac{1}{2\pi i}[G_{rs}(\omega+i\eta)-G_{rs}(\omega-i\eta)]\nonumber\\ &=&\sum_{\gamma}Q_{r\gamma}Q_{s\gamma}^{*}\delta(\omega-E_{\gamma})\nonumber\\ &=&\frac{1}{a}\sum_{\gamma}Q_{r\gamma}Q_{s\gamma}^{*}\delta(\tilde{\omega}-\xi_{\gamma})\nonumber\\ &=&\frac{1}{a}\sum_{\gamma}Q_{r\gamma}Q_{s\gamma}^{*}\sum_{k=0}^{\infty}\frac{W(\tilde{\omega})}{\nu_{k}}\phi_{k}(\tilde{\omega})\phi_{k}(\xi_{\gamma})\nonumber\\ &=&\frac{1}{a}\sum_{k=0}^{\infty}\frac{W(\tilde{\omega})}{\nu_{k}}\phi_{k}(\tilde{\omega})\phi_{k}(\tilde{M})_{rs},\nonumber\\ d(\omega)&=&\frac{1}{a}\sum_{k=0}^{\infty}\frac{W(\tilde{\omega})}{\nu_{k}}\phi_{k}(\tilde{\omega})\phi_{k}(\tilde{M}). \end{eqnarray} \tag{ 6 }$$

Here r, s, γ = 1,··· ,L with L = 4N_xN_y and N_x (N_y) being the number of lattice sites along $\hat {x}$ ($\hat {y}$) direction of the 2D lattice. a = (E^max_γ − E^min_γ)/(2 − ) ( > 0 is a small number), b = (E^max_γ + E^min_γ)/2, $\tilde {\omega }=(\omega -b)/a$, ξ_γ = (E_γ − b)/a and $\tilde {M}=(M-bI)/a$.

If we further define the L-dimensional vectors e(o) and h(o) as e(o)_γ = δ_γo and h(o)_γ = δ_γo+2Ns (N_s = N_xN_y and o = 1,··· ,2N_s), we can express the self-consistent parameters as

$$\begin{eqnarray} n_{j\beta\uparrow}&=&\langle c_{j\beta\uparrow}^{\dag}c_{j\beta\uparrow}\rangle=\int_{-\infty}^{\infty}\mathrm{d}\omega f(\omega)d_{mm}(\omega)\nonumber\\ &=&\int_{-\infty}^{\infty}\mathrm{d}\omega f(\omega)e(m)^{T}d(\omega)e(m)\nonumber\\ &=&\sum_{k=0}^{N-1}g_{k}e(m)^{T}e_{k}(m)\frac{T_{k}}{\nu_{k}},\nonumber\\ n_{j\beta\downarrow}&=&1-\sum_{k=0}^{N-1}g_{k}h(m)^{T}h_{k}(m)\frac{T_{k}}{\nu_{k}},\nonumber\\ \Delta_{ij,\alpha\beta}&=&\frac{V_{ij}\delta_{\alpha\beta}}{2}\sum_{k=0}^{N-1}g_{k}[e(n)^{T}h_{k}(m)\nonumber\\ &+&e(m)^{T}h_{k}(n)]\frac{T_{k}}{\nu_{k}}, \end{eqnarray} \tag{ 7 }$$

where m = 2(j_y + N_yj_x) + β and n = 2(i_y + N_yi_x) + α with i_x,j_x = 0,··· ,N_x − 1 and i_y,j_y = 0,··· ,N_y − 1. $e_{k}(m)=\phi _{k}(\tilde {M})e(m)$, $h_{k}(m)=\phi _{k}(\tilde {M})h(m)$ and $T_{k}= \int _{-1}^{1}\mathrm {d}\tilde {\omega }f(a\tilde {\omega }+b)W(\tilde {\omega })\phi _{k}(\tilde {\omega })$. N is the cutoff in the summation and g_k is the kernel we convolute to avoid Gibbs oscillations [28]. In addition, f(x) is the Fermi distribution function. At zero temperature we have

$$\begin{equation} T_{k}=\left\{\begin{array}{@{}l@{\qquad}l@{}} \pi-\arccos (-\frac{b}{a}), & {k=0}, \\ -\frac{\sin[k \arccos (-\frac{b}{a})]}{k}, & {k \neq 0}. \end{array}\right. \end{equation} \tag{ 8 }$$

Then we can solve the BdG equations self-consistently and the chemical potential is determined by the doping concentration. The calculation is repeated until the absolute error of the OP between two consecutive iteration steps and that of the total electron number are less than 10⁻⁴. The local density of states (LDOS) is given by

$$\begin{eqnarray} \rho_{j}(\omega)&=&\frac{1}{a}\sum_{\beta}\sum_{k=0}^{N-1}\frac{g_{k}}{\nu_{k}}[W(\tilde{\omega})\phi_{k}(\tilde{\omega})e(m)^{T}e_{k}(m)\nonumber\\ &+&W(\tilde{\omega'})\phi_{k}(\tilde{\omega'})h(m)^{T}h_{k}(m)], \end{eqnarray} \tag{ 9 }$$

where $\tilde {\omega '}=(-\omega -b)/a$. Following ref. [18], we further define

$$\begin{equation} Z_{j}(\omega)=\rho_{j}(\omega)/\rho_{j}(-\omega), \end{equation} \tag{ 10 }$$

and $Z_{\vec{t} {q}}(\omega )$ is the Fourier transform of Z_j(ω). Here $\vec{t} {q}=(q_{x},q_{y})$ is defined in the 1Fe/cell BZ in accordance with ref. [18].

The benefits of this method are threefold. First, it requires much less storage than the exact diagonalization method since the matrix M is sparse, thus we can solve large-scale BdG equations and obtain the QPI by Fourier transforming the real-space LDOS in sufficiently wide range. Second, it is applicable in parallel computation because the self-consistent parameters on each lattice site can be calculated separately. Third, the expansion scheme is very stable and efficient.

In our calculation, the magnitudes of the parameters are chosen as t₁₋₄ = 1,0.4,−2,0.04. The FS in the 2Fe/cell BZ derived by using this set of parameters is shown in fig. 1. There are hole pockets around Γ and electron ones around M, thus the FS topology as well as the band structure we use can qualitatively represent iron-based superconductors. Furthermore, magnetic unit cells are introduced where each unit cell accommodates four SC flux quanta and the linear dimension is N_x × N_y = 80 × 80, corresponding to a magnetic field B ≈ 8.32 Tesla, close to the experimental value (10 Tesla) [18]. V_ii and V_ij ($i=j\pm (\hat {x}\pm \hat {y})$) are chosen to be −2.8 and −2, respectively. Moreover we introduce 12 randomly distributed impurities with V_imp = 0.3. The ratio of N_impurity/N_vortex is chosen to be 3 in order to be consistent with the experimental observation (see fig. S2 in the supporting online material for ref. [18]). E^max_γ (E^min_γ) is chosen as 1.5 (−1.5) band width. Throughout the paper, we set the system to be 20% hole-doped. In calculating the self-consistent parameters, we use the Jackson kernel

$$\begin{equation} g_{k}=\frac{(N-k+1)\cos\frac{\pi k}{N+1}+\sin\frac{\pi k}{N+1}\cot\frac{\pi}{N+1}}{N+1}, \end{equation} \tag{ 11 }$$

with  = 0.001 and N = 500. For the LDOS we convolute the Lorentz kernel

$$\begin{equation} g_{k}=\frac{\sinh[\lambda(1-\frac{k}{N})]}{\sinh\lambda}, \end{equation} \tag{ 12 }$$

with λ = 4,  = 0.004 and N = λ/.

Zoom In Zoom Out Reset image size

First we consider the s_± case. For V_imp = 0, fig. 2(a) shows that there exists a negative-energy in-gap peak in the LDOS at the vortex core center, in agreement with the experimental observation [30] and our previous results based on the exact diagonalization method [25], indicating the reliability of the Chebyshev polynomial-expansion scheme. Figure 2(c) plots the spatial distribution of Z_j(ω), for B = 0 and at ω = Δ. The locations of the impurities can be clearly identified as the low-intensity spots. The corresponding $Z_{\vec{t} {q}}(\omega )$ is shown in fig. 2(d). We notice that there are high-intensity peaks at $\vec{t} {q}_{2}$ ((±π,0) and (0,±π)) which arise from the interpocket scattering between the hole and electron pockets where the SC OP changes sign. Furthermore the intensities at $\vec{t} {q}_{3}$ ((±π,±π)) are much weaker than those at $\vec{t} {q}_{2}$, consistent with the experimental observation [18]. This is because the intensity of the QPI is influenced by the relative sign of the SC OP between the two FS involved in the scattering process, through the coherence factor $C(\vec{t} {q})$ [18]. For iron-based superconductors, $\vec{t} {q}_{3}$ connects the electron pockets around (±π,0) and (0,±π), which preserves the sign of the SC OP if the paring symmetry is s_±. In the case of weak scalar potential scattering as considered in the present work, $C(\vec{t} {q})$ is strongly suppressed for the $\vec{t} {q}$ that preserves the sign of the SC OP [31], leading to the lower intensity at $\vec{t} {q}_{3}$ than that at $\vec{t} {q}_{2}$. On the other hand, along (±π,0) to (0,±π), there are broad high-intensity features. These should arise from the intrapocket scattering since the hole and electron pockets in our band structure are not circular, but more square-like. On the contrary, in ref. [27], these features do not show up, possibly due to the different band structure used. Upon applying the magnetic field, vortices are introduced into the system and their locations are denoted as the high-intensity spots in fig. 2(e). We notice that the vortices form a triangular lattice and none of them is pinned by the impurities since the impurity density is dilute and its strength is weak. In this work, the arrangement of the vortices as well as the position of the vortex cores are completely determined by the self-consistent calculation, and we did not assume any form of the vortex lattice. That is, we have tried various initial values and the converged results always led to the triangular vortex lattice as shown in fig. 2(e) (with or without the impurities). This means that the triangular vortex lattice is a stable state of the system. In all the calculations, the square vortex lattice as made in ref. [27] or other forms never appear as the converged results, thus they are energetically unfavorable for iron-based superconductors. Experimentally, the triangular vortex lattice has been observed in Ba_0.6K_0.4Fe₂As₂ (see fig. 5 in ref. [30]). In addition, in LiFeAs the vortices form quasihexagonal lattice at low field and become disordered at higher field. Nevertheless, they are never arranged into a square lattice [32]. Therefore in the present work, since we can only get the triangular vortex lattice self-consistently and the experiments did not observe the square vortex lattice in iron-based superconductors, we concentrate on the triangular one and do not consider other forms. In this case, from fig. 2(f) we can see that there exist additional peaks at $\vec{t} {q}_{3}$, whose intensities are enhanced by the application of the magnetic field and they are due to the interpocket scattering between different electron pockets (see fig. 1 in ref. [18]). Figure 2(b) shows the magnetic-field–induced change in the QPI intensities defined as $Z_{\vec{t} {q}}(B\neq 0)-Z_{\vec{t} {q}}(B=0)$. In the presence of the time-reversal symmetry breaking due to the magnetic field, the phase of the SC OP precesses by 2π about each vortex and its amplitude vanishes at the vortex core center. Both the phase variation and the inhomogeneity in the amplitude can scatter quasiparticles. While the former produces Doppler-shift scattering [33] that is odd under time reversal with $C(\vec{t} {q})$ being like that of magnetic impurities, the latter causes inhomogeneous Andreev scattering. Both of the scattering will enhance (suppress) the QPI intensities at those $\vec{t} {q}$ points that preserve (reverse) the sign of the SC OP [31]. Thus the sign-preserving scatterings at $\vec{t} {q}_{3}$ ($\vec{t} {q}_{3}$ connects the FS with the same sign of the SC OP) are enhanced while the sign-reversing scatterings at $\vec{t} {q}_{2}$ ($\vec{t} {q}_{2}$ connects the FS with the opposite sign of the SC OP) are suppressed. These phenomena are most pronounced around ω ≈ Δ. Away from it, the intensity at $\vec{t} {q}_{2}$ diminishes, since $\vec{t} {q}_{2}$ connects the hole and electron bands with different dispersion. On the other hand, spotty features remain at $\vec{t} {q}_{3}$ even at low and high energies, which may be ascribed to a field-induced Bragg (-like) component [18]. Both the locations and the sharpness of the QPI peaks shown in figs. 2(b), (d) and (f) are in reasonable agreement with the experimental data [18], suggesting that the experimentally observed peaks are indeed due to the QPI but not the Bragg scattering as argued in ref. [19].

Zoom In Zoom Out Reset image size

Next we consider the s-wave case. For V_imp = 0, the LDOS at B = 0 (black dotted line) and at the vortex core center (red solid line) shown in fig. 3(a) are also consistent with those obtained by exact diagonalization [34], again suggesting the validity of the current polynomial-expansion scheme. From fig. 3(c) we notice, the QPI in the presence of impurities at B = 0 exhibits no pronounced peaks at either $\vec{t} {q}_{2}$ or $\vec{t} {q}_{3}$, as compared to the clear peaks at $\vec{t} {q}_{2}$ in the s_± case as shown in fig. 2(d). This is because for s-wave pairing, $\vec{t} {q}_{2}$ also connects the FS with the same sign of the SC OP, thus the QPI intensity at this wave vector is suppressed for potential scattering as compared to the s_± case. After applying the magnetic field, the intensities at $\vec{t} {q}_{3}$ are enhanced and they form sharp peaks as shown in fig. 3(d), similar to the s_± case. At last, from the magnetic-field–induced change in the QPI intensities plotted in fig. 3(b) we can see, the intensities are enhanced at $\vec{t} {q}_{3}$ but remain almost unchanged at $\vec{t} {q}_{2}$. The lack of distinct structures at $\vec{t} {q}_{2}$ is in stark contrast to the s_± case and is inconsistent with the experimental observations [18]. Therefore, the different behaviour of the QPI intensities at $\vec{t} {q}_{2}$ in the s_± and s-wave pairing cases makes it possible to distinguish these two types of pairing symmetry since the STM experiment observed clear structures at $\vec{t} {q}_{2}$, thus excluding the possibility of s-wave pairing in iron-based superconductors.

Zoom In Zoom Out Reset image size

In summary, by using the Chebyshev polynomial expansion to directly solve large-scale BdG equations in real space, we have investigated the mixed-state effect on QPI in iron-based superconductors by assuming s_± or s-wave pairing symmetry. For the s_± pairing, the QPI intensities at $\vec{t} {q}_{2}$ which connects the FS with the opposite sign of the SC OP are suppressed by the application of the magnetic field while the situation at $\vec{t} {q}_{3}$ is reversed where $\vec{t} {q}_{3}$ connects the FS with the same sign of the SC OP. The obtained results at both B = 0 and B ≠ 0 are in qualitative agreement with experiment, suggesting that the experimentally observed peaks are indeed due to QPI. On the other hand, for the s-wave pairing, the QPI intensities at $\vec{t} {q}_{2}$ are featureless both with and without the magnetic field. Based on the available experimental data, the s-wave pairing can be excluded in iron-based superconductors.

In addition, we want to take this opportunity to compare our work with the previous ones which also studied the QPI in iron-based superconductors by using different band structures or methods. For s_± paring, without the magnetic field, our results are similar to those in refs. [20, 22, 35, 36] where strong QPI peaks appear at $\vec{t} {q}_{2}$ (see fig. 4 in refs. [20, 35, 36] and fig. 2 in ref. [22]). With the magnetic field, the QPI patterns in refs. [20–22, 35, 36] all lack clear peaks at $\vec{t} {q}_{3}$ (see fig. 4 in refs. [20, 36], fig. 8 in ref. [21], fig. 2 in ref. [22] and fig. 5 in ref. [35]). On the contrary, since we take both the phase variation and amplitude suppression induced by vortices into account, thus less approximation of the physics in involved in our calculation, we can get sharp QPI peaks at $\vec{t} {q}_{3}$ in agreement with experiment.

Furthermore we need to point out that although the method used in the present work is similar to that in ref. [27], the model and results are completely different. Our model is a multi-orbital one as compared to the single-orbital one adopted in ref. [27]. Thus the FS topology is drastically different between this two models (see fig. 1 and fig. 4 in ref. [27]). Besides, we assumed s_± or s-wave pairing symmetry in the present work while the authors in ref. [27] considered d-wave pairing symmetry. Since the QPI is closely related to the FS topology and the SC pairing symmetry, thus the QPI patterns derived in our work and those in ref. [27] share no similarity and cannot be compared to each other. Most importantly, in ref. [27] the authors concluded that only pinned vortex can lead to the correct QPI patterns. However in our work, the vortices are not pinned by any impurities and we can still get the QPI patterns in qualitative agreement with the experimental observations. This is consistent with the experimental conclusion that no appreciable correlation is observed between the location of vortices and the magnitude of field-induced change in the QPI intensity (see the supporting online material for ref. [18]), suggesting different scattering mechanism in iron-based superconductors and that studied in ref. [27].

This work was supported by NSFC (Grants No. 11204138 and No. 11175087), Natural Science Foundation of Jiangsu Province, China (Grant No. BK2012450), Natural Science Foundation of the Higher Education Institution of Jiangsu Province, China (Grant No. 12KJB140009), National Key Projects for Basic Research of China (Grant No. 2009CB929501), and by China Postdoctoral Science Foundation (Grant No. 2012M511297). PQT acknowledges the hospitality of ICTP.