\section{The Orbital Selective Spin Fluctuations Model} 
\label{Sec.Model}

We start briefly revising the distinctive features of the low-energy OSSF model originally derived in \cite{OSSFModel2}. We consider a general four-pocket model with two h-pockets at $\Gamma$, $\Gamma_\pm$ and two e-pockets at $X$ and $Y$, Figure~\ref{fs}. Mostly three orbitals contribute to the Fermi Surface, $yz$ and $xz$ for the $\Gamma_\pm$ pockets and $xy$, $yz/xz$ for the $X/Y$ pockets. The particular orbital arrangement follows from the space group of the iron-plane \cite{Vafekmodelo}.
A crucial consequence of this orbital composition is that the  $\Gamma_\pm$ hole pockets present opposite orbital nesting with the electron pockets, i.e. there is orbital mismatch between the $\Gamma_+$ and $X/Y$ and orbital match between $\Gamma_-$ and $X/Y$.

\begin{figure}[tbh]
\centering
{\includegraphics[width=0.45\linewidth]{fsmap_.pdf}}
\caption{Fermi surface for the four-pockets model. The colors represent the main orbital character of the Fermi surface. 
Notice the $yz/xz$ orbital nesting between the inner hole pocket $\Gamma_-$ and the $X/Y$ electron pockets. On the contrary the Fermi surface of the outer hole pocket $\Gamma_+$ presents orbital mismatch with the electron pockets.
The green/red arrows denote the orbital selective spin fluctuations (OSSF), connecting hole and electron pockets at different momenta, see Eqs.~\ref{sx}-\ref{hint_low}.}
\label{fs}
\end{figure}

The kinetic part of the Hamiltonian is derived adapting the low-energy model considered in~\cite{Vafekmodelo} that accounts for the orbital symmetry of the system. Each pocket is described using a spinor representation in the pseudo-orbital space 
\begin{equation}
\label{h0}
H_0^l=\sum_{{\bf k},\sigma} \psi^{l,\dagger}_{{\bf k}\sigma} \hat H_0^l \psi^l_{{\bf k}\sigma},
\end{equation}
with $\hat H^l_0= h_0^l\tau_0+\vec{h}^l\cdot\vec{\tau}$, $l=\Gamma,X,Y$ and $\tau$ matrices represent the pseudo-orbital spin. The spinors are defined as $\psi^{\Gamma}=(c_{yz},c_{xz})$ and $\psi^{X/Y}=(c_{yz/xz},c_{xy})$. Diagonalizing $\hat H_0$ we find the dispersion relations $E^{l\pm}=h_0^l\pm h^l$ with $h^l=|\vec{h}^l|$. We introduce the rotation from the orbital to the band basis, 
\begin{equation}
\label{hole_fer}
\begin{pmatrix}
 h_+ \\
 h_- \\
\end{pmatrix} = 
\begin{pmatrix}
u_{\Gamma} &-v_{\Gamma}\\
v_{\Gamma} & u_{\Gamma}\\
\end{pmatrix}
\begin{pmatrix}
 c_{yz} \\
 c_{xz} \\
\end{pmatrix} 
\end{equation}
with an analogous expression for the $X/Y$ pockets, provided that the corresponding orbital spinor is used. At  $X/Y$ only the $E^{X/Y+}$ band crosses the Fermi level, so in the following we will use $e_{X/Y}$ for the corresponding fermionic operators dropping the $+$ superscript. 
The interacting Hamiltonian in the spin channel is
\begin{equation}
\label{hint}
H_{int}=-1/2\sum_{{\bf q}}U_{\eta\eta'} \vec{S}^\eta ({\bf q}) \cdot \vec{S}^{\eta'} (-{\bf q}).
\end{equation}
$\eta,\eta'=yz,xz,xy$ are orbital indices and $U_{\eta\eta'}\sim U \delta_{\eta\eta'} + J_H (1-\delta_{\eta\eta'})$, with $U$ and $J_H$ being the usual Hubbard and Hund couplings. We consider only spin operators with intraorbital character $\vec{S}^\eta ({\bf q})=\sum_{{\bf k} ss'} \, c^{\eta \dagger}_{{\bf k} s} \vec \sigma _{s s'}c^\eta_{{\bf k}+{\bf q} s'}$ with $\sigma_{ss'}$ are the Pauli matrices for the spin operator. This choice is motivated by the general finding that intraorbital
magnetism is the dominant channel in IBS \cite{Pnictogenheight, Graser2009, Hirschfeldmultiorbital2011, MagneticReview, Nodalintra}. The relevant magnetic fluctuations occur at momenta ${\bf q}$ near ${\bf Q}_X$ or ${\bf Q}_Y$. At low energy we can project out the general interaction, Eq.~\ref{hint}, onto the fermionic excitations defined by the model Eq.~\ref{h0}. By using the rotation to the band basis, Eq.~\ref{hole_fer}, one can then establish a precise correspondence between the orbital and momentum character of the spin operators $\vec S_{X/Y}^\eta\equiv \vec S^\eta ({{\bf q}={\bf Q}_{X/Y}})$: 
\begin{eqnarray}
\label{sx}
\vec{S}_{X}^{yz} &=& \sum_{\bf k} (u_{\Gamma}h_{+}^{\dagger} 
+ v_{\Gamma} h_{-}^{\dagger} )\, \vec{\sigma} \, u_{X} e_X \\
\label{sy}
\vec{S}_{Y}^{xz} &=& \sum_{\bf k} (-v_{\Gamma}h_{+}^{\dagger} 
+ u_{\Gamma} h_{-}^{\dagger} )\, \vec{\sigma} \, u_{Y} e_Y 
\end{eqnarray}
where we drop for simplicity the momentum and spin indices of the fermionic operators. It then follows that the interacting Hamiltonian Eq.~\ref{hint} reduces to 
\begin{equation}
\label{hint_low}
H_{int}= - \frac{\tilde{U}}{2} \vec{S}^{yz/xz}_{X/Y} \cdot \vec{S}^{yz/xz}_{X/Y}
\end{equation}
where $\tilde{U}$ is the intraorbital interaction renormalized at low energy. As it is clear from the above equation, it is the projection of the generic interaction Hamiltonian Eq.~\ref{hint} onto the low-energy model Eq.~\ref{h0} that  generates OSSF. Indeed, since at low energy the $xz/yz$-fermionic states exist only around
$\mathbf{Q}_Y/\mathbf{Q}_X$, it turns out that the spin operators $\vec{S}_{X}^{\eta}$ with $\eta\neq yz$ and $\vec{S}_{Y}^{\eta}$ with $\eta\neq xz$ are absent in Eq.~\ref{hint_low}, so that there are no terms involving the Hund's coupling. 

