The ab initio study of unconventional superconductivity in CeCoIn$_{5}$ and FeSe

Electronic structure and the shape of the Fermi surface are known to be of fundamental importance for the superconducting instability in real materials. We demonstrate that such an instability may be explored by static Cooper pair susceptibility renormalized by pairing interaction and present an efficient method of its evaluation using Wannier orbitals derived from ab initio calculation. As an example, this approach is used to search for an unconventional superconducting phase of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) type in a heavy-fermion compound CeCoIn$_5$ and an iron-based superconductor FeSe. The results suggest that the FFLO superconducting phase occurs at finite magnetic field in both materials.


Introduction
In the standard theory of superconductivity developed by Bardeen, Cooper and Schrieffer (the BCS theory) introduced in 1957 [1,2] the concept of the Cooper pairs which are in a singlet state with zero total momentum plays a fundamental role. Shortly after this discovery in 1964 two independent groups proposed an unconventional superconducting state with the Cooper pairs having non-zero total momentum. First group of Fulde and Ferrell proposed a state where the Cooper pairs have only one possible momentum [3], while the second one of Larkin and Ovchinnikov assumed that pairs with two opposite momenta exist in superconducting state [4]. Thereafter, materials where such unconventional Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states could be realized have been looked for very intensively.
In the absence of an external magnetic field, the Fermi surfaces for electrons with opposite spins are similar. Then, in conventional BCS-type superconductor, every electron in state (k, ↑) is paired with an electron in state (−k, ↓) and the total momentum of this pair equals zero. Here k is the momentum of an electron with spin σ =↑, ↓. For this reason one can say that the main source of unconventional superconductivity of the FFLO-type is a shift of the Fermi surface. In such a case, an electron in state (k, ↑) can be paired with an electron in the state with shifted momentum, (−k+q, ↓). Then, the total momentum of the Cooper pairs in this case is non-zero and equals q. Indeed, the main source of the FFLO phase can be an external magnetic field which leads to the Zeeman effect and to the occurrence of the shift of the Fermi surfaces for electrons with opposite spins. However, also other effects such like the mass-imbalance in the system can lead to the shift of the Fermi surfaces [5][6][7][8], which is used in the recent experiments for realisation of the FFLO phase in ultra-cold fermion gases on optical lattices.
According to the previous theoretical studies, a material where one can expect a realization of the FFLO phase has to fulfill a few restrictions and conditions. In such a system the main factor determining the upper critical magnetic field has to be given by Zeeman (paramagnetic) effect whereas the orbital (diamagnetic) effect should be less important or negligible. Such a property can be associated with the Maki parameter, α M = √ 2H orb c2 /H P c2 [9], which describes the ratio of the critical magnetic fields at zero temperature given only by diamagnetic (H orb c2 ) and paramagnetic (H P c2 ) effects. This suggests that the interesting systems in the context of the occurrence of the FFLO phase are materials where α M ≥ 1 (so-called Pauli limited systems). Moreover, the FFLO phase can be realized at low temperature and high magnetic field only. Under these conditions, in the absence of orbital effects, one expects the first order phase transition from the superconducting FFLO state to the normal (non-superconducting) phase.

Heavy fermion systems
CeCoIn 5 is a representative heavy fermion system. This compound exhibits a layered structure (its crystal structure is presented in the left panel of figure 1). The heavy fermion systems are characterized by relatively high Maki parameter, which is estimated as α M ∼ 5 [22]. The symmetry of superconducting gap is in agreement with d-wave-type pairing [23,24].
A possible existence of the FFLO phase in this compound can be concluded from several untypical physical properties. Many experiments provide evidence in favour of the first order phase transition from the superconducting state to the normal state for the magnetic field parallel to the main crystallographic directions (H ab and H c). It has been shown in the measurements of a jump of thermal conductivity [24,25], magnetization [26], penetration depth [27], ultrasound velocity [28], thermal expansion [29], magnetostriction [30], in NMR experiments [22,31,32] and from the shape of specific heat [33][34][35][36]. Moreover, the experimental results show an existence of incommensurate spin-density wave (SDW) state in a regime of occurrence of superconducting state [37][38][39]. This cannot be explained by a simple assumption of the existence of superconducting BCS state [40][41][42][43]. The presence of the SDW also enhances a tendency of the system to stabilize the FFLO phase [44]. The chemical structure of superconducting heavy-fermion system CeCoIn 5 (left) and iron-based superconductor FeSe (right). The image was rendered using VESTA software [21].

