Exploring multichannel superconductivity in ThFeAsN

We investigate theoretically the superconducting state of the undoped Fe-based superconductor ThFeAsN. Using input from $ab~initio$ calculations, we solve the Fermi-surface based, multichannel Eliashberg equations for Cooper-pair formation mediated by spin and charge fluctuations, and by the electron-phonon interaction (EPI). Our results reveal that spin fluctuations alone, when coupling only hole-like with electron-like energy bands, can account for a critical temperature $T_c$ up to $\sim7.5\,\mathrm{K}$ with an $s_{\pm}$-wave superconducting gap symmetry, which is a comparatively low $T_c$ with respect to the experimental value $T_c^{\mathrm{exp}}=30\,\mathrm{K}$. Other combinations of interaction kernels (spin, charge, electron-phonon) lead to a suppression of $T_c$ due to phase frustration of the superconducting gap. We qualitatively argue that the missing ingredient to explain the gap magnitude and $T_c$ in this material is the first-order correction to the EPI vertex. In the noninteracting state this correction adopts a form supporting the $s_{\pm}$ gap symmetry, in contrast to EPI within Migdal's approximation, i.e., EPI without vertex correction, and therefore it enhances tendencies arising from spin fluctuations.


I. INTRODUCTION
The discovery of superconductivity below T exp c ∼ 30 K in ThFeAsN [1] has added another member to the growing family of Fe-based superconductors exhibiting large critical temperatures [2,3]. Most compounds in this class of materials show magnetic order in the ground state and allow for a Cooper pair condensate only upon sufficient doping or pressure [4,5]. ThFeAsN takes an unusual role in this respect, since the ground state has been reported to already be non-magnetic and superconducting, without the need of external parameter tuning [1,6,7]. Additionally, it has been shown that applying pressure to this system decreases T c [8], which is in contrast to other Fe-based superconductors, such as bulk FeSe [9].
Density Functional Theory (DFT) studies revealed that the electronic structure of ThFeAsN can be considered as quasi-2D, while the Fermi surface (FS) is prototypical for the Fe-based compounds [10][11][12]. It consists of hole-like bands close to the folded Brillouin zone (BZ) center, and electron-like pockets at the corners. The signchange of the superconducting gap between electron and hole pockets is commonly associated with antiferromagnetic spin fluctuations as driving force for the Cooper pair formation [13,14]. As for ThFeAsN, Barbero et al. have found a pressure independent s ± -wave symmetry of the order parameter by performing Knight shift and muon-spin rotation measurements [8]. Similarly, Wang et al. classified this compound as comparable to other Fe-based superconductors with respect to outcomes from magnetic susceptibility measurements, attributing an important role to antiferromagnetic spin fluctuations [1].
Although the discovery of ThFeAsN has triggered some interest from both theory and experiment [6-8, 10-12, 15-20], there has so far not been any attempt to theoretically explain superconductivity in this compound. This might be due to the apparent similarity of ThFeAsN to other Fe-based superconductors, which have been studied in great detail [2,3,14,21,22]. However, the aforementioned ground state properties make ThFeAsN one of very few intrinsically superconducting Fe-based compounds, along with e.g. LiFeAs [23], which makes it worthwhile to be studied in detail. It has been found that FS nesting properties play a crucial role in the family of Fe-based superconductors [24,25], which is why we employ here a FS-based multichannel Eliashberg theory (see e.g. Ref. [26]) to solve for the characteristics of the superconducting state. Our formalism allows to selfconsistently include electron-phonon interaction (EPI), spin and charge fluctuations, all on the same footing.
We find the highest critical temperature as T c ≃ 7.5 K when considering a spin-fluctuations interaction kernel, where only electron bands are coupled to hole bands. Together with an s ± -wave superconducting gap, exhibiting a maximum value of 1.43 meV, this is insufficient to explain the experimental results in superconducting ThFeAsN. Taking into account the full spin-fluctuations kernel suppresses T c to 3.15 K due to phase frustration effects of the superconducting gap, which is a generic property of Fe-based superconductors as recently reported in Ref. [27]. In additional calculations, where we combine the influence due to spin fluctuations with EPI and charge fluctuations, no superconducting state is found for temperatures down to 2 K. Again, the reason lies in a phase frustration of the superconducting gap. When considering the McMillan equation for electron-phonon mediated Cooper pairing, the T c vanishes in agreement with Ref. [12]. Since additional Eliashberg calculations for EPI only do not lead to a finite solution either, it is apparent that the FS-based description employed here is not sufficient to explain superconductivity in this compound.
However, the ratio between characteristic phonon energy scale and minimal band shallowness suggests that vertex corrections to the EPI are important for describing the superconducting and, more generally, the interacting state. We therefore consider an expression for the renormalized vertex function due to non-adiabatic contributions to the EPI, which was derived in Ref. [28]. By following Ref. [29] for a simplified version of this function in the non-interacting state, we obtain the renormalized vertex in a 'one-shot' calculation. The resulting function has, similar to the spin-fluctuation kernel, repulsive long wave vector contributions to the superconducting gap, which support the global s ± -wave symmetry. For small momenta, the vertex corrected EPI counteracts repulsive spin fluctuations contributions, which potentially minimizes the encountered phase frustration problems. We therefore argue that the superconducting state can potentially be explained by including spin fluctuations and vertex corrected EPI in a self-consistent Eliashberg formalism.
The remainder of this paper is organized as follows: Starting with electronic structure calculations in Section II, we obtain two sets of energies that serve as input for the discussion on bosonic interactions in Section III. The presentation of our ab initio calculations for the EPI, Section III A, is followed in Section III B by the description of kernels due to spin and charge fluctuations. Section IV is subject to the discussion of the superconducting state in ThFeAsN, where we solve self-consistent multichannel Eliashberg equations in Section IV A. Additionally we argue why vertex corrections to the EPI are likely to contribute cooperatively with the spin-fluctuations interaction to the gap magnitude and T c , Section IV B. Our conclusions and a brief outlook are presented in the final Section V.