The orbital selective character of the low-energy spin-excitations makes the interacting Hamiltonian for the spin-channel, Eq.~\ref{hint_low}, considerably simpler than the one obtained within a five-orbital model (see e.g. \cite{Graser2009}). As a matter of fact,  Eq.~\ref{hint_low}, while retaining the orbital dependence of the spin excitations, does not acquire a complex tensorial form and is instead formally equivalent to the spin-spin interacting Hamiltonian written in the band-basis $H^{band}_{int} \sim \ - \tilde{U} \ \vec{S}_{X/Y} \cdot \vec{S}_{X/Y}$ (see e.g. \cite{PreemptiveFernandes2012}).
This implies that one can analyze the spin-nematic phase following the same strategy of \cite{PreemptiveFernandes2012}, in which the nematic instability has been studied within an effective action derivation as a precursor effect of magnetism. 
The consequences of the inclusion of the orbital degree of freedom within a spin-nematic action  have been widely discussed especially in the analysis of the nematic phase of FeSe \cite{ OSSFModel2, OSSFModel1, OSSFNematic, Paper1}. 

In what follows we focus on the analysis of the superconductivity mediated by OSSF in the tetragonal phase for a generic four-pocket model representative of IBS. We will show that the orbital selectivity of the spin-fluctuations makes the RPA treatment extremely simple with respect to multiorbital calculations. In fact, given the intraorbital scalar character of the OSSF, the analysis turns out to be mathematically equivalent to the study of single-band systems, while retaining all the multiband and multiorbital information.

\section{Results}
\label{Sec:RPA}

Within standard RPA analysis \cite{Kuroki2008, Graser2009, GapsymmetryHirschfeld2011, Acommonthread}, the pairing is assumed to be mediated by spin and charge fluctuations. It has been shown that for IBS the charge susceptibility is more than one order of magnitude smaller than the spin susceptibility (see e.g. \cite{Graser2009}), therefore hereafter we focus on the spin channel only. 