Iron based superconductors
Superconductivity in the iron-based superconductors was discovered in 2008 by Kamichara et al. in LaFeAsO doped by F at the oxygen site below 26 K [45]. This new class of high-temperature superconductors has been intensively investigated in the last few years [46][47][48]. These materials are characterized by iron-arsenide or ironselenide layers [46,49]. As a consequence, they have characteristic Fermi surfaces (with coexisting hole and electron pockets around Γ and M point of the first Brillouin zone). Layered structure and quasi-two-dimensional character of the Fermi surfaces in many representative systems from the family of iron-based superconductors make these materials placed in the Pauli limit. The Maki parameter in this group of compounds can be estimated as α M ∼ 1 − 5 [50][51][52][53][54][55][56][57][58][59][60][61][62][63]. Thus, the first order phase transition from the superconducting to the normal state and high anisotropy of upper critical magnetic field have been reported in many experiments for different ironbased superconductors (e.g. Ba 1−x K x Fe 2 As 2 [50][51][52][53], LaO 0.9 F 0.1 FeAs 1−δ [54,55], LiFeAs [56][57][58], The simplest representative of this group of materials is FeSe (right panel of figure 1), with the critical temperature T c = 8 K [64]. It undergoes a structural transition from the tetragonal to the orthorhombic phase at T s 87 K with no magnetic order [65,66], but with emerging nematic electronic structure [67]. On the other hand, the calculations based on density functional theory (DFT) predict magnetic states both for FeSe [68] and FeSe 1−x bulk materials [69]. They indicate the proximity of FeSe to the magnetic instability, which is not observed due to electron correlations beyond the DFT [70,71]. Possible consequences of this proximity have been debated recently and exotic states with the hidden magnetic order have been proposed: antiferroquadrupolar spin order [72] and nematic quantum paramagnetic phase [73]. The absence of magnetic order in FeSe makes this material a good candidate to test the possibility of electron-pairing mechanism of superconductivity, which is still unknown in iron-based superconductors [49,74,75]. In fact, it is likely that magnetic interactions at orbital degeneracy play an important role [76]. The coexistence of superconducting and antiferromagnetic orders in some of iron-based superconductors suggests the existence of superconductivity with s ± symmetry [77] where the superconducting gap has opposite signs on the Fermi surface around the Γ and M points of the first Brillouin zone. However, for nonmagnetic FeSe materials this hypothesis is still under debate [78,79].
In the context of the present paper, there are very important observations of the phase transition in high magnetic field (above 13.5 T) and at low temperatures (below 1 K) in pure single crystals of superconducting FeSe reported by Kasahara et al. [63]. Namely, in the strongly spin-imbalanced state the zero-resistivity (superconducting) state has been observed together with an additional phase transition, what can be interpreted as the prediction of an occurrence of the FFLO phase.

Motivation
As we showed above, relatively large number of various chemical compounds manifest properties which are typical for systems where one can expect the FFLO phase occurrence. Their studies enforce a quest for new theoretical techniques. In this paper, we propose a method for studying the properties of superconducting states using a combination of the ab initio (DFT) and the Cooper pair susceptibility calculations.
The use of DFT leads us to design realistic models to capture the relevant part of the electronic structure of materials, while the Cooper pair susceptibility calculations show a tendency of the systems towards different types of superconductivity without defining its source. The theoretical background is described in detail in section 2, whereas the methods used are explained in 3. Numerical results and their discussion are presented in section 4 for two exemplary systems: a heavy fermion compound CeCoIn 5 and an iron-based superconductor FeSe. The summary and general discussion are presented in section 5.

Theoretical background
In a general case the band structure of the system can be represented by the noninteracting Hamiltonian in the diagonal form: where c † εkσ (c εkσ ) is creation (annihilation) operator (in momentum representation) of an electron in band ε with momentum k and spin σ. Similarly, E εkσ denotes the band energy for an electron from band ε with momentum k and spin σ. For nonmagnetic states considered below these energies are equal for both directions of electron spin σ =↑, ↓.
To study an unconventional superconducting state, we define the static Cooper pair (superconducting) susceptibility χ ∆ ε (q) [44,[80][81][82]: where N is the number of sites, · · · r ω is the retarded Green's function and ∆ ε (i) = j ϑ(j −i)c εi↑ c εj↓ defines the intraband Cooper pair annihilation operator in the real space at site i, while the vector q is total momentum of this pair. The factor ϑ(j − i) defines the type of the pairing interaction in the real space and corresponds to the symmetry of the order parameter in the momentum space [80]. For example, for the on-site s-wave pairing it is given as ϑ(j − i) = δ ij , whereas for the d-wave pairing ϑ(j − i) = δ i±x,j − δ i±ŷ,j , wherex andŷ are unit vectors in the x-and y-direction of the lattice.
In the momentum space, the static Cooper pairs susceptibility can be rewritten in a form where the Green's function is f (ω) = 1/{1 + exp(ω/k B T )} is the Fermi-Dirac distribution, and η(k) is the formfactor describing the symmetry of the order parameter in the momentum space -it could be equal to 1, 2(cos k x +cos k y ), 4 cos k x cos k y , 2(cos k x −cos k y ), or 4 sin k x sin k y , for s, s x 2 +y 2 , s x 2 y 2 (s ± ), d x 2 −y 2 , and d x 2 y 2 symmetry, as shown in figure 2. In a case of the cylindrical Fermi surface centered at the Γ point and d-wave pairing, as consequence, one can expect the nodal line in the gap (see 2(c) and 2(d)). This has measurable consequences as then the properties of the superconducting phase are quite different [83,84]. In the following of this paper the results presented are obtained for η(k) = 1. From equation (2) we find that the superconducting susceptibility, χ ∆ ε (q), can be associated with the effective interaction H SC describing a superconducting state which in the basis of band operators {c † εkσ , c εkσ } can be given in a form of a phenomenological BCS-like term: where U ε ∆ εk denotes an energy-gap function in a superconducting state (for a given symmetry of the order parameter in the momentum space) [80,81]. The renormalised static Cooper pair susceptibility (taking into account the effective U ε interaction) is given in the random phase approximation as Here, χ ∆ ε (q) is a superconducting susceptibility (2) at U ε = 0 (in the normal state), whileχ ∆ ε (q) is a susceptibility in the presence of the effective pairing interaction U ε = 0 (in the superconducting state) in a given band ε.
Similarly like in the studies of magnetic properties using the Lindhard spin susceptibility [85], a divergence of the susceptibilityχ ∆ ε (q) (6) suggests the most likely state, while U ε −1/χ ∆ ε (q) can be treated as a minimal value of interaction U ε needed to induce the superconductivity. As shown below, in the absence of the external magnetic field, χ ∆ ε (q) has a maximal value for q = 0 suggesting a tendency to stabilise the BCS state. In the context of the FFLO phase, it is of interest to consider different values of q to establish for which one χ ∆ ε (q) has a maximal value in the presence of the external magnetic field. In such a case, a distinct maximum for q = 0 can suggest a possibility of realisation of the FFLO phase in the system. For this reason we calculate and compare superconducting susceptibilities in the absence as well as in the presence of the external magnetic field.
Because the physical properties of superconductors are very sensitive to their electronic structure (see e.g. [86,87]), it is crucial to describe it as accurately as possible. In this paper, we present for the first time the calculations of the superconducting susceptibility (2) combined with the ab initio (DFT) band-structure approach on the examples of two systems: CeCoIn 5 and FeSe.

Conditions imposed on DFT calculations
We remark that the previous calculations of the superconducting susceptibility were performed in the tight binding models [44,[80][81][82]. It was found that a relatively small external magnetic field corresponds to a large total momentum of Cooper pairs in the FFLO state. For example, in a case of the two band model presented by Raghu et al [88], the total momentum of the Cooper pairs was estimated as ∼ 0.3/a [80,89], whereas in a case of three band model proposed by Daghofer et al [90,91] its value was given as ∼ 0.05/a [81] (where a is lattice constant). On the contrary, for a case of realistic values of model parameters such like external magnetic field in ironbased superconductors one can expect a very small total momentum of Cooper pairs, with value ∼ 0.002/a [92]. As one can see, the results for momenta of pairs in the FFLO state depend strongly on models and parameters used. For this reason, in the calculation method of the superconducting susceptibility χ ∆ ε (q) used in this work, the DFT calculation with an extremely dense k-grid mesh is employed. The choice of the grid size has a significant influence on the accuracy of numerical results (what will be shown in the next paragraph). Therefore, there are significant differences from calculations of spin [86,[93][94][95] or charge [96] susceptibilities, where the nesting vectors are of the same order of magnitude as the Fermi vectors.
For the reasons mentioned above, in this paper the main idea of calculations of the Cooper pairs susceptibilities is to connect them with a realistic description of band structures of studied materials (i.e., CeCoIn 5 and FeSe), given by equation (1). To reproduce the band structure of real materials we have performed DFT calculations (section 3.2) on a smaller k-grid. The results are used in order to construct the tight binding model in Wannier orbital basis (more details can be found in section 3.3). The use of the constructed tight binding model allows us to find the dispersion relation on the sufficiently dense k-grid, which is much denser than those used in DFT calculations. In the calculations of the Cooper pair susceptibility presented in this work we will use an effective k-grid (obtained from the tight binding model) with dense ∼ 10 4 times bigger than that used in typical DFT calculations. This helps us to increase the accuracy of calculations what is important due to the following two reasons: (i) to calculate the Cooper pairs susceptibility accurately it is necessary to use relatively dense k-grid (while to get the band structure directly from the DFT calculation on the dense enough k-grid would take a very long time); (ii) a relatively small value of the total momentum of the Cooper pairs in the FFLO phase [89,92] is expected and thus a step in k-grid as small as possible is needed what increases substantially the number of k-points used in calculations.
In calculations one has to consider also the external magnetic field which can be a source of the FFLO phase in studied materials. For layered materials, where orbital effects can be neglected, it can be simply done by using the Zeeman term in Hamiltonian (for more details see section 3.4). In the ab initio calculations it is equivalent to non-equal numbers of electrons with opposite spins.

Details of the ab initio DFT calculations
The main DFT calculations have been carried out using the Quantum-ESPRESSO software [97,98] within the generalized gradient approximation (GGA) [99]. The interactions between the core and the valence electrons are implemented through the projector augmented-wave (PAW) method [100] employing PSlibrary [101,102] pseudopotentials. Perdew and Wang (PW91) parameterized GGA functionals [103] are used to describe the exchange-correlation interactions. Before we determined the corresponding tight binding models we have calculated the band structures for both compounds according to the following scheme: All DFT calculations of band structures presented in this paper have been performed over 10 × 10 × 10 k-points Monkhorst-Pack mesh [105] and the cut-off energy for the plane waves expansion was equal to ∼ 1088 eV (80 Ry) and ∼ 950 eV (70 Ry) for CeCoIn 5 and FeSe, respectively.
Experimental (from [64,104]) and relaxed (obtained from DFT calculations presented in this paper) atomic coordinates for CeCoIn 5 and FeSe. The upper part of the table presents atomic coordinates as a function of parameter z (cf. figure 1).
In the DFT calculation, we have used lattice constants obtained from the experimental measurements (table 1). As a result of relaxation of In and Se atoms, we have found a very good compatibility with the experimental atomic coordinates for these atoms (see table 2).

The tight binding model in Wannier orbitals basis
Using the results of the DFT calculation for electron band structure one can find the tight binding model in the basis of the Wannier orbitals, which are located on selected atoms [106]. It can be performed by using the Wannier90 software [107,108] which finds the tight binding model in a base of the maximally localized Wannier functions (MLWF). As a result of this step one gets a tight binding Hamiltonian of the electrons with creation operators {c † Rµσ } on the lattice which is given by where t µν R,R are hopping elements between orbitals µ and ν localized on the atoms at sites labeled with R and R . The matrix of the normal state (i.e., non-

superconducting) Hamiltonian in the momentum space reads:
where δ = R − R is the distance (in the real space) of hopping t µν δ ≡ t µν R,R . A band structure for given k-points, i.e., E kεσ appearing in equation (1), can be found by diagonalization of the matrix H(k).
The band structures obtained from the DFT calculations (cf. section 3.2) for CeCoIn 5 and FeSe are shown by dotted red lines in figures 3(a) and 3(b), respectively). To describe accurately superconducting properties of studied compounds it is necessary to have a good description of states near the Fermi level whereas states far above (below) are irrelevant. The Wannier-based tight binding model [108][109][110] has been found from 10 × 10 × 10 full k-point DFT calculation for random projected 25 orbitals for CeCoIn 5 and all 10 d-states in Fe in a case of FeSe. For both cases the convergence tolerance of the invariant spread, n r 2 n − r 2 n , of the Wannier function in real space [108][109][110] is equal to 10 −10 (in a calculation within the defined energy window). As the result we find 25-and 10-orbital tight binding models for the description of CeCoIn 5 and FeSe, with spread smaller that 3.1Å 2 and 5.0Å 2 , respectively. The band structures obtained from tight binding models are presented in figure 3 by black solid lines.
In a general case, in order to find the tight binding model we can generate the MLWF in an appropriate energy window [107,108]. Because for FeSe its band structure obtained from DFT calculations is given by well separated bands around the Fermi level, see figure 3(b), we can set the energy window from −3.0 eV to 2.7 eV (with respect to the Fermi level). Thus the tight binding model found reproduces the DFT data very well for this iron-based superconductor. In our approach we restrict ourselves only to the separated states (the energy gap is located below and above these states). As a consequence the bands in this energy range are described by the same number of maximally localized orbitals (as number of bands). These orbitals does not give any contribution to the other band beyond the energy window (which does not have any effects on the problems studied in the paper, because they are located far away from the Fermi level). A situation is more complicated for CeCoIn 5 which has more complex band structure with larger number of bands, see figure 3(a). For a better description of states near the Fermi level we decided to find a larger number of the MLWF (containing also fully occupied bands below the Fermi level). As a consequence, the results reproduced by the tight binding model are in a good agreement with those obtained within the DFT calculation up to energy 0.5 eV above the Fermi level, whereas the unoccupied bands above it are not perfectly represented. These issues do not affect our main results which concern superconductivity in the system and thus the location of states in the neighborhood the Fermi level is the most significant.

External magnetic field
The main source of the shifted Fermi surfaces can be external magnetic field given by the Zeeman term. It relatively well describes the influence of the magnetic field on electrons for a case of Pauli limited materials where the orbital effect can be neglected. This situation occurs also for the external magnetic field applied parallel to the layers of materials, with a weak coupling between them. As a consequence, the magnetic field has been taken into account in the form of the Zeeman coupling term added to the Hamiltonian which is given by where g 2 is the gyromagnetic ratio, µ B is the Bohr magneton and H is the external magnetic field (in Teslas). Notice that Zeeman term (9) is diagonal and thus it can be easily included in equation (3). Similar approaches associated with addition of terms in the Hamiltonian have been successfully used for the description of superconductors, with e.g. impurities [111,112], orbital ordering [113,114], or spin-orbit coupling [115,116].

Heavy fermion system: CeCoIn 5
In this section we discuss the results for the heavy fermion CeCoIn 5 compound obtained within the scheme described above. The band structure obtained for this compound has been presented on the figure 3(a). By using our 25-orbital tight binding model we find the Fermi surface of CeCoIn 5 , which is built from pockets originating from only three different bands. It is presented in figure 4. Our results are in an agreement with the previous DFT studies presented for this compound [117][118][119][120][121][122]. The Fermi surface consists of several complex and anisotropic pieces. There are the edge-centered, more or less cylindrical parts of the electron-like (two-dimensional) sheets of 2nd and 3rd bands, centered at (π, π, k z )-points and hole-like sheets centered at the Γ point. In this case, the Fermi surface is made by only three out of 25 bands. As a consequence, in the following studies of CeCoIn 5 we concern only three of the Cooper pairs susceptibilities (for every Fermi pocket separately), due to the fact that the filled (empty) bands do not affect the superconducting properties of the system.
Next, we are able to calculate the Cooper pairs susceptibility χ ε (q) defined by equation (2) with or without the magnetic field using a relatively sparse 50×50×50 kgrid mesh. As we have stated above, the location of the maximal value of the Cooper pairs susceptibility in the momentum space is the most important here because it contains information about preferred momentum of the Cooper pairs realised in the system. In the absence of magnetic field maximum values of susceptibilities χ ε (q) for each band ε = 1, 2, 3 are located at the center of the first Brillouin zone (at the Γ = (0, 0, 0) point), whereas in the presence of magnetic field they are located near the Γ point and they are not zeros (at least one q α = 0, α = x, y, x).
It should be noted that the result can be visible only for relatively very large (and nonphysical) magnetic field (here the calculations were performed for H = 200 T). In this case, the maxima of χ ε (q) are located at |q| ∼ ±2π/20. It is a consequence of a relatively small density of the k-grid used in the calculations (i.e., 50 × 50 × 50 k-grid mesh in this example). To obtain more realistic values of magnetic field at which the FFLO instability can occur, i.e., these from the FFLO phase regime, H ≈ H c which equals approximately 15 T, we have to perform calculations using an extremely dense k-grid.
Thus, for more accurate description of the system in smaller magnetic field which may be accessible experimentally, we used the increased mesh of the k-grid in the following way: . Fermi surfaces for CeCoIn 5 obtained from the tight binding model. The first row presents the total Fermi surface originating from all bands (views from the top -parallel to the z axis and from the front -perpendicularly to the z axis, on left and right panel, respectively). Panels in the second and third rows (view from the top and view from the front, respectively) present the pockets given by each band, separately. Results from the tight binding model used in this paper are obtained for 50 × 50 × 50 k-grid mesh. The black boxes denote the renormalized first Brillouin zone k i ∈ (−π, π) (i = x, y, z).
• we find the band structure with a extremely dense grid (i.e., 2000 × 2000 k-point) in a xy-plane of the total momentum of the Cooper pairs q using the tight binding Hamiltonian, and then • we sparse the grid (by 10 k-points) in the direction perpendicular to q. As a result one is able to find locations of maxima of χ ∆ ε (q) with use of 2000×2000×10 k-points in the calculations. We conducted the calculations for the total momenta q of Cooper pairs which are located on xy, xz and yz planes of the momentum space. The results for a case of external magnetic field equal to 20 T are shown in figure 5. As one can see the maximum values of the χ ∆ ε are located at finite values of momenta, at least for the planes considered. Thus one can conclude that, similarly like previously at H = 200 T, the maxima of χ ∆ ε occur for non-zero q near the Γ point in each band Fermi surfaces for FeSe obtained from the tight binding model. The first row presents the total Fermi surface (views from the top -parallel to z-axis and from the front -perpendicularly to z-axis, on left and right panel, respectively). Panels in the second and third rows (view from the top and view from the front, respectively) presents the pockets given by each band, separately. The results are obtained from the tight binding model used in this paper for 50 × 50 × 50 k-grid mesh. The black boxes denotes the renormalized first Brillouin zone k i ∈ (−π, π) (i = x, y, z).
considered. Let us emphasize that the maximum of χ ∆ ε in the second band has a relatively large value when it is compared with those of the other bands.

Iron based superconductor: FeSe
In this section we present the results for the bi-atomic FeSe compound from the family of the iron-based superconductors. The band structure for the compound have been presented in figure 3(b). The four bands compose the Fermi surfaces which is obtained from a 10-band tight binding model and it is presented in figure 6. Two hole-like bands are centered at the Γ point and two electron-like bands at the M point. The calculated band structure along the Γ-Z line suggests that, because the dispersion of the bands is linear in k z -direction this system can be treated as a two dimensional one. This result agrees qualitatively with the previous calculations [67,[123][124][125][126][127][128][129][130].
Similarly like in CeCoIn 5 , in the absence of the external magnetic field χ ∆ ε (q) for FeSe clearly shows a tendency of the system to realise the homogeneous superconductivity of the BCS-type in each band. In a case of calculations with the sparse k-grid mesh (i.e., 50 × 50 × 50) results obtained in a presence of external magnetic field are similar to those for CeCoIn 5 . When one uses too small accuracy for a grid of k-points the results which are qualitatively different than the BCS (i.e., the results with maxima of χ ∆ ε (q) located not at the Γ point), are possible to find only in extremely large nonphysical magnetic fields. Thus, similarly like in calculation for the previous heavy-fermion compound, we use dense grid. In this case in the non-zero magnetic field, maxima of susceptibilities χ ∆ ε (q) in each band among these building the Fermi surface move to finite values of the Cooper pair momentum q = q max = 0. The results for H = 20 T external magnetic field are shown in figure 7 which can be treated as a good approximation of the real value of the critical magnetic field of this compound.
We would like to emphasize that the DFT predicts low-energy bands which deviate quantitatively from the ones observed, e.g. by angle-resolved photoemission (ARPES) [131,132]. The both electron and hole Fermi pocket sizes are smaller than those obtained from the DFT calculations. This corresponds to the effective shift of the Fermi level in Γ and M points [132]. However, this shift should not change qualitatively the results obtained. The correct band structure can be found by the dynamical mean field theory (DFT+DMFT) calculation [133,134], but this method gives blurred bands (i.e., the bands with finite lifetime) so the results cannot be directly used for the calculation of superconducting susceptibility.
Additional remark. As it was written in section 2, 1/χ ∆ ε (q) can be connected with a minimum value of interaction U ε which is needed to induce superconductivity with the Cooper pairs with momentum q in the system. In a case of the BCS state where q = 0, a value of χ ∆ ε (0) can be found using a sparse k-grid mesh. A situation is more complicated when one considers the FFLO state, where q = 0, but |q| is relatively small. Such conditions can be expected for experimentally available (realistic) magnetic fields. High critical magnetic fields of the order of 100 − 200 T have been found in many compounds, even for some from the group of iron-based superconductors [135]. But such large upper critical fields are given rather by the orbital effect, what excludes the existence of the FFLO phase in such systems. The second issue which should be stressed is the fact that a maximum value of χ ∆ ε in the presence of magnetic field is smaller than that obtained in the absence of the field (they occur for different momenta q: q = 0 and q = 0, respectively). As a consequence, the magnitude of the effective pairing interaction leading to the BCS phase, U ε ∼ −1/χ ∆ ε (q = 0), is bigger than the effective pairing interaction which triggers the FFLO state (because χ ∆ ε (q = 0) < χ ∆ ε (q = 0)). Thus, for U ε the FFLO phase can be induced whereas the BCS-type superconductivity cannot be realised.

Summary
In this paper we proposed a method for studying the properties of superconducting states using a combination of the ab initio density functional theory (DFT) and the Cooper pair susceptibility calculation. We have applied this method to study a tendency to stabilise the unconventional superconductivity of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) type in two compounds: CeCoIn 5 and FeSe. Using the realistic band structure determined by the DFT calculations we have derived the respective tight binding models and obtained the static Cooper pair susceptibility. We show that the conventional BCS-type superconductivity is favoured for each band in the absence of magnetic field whereas the presence of magnetic field can stabilize the inhomogeneous FFLO phase. Similarly like in the previous theoretical works [92], the total momenta of Cooper pairs are very small.
The present calculations show that in a general case the total momentum of Cooper pairs q changes nonlinearly with magnetic field [92]. Moreover, we can expect that the FFLO superconducting phase can survive for a case of the band with strong linearity near the Fermi energy because hot-spot points k and (−k + q) can be found there what causes a growth of χ ∆ ε (q) for the non-zero total momentum of Cooper pairs (q = 0). The total momenta of the Cooper pairs which can be realised in this system are in good agreement with the previous studies [92].
The calculations have also been performed for other symmetries of the superconducting gap displayed in figure 2 as well, i.e., s x 2 +y 2 , s x 2 y 2 (s ± ), d x 2 −y 2 , and d x 2 y 2 . The results are qualitatively similar to these discussed above for the swave superconductivity and give enhanced Cooper pair susceptibility which favours the FFLO phase. In the case of the d-type symmetries, the numerical results of χ ∆ ε (q) found in the presence of external magnetic field are not distinct and maxima of χ ∆ ε (q) are fuzzy. Moreover, only ground state calculations for various superconducting states can resolve the question concerning the realisation of a particular phase [80].
In summary, we would like to emphasize that the pair susceptibility calculations presented here can be a useful tool to investigate a tendency of the system to stabilize some kind of superconducting phase (in this work we discuss the BCS and FFLO phases). The calculations for different symmetries of the superconducting gap were performed and they gave results which are qualitatively similar to those obtained for the case of the s-wave symmetry. Indeed, they give very similar dependence of the superconducting susceptibility on the Cooper pair momenta and the external magnetic field. Therefore the results obtained for the s-wave symmetry are generally valid and we conclude that the electronic structures of both CeCoIn 5 and FeSe support the existence of the FFLO superconducting phase.