II. ELECTRONIC STRUCTURE CALCULATIONS
The space-group symmetry of ThFeAsN is P4/nmm, with space-group number 129, and we show the twoformula unit cell in Fig. 1. The Fe atoms in this structure are in tetragonal coordination, so that neighboring Fe-As bonds enclose an angle of 109 • . We highlight these angles in Fig. 1 in red and blue colors, and come back to this aspect below when discussing the distance between Fe and As planes.
The electronic energies ξ k,n are calculated using the plane-wave DFT computational package Quantum Espresso [30,31]. We choose the generalized gradient approximation with the exchange-correlation potential by Perdew, Burke, and Ernzerhod [32,33]. Using scalar relativistic Projector Augmented Wave pseudopotentials with non-linear corrections for the core electrons [34], we include spin-orbit coupling in our non-collinear calculation. The momentum space sampling is done on a 24×24×12 Monkhorst-Pack k-grid and we set the kinetic energy (charge density) cutoff to 110 Ry (1100 Ry).
As a first step we relax the structure by minimizing Fe (1) Fe (2) Th (1) Th (2) As (2) As (2) As (2) As (2) As (1) As (1) As (1) As (1)  the forces on the atoms. Setting a convergence threshold of 10 −8 Ry for the electronic self-consistency loop, we find relaxed structural parameters a = 4.0374Å and c = 8.4307Å. The corresponding electron energies are shown in Fig. 2(a) along high-symmetry lines of the BZ. The associated FS, drawn in Fig. 2(b), consists of three hole-like bands that cross the Fermi level around Z − Γ, and two electron bands at A − X. The relaxed structure is then used to calculate the lattice dynamics of the system, as we describe in Section III A.
Before we discuss the ionic degrees of freedom, we note that there is a non-negligible deviation of the obtained lattice parameters when comparing to experimental findings of a exp = 4.0367Å, c exp = 8.5262Å [1]. This observation in ThFeAsN was already made by other authors [10,12]. We therefore perform another electronic structure calculation by fixing a = a exp and c = c exp , the results for the band structure and FS are shown in Fig. 3.
When comparing Figs. 2(b) and Fig. 3(b) it is apparent that the relaxed structure is more dispersive along the k zdirection. This stems from the fact that two energy bands in close vicinity of the Fermi level along Z − Γ change roles, compare panel (a) in both Figures. In Fig. 2(a) we observe a dispersive band along Z − Γ crossing the FS, and a nearly constant one slightly below. This is in contrast to a nearly flat band above the Fermi level on this momentum path, and a dispersive ξ k,n below the FS in Fig. 3(a). The energy band with higher dispersion has dominant Fe-d z 2 orbital character, while the nearly flat band has mainly contributions from Fe-d xy states [10]. As mentioned before, in the ThFeAsN compound Fe is in tetragonal coordination, hence the crystal-field splitting dictates that the d z 2 state is lower in energy than the d xy state. According to this picture we conclude that Fig. 3 represents results for ξ k,n that is closer to experiment. When calculating the electronic energies with relaxed lattice parameters, Fig. 2, these two bands along Z − Γ change role, i.e. the d xy dominated band lies lower in energy than the d z 2 dominated band. This behavior indicates a structural distortion in our unit cell.
The mismatch in the unit cell height (c = 8.4307Å vs. c exp = 8.5262Å) has the effect that the Fe-As tetragons are slightly squeezed, changing the angle between Fe-As bonds from 109 • to approximately 105 • and 119 • (blue and red angles in Fig. 1). As a consequence, the system does not exhibit tetragonal crystal-field splitting, which is the reason that the d xy states fall lower in energy than the d z 2 states. The observed mismatch between relaxed and experimental lattice structure goes in line with DFT calculations on LaFeAsO [35], and has been shown to potentially influence the superconducting state significantly [36].
The results we obtain for ξ k,n (Figs. 2 and 3) are in excellent agreement with previous works [10][11][12]. Despite the structural differences that we discussed above, we employ both, the relaxed and experimental structure in the following when looking at different superconducting channels. We use the experimental lattice parameters when being concerned with purely electronic properties, such as susceptibilities and spin-fluctuation interaction kernels, see Section III. The reason for this is that such quantities require as accurate inputs as possible for yielding the correct characteristics [36]. On the other hand, for electron-phonon calculations we use the relaxed structure so as to have minimized forces on each atom which lead to reliable phonon modes.
The theory applied in later Sections III and IV is based on electronic degrees of freedom close to the Fermi level. We therefore end the current discussion by examining the density of states (DOS) N (E) in a narrow low-energy in- terval in Fig. 4, calculated from experimental lattice parameters. The red curve, representing the full DOS per formula unit, has a characteristic shape for the family of Fe-based superconductors [37,38], where for our undoped system the Fermi level falls into the prototypical dip of N (E) [39,40]. As for other members of this family, most low-energy states originate from the iron atoms, shown as blue line in Fig. 4. The remaining atom species play only a very minor role in the energy interval under consideration. At the Fermi level we find the values N (0) = 1.92 states/eV/f.u. for the total DOS, and N Fe (0) = 1.69 states/eV/f.u. for the iron DOS. These results are in reasonable agreement to related works on this material [7,[10][11][12].

III. BOSONIC INTERACTIONS
Here we discuss the EPI obtained by DFT calculations in Section III A and show how to calculate spinfluctuation kernels within the Random Phase Approximation (RPA) in Section III B. All calculations presented from here on, including Section IV, are carried out with the Uppsala Superconductivity (UppSC) code [41][42][43][44][45], partially with input from the Quantum Espresso package.

A. Electron-phonon interactions
Using the relaxed lattice parameters (compare Fig. 2) we perform Density Functional Perturbation Theory calculations with Quantum Espresso. The k-grid is chosen as in the electronic structure calculation, and we employ a 4 × 4 × 2 Monkhorst-Pack q-grid. The cutoffs for kinetic energy and charge density are fixed at 130 Ry and 1300 Ry, respectively. With an effective Coulomb repulsion of µ ⋆ = 0.136 we find an isotropic coupling strength of λ 0 = 0.40. Using the modified McMillan equation [46] for the superconducting critical temperature, with ω log = 6.09 · 10 −4 Ha, suggests that superconductivity in this system cannot simply be explained by an isotropic, weak-coupling conventional BCS approach. From our ab inito calculations we obtain the branch ν resolved phonon frequencies ω q,ν and the couplings λ (ep) q,ν . The latter are defined in terms of the electron-phonon scattering matrix elements g q,ν : By using bosonic Matsubara frequencies q l = 2πT l, l ∈ Z, we can calculate the momentum and frequency dependent electron-phonon coupling as In Fig. 5 we show λ q,ν for q z = −π in panel (a), and q z = 0 in panel (b). In both cases we observe an enhanced magnitude at (q x , q y ) = (0, 0) and (q x , q y ) = (−π, −π). The contributions at the BZ corners are comparatively smaller for q z = −π. The overall magnitude of λ (ep) q,l=0 indicates a weak-coupling situation. Before using this EPI in our multichannel Eliashberg formalism in Section IV, we first calculate the analogue of λ

B. Spin and charge fluctuations
As mentioned before, we use from here on the electronic energies as presented in Fig. 3. Assuming that the essential physics are captured by processes at the Fermi level, we can write the band resolved bare susceptibility of the system as [26] χ n,n ′ q,l = k δ(ξ k,n )δ(ξ k+q,n ′ ) n F (ξ k,n ) − n F (ξ k+q,n ′ ) ξ k+q,n ′ − ξ k,n + iq l , where we use the Fermi-Dirac function n F (·) and adopt the notation χ n,n ′ q,l = χ n,n ′ q (iq l ). A direct approach for obtaining a band-insensitive bare susceptibility is performing a double summation as χ (0) q,l = n,n ′ χ n,n ′ q,l . However, it is worthwhile splitting theses summations according to certain subsets of energy bands. As encountered in Section II, there are three hole bands and two electron bands crossing the Fermi level. Let us denote the set of such bands by h and e, respectively. We can then define where the label (e − h) means that all contributions involving one electron band and one hole band are included, and likewise for (e − e) and (h − h). Next, we use the bare response of the system to calculate the interaction kernel due to spin fluctuations within the RPA under the assumption of spin-singlet Cooper pairing. Making use of the Stoner factor U we compute with r ∈ {e − e, h − h, e − h, 0}, where the magnetic instability is marked by the condition 1 − U χ (r) q,l → 0. In a similar way the interaction kernel due to charge fluctuations is given by We show the results for λ  q,l=0 , i.e. the result for only coupling electron bands with hole bands. The full interaction kernel is drawn in Fig. 6(d).
We observe in the first two panels of Fig. 6 that λ are both peaked at q = Γ. This behavior is consistent with intuition, since the two electronic bands at the FS are located in close vicinity to each other, and likewise for the three hole bands. Therefore they are coupled via small exchange momenta, for both interband and intraband terms. Such contributions are expected to lead to a suppression of superconductivity, as noted by some of the authors in a recent work on bulk and monolayer FeSe [47]. This 'self-restraint effect' has been discussed convincingly for a more generic model system of Fe-based superconductors by Yamase and Agatsuma [27], and we come back to it in Section IV A. Together with couplings between electron and hole bands, Fig. 6(c), the full interaction has prominent peaks at q = M and a notable hump at Γ, compare Fig. 6(d). How these features influence the superconducting solution is discussed in the following Section IV.

IV. SUPERCONDUCTIVITY OF ThFeAsN
We are now in a position to solve a self-consistent set of multichannel Eliashberg equations in Section IV A, where we use the interaction kernels as introduced in Section III. As described in the following, the theory employed here is not capable of fully explaining the superconducting state in ThFeAsN, even though we consider the EPI, spin-fluctuations and charge-fluctuations channels. Therefore we argue in Section IV B that a proper inclusion of the first vertex correction to the EPI into our Eliashberg theory would likely increase the computed gap magnitude and critical temperature, so as to compare better with experimental findings.

A. Eliashberg calculations
As introduced in Section III B, we distinguish different contributions to the spin and charge-fluctuation interactions, according to the pair of bands they originate from. Consequently, the full kernel to be used in our Eliashberg theory similarly depends on r ∈ {e − e, h − h, e − h, 0}. The electron-mass renormalization Z numerically.
Here we use kernels λ q,l ), where in the following we employ different combinations of spin and/or charge fluctuations, and/or EPI, as will be specified explicitly.
We already stated before that the magnetic instability of the system is marked by choosing the Stoner parameters, such that 1 − U χ (r) q,l → 0. Therefore we can define the maximally allowed value for U as For convenience we specify U in the following as percentage value of its maximum, which depends on the choice of r, so that U = p 100 · U (r) max with p ∈ (0, 100). Further, we assume from here on a 2D system which is justified due to a very non-dispersive k z -direction, compare Fig. 3. We carefully crosschecked that our results do not change qualitatively when taking into account the full 3D system.
Let us start with the conceptually easiest case, i.e. considering only spin fluctuations and setting r = e − h. From Fig. 6(c) we know that the interaction kernel entering in the Eliashberg equations, λ (±,e−h) q,l = ±λ (sf,e−h) q,l , peaks at q = M and is otherwise small in magnitude and featureless. It is hence to be expected that the superconducting gap at the FS changes sign between electron and hole bands. In Fig. 7 we show our numerical solutions to Eqs. (11) and (12), setting T = 2 K, p = 99 and r = e − h. Indeed, we find max k ∆ (e−h) k,m=0 ∼ 1 meV > 0 around the BZ center in Fig. 7(a), while a minimum negative value min k ∆ (e−h) k,m=0 ∼ −1.16 meV < 0 occurs around k = M . This is the prototypical s ± -wave state that is shown here, which is the only stable symmetry for the superconducting gap in ThFeAsN, which is true for any choice of parameters we explored in this work.
The electron-mass enhancement is shown in Fig. 7(b), and takes values between unity and ∼ 2.7. As apparent, the hole bands exhibit clearly smaller values of Z k,m=0 than the electron bands. Consequently, we predict that electron masses at the BZ corners are noticeably higher than at the center. However, it is under question whether such a behavior can be seen in experiment, since our results change as function of U and, as we discuss below, the critical temperature and maximum gap magnitude found here are lower than observed experimentally. The next step is to perform a parameter variation in T and U to see how these quantities affect the solutions of the Eliashberg equations.
In Fig. 8(a) we plot our results for ∆ (e−h) = max k |∆ (e−h) k,m=0 |, found from numerically solving Eqs. (11) and (12), as function of p for various temperatures as indicated in the legend. It is easily observed that the maximum superconducting gap grows nearly linear with p, while the onset of superconductivity is strongly temperature dependent. For all solutions found in the shown parameter range, the gap function has s ± -wave symmetry. The reasons for a growing ∆ (e−h) with p is that the interaction kernel is approaching a magnetic instability as U → U max . As was shown in Ref. [47], the spin fluctuation kernel at the instability wave vector scales like hence an increase in coupling strength is responsible for the enhancement of the gap magnitude with p.
It is also apparent that our theory inherits an upper limit of ∆ (e−h) , as the gap size does not diverge for p → 100. The reason lies again in Eq. (14): The mass renormalization scales approximately like Z k,m=0 ∝ 1 + max q λ (+,e−h) q,l=0 , while the superconducting order parameter has a similar behavior, φ q,l=0 |. Since the gap function is found from ∆ Next, let us turn to the maximum superconducting gap as function of T , which we show in Fig. 8(b), still for the choice r = e − h. Different choices for p are colored as written in the legend. Note, that the values for p are very critical, i.e. we have to consider a Stoner factor very close to U (e−h) max , similar to results shown in panel (a) of the same figure. Here we find that both the maximum superconducting gap and critical temperature increase with p, where the largest possible values are ∆ (e−h) ∼ 1.43 meV and T c ∼ 7.5 K. Naively one might expect that these quantities can be arbitrarily increased by choosing U closer to U (e−h) max , but as discussed in connection to Fig. 8(a) this is not possible due to an upper bound on the gap magnitude. The highest critical temperature found here is significantly lower than the experimental value, T exp c = 30 K [1,16], which indicates that some important ingredients relevant for ThFeAsN are missing in our theory.
We now turn to a different interaction kernel used for solving the Eliashberg equations, namely we include all contributions due to spin fluctuations, setting λ (±,0) q,l = ±λ (sf,0) q,l . When solving Eqs. (11) and (12) as function of U we find that a finite superconducting gap (still s ±wave) can only be found if U is at least p = 99.98% of U (0) max . The results ∆ (0) from three of such very critical choices are shown in Fig. 8(c) as function of T . It is easily observed that, compared to panel (b), both gap magnitude (max. ∼ 0.3 meV) and T c (max. ∼ 3.15 K) are reduced by more than a factor of two. This behavior is, however, not very surprising when considering the shape of the interaction kernel.
As shown earlier, the large wave vector contributions to the spin fluctuations kernel, compare Fig. 6(d), promote an s ± -wave symmetry of ∆ k,m=0 , because they enter the equation for the superconducting gap with an overall minus sign. On the other hand, finite values of λ (sf,0) q,l around the BZ center, which originate from electron-electron and hole-hole band couplings, can numerically induce a different symmetry of ∆ k,m=0 , e.g. a nodal s ± or d-wave state [21]. If the Stoner instability is not approached sufficiently close, the competition between symmetries promoted by λ (sf,0) q∼Γ,l and λ (sf,0) q∼M,l leads to a phase oscillation of the superconducting gap, which does not represent a valid solution, i.e., ∆ = 0 meV. For more details on this self-restraint effect in Fe-based superconductors we refer to Ref. [27].
When including EPI and/or charge fluctuations in the Eliashberg equations we do not find any superconductivity down to T = 2 K. Choosing λ k,m = 0, too, due to the aforementioned reasons. We therefore conclude that results from our self-consistent Eliashberg theory are underestimating T c and the superconducting gap magnitude. Potential cures to this behavior are discussed in the following Sections IV B and V.