The spin susceptibility for a generic multiorbital system is a four orbital indices tensor, $\chi^{\eta \eta'}_{\delta \delta'} (\mathbf{q},\Omega)$. This is obtained from the analytical continuation  $i\Omega_{m} \rightarrow \Omega + i0^{+}$ of the Matsubara spin-spin correlation function  
\begin{equation}
\chi^{\eta \eta'}_{\delta \delta'}(\mathbf{q}, i\Omega_{m}) = \int_{0}^{\beta} d\tau \, e^{i \Omega_{m} \tau} \,  \bigg\langle T_{\tau}  \vec{S}^{\eta\delta} (\mathbf{q}, \tau) \vec{S}^{\eta' \delta'} (\mathbf{-q}, 0)  \bigg\rangle 
\label{SFMulti1}
\end{equation}
where $\mathbf{q}$ is the momentum vector, $\beta = 1/k_BT$ is the inverse temperature, $\tau$ is the imaginary time and $\Omega_m = 2m\pi k_BT$ is the bosonic Matsubara frequency. The spin operator in the orbital space for the $\eta,\delta$ orbitals is defined as
\begin{equation}
\vec{S}^{\eta\delta} (\mathbf{q}, \tau)=\sum_{\mathbf{k} ss'} \, d^{\eta \dagger}_{\mathbf{k} s} \vec \sigma _{s s'}d^\delta_{\mathbf{k}+\mathbf{q} s'}
\label{SFMulti2}
\end{equation}
with $ \vec \sigma _{s s'}$ the Pauli matrices for the spin operator being $s,s'$ spin indices. Using this explicit definition given by Eq.~\ref{SFMulti2} and applying the Wick's theorem to Eq.~\ref{SFMulti1} the non-interacting spin susceptibility can be rewritten as
\begin{equation} 
\chi^{\eta \eta'}_{\delta \delta'}(\mathbf{q}, \Omega_m) = - \frac{1}{\beta} \sum_{\mathbf{k},n} Tr \left[G^{\delta \eta}  (\mathbf{k},i\omega_n)
 G^{\delta' \eta'} (\mathbf{k}+\mathbf{q},i\omega_n+i\Omega_m)  \right] 
\label{SFMulti3}
\end{equation}
where the spectral representation of the Green's function is given by the rotation to the orbital basis of the non-interacting Green's function in the band basis 
\begin{equation}
G^{\delta \eta}  (\mathbf{k},i\omega_n) = \sum_m \frac{a_m^\delta(\mathbf{k})a_m^{\eta*}(\mathbf{k})}{i\omega_n -E_m(\mathbf{k})}
\label{SFMulti4}
\end{equation}
where $\omega_n = (2n+1)\pi k_BT$ is the fermionic Matsubara frequency and $a_m^\eta(\mathbf{k})$ the matrix elements connecting the orbital ($\eta$) and the band space ($m$) determined by diagonalization of the tight-binding Hamiltonian.
Performing the Matsubara frequency summation and setting $\Omega \rightarrow 0$, the static spin susceptibility reads 
\begin{equation}
\chi^{\eta \eta'}_{\delta \delta'}(\mathbf{q}) =  - \sum_{\mathbf{k},m n} \frac{a_m^\delta(\mathbf{k})a_m^{\eta*}(\mathbf{k})a_n^{\eta'}(\mathbf{k}+\mathbf{q})a_n^{\delta'*}(\mathbf{k}+\mathbf{q})}{E_{n}(\mathbf{k}+\mathbf{q})-E_{m}(\mathbf{k})} {f(E_{n}(\mathbf{k}+\mathbf{q}))-f(E_m(\mathbf{k}))}
\label{SFMulti5}
\end{equation}
with $f(E_m(\mathbf{k}))$ the Fermi distribution function. The RPA spin-fluctuation is given in the form of a Dyson-type equation with the spin interaction $\hat{U}_S$ defined in terms of the multiorbital interaction parameters $U,U',J_H$ \cite{Kuroki2008, Graser2009}.
Analogously, the singlet pairing vertex driven by spin-fluctuations can be computed on the low-energy sector in terms of the RPA spin susceptibility \cite{Graser2009, Roleofvertex}. 
The variety of diagrams contributing to the SC vertex is large given the number of orbitals included making unfeasible to draw the possible Feynman diagrams up to orders larger than one.
The gap equation for the multiorbital model can be computed numerically by taking into account the singlet pairing vertex as an eigenvalue problem in which the largest eigenvalue leads to the highest transition temperature and its eigenfunction determines the symmetry of the gap (see e.g. \cite{Graser2009, GapsymmetryHirschfeld2011, Acommonthread, Kemper2010}). The main result is that an anisotropic sign changing s-wave, $s_\pm$, is found as the dominant symmetry for system parameters compatible with moderately doped IBS, in agreement with experiments, e.g. \cite{NMRgap, NMRgap2, ResonancepeakNSIBS,ARPESreview}.
\subsection{Magnetic excitations in the OSSF model: RPA analysis} \label{Sec.SFRPA}

Within the OSSF model, the situation is substantially simplified as compared with the five-orbital RPA approach due to the orbital-selectivity of the spin-fluctuations.
By assuming the spin-operator to be intraorbital, the spin-susceptibility of Eq.~\ref{SFMulti1} can be reduced to a two-orbital indices matrix
\begin{equation}
\chi^{\eta \eta'}(\mathbf{q}, i\Omega_{m}) = \int_{0}^{\beta} d\tau \, e^{i \Omega_{m} \tau} \,  \bigg\langle T_{\tau}  \vec{S}^{\eta} (\mathbf{q}, \tau) \vec{S}^{\eta'} (\mathbf{-q}, 0)  \bigg\rangle 
\label{SFOSSF1}
\end{equation}
In addition, the low-energy projection further simplified the spin-susceptibility structure as the low-energy states are defined only around high symmetry points and have a well defined orbital character described by Eq.~\ref{h0}. 
As a consequence, also the Green's functions are defined only for  $l= \Gamma, X, Y$ as $G^{l}(\mathbf{k},i\omega_n)={\hat{\cal U}^l}(\mathbf{k},i\omega_n)diag (g^{l_+}(\mathbf{k},i\omega_{n}),g^{l_-}(\mathbf{k},i\omega_{n})){\hat{\cal U}^{l^{-1}}}(\mathbf{k},i\omega_n)$. Here $\hat{\cal U}^l$ are the matrices that diagonalize the $l$-Hamiltonian and $g^{l\pm} (\mathbf{k}, i\omega_n)= (i\omega_n-E^{l\pm}_\mathbf{k})^{-1}$ the Green's functions in the band basis. Substituting the intraorbital spin operator $\vec{S}^{\eta} (\mathbf{q}, \tau)$ and applying Wick’s theorem to Eq.~\ref{SFOSSF1}, the intraorbital spin susceptibility in the low-energy projection can be read as
\begin{equation} 
\chi^{l l'}(\mathbf{q}, i\Omega_m) = - \frac{1}{\beta} \sum_{\mathbf{k},n} Tr \left[\hat{G}^l(\mathbf{k},i\omega_n)
 \hat G^{l'} (\mathbf{k}+\mathbf{q},i\omega_n+i\Omega_m)  \right] \,.
 \label{SFOSSF2}
\end{equation}
Eq.~\ref{SFOSSF2} represents the spin susceptibility between two pockets $l$ and $l^{'}$ and depends on the transferred momentum $\mathbf{q}=\mathbf{k'}-\mathbf{k}$ and the external frequency $\Omega$. 
Performing the Matsubara frequency summation and setting $\Omega \rightarrow 0$ in Eq.~\ref{SFOSSF2} 
we find the static susceptibility for the $l_{\pm}$ $l^{'}_{\pm}$ pockets in terms of the Fermi distribution function $f(E^{l_\pm}_{\mathbf{k}})$ and  $|(u/v)^{l}(\mathbf{k})|$ that are the elements of the rotational matrix $\hat{\cal U}^l$ connecting the orbital and the band space 
\begin{equation}
\chi^{l_\pm l'_\pm}_{\eta}(\mathbf{q}) = \sum_{\mathbf{k}} \frac{f(\epsilon^{l'_\pm}_{\mathbf{k}+\mathbf{q}})-f(\epsilon^{l_\pm}_{\mathbf{k}})}{\epsilon^{l'_\pm}_{\mathbf{k}+\mathbf{q}}-\epsilon^{l_\pm}_{\mathbf{k}}} |(u/v)^l_{\mathbf{k}}|^2 |(u/v)^{l'}_{\mathbf{k}+\mathbf{q}}|^2 \,.
\label{SFOSSF3}
\end{equation}
The resulting static susceptibility Eq.~\ref{SFOSSF3} acquires a much simpler expression with respect to the multiorbital spin-susceptibility in Eq.~\ref{SFMulti5}. 
In particular, within the OSSF model, we only consider magnetic fluctuations close to momenta $\mathbf{Q}_X=(\pi,0)$ and $\mathbf{Q}_Y=(0,\pi)$, so the two most relevant spin susceptibilities for a four-pocket model will involve the $yz$ orbital coming from the interaction between the $\Gamma_\pm$ holes with the $X$ electron pockets, and the $xz$ orbital coming from the $\Gamma_\pm$ holes with the $Y$ electron pockets near $\mathbf{Q}_X$ and $\mathbf{Q}_Y$ respectively. To better compare with the results obtained within five-orbital model calculations, in what follows we account also for the the spin susceptibility centered around $\mathbf Q_M=(\pi,\pi)$ that involves the $xy$ orbital only and describes the spin exchange between the $X$ and $Y$ electron pockets. Notice that the inclusion of the $\mathbf Q_M=(\pi,\pi)$ susceptibility within the OSSF is straightforward due to the orbital selectivity of the spin fluctuation involved in the exchange process. 

The RPA spin susceptibilities are obtained in the form of Dyson-type equations. The results of the resummation reads  
\begin{equation}
\chi^{l_\pm l'_\pm}_{\eta \rm{RPA}}(\mathbf{q})  = \frac{\chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})}{1- \tilde{U}\chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})}
\label{SFOSSF4}
\end{equation}
with $\chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})$ the non-interacting spin susceptibility given by Eq.~\ref{SFOSSF3}. Notice that, within our model Eq.~\ref{hint_low}, both the low-energy effective coupling $\tilde{U}$  and the OSSF $\chi_{\eta}^{l_\pm l'_\pm}$ are intraorbital and have scalar character. As a consequence, the RPA spin susceptibility given by Eq.~\ref{SFOSSF4} is straightforward and inherits the orbital selectivity and scalar form of the OSSF. As we will discuss in the following Sections~\ref{Sec.VertexRPAandGAP} and ~\ref{SubSec.GAP} this aspect allows us to derive analytical expressions for the pairing vertex and SC gaps, making the computations much more manageable than for the multiorbital models. 

To get insight into the previous result, we perform a numerical estimate for the RPA spin susceptibility given by Eq.~\ref{SFOSSF4} for the four-pocket model in Figure~\ref{fs}. 
We show  $\chi^{\Gamma_+ X}_{yz \rm{RPA}}(\mathbf{q})$, $\chi^{\Gamma_- X}_{yz \rm{RPA}}(\mathbf{q})$ and $\chi^{XY}_{xy \rm{RPA}}(\mathbf{q})$ n Figure~\ref{fig:122RPA}. 
In the upper panel we show 3D color maps in $(q_x,q_y)$. In the bottom panel we show the $q_x/q_y$ 2D cuts. The middle point for all the 2D plots is $\mathbf{Q}_X=(\pi,0)$ for the electron-hole spin susceptibilities and $\mathbf{Q}_M=(\pi,\pi)$ for the electron-electron one. The right part of each graphic represents the $\chi^{l_\pm l'_\pm}_{\eta \rm{RPA}}$ by moving from the center of the 3D plots in the $q_x$ momentum, maintaining $q_y$ constant, and analogously on the other side. This representation allows to easily compare the relative weight of the different susceptibilities.

Notice that the contributions of the $Y$ pocket (not shown) are equivalent to those for the $X$ pockets with a $\pi/2$ rotation, since in the tetragonal phase the susceptibility is isotropic in both directions. 
\begin{figure}[tbh]
\centering
{\includegraphics[width=0.95\linewidth]{Fig2_122_5.pdf}}
\caption{RPA spin susceptibility for a four-pocket model in the tetragonal phase in 3D (upper panel) and 2D cuts (bottom panel) around the high symmetry points $X$ for the hole-electron sector and $M$ for the electron-electron sector. The right-hand part of the 2D cuts represents the $q_x$ direction and the left-hand part the $q_y$ direction starting from $\mathbf{Q}_X=(\pi,0)$ or $\mathbf{Q}_M=(\pi,\pi)$ points.}
\label{fig:122RPA}
\end{figure}

From Figure~\ref{fig:122RPA} we can highlight some main results:

(i) The orbital-selective RPA spin susceptibilities around $\mathbf{Q}_X$ and $\mathbf{Q}_M$ show a clear momentum-dependent structure of the peaks. This can be explained due to the degree of orbital nesting between pockets. The orbital nesting indicates the relative orbital composition between the two pockets involved in the spin-exchange mechanism. 
In Figure~\ref{fig:122RPA} we can see that when there is an orbital mismatch, as is the case of $\Gamma_+$ and $X$ pockets, the spin susceptibility develops two incommensurate peaks around $\mathbf{Q}_X$. In contrast, if there is an orbital match between pockets, i.e. the case of $\Gamma_-$ and $X$, the spin susceptibility develops a single commensurate peak at the $\mathbf{Q}_X=(\pi,0)$. For the $\chi^{XY}_{xy \rm{RPA}}(\mathbf{q})$ susceptibility there is a total mismatch between the $xy$ orbital of the electrons pockets. Thus, the spin susceptibility is totally incommensurate and develops four symmetric peaks that correspond to the overlaps of the $xy$ orbital contribution around the $(\pi,\pi)$ point. 

(ii) The main contribution to the spin susceptibility comes from the $\mathbf{Q}_{X/Y}$ spin-mode, between the $\Gamma_-$ and the $X/Y$ pockets. This is a consequence of the orbital composition of the nested bands. The Fermi surface of $\Gamma_-$ and and $X/Y$ are, in fact, characterized by fully matching orbitals.

We compare our results with a numerical calculation of the RPA spin susceptibility for a five-orbital model.
By tuning the filling and the crystal field we can consider separately two different cases: the first one corresponds to a four-pocket model with better band nesting between the $\Gamma_-$ and the electron pockets, and the second case to a four-pocket model with better band nesting between the $\Gamma_+$ and the $X/Y$ pockets. We obtain for the first case a commensurability of the RPA spin susceptibility at the $\mathbf{Q}_X=(\pi,0)$, whereas in the second case we get non-commensurate peaks around $\mathbf Q_X=(\pi,0)$. Therefore, the same orbital modulation for the momentum dependence of the RPA spin susceptibility is obtained within the OSSF and the multiorbital models. 
We also compute the RPA spin susceptibility coming from the electron-electron sector within the five-orbital model. As expected, we find that the contribution from this sector is negligible in comparison with the one for the hole-electron sector \cite{ReviewHirschfeld}. 
Notice that the orbital modulation of the RPA spin susceptibility is an effect that comes directly from taking into account the orbital degree of freedom of IBS. Within a low-energy multiband model, which is blind to the orbital content of the system, the spin susceptibility gives rise to two equivalents peaks with the same commensurable narrow structure. Moreover, it is worth noticing that, while a generic correlation between the orbital make-up of the Fermi surface and the momentum-dependent structure of the RPA spin excitation has already been highlighted within multiorbital models (e.g. \cite{GapanisotropycomparitionRPA}), the explicit link between orbital nesting and momentum dependence of the spin-susceptibility is a new result of the RPA analysis within the OSSF model.
\begin{figure}[tbh]
\centering     
\subfloat[OSSF model]{\includegraphics[width=75mm]{CaminoRaquel.pdf}} 
\subfloat[Five-orbital model]{\includegraphics[width=75mm]{CaminoMJ.pdf}}
\caption{Cuts of the total RPA spin susceptibility along the high symmetry directions $\Gamma X M \Gamma$ in the Brillouin zone within the (a) OSSF model and the (b) five-orbital model. The interaction parameters have been chosen for (a) as $\tilde{U}$ = 1~eV and (b) as $U$ = 1.2~eV, $J_H=0.25U$ and $U'=U-2J_H$.  The discontinuity of the susceptibility for the OSSF model is due to the fact that this model is limited to low energies and we only get the contribution near the Fermi surface around the high symmetry points $\Gamma, X, M$.}
\label{fig:CaminoGXMG}
\end{figure}

In Figure~\ref{fig:CaminoGXMG} we compare the cuts along the main symmetry directions of the RPA spin susceptibility (including the intraband ones) computed within the OSSF model and the multiorbital model. The calculation performed within the OSSF model reproduces remarkably well the overall momentum-dependence of the spin-spectrum as well as the relative height and width of the various peaks.
This comparative analysis proves that we can obtain a reliable description of the spin spectrum within the OSSF model of Eq.~\ref{SFOSSF3}, without dealing with the tensorial form of the spin susceptibility given in the five-orbital model, Eq.~\ref{SFMulti5}. 

\subsection{Superconductivity mediated by OSSF} \label{Sec.VertexRPAandGAP}

In Figure~\ref{fig:Vertex} we draw the Feynman's diagrams that contain only the (scalar) spin susceptibilities up to finite order in perturbation of $\tilde{U}$ for electrons of opposite spin and momentum. 
\begin{figure}[tbh]
\centering
{\includegraphics[width=0.99\linewidth]{Pairing_vertex.pdf}}
\caption{Pairing vertex in the random phase approximation up to fourth order within the OSSF model.}
\label{fig:Vertex}
\end{figure}

From this diagrammatic expansion we learn that:

(i) It is practically unfeasible to draw analytically the Feynman's diagrams involved in the pairing interaction vertex within the five-orbital model. This is due to the complex tensorial structure of the pairing vertex. On the contrary, given the intraorbital scalar character of the low-energy orbital-selective spin-susceptibility we can easily perform the diagrammatic expansion explicitly.

(ii) As a consequence of the projection to a constrained orbital space within the OSSF model, the diagrammatic expansion for the RPA pairing vertex is formally equivalent to the one for the band-basis model which does not contain the orbital information of the spin-fluctuations exchange. Within the multiband model, the RPA pairing vertex is composed by the exchange of spin fluctuations connecting an electron pocket with a hole pocket. On the other hand, within the OSSF model the essential difference is that the low energy exchanged spin fluctuations connect a hole with an electron pocket with the {\it same orbital content} $yz$ or $xz$. In this way, we retain the simplicity of the analysis of the Feynman's diagrams within the multiband model and, at the same time, we account for the orbital degree of freedom of the system. 

The leading RPA diagrams for the vertex, Figure~\ref {fig:Vertex}, can be written as
\begin{equation}
\Gamma^{l_\pm l'_\pm}_{\eta}(\mathbf{q}) =  
\tilde{U} 
+   \frac{ \tilde{U}^2  \chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})}{1-\tilde{U} \chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})}
+ \frac{\tilde{U}^3  \chi_{\eta}^{{l_\pm l'_\pm}^2}(\mathbf{q})}{1-\tilde{U}^{2} \chi_{\eta}^{{l_\pm l'_\pm}^2}(\mathbf{q})} 
\label{SFOSSF5}
\end{equation}
where $\mathbf{q}=\mathbf{k'}-\mathbf{k}$ is the transferred momentum, $\tilde{U}$ is the intraorbital effective coupling and $\chi_{\eta}^{l_\pm l'_\pm}(\mathbf{q})$ the intraorbital susceptibility given in formula Eq.~\ref{SFOSSF3}. 
As we are considering the spin channel only, the SC vertex is proportional to the RPA spin susceptibility and preserves the orbital and momentum dependencies and the physical properties discussed in the previous section. For instance, we get the same criterion of commensurability or incommensurability depending on the orbital matching between pockets. In the same way, we also obtain that the dominant contribution to the RPA pairing vertex is given by the spin-fluctuations exchange between hole and electron pockets, being the largest between pockets having matching orbital, i.e. $\Gamma_- - X/Y$.

\subsection{Superconducting gaps} \label{SubSec.GAP}
We solve the BCS gap equations mediated by the OSSF for a four pocket model for IBS in the tetragonal phase. In order to better highlight the effect of the orbital nesting between the hole and electron nested pockets we consider a case in which the two hole pockets are almost equivalent in size, see Figure \ref{fig:gapnodo}. In this way the degree of band nesting between the hole and electron pockets is the same for the inner and outer hole pockets, but the orbital matching condition determined by the space group of the iron-plane \cite{Vafekmodelo} is different for $\Gamma_\pm$.
For simplicity we show here the equations for the $yz$ and $xz$ components only. However, we have numerically solved the full orbital system in which we also consider the SC vertex mediated by the exchange of spin-fluctuations between the $X-Y$ pockets.
The pairing Hamiltonian for the $yz$ and $xz$ orbital contributions reads
\begin{eqnarray}
H^{pairing}_{yz,xz} &=& - \sum_{\mathbf{k},\mathbf{k'}} \Gamma^{\Gamma_+ X}_{yz \ \mathbf{k}\mathbf{k'}} \left[(u^\Gamma_{\mathbf{k}})^2 h^{+\dagger}_{\mathbf{k}} h^{+\dagger}_{\mathbf{-k}} (u^X_{\mathbf{k}})^2 e^{X}_{\mathbf{k'}} e^{X}_{\mathbf{-k'}} \right] - \sum_{\mathbf{k},\mathbf{k'}} \Gamma^{\Gamma_- X}_{yz \ \mathbf{k}\mathbf{k'}} \left[(v^{\Gamma}_{\mathbf{k}})^2 h^{-\dagger}_{\mathbf{k}} h^{-\dagger}_{\mathbf{-k}} (u^X_{\mathbf{k}})^2 e^{X}_{\mathbf{k'}} e^{X}_{\mathbf{-k'}} \right] \nonumber \\
&-& \sum_{\mathbf{k},\mathbf{k'}} \Gamma^{\Gamma_+ Y}_{xz \ \mathbf{k}\mathbf{k'}} \left[(v^\Gamma_{\mathbf{k}})^2 h^{+\dagger}_{\mathbf{k}} h^{+\dagger}_{\mathbf{-k}} (u^Y_{\mathbf{k}})^2 e^{Y}_{\mathbf{k'}} e^{Y}_{\mathbf{-k'}} \right] - \sum_{\mathbf{k},\mathbf{k'}} \Gamma^{\Gamma_- Y}_{xz \ \mathbf{k}\mathbf{k'}} \left[(u^{\Gamma}_{\mathbf{k}})^2 h^{-\dagger}_{\mathbf{k}} h^{-\dagger}_{\mathbf{-k}} (u^Y_{\mathbf{k}})^2 e^{Y}_{\mathbf{k'}} e^{Y}_{\mathbf{-k'}}  \right] + h.c. \nonumber \\
\label{BCS1}
\end{eqnarray}
where $\Gamma_{\eta \ \mathbf{k}\mathbf{k'}}^{l_\pm l'_\pm}$ with $\eta=yz/xz$ is the RPA pairing vertex given by Eq.~\ref{SFOSSF5} for the different pockets ${l_\pm l'_\pm}$ and $(u^l_\mathbf{k})^2$ and $(v^l_\mathbf{k})^2$ are the coherence factors that connect the orbital and the band basis and account for the pockets orbital character.

The pairing Hamiltonian given by Eq.~\ref{BCS1} is solved in the mean field approximation by defining the orbital dependent SC order parameters for the hole sector
($\Delta_{yz}^{\Gamma_+ X},\Delta_{yz}^{\Gamma_-X},\Delta_{xz}^{\Gamma_+Y},\Delta_{xz}^{\Gamma_- Y}$) 
and for the electron sector
($\Delta_{yz}^{X\Gamma_+},\Delta_{yz}^{X\Gamma_-},\Delta_{xz}^{Y \Gamma_+},\Delta_{xz}^{Y\Gamma_-}$).
The resulting linearized gap equations at $T=0$ read as
\begin{eqnarray}
\Delta^{\Gamma_+X}_{yz\mathbf{k'}}&=&-\sum_{k} \frac{\Gamma^{\Gamma_+ X}_{yz\mathbf{k}\mathbf{k'}} }{|v^X_{F\mathbf{k}}|} (u^{X}_{\mathbf{k}})^{4} \left[
\Delta^{X\Gamma_+}_{yz\mathbf{k}} +
\Delta^{X\Gamma_-}_{yz\mathbf{k}} \right]  \label{BCS2.1}  \\
\Delta^{\Gamma_- X}_{yz\mathbf{k'}} &=& - \sum_{k} \frac{\Gamma^{\Gamma_- X}_{yz\mathbf{k}\mathbf{k'}} }{|v^X_{F\mathbf{k}}|} (u^{X}_{\mathbf{k}})^{4} \left[
\Delta^{X \Gamma_+}_{yz\mathbf{k}} +
\Delta^{X \Gamma_-}_{yz\mathbf{k}} \right]  \label{BCS2.2}  \\
\Delta^{\Gamma_+ Y}_{xz\mathbf{k'}} &=& - \sum_{k} \frac{\Gamma^{\Gamma_+ Y}_{xz\mathbf{k}\mathbf{k'}} }{|v^Y_{F\mathbf{k}}|} (u^{Y}_{\mathbf{k}})^{4} \left[
\Delta^{Y \Gamma_+}_{xz\mathbf{k}} +
\Delta^{Y \Gamma_-}_{xz\mathbf{k}} \right]  \label{BCS2.3}  \\
\Delta^{\Gamma_- Y}_{xz\mathbf{k'}} &=&  - \sum_{k} \frac{\Gamma^{\Gamma_- Y}_{xz\mathbf{k}\mathbf{k'}} }{|v^Y_{F\mathbf{k}}|} (u^{Y}_{\mathbf{k}})^{4}  \left
\Delta^{Y \Gamma_+}_{xz\mathbf{k}} +
\Delta^{Y \Gamma_-}_{xz\mathbf{k}} \right]  \label{BCS2.4} \\
\Delta^{X \Gamma_+}_{yz\mathbf{k'}} &=& - \sum_{k} \frac{\Gamma^{\Gamma_+ X}_{yz\mathbf{k}\mathbf{k'}} }{|v^{h_+}_{F\mathbf{k}}|} (u^{\Gamma}_{\mathbf{k}})^2 \left[ (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ X}_{yz\mathbf{k}} + (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ Y}_{xz\mathbf{k}} \right]  \label{BCS2.5}  \\
\Delta^{X \Gamma_-}_{yz\mathbf{k'}} &=&  - \sum_{k} \frac{\Gamma^{\Gamma_- X}_{yz\mathbf{k}\mathbf{k'}} }{|v^{h_-}_{F\mathbf{k}}|} (v^{\Gamma}_{\mathbf{k}})^2 \left[ (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_- X}_{yz\mathbf{k}} + (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_-Y}_{xz\mathbf{k}} \right] \label{BCS2.6}  \\
\Delta^{Y \Gamma_+}_{xz\mathbf{k'}} &=&  -\sum_{k} \frac{\Gamma^{\Gamma_+ Y}_{xz\mathbf{k}\mathbf{k'}} }{|v^{h_+}_{F\mathbf{k}}|} (v^{\Gamma}_{\mathbf{k}})^2 \left[ (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ X}_{yz\mathbf{k}} + (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+Y}_{xz\mathbf{k}} \right] \label{BCS2.7}   \\
\Delta^{Y \Gamma_-}_{xz\mathbf{k'}} &=&   - \sum_{k} \frac{\Gamma^{\Gamma_- Y}_{xz\mathbf{k}\mathbf{k'}} }{|v^{h_-}_{F\mathbf{k}}|} (u^{\Gamma}_{\mathbf{k}})^2 \left[ (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_-X}_{yz\mathbf{k}} + (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_-Y}_{xz\mathbf{k}} \right] \label{BCS2.8}  
\end{eqnarray}
with $v^l_{F\mathbf{k}}=\partial_\mathbf{k}(\epsilon_{\mathbf{k}}^{l\pm })$ the Fermi velocity for the pocket $l$. Eqs.~\ref{BCS2.1}\,-\,\ref{BCS2.8} represent the orbital components for each gap. Then, we define the total low-energy band gaps $\Delta_\mathbf{k}^{l}$ as
\begin{eqnarray}
\Delta_\mathbf{k}^{\Gamma_+} &=& (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ X}_{yz\mathbf{k}} + (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ Y}_{xz\mathbf{k}} \label{BCS3.1} \\
\Delta_\mathbf{k}^{\Gamma_-} &=& (v^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_- X}_{yz\mathbf{k}} + (u^{\Gamma}_{\mathbf{k}})^2 \Delta^{\Gamma_+ Y}_{xz\mathbf{k}}  \label{BCS3.2}\\
\Delta_\mathbf{k}^{X} &=& (u^{X}_{\mathbf{k}})^2 ( \Delta^{X\Gamma_+ }_{yz\mathbf{k}} + \Delta^{X\Gamma_-}_{yz\mathbf{k}} ) \label{BCS3.3} \\
\Delta_\mathbf{k}^{Y} &=& (u^{Y}_{\mathbf{k}})^2 ( \Delta^{Y \Gamma_+}_{xz\mathbf{k}} + \Delta^{Y\Gamma_ -}_{xz\mathbf{k}})  \label{BCS3.4}
\end{eqnarray}
where each low-energy band gap involves the sum of the different orbital contributions weighted by the correspondent coherent factors of the pocket. 
The gap functions $\Delta_\mathbf{k}^{l}$ given by the set of coupled equations in \ref{BCS3.1}\,-\,\ref{BCS3.4}, contain information on the spatial and orbital structures of the pairs. 

We solve numerically the linearized gap equations by searching for the largest eigenvalue that corresponds to the leading instability of the system. Then, we calculate its corresponding eigenfunction which  determines the symmetry and the structure of the gap function. 

\begin{figure}[tbh]
\centering  
\subfloat[]{\includegraphics[scale=0.52]{FSnodos.pdf}} 
\subfloat[]{\includegraphics[scale=0.4]{Gapnodos.pdf}}
\caption{RPA gap functions within the OSSF model (a) plotted on the Fermi surface pockets (red circles for positive and black circles for negative gaps), (b) plotted as a function of angle from $0$ to $2\pi$. The SC gap functions present $s_\pm$ symmetry with angular modulations controlled by the orbital composition of the nested bands. The hierarchy of the band gaps follows from the orbital nesting properties of the matching Fermi surface pockets.}
\label{fig:gapnodo}
\end{figure}

The results for the gap functions are summarized in Figures~\ref{fig:gapnodo}. We find an anisotropic $s_\pm$ symmetry, in agreement with multiorbital RPA results \cite{Graser2009, GapanisotropycomparitionRPA}. 
Given the simplicity of our treatment, it can be easily explained how the overall momentum dependence of the band gaps is determined by the momentum-dependence of the pairing vertex (i.e of the spin-susceptibility) and the orbital make-up of the Fermi surface. As already mentioned, the key parameter that controls the angular modulation of the gap is the orbital matching between the low-energy states connected by the spin-fluctuations that determine the structure of the spin-susceptibilities (and thus of the pairing vertex).