B. Possibility of non-adiabatic effects
In the previous Section IV A we showed that the superconducting state in ThFeAsN can neither be explained by spin fluctuations alone, nor does the inclusion of charge fluctuations and EPI lead to a complete picture. This can be explained by phase frustration (i.e., phase oscillation) of the superconducting gap, caused by repulsive spin fluctuation kernel contributions at Γ and M , which both do not support the same global symmetry. We therefore want to discuss here the possibility of vertex corrections to the EPI, that potentially can contribute positively to the superconducting gap size and eventually T c . In a recent work by some of the authors it was shown for a model system that unconventional gap symmetries, such as s ± -wave, can arise from isotropic EPI when taking vertex corrections into account [29]. We therefore want to obtain here an estimate of the vertex function when considering the non-interacting state of the system. Let us assume that the EPI is isotropic to first-order approximation, so that scattering matrix elements are given by a single number g 0 . The characteristic energy scale of the phonon spectrum is ω log ≃ 17 meV (value taken from DFT calculation), from which we can define the isotropic EPI kernel V m−m ′ = 2g 2 0 ω log /(ω 2 log + q 2 m−m ′ ). Taking into consideration the shallowness of our system, ǫ F ≃ 52 meV, we get a non-adiabaticity ratio of α = ω log /ǫ F ∼ 0.33, which is an indicator for the nonnegligible relevance of vertex corrections [29]. For simplicity let us further assume that the system is in the non-interacting state, which translates into an electron Green's function defined in Nambu space that is spanned by a Pauli matrix basisρ i (i = 0, 1, 2, 3). As was shown in Refs. [28,29], the electron self-energy including vertex corrections can be written aŝ We use here the label (0) to stress that we are in a noninteracting situation. By making use of the relation q = k − k ′ , the vertex function can be expressed as with the definitions γ k,m = n ξ k,n /θ k,n,m and θ k,n,m = (iω m ) 2 − ξ 2 k,n . We perform the Matsubara summation over index m ′′ analytically and set m = m ′ , allowing us to visualize the renormalized vertex function below.
Considering that the non-adiabatic equation for the superconducting order parameter contains the interaction kernel −(1 + g 2 0 Γ k,k ′ ,m,m ′ ) [29], we want to find a rough estimate of the scattering strength g 0 by following the approach of Ref. [15]. According to McMillan, the electron-phonon coupling constant can to first order approximation be written as where Θ D is the Debye temperature. Inserting µ ⋆ = 0.136, Θ D = 332 K [7] and T c = T exp c = 30 K [1] gives λ ≃ 1.6. From here we can solve for the scattering strength via g 0 = λω log /2N (0).
We show the two dimensional zero-frequency result for (1 + g 2 0 Γ q,m=0 ) (0) in Fig. 9, obtained at T = 10 K < T exp c . At q = M the vertex correction is negative, while in the remaining BZ our result is strictly positive. These contributions enter the equation for the order parameter in a repulsive and attractive way, respectively. Therefore, in combination with the spin-fluctuations kernel, see Section III B, we get an enhanced repulsive interaction at q = M , which supports the s ± -wave symmetry of the gap. On the other hand, the repulsive small-q contribution of Fig. 6(d) might likely be compensated by an attractive contribution from Γ (0) q,0 around the BZ center, such that the 'self-restraint' effect is minimized [27].