The orbital degree of freedom is also responsible for the hierarchy of the band gaps. In fact, given a similar condition of band nesting between the two hole pockets and the electrons, one could naively predict similar gaps opening along $\Gamma_\pm$. However, from Figure \ref{fig:gapnodo}(b), it can be seen that the gap in  $\Gamma_-$ is larger than the gap in $\Gamma_+$. This is a direct consequence of the orbital composition of the nested pockets: in fact, the exchange of intraorbital spin-fluctuations between nested pockets having matching orbital, as $\Gamma_- - X/Y$ is stronger than the spin exchange between the pockets presenting orbital mismatch, such as $\Gamma_+ - X/Y$. This reflects in a larger spin-susceptibility and SC vertex for the pockets presenting a better degree of orbital matching. 

The results we obtain can be found within more complete RPA five-orbital model calculations. The main advantage of our derivation is that it gives us a new physical insight on the role of the orbital nesting in controlling the features of the SC gaps. We can compare our results, for example, with the ones obtained in \cite{GapanisotropycomparitionRPA}, where the anisotropy of the gaps for a five-orbital model is discussed in detail by using RPA calculations for the exchange of spin and charge fluctuations. The authors conclude that the anisotropy of the gap on the different Fermi surface pockets arises from an interplay of the orbital make-up of the states on the Fermi surface together with the momentum dependence of the fluctuation-exchange pairing interaction. In addition, to minimize the repulsion between the electron pockets, the gap can present accidental nodes. 
By comparing the results obtained within the OSSF model, Figures~\ref{fig:gapnodo}, and the five-orbital calculation, Figure 4b in  \cite{GapanisotropycomparitionRPA}, one can see that the OSSF calculations qualitatively reproduce the main feature of the multiorbital ones i.e. the modulated $s_\pm$ symmetry for the gap, with accidental nodes. We also reproduce the correct hierarchy of the band gaps with $|\Delta_\mathbf{k}^{\Gamma_-}| > |\Delta_\mathbf{k}^{\Gamma_+}|$. The main novelty here is that our analysis shows how the band gap hierarchy results from the degree of orbital matching between nested Fermi surfaces.  

\section{Conclusions}

In this work we show that the OSSF model provides a reliable description of the magnetic excitations in IBS and represents the minimal model to study spin-mediated superconductivity in this system. 

We compute the spin susceptibility for a four pocket model in the tetragonal phase and compare the results with the ones obtained for the five-orbital model. Depending on the degree of orbital nesting between pockets, we get an orbital modulation of the spin susceptibility that gives rise to commensurate or incommensurate peaks in the spin susceptibility when there is an orbital match or mismatch between the hole and the electron pockets respectively. The spin-exchange between pockets with matching orbitals is the larger contribution to the total spin-susceptibility.
By comparing the total spin susceptibility of the OSSF model with the one obtained within the five-orbital model we show that the OSSF reproduces qualitatively well the overall momentum dependence and the relative heights and widths of the peaks located at different momenta.
This is a remarkable result considering that the OSSF model is a low-energy effective model that only considers the $yz,xz,xy$ orbitals.

Within a spin-driven scenario, we compute the pairing vertex and the corresponding gap equations. 
Notably, within the OSFF model it is possible to draw analytically the Feynman's diagrams involved in the pairing vertex, something which is almost unfeasible within the five-orbital model due to the large number of different possible diagrams. 
We solve the corresponding BCS gap equations, analyse the gap symmetry and structure of the band gaps for the different pockets. We find anisotropic $s_\pm$ gap functions in agreement with multiorbital calculations and explain the angular modulation of the gaps as due to the orbital composition of the pockets connected by the OSSF mediating the pairing. 