V. DISCUSSION
In summary, a self-consistent FS-based Eliashberg theory, including spin and charge fluctuations, and EPI, can not explain the critical temperature in ThFeAsN as it is observed experimentally. In Section IV A we show that this is due to phase frustration in the superconducting gap due to competing contributions to the interaction kernel. From the strong similarity to other Febased superconductors, especially concerning FS properties, we conclude that such frustration behavior is likely to be generic for this family of compounds. However, even in the most simplified case of solely including spin-fluctuation couplings between electron and hole bands, where no phase frustration occurs, we find a maximum T c = 7.5 K, which is still significantly smaller than T max c = 30 K. As mentioned above, our Eliashberg calculations lead to an underestimation of both T c and superconducting gap magnitude, due to a phase frustration in the gap, caused by competing short and long range wave vector contributions to the spin-fluctuation interaction. It is possible that an inclusion of orbital dependent matrix elements in the calculation of bare susceptibilities would improve this aspect significantly [48,49]. Such matrix elements can alter the momentum structure of the susceptibility, and therefore also of the charge and spin kernels, such that the problematic features at the BZ center might be suppressed. However, at this point we can only speculate about the influence of orbital content on the level of susceptibilities, because the self-restraint effect in Fe-based superconductors was only demonstrated in calculations without such matrix elements [27].
Another potentially important aspect is the level on which the RPA formalism is introduced. Here we use a global Stoner parameter, i.e. only a single quantity to tune in order to find a valid solution for the superconducting state. More general Hubbard Hamiltonians can include additional physics, such as Hundness and intra-/inter-orbital hopping, as e.g. employed in Refs. [50][51][52][53]. Not only can such modeling provide more insights into the important physics of the system, but, additionally, a larger phase space for parameter exploration is likely to allow the correct solution of the interacting state. For a realistic anisotropic, full bandwidth and multi-band/orbital description one might employ the formalism of Ref. [47]. However, a faithful tight-binding fit to the electronic dispersion is crucial for this approach, which is why we did not attempt it here.
The potential boost for T c due to vertex corrected EPI as discussed in Section IV B has its root in the quasi-2D character of the electronic energies, together with a degree of non-adiabaticity, that suggests a treatment of EPI beyond Migdals approximation as important. For representative model systems of high-T c materials, including Fe-based superconductors, it has been shown that a wellnested FS leads to self-consistent unconventional superconductivity, mediated by isotropic electron-phonon scattering [29]. Although only a fully self-consistent calculation in ThFeAsN can test our conjecture, the renormalized vertex function obtained here resembles to a good degree recent model systems [29]. We are hence led to conclude that non-adiabatic corrections are likely effective to enhance spin-fluctuations mediated Cooper pairing in ThFeAsN, while supporting the s ± symmetry of the gap.