The analysis of the bands gap structure shows that the angular dependence and the magnitudes of the different gaps depend directly on the degree of orbital matching between the hole and the electron pockets. In particular, the orbital nesting between the different Fermi surfaces is the parameter that controls the hierarchy of the band gaps, in contrast with the naive expectation of having a gap proportional to the degree of band-nesting. 
This new insight on the role of the orbital composition of the nested pockets breaks the usual paradigm based on the Fermi surface matching as the parameter controlling the instabilities and establish the orbital nesting as the crucial control parameter. It is worth noticing that this new scenario is more robust since it does not depend on the fine tuning matching of size and shape of the Fermi surfaces, but relies instead on the orbital composition of the bands determined by the space group of the iron-plane of IBS.\\

In conclusion, in this work we analyze the magnetic excitations and SC gaps for a generic IBS system in the tetragonal phase using the OSSF model. Despite its simplicity, the OSSF model is able to reproduce well all the relevant qualitative features of the spin-excitations and pairing interactions of the multiorbital description. It allows for analytical treatment making the interpretation of the results straightforward. Due to the orbital-selectivity of the spin-fluctuation, it allows to include in a very simple way additional interacting channels, as we showed explicitly by considering the electron-electron interaction besides the hole-electron spin-exchange. The simplified frame of the OSSF allows us to gain physical insight on the role of orbital nesting in determining the angular modulations and the hierarchy of the band gaps. 
This result proves that the OSSF model is the minimal low-energy model to study spin-mediated superconductivity in IBS as it correctly incorporates in the low-energy description the orbital composition of the bands determined by the space group of the iron-plane.  

L.~F. acknowledges financial support from the European Unions Horizon 2020 research and innovation programme under the Marie Sk\lambda odowska-Curie grant SuperCoop (Grant No 838526). We acknowledge funding from Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación and FEDER funds (EU) via Grant No. FIS2015-64654-P and  PGC2018-099199-B-I00.

